Parallelized Simple Random Constrained Portfolio Generation

Generate an efficient frontier of random portfolios with custom constraints in R. Accelerate your project with the power of cloud computing using Techila Technologies’ solution in combination with R’s plyr functions. Source code included.

The embedded source code below is available due to the generous sponsorship of Techila Technologies.

  1. #############################
  2. # Set Environment Variables #
  3. #############################
  4.  
  5. setwd("C:/PathToWorkingDirectory")
  6. Sys.setenv("JAVA_HOME"="C:\\Program Files\\Java\\PathToJava")
  7. Sys.setenv("TECHILA_SDKROOT" = "C:/PathToTechila")
  8.  
  9. ##################
  10. # Load Libraries #
  11. ##################
  12.  
  13. library("rportfolios") # Random portfolio generation
  14. library("plyr") # Merging and applying functions
  15. library("quantmod") # Downloading returns
  16. library("PerformanceAnalytics") # Charting and calculating return metrics
  17. library("ggplot2") # Visualizing data
  18. library("reshape2") # Melting objects, passed into ggplot
  19.  
  20. ########################
  21. # Techila Installation #
  22. ########################
  23.  
  24. library("rJava")
  25. library("R.utils")
  26. library("techila")
  27. registerDoTechila()
  28.  
  29. # Set TRUE to run in parallel with Techila, FALSE to run locally
  30. isParallel <- T
  31.  
  32. # Verify Techila is working
  33. foreach(b=1:4, .combine="cbind") %:%
  34.   foreach(a=1:3) %dopar% {
  35.     a* b
  36.   }
  37.  
  38. ####################
  39. # Data Preparation #
  40. ####################
  41.  
  42. # Now use sector ETFs as tickers
  43. Tickers <- sort(c("XLU", "XLV", "XLI", "XLY", "XLK", "XLF", "XLP", "XLE", "XLB"))
  44. numAssets <- length(Tickers) # How many assets do we have?
  45.  
  46. # Assign dates to set range for stock quotes
  47. sDate <- as.Date("2000-01-01")
  48. eDate <- Sys.Date()
  49.  
  50. # Download the data for each ticker
  51. # library("quantmod")
  52. # library("plyr")
  53. tickerData <- llply(Tickers, function(ticker, sDate, eDate){
  54.   getSymbols(ticker, src = "yahoo", from=sDate, to=eDate, auto.assign = F)
  55. },.parallel = isParallel, sDate=sDate, eDate=eDate, .paropts=list(.options.packages=list("plyr","quantmod"))) 
  56.  
  57. # Merge and calculate Returns
  58. # library("PerformanceAnalytics")
  59. stockMerge <- do.call(merge, lapply(tickerData, Ad)) # Take adjusted close prices
  60. stockReturns <- ROC(na.omit(stockMerge), type="discrete")[-1,] # Get rid of first NA return
  61. colnames(stockReturns) <- Tickers # Rename the columns
  62.  
  63. # Chart cumulative returns of our assets
  64. charts.PerformanceSummary(stockReturns, main="Individual Stock Cumulative Performance", geometric = T)
  65.  
  66. # How many simulations should we run?
  67. numSims <- 1000 # 100 max for non parallel
  68.  
  69. ### One Variable ###
  70.  
  71. # library("rportfolios")
  72.  
  73. # Generate a random portfolio return timeseries given parameters
  74. randomSim <- function(seed, numAssets, probLong, maxWeight, stockReturns){
  75.   weight <- random.general.test( n=numAssets, p=probLong, x.u = maxWeight)$x
  76.   portRet <- as.xts(apply(sweep(stockReturns, 2, weight, FUN="*"),1,sum))
  77.   return(portRet)
  78. }
  79.  
  80. # Dynamically plot the return timeseries as the code is runinng
  81. plotPrices <- function(preresult){
  82.   invisible(lines(as.zoo(cumprod(1+preresult[[1]])-1), col=adjustcolor("black", alpha.f=0.1))) # Realtime Charting
  83.   return(preresult)
  84. }
  85.  
  86. # Verify function works correctly
  87. test <- randomSim(1, numAssets, 1, 0.5, stockReturns)
  88. head(test)
  89.  
  90. # Generate base plot
  91. plot(as.zoo(cumprod(1+test)-1), ylim=c(-1.0,10), ylab="Cumulative Return", col=adjustcolor("black", alpha.f=0.0))
  92.  
  93. # Now generate portfolios and dynamically generate the chart, capture the pre-results
  94. result <- maply(.data=1:numSims,
  95.                 .fun=randomSim,
  96.                 .parallel=isParallel, 
  97.                 probLong=1, 
  98.                 numAssets=numAssets, 
  99.                 maxWeight=0.50,
  100.                 stockReturns=stockReturns,
  101.                 .paropts=list(.options.callback=plotPrices,
  102.                               .options.packages=list("rportfolios","xts"),
  103.                               .options.steps=10#,.options.projectid=206
  104.                 ))
  105.  
  106. # Combine results into a zoo object and re-index
  107. result.zoo <- as.zoo(t(result))
  108. index(result.zoo) <- index(stockReturns)
  109.  
  110. # Visualize results some more
  111. charts.PerformanceSummary(result.zoo[,1:10], main="Cumulative Performance of Some Simulations", geometric = F)
  112.  
  113. ### Two Variables / No Live Charting ###
  114.  
  115. # Vector of maxweight parameters to iterate through
  116. maxWeights <- c(1/numAssets, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50)
  117.  
  118. # All possible combinations of parameters  
  119. paramGrid <- expand.grid(seed=1:numSims, maxWeight=maxWeights)
  120. dim(paramGrid)
  121.  
  122. # Run simulation in Techila
  123. result2d <- maply(.data=paramGrid,
  124.                   .fun=randomSim,
  125.                   .parallel=isParallel, 
  126.                   probLong=1, 
  127.                   numAssets=numAssets,
  128.                   stockReturns=stockReturns,
  129.                   .paropts=list(.options.packages=list("rportfolios","xts")
  130.                                 #,.options.steps=10
  131.                   ))
  132.  
  133. # Melt objects before charting with ggplot
  134. #library(ggplot2)
  135. #library(reshape2)
  136.  
  137. # Pre-allocate
  138. meltedPath <- NA
  139.  
  140. # Loop through each weight parameter
  141. for(i in 1:length(maxWeights)){
  142.  
  143.   # Extract the simulation returns
  144.   pathrets <- data.frame(t(result2d[,i,]))
  145.  
  146.   # Calculate cumulative returns and reindex the path
  147.   path <- as.zoo(apply(pathrets+1,2,cumprod)-1)
  148.   index(path) <- index(stockReturns)
  149.  
  150.   # Change frequency to monthly and re-index
  151.   path.monthly <- to.monthly(path, OHLC=F)
  152.   path.monthly.df <- data.frame(path.monthly)
  153.   path.monthly.df$Date <- index(path.monthly) # Add date column for indexing
  154.   path.monthly.df$MaxWeight <- maxWeights[i] # Fill in weight parameter 
  155.  
  156.   # Calculate the average path across each simulation for this weight parameter
  157.   path.monthly.df$Average <- apply(path.monthly.df[,1:numSims], 1, mean)
  158.  
  159.   # Return all values to the melted object
  160.   meltedPath <- rbind(meltedPath, melt(path.monthly.df, id.vars = c("Date","MaxWeight","Average"), variable.name = "Seed", value.name="Return"))
  161. }
  162. meltedPath <- meltedPath[-1,]
  163. dim(meltedPath)
  164.  
  165. # Chart Averages using GGPlot
  166. p <- ggplot(meltedPath, aes(x=Date, y=Average, group=MaxWeight))
  167. p + geom_line(aes(colour = MaxWeight)) + scale_colour_gradient(low="red") + scale_x_yearmon()
  168.  
  169. # Chart all results in ggplot2
  170. p <- ggplot(meltedPath, aes(x=Date, y=Return, group=Seed))
  171. p + geom_line(aes(colour = MaxWeight)) + scale_colour_gradient(low="red") + scale_x_yearmon()
  172.  
  173. # Melt for Performance Metrics
  174. perfMetrics <- NA # Pre-allocate
  175. for(i in 1:length(maxWeights)){
  176.   # Gather and re-index simulation results
  177.   pathrets <- data.frame(t(result2d[,i,]))
  178.   path <- as.zoo(apply(pathrets+1,2,cumprod)-1)
  179.   index(path) <- index(stockReturns)
  180.  
  181.   # Calculate annualized return, risk, and sharpe ratios
  182.   Return <- as.numeric(tail(path+1,1)^(252/nrow(path))-1)
  183.   Risk <- as.numeric(apply(pathrets, 2, StdDev.annualized, scale=252))
  184.   Sharpe <- as.numeric(Return/Risk)
  185.  
  186.   # Gather results into a data frame and return
  187.   metricsDF <- data.frame(Return=Return, Risk=Risk, Sharpe=Sharpe, MaxWeight=maxWeights[i])
  188.   perfMetrics <- rbind(perfMetrics, metricsDF)
  189. }
  190. perfMetrics <- perfMetrics[-1,]
  191.  
  192. # Analyze performance metrics
  193. p <- ggplot(perfMetrics, aes(x=Risk, y=Return, color=MaxWeight))
  194. p + geom_point(alpha=0.5)+ scale_colour_gradient(low="red")
  195.  
  196. p <- ggplot(perfMetrics, aes(x=Risk, y=Return, color=Sharpe))
  197. p + geom_point(alpha=0.5)+ scale_colour_gradient(low="red",high="green")
  198.  
  199. p <- ggplot(perfMetrics, aes(x=Risk, y=Sharpe, color=MaxWeight))
  200. p + geom_point(alpha=0.5)+ scale_colour_gradient(low="red")

Leave a Reply