Introduction

I am very excited to finally share some of my research exploring Meucci’s (Meucci (2005)) portfolio optimization methods, and how the resulting portfolios compare to the use of historical data. For those unfamiliar with Attilio Meucci, he runs an annual Advanced Risk and Portfolio Managment Bootcamp in New York City every summer. The bootcamp attracts academics and professionals within the industry, and over 6 intense days, topics and techniques in Risk Management and Portfolio Management are discussed in depth. I was fortunate to attend the bootcamp during the summer of 2014. In this post, I will be exploring portfolio optimization of the 9 Sector SPDR ETFs, using Meucci’s methods, specifically his framework outlined in “The Checklist” (Meucci (2011)), and I will be verifying if Meucci’s methods add value over the use of historical data. Lastly, I will be using Markowitz’s (Markowitz (1952)) quadratic utility framework for portfolio optimization, and will be comparing results using Minimum Variance and Mean Variance optimizations.

Background

Meucci’s core methods differ from the standard practice insofar as he prefers the use of Profit-and-Loss [P&L] over returns, and he prefers Monte Carlo Simulation over the direct use of historical data. Details of his framework can be found in “The Checklist” (Meucci (2011)), however I will provide a brief summary here:

1. For each instrument, one needs to identify the “invariants”, or the drivers of the changes in price that are “Identically and Independently Distributed” [IID]. For equities, this is usually log-returns.
2. Once the invariants are identified (e.g. log-returns for equities), this yields a matrix of invariants. To reduce estimation error, dimension reduction is performed by computing the eigenvalues and eigenvectors. Then, a subset of eigenvalues/eigenvectors are identified as the “signal”, whereas the remaining ones are classified as “noise”.
3. For the “signal” eigenvalues and eigenvectors, these are converted to their Principal Components [PCs] (i.e. the product of the invariants with the signal-eigenvectors). The distribution for each of the PCs is estimated, as well as the dependence between them (the copula). For the noise components, these are assumed to be distributed with mean 0, and standard deviation equal to the corresponding noise-eigenvalue, and are assumed to be independent/uncorrelated with each other.
4. Next, the Monte Carlo Simulation is conducted by drawing N trials, for P periods, from the meta-distribution estimated in the previous step. The idea is that the optimization should occur at some future investment horizon, and thus prices should be projected to where they likely will be.
5. Finally, these returns need to be converted to P&Ls, thus the distribution of N invariants at the investment horizon is multiplied by the latest security prices to obtain a distribution of P&Ls.
6. Lastly, the portfolio optimization is conducted on the P&Ls, instead of percent returns. Thus, the solution yields holdings (i.e. number of shares for each security) instead of the traditional asset weights.

In order to simplify this process, a number of assumptions were made for this study.

Assumptions

• Only equities are considered, whereas Meucci’s framework includes Fixed Income and Derivatives
• All signal-PCs are assumed to have a skew-student-t distribution
• All noise-PCs are assumed to either be normally distributed, or have student-t distribution, and are selected based on the Anderson-Darling goodness of fit test
• Given that there are 9 sector spdr ETFs, it is assumed that there are always 3 signal PCs, and 6 noise PCs
• The dependence between the 3 signal PCs is assumed to have a Vine Copula (Kramer and Schepsmeier (2011)) structure
• The dependence between the 6 noise PCs is assumed to be uncorrelated within the 6 PCs, as well as uncorrelated with the 3 signal PCs
• Transaction costs are ignored for rebalancing
• The optimization is long-only, fully-invested with a max-weight constraint on the 2 highest volatility ETFs
  • For each period, the historical volatility is computed, and for the 2 highest volatility ETFs, their weights are constrained by the inverse-volatility weight calculated for that period
  • All other “non-high-volatility” ETFs can be fully-invested (100% of wealth)
• For the Auto-Regressive Moving-Average [ARMA] trials, the ARMA model is fit if the Ljung-Box test indicates autocorrelation with 2.5% statistical significance
  • I am aware of criticisms of the Ljung-Box test for autocorrelation
• Weekly data is used to calculate the invariants, eigenvalues/eigenvector and to ultimately solve for holdings
  • The rationale is that this eliminates some of the noise from the daily prices, yet provides enough data to draw conclusions
• Performance is evaluated using daily data
• Covariance matrices are estimated using “robust” methodologies
  • robust::covRob() was used to estimate covariance matrices
• 3 Windows are selected to evaluate historical horizons
  • 18 months is assumed to be a short-horizon window
  • 36 months is assumed to be a medium-horizon window
  • 60 months is assumed to be a long-horizon window
• Portfolios are rebalanced monthly, thus the Meucci investment horizon is 4-weeks ahead
  • The invariants are projected 4 periods ahead, and then re-priced and optimized at this 4-period-ahead horizon

It should also be noted that, I have made many assumptions above. If anyone happened to read my first post on the Low-Volatilty Trading Strategy, then they may wonder why I did not take the hypothesis-driven approach outlined by Peterson (2014). Work on this post started in Q1 2015, months after I had attended Meucci’s bootcamp. Much of the work, assumptions, and results were obtained throughout 2015. I had taken Brian Peterson’s CFRM561 course in the spring/summer of 2016, after much of the work was completed. Due to the number of assumptions, all of which can arguably be converted to testable hypotheses, this would significantly delay the release of this information. I believe I have some valuable conclusions to share based on the work that has already been completed, and thus decided to write this post. Ideally though, all assumptions I have made above would be tested using the methods outlined in Peterson (2014), to ensure that the distributional and dependence assumptions are sound, and to reduce the risk of overfitting.

Result Replication

To replicate the results, two of my packages, titled propfolio and tsconv need to be installed. The following lines of code will install the packages, provided that the devtools library is already installed. Warning: propfolio has a lot of dependencies, so installing it will install many additional packages. Please review the description of the package prior to installing. Code to replicate the results can be found in the appendix.

devtools::install_github("erolbicero/tsconv")
devtools::install_github("erolbicero/propfolio")

ARMA Approach

Throughout 2015, at times during the analysis, I noticed that some of the sectors exhibited autocorrelation (e.g. XLU at times appears to have significant autocorrelation). During the “quest for invariance” step (Meucci (2011)), the result should yield no autocorrelation once the invariant has been identified. Given that I observed autocorrelation, this implied that there was information in the return series that was not being captured. In an effort to obtain uncorrelated residuals (invariants), I decided to verify if there is statistically significant autocorrelation using the Ljung-Box test at 2.5% significance. If the null hypothesis was rejected, then I used forecast::auto.arima() to fit an ARMA model. The only time this additional step is conducted is for the 18 month window. It is an extension of the Monte Carlo approach, and technically more accurate as it ensures that the invariants are IID. I did not do this for the 36 month and 60 month windows, since I did test this previously in 2015, and found that using forecast::auto.arima() did not make performance improvements; if anything it negatively impacted performance.

Optimization Criteria, Objectives and Choice of Benchmark

Overall Structure and Process

First, the Minimum Variance Optimization will be evaluated for each time window, where Meucci’s Monte Carlo approach will be compared to the historical method. The stronger strategy, i.e. Historical vs. Monte Carlo, will be identified for each window length.

Next, the same comparison (Historical vs. Monte Carlo) will be evaluated for the Optimal Sharpe criterion, within each window length.

Lastly, a summary of outperforming strategies relative to the benchmark will be identified, and then it can be assessed whether or not the Monte Carlo [MC] approach adds value over the benchmark and the historical method.

Historical vs. Monte Carlo Performance for Minimum Variance Optimization

## [1] "Paired t-test (mean difference) for 'Minimum Variance 60mo Hist' vs. 'Minimum Variance 60mo MC'"
## [1] "p-value = 0.94"
## [1] "F-test (ratio of variances) for 'Minimum Variance 60mo Hist' vs. 'Minimum Variance 60mo MC'"
## [1] "p-value = 0.92"
## [1] "Paired t-test (mean difference) for 'Minimum Variance 60mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0.82"
## [1] "F-test (ratio of variances) for 'Minimum Variance 60mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0"

## [1] "Paired t-test (mean difference) for 'Minimum Variance 36mo Hist' vs. 'Minimum Variance 36mo MC'"
## [1] "p-value = 0.68"
## [1] "F-test (ratio of variances) for 'Minimum Variance 36mo Hist' vs. 'Minimum Variance 36mo MC'"
## [1] "p-value = 0.67"
## [1] "Paired t-test (mean difference) for 'Minimum Variance 36mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0.41"
## [1] "F-test (ratio of variances) for 'Minimum Variance 36mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0"

## [1] "Paired t-test (mean difference) for 'Minimum Variance 18mo Hist' vs. 'Minimum Variance 18mo ARMA'"
## [1] "p-value = 0.85"
## [1] "F-test (ratio of variances) for 'Minimum Variance 18mo Hist' vs. 'Minimum Variance 18mo ARMA'"
## [1] "p-value = 0.8"
## [1] "Paired t-test (mean difference) for 'Minimum Variance 18mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0.79"
## [1] "F-test (ratio of variances) for 'Minimum Variance 18mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0"
## [1] "Paired t-test (mean difference) for 'Minimum Variance 18mo ARMA' vs. 'Benchmark'"
## [1] "p-value = 0.86"
## [1] "F-test (ratio of variances) for 'Minimum Variance 18mo ARMA' vs. 'Benchmark'"
## [1] "p-value = 0"

Historical vs. Monte Carlo Performance for Optimal Sharpe Ratio Optimization

## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 60mo Hist' vs. 'Optimal Sharpe 60mo MC'"
## [1] "p-value = 0.13"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 60mo Hist' vs. 'Optimal Sharpe 60mo MC'"
## [1] "p-value = 0.2"
## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 60mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0.87"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 60mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0"

## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 36mo Hist' vs. 'Optimal Sharpe 36mo MC'"
## [1] "p-value = 0.89"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 36mo Hist' vs. 'Optimal Sharpe 36mo MC'"
## [1] "p-value = 0.04"
## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 36mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0.21"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 36mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0"

## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 18mo Hist' vs. 'Optimal Sharpe 18mo ARMA'"
## [1] "p-value = 0.59"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 18mo Hist' vs. 'Optimal Sharpe 18mo ARMA'"
## [1] "p-value = 0.01"
## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 18mo Hist' vs. 'Optimal Sharpe 18mo MC'"
## [1] "p-value = 0.96"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 18mo Hist' vs. 'Optimal Sharpe 18mo MC'"
## [1] "p-value = 0"
## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 18mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0.98"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 18mo Hist' vs. 'Benchmark'"
## [1] "p-value = 0"
## [1] "Paired t-test (mean difference) for 'Optimal Sharpe 18mo ARMA' vs. 'Benchmark'"
## [1] "p-value = 0.67"
## [1] "F-test (ratio of variances) for 'Optimal Sharpe 18mo ARMA' vs. 'Benchmark'"
## [1] "p-value = 0"

Overall Comparison and Identifying the Optimal Strategy

Upon comparing the results above, the clear performers are:

• Historical Minimum Variance (60 Months)
• Monte Carlo Optimal Sharpe (18 Months)
• ARMA Optimal Sharpe (18 Months)

Which is summarized in the table below.

Where:

• H: Historical
• M: Monte Carlo
• A: ARMA
• OS: Optimal Sharpe
• MV: Minimum Variance

The 36 month window is excluded since it could not outperform the benchmark in all tests.

## [1] "Paired t-test (mean difference) for 'Minimum Variance 60mo Hist' vs. 'Optimal Sharpe 18mo ARMA'"
## [1] "p-value = 0.12"
## [1] "F-test (ratio of variances) for 'Minimum Variance 60mo Hist' vs. 'Optimal Sharpe 18mo ARMA'"
## [1] "p-value = 0"
## [1] "Paired t-test (mean difference) for 'Minimum Variance 60mo Hist' vs. 'Optimal Sharpe 18mo MC'"
## [1] "p-value = 0.46"
## [1] "F-test (ratio of variances) for 'Minimum Variance 60mo Hist' vs. 'Optimal Sharpe 18mo MC'"
## [1] "p-value = 0"

Final Thoughts and Overall Conclusion

There is one additional consideration that needs attention. I assumed no transaction costs, and given the shorter time horizon of the 18 month strategies, as well as the optimization method (Optimal Sharpe vs. Minimum Variance), when the weights are plotted for each of the strategies, it is clear that the 18 month ARMA and MC strategies are not really portfolio management strategies; they exhibit characteristics of tactical trading strategies due to the aggressively changing portfolio weights.

Thus, if the objective is for stable portfolio weights, then a long horizon window and minimum variance optimization appears to be a better choice under this particular study, due to the relatively stable weights. However, if there is appetite for complexity, and individual transaction costs are low, then there is potential to add more value using Attilio Meucci’s methods.

It may be possible to stabilize the weights under the 18 month scenarios, however it is out of scope for this study, and may be something valuable to investigate at a later date.

Rolling Weights for Minimum Variance Optimal Strategy

Rolling Weights for Optimal Sharpe Optimal Strategy

Rolling Weights for Optimal Sharpe Optimal Strategy

Opportunities for Future Research

Based on the work in this study, I think there is a lot of opportunity to explore the power of Meucci’s methods. I do think there is a lot of potential for horizons

• Start hypothesis testing from first-principles to reduce the number of assumptions (e.g. signal-PC distribution)
  • Add additional well-thought-out bias to assist with out-of-sample performance
  • Use Cross Validation to evaluate distribution and copula fitting
• Verify predictive ability of the ARMA process
• Add concentration constraints to Optimal Sharpe optimization to stabilize asset weights, especially in the 18-month window
• Identify different regimes (e.g. High vs. Low volatility) and apply different models/hypotheses/assumptions to each regime

Please contact me if you have any questions/comments/feedback or you have any ideas. I hope you found this interesting and useful.

© 2016 Erol Biceroglu

Appendix: Code to Replicate Results

Historical strategy code:

library(propfolio)
library(robust)
library(PerformanceAnalytics)
library(foreach)
library(doParallel)

nYears <- c(5,3,1.5) #set this to the 1st index for 60 months, 2nd index for 36 months or 3rd for 18 months
nYears <- nYears[1]
freq <- "weekly"
freqScale <- switch(freq, daily = 252, weekly = 52, monthly = 12, quarterly = 4, yearly = 1, stop("undefined"))

#########
#Start from scratch
symb <- c("XLE"
          ,"XLU"
          ,"XLK"
          ,"XLB"
          ,"XLP"
          ,"XLY"
          ,"XLI"
          ,"XLV"
          ,"XLF"
)

optTestPriceXtsMatrixALLweekly <- importEquityPricesToXtsMatrix(symbolVector = symb
                                , startDate = "1998-01-01"
                                , endDate = "2016-12-30"
                                , priceType = "Adjusted"
                                , freq = freq
                                , keep = FALSE
)

optTestPriceXtsMatrixALLdaily <- importEquityPricesToXtsMatrix(symbolVector = symb
                                                                , startDate = "1998-01-01"
                                                                , endDate = "2016-12-30"
                                                                , priceType = "Adjusted"
                                                                , freq = "daily"
                                                                , keep = FALSE
)

SPY <- importEquityPricesToXtsMatrix(symbolVector = "SPY"
                                                                , startDate = "1998-01-01"
                                                                , endDate = "2016-12-30"
                                                                , priceType = "Adjusted"
                                                                , freq = freq
                                                                , keep = FALSE
                                    )

SPYd <- importEquityPricesToXtsMatrix(symbolVector = "SPY"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = "daily"
                                     , keep = FALSE
)

RSP <- importEquityPricesToXtsMatrix(symbolVector = "RSP"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = freq
                                     , keep = FALSE
)

RSPd <- importEquityPricesToXtsMatrix(symbolVector = "RSP"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = "daily"
                                     , keep = FALSE
)

  cl <- makeForkCluster(detectCores())  

  registerDoParallel(cl)

  loopStart <- Sys.time()

 rollingOptimizationResults<-

foreach(dateIndex = toRowList(createRollingTimeIndices(index(optTestPriceXtsMatrixALLweekly), estimationLength = 52*nYears, timeUnit = "weeks", investmentHorizonLength = 4))
  , .errorhandling = "pass"
) %dopar%
{
  #allocate date indices
  OOSIndex <- as.vector(dateIndex[2])
  pseudoISPeriodPlot <- as.vector(dateIndex[3])
  dateIndex <- as.vector(dateIndex[1])

  #filter out time series object
  optTestPriceXtsMatrix <- optTestPriceXtsMatrixALLweekly[dateIndex]

  #extract "horizon" prices
  optTestLatestPrices <- as.vector(last(optTestPriceXtsMatrix))

  #compute linear returns
  optTestLRetXtsMatrix <- retMatrixL(optTestPriceXtsMatrix)

  #convert historical prices to change in prices (using horizon prices)
  optTestLRetXtsMatrixPnL <- coredata(do.call(merge,lapply(1:length(optTestLatestPrices)
                                                           ,function(x){optTestLRetXtsMatrix[,x] * optTestLatestPrices[x]}
  )))

  #compute moments
  covRobResults <- robust::covRob(data = optTestLRetXtsMatrixPnL)
  optTestRobMean <- as.vector(covRobResults$center)
  optTestRobCovPnL <- covRobResults$cov

  #compute historical inverse volatility weights (for constraints)
  sdWeightVec <- 1/apply(optTestLRetXtsMatrix,2,sd)/sum(1/apply(optTestLRetXtsMatrix,2,sd))

  #extract 7 lowest volatility names (highest 2 will be constrained by inverse volatility weight)
  constraintNames <- names(which(rank(sdWeightVec)>=7))

  sdWeightVecConstraint <- sdWeightVec[constraintNames]

  #compute efficient frontier
  LongOnly <- TRUE
  optTestOptimizedPortfolio <-
    optimizeMeanVariancePortfolio( meanVector = optTestRobMean
                                   , covarianceMatrix = optTestRobCovPnL
                                   , horizonPriceVector = optTestLatestPrices
                                   , horizonPnLMatrix = optTestLRetXtsMatrixPnL
                                   , wealthBalance = 10000
                                   , maxWeightConstraint = list(    names = constraintNames #safe weights
                                                                    , constraintWeight = sdWeightVec #all weights
                                                                    , maxWeight = ifelse(LongOnly,1,100)
                                   )
                                   , colNames = names(optTestPriceXtsMatrix)
                                   , longOnly = LongOnly
    )

  return(
    list(
      estimationPeriod = dateIndex
      , OOSPeriod = OOSIndex
      , pseudoISPeriodPlot = pseudoISPeriodPlot
      , optimalResult = optTestOptimizedPortfolio$optimalResult
      , minVarResult = optTestOptimizedPortfolio$minVarResult
      , expectedReturn = optTestOptimizedPortfolio$expectedReturn
      , standardDeviation = optTestOptimizedPortfolio$standardDeviation
    )

  )

}

 stopCluster(cl=cl)
 closeAllConnections()
 gc()
 cl <- NULL

 loopEnd <- Sys.time()
 loopEnd - loopStart
 # Time difference of 52.48464 mins on 8 cores

#create final result list object
rollingOptimizationResults <- list( portfolioNames = names(optTestPriceXtsMatrixALLweekly)
                  , portfolioFreq = "weeks"
                  , portfolioStats = rollingOptimizationResults
)

Monte Carlo Simulation strategy code

library(propfolio)
library(robust)
library(PerformanceAnalytics)
library(foreach)
library(doParallel)

nYears <- c(5,3,1.5) #set this to the 1st index for 60 months, 2nd index for 36 months or 3rd for 18 months
nYears <- nYears[1]
freq <- "weekly"
freqScale <- switch(freq, daily = 252, weekly = 52, monthly = 12, quarterly = 4, yearly = 1, stop("undefined"))

#########
#Start from scratch
symb <- c("XLE"
          ,"XLU"
          ,"XLK"
          ,"XLB"
          ,"XLP"
          ,"XLY"
          ,"XLI"
          ,"XLV"
          ,"XLF"
)

optTestPriceXtsMatrixALLweekly <- importEquityPricesToXtsMatrix(symbolVector = symb
                                                                , startDate = "1998-01-01"
                                                                , endDate = "2016-12-30"
                                                                , priceType = "Adjusted"
                                                                , freq = freq
                                                                , keep = FALSE
)

optTestPriceXtsMatrixALLdaily <- importEquityPricesToXtsMatrix(symbolVector = symb
                                                               , startDate = "1998-01-01"
                                                               , endDate = "2016-12-30"
                                                               , priceType = "Adjusted"
                                                               , freq = "daily"
                                                               , keep = FALSE
)

SPY <- importEquityPricesToXtsMatrix(symbolVector = "SPY"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = freq
                                     , keep = FALSE
)

SPYd <- importEquityPricesToXtsMatrix(symbolVector = "SPY"
                                      , startDate = "1998-01-01"
                                      , endDate = "2016-12-30"
                                      , priceType = "Adjusted"
                                      , freq = "daily"
                                      , keep = FALSE
)

RSP <- importEquityPricesToXtsMatrix(symbolVector = "RSP"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = freq
                                     , keep = FALSE
)

RSPd <- importEquityPricesToXtsMatrix(symbolVector = "RSP"
                                      , startDate = "1998-01-01"
                                      , endDate = "2016-12-30"
                                      , priceType = "Adjusted"
                                      , freq = "daily"
                                      , keep = FALSE
)

  cl <- makeForkCluster(detectCores())
  loopStart <- Sys.time()

   rollingOptimizationResults<-

parLapply( cl = cl
           ,  toRowList(createRollingTimeIndices(index(optTestPriceXtsMatrixALLweekly), estimationLength = 52*nYears, timeUnit = "weeks", investmentHorizonLength = 4))
  ,function(dateIndex)
{

    #allocate date indices
OOSIndex <- as.vector(dateIndex[2])
pseudoISPeriodPlot <- as.vector(dateIndex[3])
dateIndex <- as.vector(dateIndex[1])

#filter out time series object
optTestPriceXtsMatrix <- optTestPriceXtsMatrixALLweekly[dateIndex]

#calculate log returns
optTestCRetXtsMatrix <- retMatrixC(optTestPriceXtsMatrix)

#eigen/dimension reduction
optTestRobCov <- robust::covRob(coredata(optTestCRetXtsMatrix))$cov

optTesteig <- eigenDimensionReduction(optTestRobCov)

#Identify dominant eigenvectors
optTestSigNoise <- dominantEigenvectors(optTesteig)

#fitting signal and noise marginals
optTestSigMarginalFits <- fitMarginalSignal(optTestCRetXtsMatrix, optTestSigNoise$signal$eigenValues, optTestSigNoise$signal$eigenVectors)

optTestNoiseMarginalFits <- fitMarginalResidual(optTestCRetXtsMatrix, optTestSigNoise$noise$eigenValues, optTestSigNoise$noise$eigenVectors)

#fit copula
optTestCopulaFit <- fitCopula( coredata(tsEigenProduct(optTestCRetXtsMatrix, optTestSigNoise$signal$eigenVectors)))

#simulate data to horizon
optTestSimulatedInvariantsToHorizon <- simulateDataToHorizon(     numSimulations = 52*200
                         , numPeriodsForward = 4 #22
                         , vineCopulaFit = optTestCopulaFit
                         , signalMarginalFit = optTestSigMarginalFits$fits
                         , residualMarginalFit = optTestNoiseMarginalFits$fits
                         , eigenVectors = optTesteig$eigenVectors #all eigenvectors
                      )
#dim(optTestSimulatedInvariantsToHorizon)

#pricing
optTestHorizonPriceVector <- as.vector(coredata(optTestPriceXtsMatrix[nrow(optTestPriceXtsMatrix),]))

optTestHorizonPrices <- horizonPrices(     simulatedHorizonInvariants = optTestSimulatedInvariantsToHorizon
                 , horizonPriceVector = optTestHorizonPriceVector
)

#extract horizon profit-and-loss
optTestHorizonPnLs <- horizonPnLs( simulatedHorizonPrices = optTestHorizonPrices
               , horizonPriceVector = optTestHorizonPriceVector
  )

#compute moments
covRobResults <- robust::covRob(data = optTestHorizonPnLs)
optTestHorizonRobMean <- as.vector(covRobResults$center)
optTestHorizonRobCovPnL <- covRobResults$cov

#compute historical inverse volatility weights (for constraints)
sdWeightVec <- 1/apply(optTestCRetXtsMatrix,2,sd)/sum(1/apply(optTestCRetXtsMatrix,2,sd))

#extract 7 lowest volatility names (highest 2 will be constrained by inverse volatility weight)
constraintNames <- names(which(rank(sdWeightVec)>=7))

sdWeightVecConstraint <- sdWeightVec[constraintNames]

#compute efficient frontier
LongOnly <- TRUE
optTestOptimizedPortfolio <-
    optimizeMeanVariancePortfolio( meanVector = optTestHorizonRobMean
                                 , covarianceMatrix = optTestHorizonRobCovPnL
                                 , horizonPriceVector = optTestHorizonPriceVector
                                 , horizonPnLMatrix = optTestHorizonPnLs
                                 , wealthBalance = 10000
                                 , maxWeightConstraint = list(    names = constraintNames #safe weights
                                                                , constraintWeight = sdWeightVec #all weights
                                                                , maxWeight = ifelse(LongOnly,1,100)
                                                              )
                                 , colNames = names(optTestPriceXtsMatrix)
                                 , longOnly = LongOnly
                                )

return(
      list(
           estimationPeriod = dateIndex
         , OOSPeriod = OOSIndex
         , pseudoISPeriodPlot = pseudoISPeriodPlot
         , optimalResult = optTestOptimizedPortfolio$optimalResult
         , minVarResult = optTestOptimizedPortfolio$minVarResult
         , expectedReturn = optTestOptimizedPortfolio$expectedReturn
         , standardDeviation = optTestOptimizedPortfolio$standardDeviation
         , fitData = list(signalMarginals = optTestSigMarginalFits$evaluation
                        , noiseMarginals = optTestNoiseMarginalFits$evaluation
                        , copula = -999
                          )

       )

  )

})

stopCluster(cl=cl)
closeAllConnections()
gc()
cl <- NULL

loopEnd <- Sys.time()
loopEnd - loopStart
# Time difference of 11.50882 mins  #aws 64 cores

#create final result list object
rollingOptimizationResults <- list( portfolioNames = names(optTestPriceXtsMatrixALLweekly)
                  , portfolioFreq = "weeks"
                  , portfolioStats = rollingOptimizationResults
)

Monte Carlo with ARMA strategy code:

library(propfolio)
library(robust)
library(PerformanceAnalytics)
library(foreach)
library(doParallel)

nYears <- c(5,3,1.5) #set this to the 1st index for 60 months, 2nd index for 36 months or 3rd for 18 months
nYears <- nYears[1]
freq <- "weekly"
freqScale <- switch(freq, daily = 252, weekly = 52, monthly = 12, quarterly = 4, yearly = 1, stop("undefined"))

#########
#Start from scratch
symb <- c("XLE"
          ,"XLU"
          ,"XLK"
          ,"XLB"
          ,"XLP"
          ,"XLY"
          ,"XLI"
          ,"XLV"
          ,"XLF"
)

optTestPriceXtsMatrixALLweekly <- importEquityPricesToXtsMatrix(symbolVector = symb
                                                                , startDate = "1998-01-01"
                                                                , endDate = "2016-12-30"
                                                                , priceType = "Adjusted"
                                                                , freq = freq
                                                                , keep = FALSE
)

optTestPriceXtsMatrixALLdaily <- importEquityPricesToXtsMatrix(symbolVector = symb
                                                               , startDate = "1998-01-01"
                                                               , endDate = "2016-12-30"
                                                               , priceType = "Adjusted"
                                                               , freq = "daily"
                                                               , keep = FALSE
)

SPY <- importEquityPricesToXtsMatrix(symbolVector = "SPY"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = freq
                                     , keep = FALSE
)

SPYd <- importEquityPricesToXtsMatrix(symbolVector = "SPY"
                                      , startDate = "1998-01-01"
                                      , endDate = "2016-12-30"
                                      , priceType = "Adjusted"
                                      , freq = "daily"
                                      , keep = FALSE
)

RSP <- importEquityPricesToXtsMatrix(symbolVector = "RSP"
                                     , startDate = "1998-01-01"
                                     , endDate = "2016-12-30"
                                     , priceType = "Adjusted"
                                     , freq = freq
                                     , keep = FALSE
)

RSPd <- importEquityPricesToXtsMatrix(symbolVector = "RSP"
                                      , startDate = "1998-01-01"
                                      , endDate = "2016-12-30"
                                      , priceType = "Adjusted"
                                      , freq = "daily"
                                      , keep = FALSE
)

  cl <- makeForkCluster(detectCores())  

  registerDoParallel(cl)

  loopStart <- Sys.time()

 rollingOptimizationResults<-

foreach(dateIndex = toRowList(createRollingTimeIndices(index(optTestPriceXtsMatrixALLweekly), estimationLength = 52*nYears, timeUnit = "weeks", investmentHorizonLength = 4))
  , .errorhandling = "pass"
) %dopar%
{

    #allocate date indices
OOSIndex <- as.vector(dateIndex[2])
pseudoISPeriodPlot <- as.vector(dateIndex[3])
dateIndex <- as.vector(dateIndex[1])

#filter out time series object
optTestPriceXtsMatrix <- optTestPriceXtsMatrixALLweekly[dateIndex]

#calculate log returns
optTestCRetXtsMatrix <- retMatrixC(optTestPriceXtsMatrix)

#verify autocorrelation
corrTheshold <- 0.025
#stats::Box.test(x = , type = "Ljung-Box")
#check if correlation is above threshold, for lag = 1
corrVector<-
sapply(optTestCRetXtsMatrix
       , function(x){

         LBtest <- stats::Box.test(x = x, type = "Ljung-Box")
         LBpvalue <- LBtest$p.value
         return(LBpvalue)
       }
      )

corrIDX <- which(corrVector <= corrTheshold)

#if autocorrelation, fit ARMA model
arimaModels <-
  lapply(corrIDX
         , function(x){
           vectorToFit <- optTestCRetXtsMatrix[,x]
           arimaFit <- forecast::auto.arima(x = vectorToFit)
           return(arimaFit)
         }
  )

residualData <- do.call(cbind,lapply(arimaModels
                       , function(x){
                         x$residuals
                       }
                       ))

combinedInvariants <- NULL

for(idx in 1:ncol(optTestCRetXtsMatrix)){
  if(idx %in% corrIDX){
  combinedInvariants <- cbind(combinedInvariants
                              , if(length(corrIDX)==1){residualData} else {residualData[,which(idx == corrIDX)]}
                              )
  } else {
    combinedInvariants <- cbind(combinedInvariants, coredata(optTestCRetXtsMatrix[,idx]))
  }

}

colnames(combinedInvariants) <- symb

#eigen/dimension reduction
optTestRobCov <- robust::covRob(coredata(combinedInvariants))$cov

optTesteig <- eigenDimensionReduction(optTestRobCov)

#Identify dominant eigenvectors
optTestSigNoise <- dominantEigenvectors(optTesteig)

#fitting signal and noise marginals
optTestSigMarginalFits <- fitMarginalSignal(optTestCRetXtsMatrix, optTestSigNoise$signal$eigenValues, optTestSigNoise$signal$eigenVectors)

optTestNoiseMarginalFits <- fitMarginalResidual(optTestCRetXtsMatrix, optTestSigNoise$noise$eigenValues, optTestSigNoise$noise$eigenVectors)

#fit copula
optTestCopulaFit <- fitCopula( coredata(tsEigenProduct(optTestCRetXtsMatrix, optTestSigNoise$signal$eigenVectors)))

#simulate data to horizon
optTestSimulatedInvariantsToHorizon <- simulateDataToHorizon3(     numSimulations = 52*200
                         , numPeriodsForward = 4 #22
                         , vineCopulaFit = optTestCopulaFit
                         , signalMarginalFit = optTestSigMarginalFits$fits
                         , residualMarginalFit = optTestNoiseMarginalFits$fits
                         , eigenVectors = optTesteig$eigenVectors #all eigenvectors
                         , corrIDX = corrIDX
                         , arimaModel = arimaModels
                      )

#pricing
optTestHorizonPriceVector <- as.vector(coredata(optTestPriceXtsMatrix[nrow(optTestPriceXtsMatrix),]))

optTestHorizonPrices <- horizonPrices(     simulatedHorizonInvariants = optTestSimulatedInvariantsToHorizon
                 , horizonPriceVector = optTestHorizonPriceVector
)

#extract horizon profit-and-loss
optTestHorizonPnLs <- horizonPnLs( simulatedHorizonPrices = optTestHorizonPrices
               , horizonPriceVector = optTestHorizonPriceVector
  )

#compute moments
covRobResults <- robust::covRob(data = optTestHorizonPnLs)
optTestHorizonRobMean <- as.vector(covRobResults$center)
optTestHorizonRobCovPnL <- covRobResults$cov

#compute historical inverse volatility weights (for constraints)
sdWeightVec <- 1/apply(optTestCRetXtsMatrix,2,sd)/sum(1/apply(optTestCRetXtsMatrix,2,sd))

#extract 7 lowest volatility names (highest 2 will be constrained by inverse volatility weight)
constraintNames <- names(which(rank(sdWeightVec)>=7))

sdWeightVecConstraint <- sdWeightVec[constraintNames]

#compute efficient frontier
LongOnly <- TRUE
optTestOptimizedPortfolio <-
    optimizeMeanVariancePortfolio( meanVector = optTestHorizonRobMean
                                 , covarianceMatrix = optTestHorizonRobCovPnL
                                 , horizonPriceVector = optTestHorizonPriceVector
                                 , horizonPnLMatrix = optTestHorizonPnLs
                                 , wealthBalance = 10000
                                 , maxWeightConstraint = list(    names = constraintNames #safe weights
                                                                , constraintWeight = sdWeightVec #all weights
                                                                , maxWeight = ifelse(LongOnly,1,100)
                                                              )
                                 , colNames = names(optTestPriceXtsMatrix)
                                 , longOnly = LongOnly
                                )

return(
      list(
           estimationPeriod = dateIndex
         , OOSPeriod = OOSIndex
         , pseudoISPeriodPlot = pseudoISPeriodPlot
         , optimalResult = optTestOptimizedPortfolio$optimalResult
         , minVarResult = optTestOptimizedPortfolio$minVarResult
         , expectedReturn = optTestOptimizedPortfolio$expectedReturn
         , standardDeviation = optTestOptimizedPortfolio$standardDeviation
         , fitData = list(signalMarginals = optTestSigMarginalFits$evaluation
                        , noiseMarginals = optTestNoiseMarginalFits$evaluation
                        , copula = -999
                          )

       )

  )

}

stopCluster(cl=cl)
closeAllConnections()
gc()
cl <- NULL

loopEnd <- Sys.time()
loopEnd - loopStart

# Time difference of 1.595575 hours
# Time difference of 13.83757 mins #ON AWS 64 cores

rollingOptimizationResults <- list( portfolioNames = names(optTestPriceXtsMatrixALLweekly)
                  , portfolioFreq = "weeks"
                  , portfolioStats = rollingOptimizationResults
)

Calculating results after running one of 3 scenarios above

#Optimal Sharpe
realizedOPPortfolioMatrix<-
  Return.portfolio(R = retMatrixL(optTestPriceXtsMatrixALLdaily)
                   , weights = createOptimalStatisticTimeSeries(optimizationResultsObject =  rollingOptimizationResults
                                                                , statisticType =  "opWeights"
                                                                , datePeriodType = "test"
                                                                , indexTimeSeriesFreq = "monthly")
  )
colnames(realizedOPPortfolioMatrix) <- "Optimal Sharpe"
realizedOPPortfolioMatrix <- merge(realizedOPPortfolioMatrix, retVectorL(RSPd), retMatrixL(optTestPriceXtsMatrixALLdaily))
colnames(realizedOPPortfolioMatrix)[2] <- "Benchmark"

#Minimum Variance
realizedMVPortfolioMatrix<-
  Return.portfolio(R = retMatrixL(optTestPriceXtsMatrixALLdaily)
                   , weights = createOptimalStatisticTimeSeries(optimizationResultsObject =  rollingOptimizationResults
                                                                , statisticType =  "mvWeights"
                                                                , datePeriodType = "test"
                                                                , indexTimeSeriesFreq = "monthly")
  )
colnames(realizedMVPortfolioMatrix) <- "Minimum Variance"
realizedMVPortfolioMatrix <- merge(realizedMVPortfolioMatrix, retVectorL(RSPd), retMatrixL(optTestPriceXtsMatrixALLdaily))
colnames(realizedMVPortfolioMatrix)[2] <- "Benchmark"

References