We will load some stock market data and divide it into a training part (for the estimation of the expected return \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\), and subsequent portfolio design) and a test part (for the out-of-sample performance evaluation).

In particular, we will start by loading some stock data from three different sectors:

• sector technology: AAPL, AMD, ADI.
• sector healthcare: ABBV, AEZS, A.
• sector basic materials: APD, AA, CF.

library(xts)  # to manipulate time series of stock data
library(quantmod)  # to download stock data
library(PerformanceAnalytics)  # to compute performance measures

# download data from YahooFinance
stock_namelist <- c("AAPL", "AMD", "ADI",  "ABBV", "AEZS", "A",  "APD", "AA","CF")
prices <- xts()
for (i in 1:length(stock_namelist)) {
  tmp <- Ad(getSymbols(stock_namelist[i], from = "2013-01-01", to = "2016-12-31", auto.assign = FALSE))
  tmp <- na.approx(tmp, na.rm = FALSE)  # interpolate NAs
  prices <- cbind(prices, tmp)
}
colnames(prices) <- stock_namelist
indexClass(prices) <- "Date"

str(prices)

R>> An 'xts' object on 2013-01-02/2016-12-30 containing:
R>>   Data: num [1:1008, 1:9] 55.5 54.8 53.2 52.9 53.1 ...
R>>  - attr(*, "dimnames")=List of 2
R>>   ..$ : NULL
R>>   ..$ : chr [1:9] "AAPL" "AMD" "ADI" "ABBV" ...
R>>   Indexed by objects of class: [Date] TZ: UTC
R>>   xts Attributes:  
R>>  NULL

head(prices)

R>>                AAPL  AMD      ADI     ABBV AEZS        A      APD       AA       CF
R>> 2013-01-02 55.47129 2.53 37.50844 27.60817  253 27.97831 66.43760 20.62187 29.44429
R>> 2013-01-03 54.77112 2.49 36.90318 27.38020  254 28.07852 66.20545 20.80537 29.30506
R>> 2013-01-04 53.24548 2.59 36.24679 27.03431  257 28.63300 67.09529 21.24121 29.96146
R>> 2013-01-07 52.93227 2.67 36.35761 27.08934  259 28.42591 67.03341 20.87419 29.84922
R>> 2013-01-08 53.07474 2.67 35.98254 26.49976  255 28.19877 67.15720 20.87419 29.41162
R>> 2013-01-09 52.24524 2.63 35.88875 26.64913  258 28.96036 68.06254 20.82831 30.44025

tail(prices)

R>>                AAPL   AMD      ADI     ABBV AEZS        A      APD    AA       CF
R>> 2016-12-22 111.8445 11.60 70.07727 55.82264 4.10 44.91877 136.0360 29.75 27.82244
R>> 2016-12-23 112.0657 11.58 70.44936 56.43826 4.10 45.14350 136.3992 29.71 28.37686
R>> 2016-12-27 112.7774 12.07 70.89778 56.58311 4.05 45.44642 137.2373 29.65 29.52265
R>> 2016-12-28 112.2965 11.55 70.18222 56.37488 3.55 44.67448 135.1328 29.43 29.22696
R>> 2016-12-29 112.2676 11.59 70.20129 56.79134 3.60 44.72544 135.2359 28.89 29.47645
R>> 2016-12-30 111.3924 11.34 69.28539 56.69175 3.60 44.64704 134.7206 28.08 29.08836

# compute log-returns and linear returns
X_log <- diff(log(prices))[-1]
X_lin <- (prices/lag(prices) - 1)[-1]

# or alternatively...
X_log <- CalculateReturns(prices, "log")[-1]
X_lin <- CalculateReturns(prices)[-1]

N <- ncol(X_log)  # number of stocks
T <- nrow(X_log)  # number of days

We can take a look at the prices of the stocks:

plot(prices/rep(prices[1, ], each = nrow(prices)), col = rainbow10equal, legend.loc = "topleft",
     main = "Normalized prices")

We now divide the data into a training set and test set:

T_trn <- round(0.7*T)  # 70% of data
X_log_trn <- X_log[1:T_trn, ]
X_log_tst <- X_log[(T_trn+1):T, ]
X_lin_trn <- X_lin[1:T_trn, ]
X_lin_tst <- X_lin[(T_trn+1):T, ]

We can now use the training set to obtain the sample estimates from the returns \(\mathbf{x}_t\) (i.e., sample mean and sample covariance matrix) as \[ \begin{aligned} \hat{\boldsymbol{\mu}} & = \frac{1}{T}\sum_{t=1}^T \mathbf{x}_t\\ \hat{\boldsymbol{\Sigma}} & = \frac{1}{T-1}\sum_{t=1}^T (\mathbf{x}_t - \hat{\boldsymbol{\mu}})(\mathbf{x}_t - \hat{\boldsymbol{\mu}})^T \end{aligned} \]

However, it is not totally clear whether we should use linear returns or log returns to estimate \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\). Clearly, for the portfolio design we need the expected return \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\) of the linear returns. There are three different main philosophies in the estimation procedure for \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\):

1. estimate them directly from the linear returns (even though linear returns are not supposed to be easily modeled): \(\hat{\boldsymbol{\mu}} = \hat{\boldsymbol{\mu}}^\textsf{lin}\) and \(\hat{\boldsymbol{\Sigma}} = \hat{\boldsymbol{\Sigma}}^\textsf{lin}\)
2. estimate them from the log returns (and ignore the approximation error): \(\hat{\boldsymbol{\mu}} = \hat{\boldsymbol{\mu}}^\textsf{log}\) and \(\hat{\boldsymbol{\Sigma}} = \hat{\boldsymbol{\Sigma}}^\textsf{log}\)
3. estimate them from the log returns but properly transforming them to linear: \[ \begin{align} \hat{\boldsymbol{\mu}} & = \exp\left( \hat{\boldsymbol{\mu}}^\textsf{log} + \frac{1}{2}\textsf{diag}(\hat{\boldsymbol{\Sigma}}^\textsf{log})\right) - 1\\ \hat{\Sigma}_{ij} & = \exp\left( \hat{\mu}_i^\textsf{log} + \hat{\mu}_j^\textsf{log} + \frac{1}{2}(\hat{\Sigma}_{ii}^\textsf{log} + \hat{\Sigma}_{jj}^\textsf{log}) \right) \times (\exp(\hat{\Sigma}_{ij}^\textsf{log}) - 1) \end{align} \]

In the following, we will explore these three different ways for daily returns. For the comparison we will use the simple Global Minimum Variance Portfolio (GMVP), \[ \mathbf{w}_\textsf{GMVP} = \frac{1}{\mathbf{1}^T\boldsymbol{\Sigma}^{-1}\mathbf{1}}\boldsymbol{\Sigma}^{-1}\mathbf{1}, \] as well as the mean-variance portfolio (MVP) as the solution to the problem: \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \mathbf{w}^{T}\boldsymbol{\mu}-\lambda\mathbf{w}^{T}\boldsymbol{\Sigma}\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^T\mathbf{w}=1. \end{array}\]

# Method 1: estimate directly from linear returns
mu_lin <- colMeans(X_lin_trn)
Xc <- X_lin_trn - matrix(mu_lin, T_trn, N, byrow = TRUE)  # remove mean (center data)
Sigma_lin <- 1/(T_trn-1) * t(Xc) %*% Xc
# or more conveniently:
Sigma_lin_ <- cov(X_lin_trn)
norm(Sigma_lin - Sigma_lin_, "F")  # sanity check

R>> [1] 4.019848e-18

# Method 2: estimate directly from log returns
mu_log <- colMeans(X_log_trn)
Sigma_log <- cov(X_log_trn)

# Method 3: estimate from log returns plus transformation
momentsReturnLog2Lin <- function(mu, Sigma) {
  N <- ncol(Sigma)
  mu_ <- exp(mu + 0.5*diag(Sigma)) - 1
  Sigma_ <- matrix(NA, nrow = N, ncol = N)
  for(ii in 1:N)
    for(jj in 1:N)
      Sigma_[ii, jj] <- exp(mu[ii] + mu[jj] + 0.5*(Sigma[ii, ii]+Sigma[jj, jj])) * (exp(Sigma[ii, jj])-1)
  return( list(mu=mu_, Sigma=Sigma_) )
}

tmp <- momentsReturnLog2Lin(mu_log, Sigma_log)
mu_log_trans <- tmp$mu
Sigma_log_trans <- tmp$Sigma

Now let’s compute and compare the three GMVP portfolios:

# create function for GMVP
GMVP <- function(Sigma) {
  ones <- rep(1, nrow(Sigma))
  Sigma_inv_1 <- solve(Sigma, ones)  #same as: inv(Sigma) %*% ones
  w <- (1/as.numeric(ones %*% Sigma_inv_1)) * Sigma_inv_1
  return(w)
}

# compute the three versions of GMVP
w_GMVP_lin <- GMVP(Sigma_lin)
w_GMVP_log <- GMVP(Sigma_log)
w_GMVP_log_trans <- GMVP(Sigma_log_trans)
w_GMVP_all <- cbind("GMVP_lin"       = w_GMVP_lin, 
                    "GMVP_log"       = w_GMVP_log, 
                    "GMVP_log_trans" = w_GMVP_log_trans)
w_GMVP_all

R>>        GMVP_lin   GMVP_log GMVP_log_trans
R>> AAPL 0.18865586 0.18725916     0.18703112
R>> AMD  0.01608760 0.01571003     0.01575742
R>> ADI  0.08229133 0.08404712     0.08407072
R>> ABBV 0.13928582 0.13836379     0.13818423
R>> AEZS 0.01377123 0.01452367     0.01463238
R>> A    0.14917832 0.14789974     0.14821945
R>> APD  0.27735998 0.28026594     0.28007130
R>> AA   0.03009397 0.03032780     0.03047517
R>> CF   0.10327589 0.10160275     0.10155822

# plot to compare the allocations
barplot(t(w_GMVP_all), col = rainbow8equal[1:3], legend = colnames(w_GMVP_all), 
        main = "GMVP allocation", xlab = "stocks", ylab = "dollars", beside = TRUE)

And now the three MVP portfolios:

library(CVXR)

# create function for MVP
MVP <- function(mu, Sigma, lmd = 0.5) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Maximize(t(mu) %*% w - lmd*quad_form(w, Sigma)),
                  constraints = list(sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# compute the three versions of MVP
w_MVP_lin <- MVP(mu_lin, Sigma_lin)
w_MVP_log <- MVP(mu_log, Sigma_log)
w_MVP_log_trans <- MVP(mu_log_trans, Sigma_log_trans)
w_MVP_all <- cbind("MVP_lin"       = w_MVP_lin, 
                   "MVP_log"       = w_MVP_log, 
                   "MVP_log_trans" = w_MVP_log_trans)
w_MVP_all

R>>           MVP_lin     MVP_log MVP_log_trans
R>>  [1,]  1.82349683  1.78880418     1.8028310
R>>  [2,] -0.44019876 -0.80020154    -0.4097064
R>>  [3,] -0.08553791  0.03403578    -0.0881237
R>>  [4,]  1.81463970  1.96492811     1.8217410
R>>  [5,] -0.64990649 -1.30141524    -0.8394611
R>>  [6,] -2.53548501 -2.23254728    -2.5396476
R>>  [7,]  1.86916151  2.40358620     1.9655884
R>>  [8,] -1.19378174 -1.20062371    -1.1803529
R>>  [9,]  0.39761187  0.34343349     0.4671312

# plot to compare the allocations
barplot(t(w_MVP_all), col = rainbow8equal[1:3], 
        legend = colnames(w_MVP_all), args.legend = list(x = "bottomright", inset = 0.04), 
        main = "GMVP allocation", xlab = "stocks", ylab = "dollars", beside = TRUE)

Finally, we can do an in-sample evaluation (no need to use the test sample for this):

# compute returns of the portfolios
ret_GMVP_all_trn <- xts(X_lin_trn %*% w_GMVP_all, index(X_lin_trn))
ret_MVP_all_trn <- xts(X_lin_trn %*% w_MVP_all, index(X_lin_trn))
# compare them
StdDev.annualized(ret_GMVP_all_trn)

R>>                                GMVP_lin  GMVP_log GMVP_log_trans
R>> Annualized Standard Deviation 0.1575621 0.1575663      0.1575672

SharpeRatio.annualized(ret_GMVP_all_trn)

R>>                                 GMVP_lin GMVP_log GMVP_log_trans
R>> Annualized Sharpe Ratio (Rf=0%) 1.193824 1.188617       1.187221

SharpeRatio.annualized(ret_MVP_all_trn)

R>>                                   MVP_lin MVP_log MVP_log_trans
R>> Annualized Sharpe Ratio (Rf=0%) 0.5270622     NaN           NaN

For the GMVP, which only uses the covariance matrix, clearly there is no apparent difference. This is not unexpected since the daily returns of the stocks are roughly within the interval [-0.04, 0.04] and the approximation \(r_t = \log(1 + R_t) \approx R_t\) is very accurate:

apply(X_lin_trn, 2, sd)

R>>       AAPL        AMD        ADI       ABBV       AEZS          A        APD         AA         CF 
R>> 0.01661471 0.03105342 0.01523081 0.01617068 0.07467274 0.01398818 0.01296057 0.01918079 0.01787465

However, for the mean-variance portfolio, which also uses the expected return, there is some small difference and the methods can be easily ranked:

1. the best method is to properly model the log-returns and then transform the moments to the linear return case
2. the second best method is to use directly the linear returns to estimate the moments
3. the worst method is to use the log-returns to estimate the moment and then use then as if they were the linear return moments.

Now, another question is whether these differences will become more evident for weekly or monthly rebalancing where the approximation \(r_t = \log(1 + R_t) \approx R_t\) may not be so accurate.

Let’s repeat the previous analysis but using weekly returns:

periodicity(prices)

R>> Daily periodicity from 2013-01-02 to 2016-12-30

# change periodicity of prices
prices_weekly <- xts()
for (i in 1:ncol(prices))
  prices_weekly <- cbind(prices_weekly, Cl(to.weekly(prices[, i])))
colnames(prices_weekly) <- colnames(prices)
indexClass(prices_weekly) <- "Date"

periodicity(prices_weekly)

R>> Weekly periodicity from 2013-01-04 to 2016-12-30

head(prices_weekly)

R>>                AAPL  AMD      ADI     ABBV AEZS        A      APD       AA       CF
R>> 2013-01-04 53.24548 2.59 36.24679 27.03431  257 28.63300 67.09529 21.24121 29.96146
R>> 2013-01-11 52.56855 2.67 36.02515 26.92660  274 29.01381 68.37980 20.50717 30.99721
R>> 2013-01-18 50.51755 2.46 36.70712 29.68687  259 29.49481 68.62740 20.64481 31.12366
R>> 2013-01-25 44.44331 2.85 37.33794 29.90961  259 30.25639 68.37205 20.71362 32.71777
R>> 2013-02-01 45.83153 2.60 38.07959 29.57551  289 30.25639 68.65834 20.64481 32.67231
R>> 2013-02-08 50.01838 2.59 38.95763 28.83573  302 30.11611 68.36431 20.57622 32.13240

# recompute returns
X_weekly_log <- CalculateReturns(prices_weekly, "log")[-1]
X_weekly_lin <- CalculateReturns(prices_weekly)[-1]
T_weekly <- nrow(X_weekly_log)  # number of weeks

# split data into training and set data
T_weekly_trn <- round(0.7*T_weekly)  # 70% of data
X_weekly_log_trn <- X_weekly_log[1:T_weekly_trn, ]
X_weekly_log_tst <- X_weekly_log[(T_weekly_trn+1):T_weekly, ]
X_weekly_lin_trn <- X_weekly_lin[1:T_weekly_trn, ]
X_weekly_lin_tst <- X_weekly_lin[(T_weekly_trn+1):T_weekly, ]

# estimate mu and Sigma
mu_weekly_lin <- colMeans(X_weekly_lin_trn)
Sigma_weekly_lin <- cov(X_weekly_lin_trn)

mu_weekly_log <- colMeans(X_weekly_log_trn)
Sigma_weekly_log <- cov(X_weekly_log_trn)

tmp <- momentsReturnLog2Lin(mu_weekly_log, Sigma_weekly_log)
mu_weekly_log_trans <- tmp$mu
Sigma_weekly_log_trans <- tmp$Sigma

# compute GMVPs
w_GMVP_weekly_lin <- GMVP(Sigma_weekly_lin)
w_GMVP_weekly_log <- GMVP(Sigma_weekly_log)
w_GMVP_weekly_log_trans <- GMVP(Sigma_weekly_log_trans)
w_GMVP_weekly_all <- cbind("GMVP_lin"       = w_GMVP_weekly_lin, 
                           "GMVP_log"       = w_GMVP_weekly_log, 
                           "GMVP_log_trans" = w_GMVP_weekly_log_trans)
barplot(t(w_GMVP_weekly_all), col = rainbow8equal[1:3], legend = colnames(w_GMVP_weekly_all), 
        main = "GMVP allocation for weekly returns", xlab = "stocks", ylab = "dollars", beside = TRUE)

# compute MVPs
w_MVP_weekly_lin <- MVP(mu_weekly_lin, Sigma_weekly_lin)
w_MVP_weekly_log <- MVP(mu_weekly_log, Sigma_weekly_log)
w_MVP_weekly_log_trans <- MVP(mu_weekly_log_trans, Sigma_weekly_log_trans)
w_MVP_weekly_all <- cbind("MVP_lin"       = w_MVP_weekly_lin, 
                          "MVP_log"       = w_MVP_weekly_log, 
                          "MVP_log_trans" = w_MVP_weekly_log_trans)
barplot(t(w_GMVP_weekly_all), col = rainbow8equal[1:3], legend = colnames(w_GMVP_weekly_all), 
        main = "MVP allocation for weekly returns", xlab = "stocks", ylab = "dollars", beside = TRUE)

The allocations seem to be very similar.

Let’s look at the portfolios performance:

# compute weekly returns of portfolios
ret_GMVP_weekly_all_trn <- xts(X_weekly_lin_trn %*% w_GMVP_weekly_all, index(X_weekly_lin_trn))
ret_MVP_weekly_all_trn <- xts(X_weekly_lin_trn %*% w_MVP_weekly_all, index(X_weekly_lin_trn))

# compare them
StdDev.annualized(ret_GMVP_weekly_all_trn)

R>>                                GMVP_lin  GMVP_log GMVP_log_trans
R>> Annualized Standard Deviation 0.1495471 0.1495629      0.1495584

SharpeRatio(ret_MVP_weekly_all_trn, FUN = "StdDev")

R>>                                 MVP_lin   MVP_log MVP_log_trans
R>> StdDev Sharpe (Rf=0%, p=95%): 0.2223933 0.2138044       0.22141

Again, there is no significant difference for the GMVP. However, for the MVP (which involves also the expected return) there is some small difference.

Let’s repeat the previous analysis but using monthly returns:

# change periodicity of prices
prices_monthly <- xts()
for (i in 1:ncol(prices))
  prices_monthly <- cbind(prices_monthly, Cl(to.monthly(prices[, i])))
colnames(prices_monthly) <- colnames(prices)
indexClass(prices_monthly) <- "Date"

# recompute returns
X_monthly_log <- CalculateReturns(prices_monthly, "log")[-1]
X_monthly_lin <- CalculateReturns(prices_monthly)[-1]
T_monthly <- nrow(X_monthly_log)  # number of months

# split data into training and set data
T_monthly_trn <- round(0.7*T_monthly)  # 70% of data
X_monthly_log_trn <- X_monthly_log[1:T_monthly_trn, ]
X_monthly_log_tst <- X_monthly_log[(T_monthly_trn+1):T_monthly, ]
X_monthly_lin_trn <- X_monthly_lin[1:T_monthly_trn, ]
X_monthly_lin_tst <- X_monthly_lin[(T_monthly_trn+1):T_monthly, ]

# estimate mu and Sigma
mu_monthly_lin <- colMeans(X_monthly_lin_trn)
Sigma_monthly_lin <- cov(X_monthly_lin_trn)

mu_monthly_log <- colMeans(X_monthly_log_trn)
Sigma_monthly_log <- cov(X_monthly_log_trn)

tmp <- momentsReturnLog2Lin(mu_monthly_log, Sigma_monthly_log)
mu_monthly_log_trans <- tmp$mu
Sigma_monthly_log_trans <- tmp$Sigma

# compute GMVPs
w_GMVP_monthly_lin <- GMVP(Sigma_monthly_lin)
w_GMVP_monthly_log <- GMVP(Sigma_monthly_log)
w_GMVP_monthly_log_trans <- GMVP(Sigma_monthly_log_trans)
w_GMVP_monthly_all <- cbind("GMVP_lin"       = w_GMVP_monthly_lin, 
                            "GMVP_log"       = w_GMVP_monthly_log, 
                            "GMVP_log_trans" = w_GMVP_monthly_log_trans)
barplot(t(w_GMVP_monthly_all), col = rainbow8equal[1:3], legend = colnames(w_GMVP_monthly_all), 
        main = "GMVP allocation for monthly returns", xlab = "stocks", ylab = "dollars", beside = TRUE)

# compute MVPs
w_MVP_monthly_lin <- MVP(mu_monthly_lin, Sigma_monthly_lin)
w_MVP_monthly_log <- MVP(mu_monthly_log, Sigma_monthly_log)
w_MVP_monthly_log_trans <- MVP(mu_monthly_log_trans, Sigma_monthly_log_trans)
w_MVP_monthly_all <- cbind("MVP_lin"       = w_MVP_monthly_lin, 
                           "MVP_log"       = w_MVP_monthly_log, 
                           "MVP_log_trans" = w_MVP_monthly_log_trans)
barplot(t(w_MVP_monthly_all), col = rainbow8equal[1:3], legend = colnames(w_MVP_monthly_all), 
        main = "MVP allocation for monthly returns", xlab = "stocks", ylab = "dollars", beside = TRUE)

The allocations seem to be very similar.

Let’s look at the portfolios performance:

# compute weekly returns of portfolios
ret_GMVP_monthly_all_trn <- xts(X_monthly_lin_trn %*% w_GMVP_monthly_all, index(X_monthly_lin_trn))
ret_MVP_monthly_all_trn <- xts(X_monthly_lin_trn %*% w_MVP_monthly_all, index(X_monthly_lin_trn))

# compare them
StdDev.annualized(ret_GMVP_monthly_all_trn)

R>>                                GMVP_lin  GMVP_log GMVP_log_trans
R>> Annualized Standard Deviation 0.1363082 0.1365655      0.1364513

SharpeRatio(ret_MVP_monthly_all_trn, FUN = "StdDev")

R>>                                 MVP_lin   MVP_log MVP_log_trans
R>> StdDev Sharpe (Rf=0%, p=95%): 0.6103515 0.6000131     0.5985821

Again, there is no significant difference.

We will estimate \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) from the in-sample log-returns:

mu <- colMeans(X_log_trn)
Sigma <- cov(X_log_trn)

Buy & Hold simply means that we allocate the whole budget to one stock and we stick to it. Since we have \(N=9\) stocks in our universe, we can define \(N=9\) different B&H portfolios, which we will store as column vectors.

# a B&H portfolio is trivially the zero vector with a one on the stock held
w_BnH <- diag(N)
rownames(w_BnH) <- colnames(X_lin)
colnames(w_BnH) <- paste0("B&H - ", colnames(X_lin))
w_BnH

R>>      B&H - AAPL B&H - AMD B&H - ADI B&H - ABBV B&H - AEZS B&H - A B&H - APD B&H - AA B&H - CF
R>> AAPL          1         0         0          0          0       0         0        0        0
R>> AMD           0         1         0          0          0       0         0        0        0
R>> ADI           0         0         1          0          0       0         0        0        0
R>> ABBV          0         0         0          1          0       0         0        0        0
R>> AEZS          0         0         0          0          1       0         0        0        0
R>> A             0         0         0          0          0       1         0        0        0
R>> APD           0         0         0          0          0       0         1        0        0
R>> AA            0         0         0          0          0       0         0        1        0
R>> CF            0         0         0          0          0       0         0        0        1

We can now compute the performance of those \(N=9\) portfolios in the training data with the package PerformanceAnalytics:

# compute returns of all B&H portfolios
ret_BnH <- xts(X_lin %*% w_BnH, index(X_lin))
ret_BnH_trn <- ret_BnH[1:T_trn, ]
ret_BnH_tst <- ret_BnH[-c(1:T_trn), ]
head(ret_BnH)

R>>              B&H - AAPL    B&H - AMD    B&H - ADI   B&H - ABBV   B&H - AEZS      B&H - A
R>> 2013-01-03 -0.012622240 -0.015810277 -0.016136502 -0.008257409  0.003952569  0.003581918
R>> 2013-01-04 -0.027854806  0.040160643 -0.017786948 -0.012632777  0.011811024  0.019747441
R>> 2013-01-07 -0.005882358  0.030888031  0.003057512  0.002035524  0.007782101 -0.007232493
R>> 2013-01-08  0.002691496  0.000000000 -0.010316326 -0.021764499 -0.015444015 -0.007990737
R>> 2013-01-09 -0.015628866 -0.014981273 -0.002606487  0.005636580  0.011764706  0.027007986
R>> 2013-01-10  0.012395942 -0.003802281  0.012114243  0.002949440 -0.003875969  0.007381675
R>>                B&H - APD     B&H - AA     B&H - CF
R>> 2013-01-03 -0.0034942112  0.008898612 -0.004728862
R>> 2013-01-04  0.0134405698  0.020948193  0.022398933
R>> 2013-01-07 -0.0009223002 -0.017278772 -0.003746079
R>> 2013-01-08  0.0018466463  0.000000000 -0.014660383
R>> 2013-01-09  0.0134809381 -0.002197690  0.034973559
R>> 2013-01-10 -0.0009096340 -0.012114520  0.015029411

# performance measures
library(PerformanceAnalytics)
t(table.AnnualizedReturns(ret_BnH_trn))

R>>            Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> B&H - AAPL            0.2629             0.2638                    0.9969
R>> B&H - AMD            -0.0773             0.4930                   -0.1568
R>> B&H - ADI             0.1517             0.2418                    0.6276
R>> B&H - ABBV            0.2167             0.2567                    0.8442
R>> B&H - AEZS           -0.7375             1.1854                   -0.6221
R>> B&H - A               0.0846             0.2221                    0.3808
R>> B&H - APD             0.2150             0.2057                    1.0448
R>> B&H - AA              0.0271             0.3045                    0.0891
R>> B&H - CF              0.1674             0.2838                    0.5901

t(table.AnnualizedReturns(ret_BnH_tst))

R>>            Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> B&H - AAPL            0.0375             0.2386                    0.1572
R>> B&H - AMD             3.2190             0.8017                    4.0153
R>> B&H - ADI             0.2000             0.2453                    0.8151
R>> B&H - ABBV            0.1532             0.3080                    0.4973
R>> B&H - AEZS           -0.3470             1.3677                   -0.2538
R>> B&H - A               0.2220             0.2329                    0.9532
R>> B&H - APD             0.1450             0.1917                    0.7561
R>> B&H - AA              0.2154             0.4755                    0.4531
R>> B&H - CF             -0.3103             0.5015                   -0.6188

Note how the in-sample (ex ante) performance is not maintained out-of-sample (ex post).

We can compute many other performance measures:

table.DownsideRisk(ret_BnH_trn)

R>> VaR calculation produces unreliable result (inverse risk) for column: 1 : -0.434364986114259

R>>                               B&H - AAPL B&H - AMD B&H - ADI B&H - ABBV B&H - AEZS B&H - A
R>> Semi Deviation                    0.0120    0.0221    0.0104     0.0120     0.0423  0.0099
R>> Gain Deviation                    0.0111    0.0231    0.0118     0.0095     0.0936  0.0091
R>> Loss Deviation                    0.0123    0.0237    0.0104     0.0117     0.0528  0.0093
R>> Downside Deviation (MAR=210%)     0.0162    0.0263    0.0149     0.0163     0.0471  0.0148
R>> Downside Deviation (Rf=0%)        0.0115    0.0220    0.0100     0.0115     0.0434  0.0097
R>> Downside Deviation (0%)           0.0115    0.0220    0.0100     0.0115     0.0434  0.0097
R>> Maximum Drawdown                  0.2586    0.6524    0.2639     0.2630     0.9845  0.2296
R>> Historical VaR (95%)             -0.0246   -0.0441   -0.0231    -0.0272    -0.0667 -0.0227
R>> Historical ES (95%)              -0.0370   -0.0744   -0.0324    -0.0377    -0.1512 -0.0304
R>> Modified VaR (95%)               -0.0264   -0.0474   -0.0192    -0.0273         NA -0.0217
R>> Modified ES (95%)                -0.0531   -0.0756   -0.0192    -0.0404    -0.0038 -0.0332
R>>                               B&H - APD B&H - AA B&H - CF
R>> Semi Deviation                   0.0087   0.0130   0.0120
R>> Gain Deviation                   0.0099   0.0145   0.0140
R>> Loss Deviation                   0.0081   0.0122   0.0117
R>> Downside Deviation (MAR=210%)    0.0134   0.0178   0.0165
R>> Downside Deviation (Rf=0%)       0.0083   0.0128   0.0116
R>> Downside Deviation (0%)          0.0083   0.0128   0.0116
R>> Maximum Drawdown                 0.1984   0.5359   0.3573
R>> Historical VaR (95%)            -0.0196  -0.0267  -0.0263
R>> Historical ES (95%)             -0.0268  -0.0417  -0.0386
R>> Modified VaR (95%)              -0.0170  -0.0276  -0.0222
R>> Modified ES (95%)               -0.0191  -0.0347  -0.0222

To compute the wealth or cumulative P&L, we have two options: one assumes the same quantity is repeateadly invested, whereas the other assumes reinvesting (compounding):

# compute cumulative wealth
wealth_arith_BnH_trn <- 1 + cumsum(ret_BnH_trn)  # initial budget of 1$
wealth_geom_BnH_trn <- cumprod(1 + ret_BnH_trn)  # initial budget of 1$

# plots
# same as: 
#   plot(wealth_arith_BnH_trn[, 1], main = "Buy & Hold performance (not compounded)", ylab = "wealth")
chart.CumReturns(ret_BnH_trn[, 1], main = "Buy & Hold performance (not compounded)", 
                 geometric = FALSE, wealth.index = TRUE)

# same as: 
#   plot(wealth_geom_BnH_trn[, 1], main = "Buy & Hold performance (compounded)", ylab = "wealth")
chart.CumReturns(ret_BnH_trn[, 1], main = "Buy & Hold performance (compounded)", 
                 geometric = TRUE, wealth.index = TRUE)

We can plot many more plots with the package PerformanceAnalytics:

# more plots
chart.CumReturns(ret_BnH, main = "Buy & Hold performance", 
                 wealth.index = TRUE, legend.loc = "topleft", colorset = rich10equal)

charts.PerformanceSummary(ret_BnH_trn, main = "Buy & Hold performance", 
                          wealth.index = TRUE, colorset = rich10equal)

chart.Boxplot(ret_BnH_trn)

chart.RiskReturnScatter(ret_BnH_trn, symbolset = 21, bg = "red")

We are now ready to consider several portfolio designs and compare their performance (the underlying constraint is simply \(\mathbf{w}\ge\mathbf{0}\) and \(\mathbf{1}^T\mathbf{w}=1\)).

The EWP or \(1/N\) portfolio allocates equal dollar weight to each stock: \(\mathbf{w} = \frac{1}{N}\mathbf{1}\).

w_EWP <- rep(1/N, N)
names(w_EWP) <- colnames(X_lin)
w_EWP

R>>      AAPL       AMD       ADI      ABBV      AEZS         A       APD        AA        CF 
R>> 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111

Quintile portfolios are widely used by practitioners. The idea is to 1) rank the \(N\) stocks, 2) divide them into five parts, and 3) long the top part (and possibly short the bottom part). One can rank the stocks in a multitude of ways (typically based on expensive factors that investment funds buy at a premium price). For our experiments, we will consider three possible rankings according to:

1. \(\boldsymbol{\mu}\)
2. \(\boldsymbol{\mu}/\textsf{diag}(\boldsymbol{\Sigma})\)
3. \(\boldsymbol{\mu}/\sqrt{(\textsf{diag}(\boldsymbol{\Sigma}))}\)

# find indices of sorted stocks
i1 <- sort(mu, decreasing = TRUE, index.return = TRUE)$ix
i2 <- sort(mu/diag(Sigma), decreasing = TRUE, index.return = TRUE)$ix
i3 <- sort(mu/sqrt(diag(Sigma)), decreasing = TRUE, index.return = TRUE)$ix

# create portfolios
w_QuintP_1 <- w_QuintP_2 <- w_QuintP_3 <- rep(0, N)
w_QuintP_1[i1[1:round(N/5)]] <- 1/round(N/5)
w_QuintP_2[i2[1:round(N/5)]] <- 1/round(N/5)
w_QuintP_3[i3[1:round(N/5)]] <- 1/round(N/5)
w_QuintP <- cbind("QuintP (mu)"        = w_QuintP_1, 
                  "QuintP (mu/sigma2)" = w_QuintP_2, 
                  "QuintP (mu/sigma)"  = w_QuintP_3)
rownames(w_QuintP) <- colnames(X_lin)
w_QuintP

R>>      QuintP (mu) QuintP (mu/sigma2) QuintP (mu/sigma)
R>> AAPL         0.5                0.5               0.5
R>> AMD          0.0                0.0               0.0
R>> ADI          0.0                0.0               0.0
R>> ABBV         0.5                0.0               0.0
R>> AEZS         0.0                0.0               0.0
R>> A            0.0                0.0               0.0
R>> APD          0.0                0.5               0.5
R>> AA           0.0                0.0               0.0
R>> CF           0.0                0.0               0.0

The global maximum return portfolio (GMRP) chooses the stock with the hightest return during the in-sample period:

i_max <- which.max(mu)
w_GMRP <- rep(0, N)
w_GMRP[i_max] <- 1
names(w_GMRP) <- colnames(X_lin)
w_GMRP

R>> AAPL  AMD  ADI ABBV AEZS    A  APD   AA   CF 
R>>    1    0    0    0    0    0    0    0    0

We can now compare the allocations of the portfolios:

# put together all portfolios
w_heuristic <- cbind("EWP" = w_EWP, w_QuintP, "GMRP" = w_GMRP)
round(w_heuristic, digits = 2)

R>>       EWP QuintP (mu) QuintP (mu/sigma2) QuintP (mu/sigma) GMRP
R>> AAPL 0.11         0.5                0.5               0.5    1
R>> AMD  0.11         0.0                0.0               0.0    0
R>> ADI  0.11         0.0                0.0               0.0    0
R>> ABBV 0.11         0.5                0.0               0.0    0
R>> AEZS 0.11         0.0                0.0               0.0    0
R>> A    0.11         0.0                0.0               0.0    0
R>> APD  0.11         0.0                0.5               0.5    0
R>> AA   0.11         0.0                0.0               0.0    0
R>> CF   0.11         0.0                0.0               0.0    0

barplot(t(w_heuristic), col = rainbow8equal[1:5], legend = colnames(w_heuristic), beside = TRUE,
        main = "Portfolio allocation of heuristic portfolios", xlab = "stocks", ylab = "dollars")

Then we can compare the performance (in-sample vs out-of-sample):

# compute returns of all portfolios
ret_heuristic <- xts(X_lin %*% w_heuristic, index(X_lin))
ret_heuristic$`QuintP (mu/sigma2)` <- NULL  # remove since it coincides with "QuintP (mu/sigma)"
ret_heuristic_trn <- ret_heuristic[1:T_trn, ]
ret_heuristic_tst <- ret_heuristic[-c(1:T_trn), ]

# performance
t(table.AnnualizedReturns(ret_heuristic_trn))

R>>                   Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> EWP                          0.0405             0.2059                    0.1967
R>> QuintP (mu)                  0.2573             0.1988                    1.2945
R>> QuintP (mu/sigma)            0.2518             0.1870                    1.3461
R>> GMRP                         0.2629             0.2638                    0.9969

t(table.AnnualizedReturns(ret_heuristic_tst))

R>>                   Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> EWP                          0.3397             0.2770                    1.2265
R>> QuintP (mu)                  0.1103             0.2142                    0.5151
R>> QuintP (mu/sigma)            0.0988             0.1748                    0.5652
R>> GMRP                         0.0375             0.2386                    0.1572

Let’s plot the wealth evolution (cumulative PnL) over time:

{ chart.CumReturns(ret_heuristic, main = "Cumulative return of heuristic portfolios", 
                   wealth.index = TRUE, legend.loc = "topleft", colorset = rich8equal)
  addEventLines(xts("training", index(X_lin[T_trn])), srt=90, pos=2, lwd = 2, col = "darkblue") }

charts.PerformanceSummary(ret_heuristic, main = "Performance of heuristic portfolios", 
                          wealth.index = TRUE, colorset = rich8equal)

Finally, we can plot the risk-return scatter plot:

chart.RiskReturnScatter(ret_heuristic_trn, symbolset = 21, bg = "red",
                        main = "Annualized Return and Risk (in-sample)")

Markowitz’s mean-variance portfolio (MVP) with no shorting is formulated as \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \boldsymbol{\mu}^T\mathbf{w} -\lambda\mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

This problem does not have a closed-form solution and we need to resort to using a solver. It is very convenient to use the package CVXR (although the computational cost will be high and the solution not totally robust, if necessary use a QP solver like quadprog):

library(CVXR)

# create function for MVP
MVP <- function(mu, Sigma, lmd = 0.5) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Maximize(t(mu) %*% w - lmd*quad_form(w, Sigma)),
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  names(w) <- colnames(Sigma)
  return(w)
}

# this function can now be used as
w_MVP <- MVP(mu, Sigma, lmd = 2)

The Global Minimum Variance Portfolio (GMVP) with no shorting is formulated as \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0} \end{array} \]

Since a closed-form solution does not exist with the constraint \(\mathbf{w}\ge\mathbf{0}\), we need to resort to a solver. We can convenientky use the package CVXR (although the computational cost will be high and the solution not totally robust, if necessary use a QP solver like quadprog):

library(CVXR)

# create function for GMVP
GMVP <- function(Sigma) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Minimize(quad_form(w, Sigma)), 
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  names(w) <- colnames(Sigma)
  return(w)
}

# this function can now be used as
w_GMVP <- GMVP(Sigma)

The maximum Sharpe ratio portfolio (MSRP) is given by the solution to the nonconvex problem \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \dfrac{\mathbf{w}^{T}\boldsymbol{\mu}}{\sqrt{\mathbf{w}^{T}\boldsymbol{\Sigma}\mathbf{w}}}\\ \textsf{subject to} & \mathbf{1}^{T}\mathbf{w}=1, \quad\mathbf{w}\geq\mathbf{0}. \end{array}\]

We will solve this problem with the general-purpose nonlinear solver nloptr in R:

library(nloptr)

# define the nonconvex objective function
fn_SR <- function(w) {
  return(as.numeric(t(w) %*% mu / sqrt(t(w) %*% Sigma %*% w)))
}
  
# initial point
w0 <- rep(1/N, N)

res <- nloptr::slsqp(w0, fn_SR,
                     lower = rep(0, N),  # w >= 0
                     heq = function(w) return(sum(w) - 1))    # sum(w) = 1
w_nonlinear_solver <- res$par
res

R>> $par
R>> [1] 1.126793e-19 1.611956e-01 6.598538e-18 1.156482e-18 8.388044e-01 0.000000e+00 1.909453e-18
R>> [8] 5.671798e-17 0.000000e+00
R>> 
R>> $value
R>> [1] -0.07987582
R>> 
R>> $iter
R>> [1] 68
R>> 
R>> $convergence
R>> [1] 4
R>> 
R>> $message
R>> [1] "NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached."

# res <- alabama::constrOptim.nl(w0, fn_SR,
#                                hin = function(w) return(w),  # w >= 0
#                                heq = function(w) return(sum(w) - 1),    # sum(w) = 1
#                                control.outer = list(trace = FALSE))
# w_nonlinear_solver <- res$par
# res

We are going to solve the nonconvex problem \[\begin{array}{ll} \underset{\mathbf{w},t}{\textsf{maximize}} & t\\ \textsf{subject to} & t \leq \dfrac{\mathbf{w}^{T}\boldsymbol{\mu}}{\sqrt{\mathbf{w}^{T}\boldsymbol{\Sigma}\mathbf{w}}}\\ & \mathbf{1}^{T}\mathbf{w}=1, \quad\left(\mathbf{w}\geq\mathbf{0}\right). \end{array}\] via bisection on \(t\) with the following (convex) SOCP problem for a given \(t\): \[\begin{array}{ll} \underset{\;}{\textsf{find}} & \mathbf{w}\\ \textsf{subject to} & t \left\Vert \boldsymbol{\Sigma}^{1/2}\mathbf{w}\right\Vert_{2}\leq\mathbf{w}^{T}\boldsymbol{\mu}\\ & \mathbf{1}^{T}\mathbf{w}=1, \quad\left(\mathbf{w}\geq\mathbf{0}\right). \end{array}\]

# define the inner solver based on an SOCP solver 
# (we will simply use CVXR for convenience, see: https://cvxr.rbind.io/cvxr_functions/)
library(CVXR)

# square-root of matrix Sigma
Sigma_12 <- chol(Sigma)
max(abs(t(Sigma_12) %*% Sigma_12 - Sigma))  # sanity check

R>> [1] 5.421011e-20

# create function for MVP
SOCP_bisection <- function(t) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Maximize(0),
                  constraints = list(t*cvxr_norm(Sigma_12 %*% w, 2) <= t(mu) %*% w,
                                     sum(w) == 1,
                                     w >= 0))
  result <- solve(prob)
  return(list("status" = result$status, "w" = as.vector(result$getValue(w))))
}

# now run the bisection algorithm
t_lb <- 0   # for sure the problem is feasible in this case
t_ub <- 10  # a tighter upper bound coud be chose, but a Sharpe ratio of 10 surely cannot be achieved
while(t_ub - t_lb > 1e-6) {
  t <- (t_ub + t_lb)/2  # midpoint
  if(SOCP_bisection(t)$status == "infeasible")
    t_ub <- t
  else
    t_lb <- t
}
w_bisection <- SOCP_bisection(t_lb)$w

# comparison between two solutions
round(cbind(w_nonlinear_solver, w_bisection), digits = 3)

R>>       w_nonlinear_solver w_bisection
R>>  [1,]              0.000       0.325
R>>  [2,]              0.161       0.000
R>>  [3,]              0.000       0.000
R>>  [4,]              0.000       0.223
R>>  [5,]              0.839       0.000
R>>  [6,]              0.000       0.000
R>>  [7,]              0.000       0.399
R>>  [8,]              0.000       0.000
R>>  [9,]              0.000       0.053

# Sharpe ratio of two solutions
c("nonlinear_solver" = fn_SR(w_nonlinear_solver), 
  "bisection"        = fn_SR(w_bisection))

R>> nonlinear_solver        bisection 
R>>      -0.07987582       0.07747009

We can see that the nonlinear solver could not properly solve the nonconvex problem. However, using bisection we can solve the problem optimally.

We are going to solve the nonconvex problem \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \dfrac{\mathbf{w}^{T}\boldsymbol{\mu}}{\sqrt{\mathbf{w}^{T}\boldsymbol{\Sigma}\mathbf{w}}}\\ \textsf{subject to} & \mathbf{1}^{T}\mathbf{w}=1, \quad \mathbf{w}\geq\mathbf{0} \end{array}\] by iteratively solving the (convex) SOCP problem for a given \(y^{(k)}\): \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \mathbf{w}^{T}\boldsymbol{\mu} - y^{(k)} \left\Vert \boldsymbol{\Sigma}^{1/2}\mathbf{w}\right\Vert_{2}\\ \textsf{subject to} & \mathbf{1}^{T}\mathbf{w}=1, \quad \mathbf{w}\geq\mathbf{0}. \end{array}\] where \[y^{(k)} = \frac{\mathbf{w}^{(k)T}\boldsymbol{\mu}}{\sqrt{\mathbf{w}^{(k)T}\boldsymbol{\Sigma}\mathbf{w}^{(k)}}}.\]

# define the inner solver based on an SOCP solver 
# (we will simply use CVXR for convenience, see: https://cvxr.rbind.io/cvxr_functions/)
library(CVXR)

# square-root of matrix Sigma
Sigma_12 <- chol(Sigma)
max(abs(t(Sigma_12) %*% Sigma_12 - Sigma))  # sanity check

R>> [1] 5.421011e-20

# create function for MVP
SOCP_Dinkelbach <- function(y) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Maximize(t(mu) %*% w - y*cvxr_norm(Sigma_12 %*% w, 2)),
                  constraints = list(sum(w) == 1,
                                     w >= 0))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# initial point (has to satisfy t(w_k) %*% mu>=0)
i_max <- which.max(mu)
w_k <- rep(0, N)  
w_k[i_max] <- 1

# now the iterative Dinkelbach algorithm
k <- 1
while(k == 1 || max(abs(w_k - w_prev)) > 1e-6) {
  w_prev <- w_k
  y_k <- as.numeric(t(w_k) %*% mu / sqrt(t(w_k) %*% Sigma %*% w_k))
  w_k <- SOCP_Dinkelbach(y_k)
  k <- k + 1
}
w_Dinkelbach <- w_k
cat("Number of iterarions:", k-1)

R>> Number of iterarions: 4

# comparison among three solutions
round(cbind(w_nonlinear_solver, w_bisection, w_Dinkelbach), digits = 3)

R>>       w_nonlinear_solver w_bisection w_Dinkelbach
R>>  [1,]              0.000       0.325        0.325
R>>  [2,]              0.161       0.000        0.000
R>>  [3,]              0.000       0.000        0.000
R>>  [4,]              0.000       0.223        0.223
R>>  [5,]              0.839       0.000        0.000
R>>  [6,]              0.000       0.000        0.000
R>>  [7,]              0.000       0.399        0.399
R>>  [8,]              0.000       0.000        0.000
R>>  [9,]              0.000       0.053        0.052

# Sharpe ratio of three solutions
c("nonlinear_solver" = fn_SR(w_nonlinear_solver), 
  "bisection"        = fn_SR(w_bisection),
  "Dinkelbach"       = fn_SR(w_Dinkelbach))

R>> nonlinear_solver        bisection       Dinkelbach 
R>>      -0.07987582       0.07747009       0.07747018

As expected, the Dinkelbach method converges to the optimal solution, like the bisection method (but unlike the general-purpose nonlinear solver).

The maximum Sharpe ratio portfolio (MSRP) is the nonconvex problem \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \dfrac{\mathbf{w}^{T}\boldsymbol{\mu}}{\sqrt{\mathbf{w}^{T}\boldsymbol{\Sigma}\mathbf{w}}}\\ \textsf{subject to} & \mathbf{1}^{T}\mathbf{w}=1, \quad \mathbf{w}\geq\mathbf{0} \end{array}\] that can be rewritten in convex form as \[ \begin{array}{ll} \underset{\tilde{\mathbf{w}}}{\textsf{minimize}} & \tilde{\mathbf{w}}^T\mathbf{\Sigma}\tilde{\mathbf{w}}\\ {\textsf{subject to}} & \tilde{\mathbf{w}}^T\boldsymbol{\mu} = 1\\ & \tilde{\mathbf{w}}\ge\mathbf{0} \end{array} \] and then \(\mathbf{w} = \tilde{\mathbf{w}}/(\mathbf{1}^T\tilde{\mathbf{w}})\).

This is a quadratic problem (QP) and we can conveniently use CVXR (although one is advised to use a specific QP solver like quadprog for speed and stability):

# create function for MSRP
MSRP <- function(mu, Sigma) {
  w_ <- Variable(nrow(Sigma))
  prob <- Problem(Minimize(quad_form(w_, Sigma)),
                  constraints = list(w_ >= 0, t(mu) %*% w_ == 1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w_)/sum(result$getValue(w_)))
  names(w) <- colnames(Sigma)
  return(w)
}

# this function can now be used as
w_MSRP <- MSRP(mu, Sigma)

# comparison among solutions
round(cbind(w_nonlinear_solver, w_bisection, w_Dinkelbach, w_MSRP), digits = 3)

R>>      w_nonlinear_solver w_bisection w_Dinkelbach w_MSRP
R>> AAPL              0.000       0.325        0.325  0.325
R>> AMD               0.161       0.000        0.000  0.000
R>> ADI               0.000       0.000        0.000  0.000
R>> ABBV              0.000       0.223        0.223  0.224
R>> AEZS              0.839       0.000        0.000  0.000
R>> A                 0.000       0.000        0.000  0.000
R>> APD               0.000       0.399        0.399  0.399
R>> AA                0.000       0.000        0.000  0.000
R>> CF                0.000       0.053        0.052  0.052

# Sharpe ratio of different solutions
c("nonlinear_solver" = fn_SR(w_nonlinear_solver), 
  "bisection"        = fn_SR(w_bisection),
  "Dinkelbach"       = fn_SR(w_Dinkelbach),
  "Schaible"         = fn_SR(w_MSRP))

R>> nonlinear_solver        bisection       Dinkelbach         Schaible 
R>>      -0.07987582       0.07747009       0.07747018       0.07747018

As expected, the bisection method, Dinkelbach method, and Schaible method give the optimal solution, unlike the general-purpose nonlinear solver which solves a nonconvex problem.

The Inverse volatility portfolio (IVP) has the simple closed-form solution: \[\mathbf{w} = \frac{\boldsymbol{\sigma}^{-1}}{\mathbf{1}^T\boldsymbol{\sigma}^{-1}}.\]

Its implementation in R is trivial:

# create function for IVP
IVP <- function(Sigma) {
  sigma <- sqrt(diag(Sigma))
  w <- 1/sigma
  w <- w/sum(w)
  return(w)
}

# this function can now be used as
w_IVP <- IVP(Sigma)

The Most diversified portfolio (MDP) is formulated exactly as the maximum Sharpe ratio portfolio (MSRP) but using \(\boldsymbol{\sigma}\) in lieu of \(\boldsymbol{\mu}\): \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \dfrac{\mathbf{w}^{T}\boldsymbol{\sigma}}{\sqrt{\mathbf{w}^{T}\boldsymbol{\Sigma}\mathbf{w}}}\\ \textsf{subject to} & \mathbf{1}^{T}\mathbf{w}=1\\ & \tilde{\mathbf{w}}\ge\mathbf{0}. \end{array}\]

Therefore, we can use the same R function created to solve the MSRP:

w_MDP <- MSRP(mu = sqrt(diag(Sigma)), Sigma)

The maximum decorrelation portfolio (MDCP) is formally formulated as the GMVP but using \(\mathbf{C}\) in lieu of \(\boldsymbol{\Sigma}\): \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \mathbf{w}^{T}\mathbf{C}\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^{T}\mathbf{w}=1 \end{array}\] where \(\mathbf{C} \triangleq \mathsf{Diag}(\boldsymbol{\Sigma})^{-1/2}\boldsymbol{\Sigma}\mathsf{Diag}(\boldsymbol{\Sigma})^{-1/2}\) is the correlation matrix.

Therefore, we can use the same R function created to solve the GMVP:

# create function for MDCP based on GMVP()
MDCP <- function(Sigma) {
  C <- diag(1/sqrt(diag(Sigma))) %*% Sigma %*% diag(1/sqrt(diag(Sigma)))
  colnames(C) <- colnames(Sigma)
  return(GMVP(Sigma = C))
}

# this function can now be used as
w_MDCP <- MDCP(Sigma)

We are now ready to consider several portfolio designs and compare their performance (the underlying constraint is simply \(\mathbf{w}\ge\mathbf{0}\) and \(\mathbf{1}^T\mathbf{w}=1\)).

Recall the three heuristic portfolios we want to compare:

w_heuristic <- cbind("EWP"    = w_EWP, 
                     "QuintP" = w_QuintP[, "QuintP (mu/sigma)"], 
                     "GMRP"   = w_GMRP)

Now we stack the Markowitz-based portfolios we want to compare:

w_Markowitz <- cbind("MVP"  = w_MVP, 
                     "GMVP" = w_GMVP, 
                     "IVP"  = w_IVP, 
                     "MSRP" = w_MSRP, 
                     "MDP"  = w_MDP, 
                     "MDCP" = w_MDCP)

We can now compare the allocations of the portfolios:

w_all <- cbind(w_heuristic, w_Markowitz)
barplot(t(w_all), col = rainbow10equal[1:9], legend = colnames(w_all), beside = TRUE,
        main = "Portfolio allocation", xlab = "stocks", ylab = "dollars")

The performance is (in-sample vs out-of-sample):

# compute returns of all portfolios
ret_all <- xts(X_lin %*% w_all, index(X_lin))
ret_all_trn <- ret_all[1:T_trn, ]
ret_all_tst <- ret_all[-c(1:T_trn), ]

# performance
t(table.AnnualizedReturns(ret_all_trn))

R>>        Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> EWP               0.0405             0.2059                    0.1967
R>> QuintP            0.2518             0.1870                    1.3461
R>> GMRP              0.2629             0.2638                    0.9969
R>> MVP               0.2518             0.1719                    1.4649
R>> GMVP              0.1873             0.1576                    1.1886
R>> IVP               0.1386             0.1644                    0.8434
R>> MSRP              0.2471             0.1669                    1.4807
R>> MDP               0.0935             0.1886                    0.4955
R>> MDCP             -0.1027             0.3287                   -0.3123

t(table.AnnualizedReturns(ret_all_tst))

R>>        Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> EWP               0.3397             0.2770                    1.2265
R>> QuintP            0.0988             0.1748                    0.5652
R>> GMRP              0.0375             0.2386                    0.1572
R>> MVP               0.1196             0.1676                    0.7139
R>> GMVP              0.1430             0.1779                    0.8038
R>> IVP               0.2494             0.2110                    1.1821
R>> MSRP              0.1033             0.1664                    0.6207
R>> MDP               0.3146             0.2573                    1.2229
R>> MDCP              0.4050             0.4246                    0.9538

We can observe that:

• as expected, the MSRP achieves the maximum ex ante (in-sample) Sharpe ratio; however, this is not maintained ex post (out-of-sample);
• as expected, the EWP achieves the best ex post performance;
• the MVP performs poorly;
• other top performers include the IVP and MDP.

Let’s plot the wealth evolution (cumulative PnL) over the whole time:

{ chart.CumReturns(ret_all, main = "Cumulative return of portfolios", 
                   wealth.index = TRUE, legend.loc = "topleft", colorset = rich10equal)
  addEventLines(xts("training", index(X_lin[T_trn])), srt=90, pos=2, lwd = 2, col = "darkblue") }

and let’s zoom in the out-of-sample period:

chart.CumReturns(ret_all_tst, main = "Cumulative return of portfolios (out-of-sample)",
                   wealth.index = TRUE, legend.loc = "topleft", colorset = rich10equal)

Let’s look at the drawdown:

chart.Drawdown(ret_all_tst, main = "Drawdown of portfolios (out-of-sample)", 
               legend.loc = "bottomleft", colorset = rich10equal)

Clearly the MDCP and the GMRP have the worst drawdown and are unacceptable.

Finally, we can plot the (in-sample) return-risk scatter plot along with the efficient frontier:

# first, compute the efficient frontier
w_frontier_trn <- NULL
lmd_sweep <- exp(seq(-6, 6, by = 0.5))
for (lmd in lmd_sweep)
  w_frontier_trn <- cbind(w_frontier_trn, MVP(mu, Sigma, lmd))
ret_frontier_trn <- xts(X_lin_trn %*% w_frontier_trn, index(X_lin_trn))
mu_frontier_trn <- table.AnnualizedReturns(ret_frontier_trn)[1, ]
sd_frontier_trn <- table.AnnualizedReturns(ret_frontier_trn)[2, ]

# plot in-sample sd-mu scatter plot
maxSR <- table.AnnualizedReturns(ret_all_trn[, "MSRP"])[3, ]
chart.RiskReturnScatter(ret_all_trn,
                        main = "Annualized Return and Risk (in-sample)",
                        symbolset = c(rep(21, 3), rep(22, 6)), 
                        colorset = c(rep("darkred", 3), rep("darkblue", 6)),
                        bg = "black",
                        add.sharpe = maxSR)
lines(sd_frontier_trn, mu_frontier_trn)

Observe that this nice return-risk scatter plot totally deforms in an unpredictable way when we evaluate it in the out-of-sample set as we can see:

# compute the efficient frontier again but based on the test data
mu_tst <- colMeans(X_log_tst)
Sigma_tst <- cov(X_log_tst)

w_frontier_tst <- NULL
lmd_sweep <- exp(seq(-6, 6, by = 0.5))
for (lmd in lmd_sweep)
  w_frontier_tst <- cbind(w_frontier_tst, MVP(mu_tst, Sigma_tst, lmd))
ret_frontier_tst <- xts(X_lin_tst %*% w_frontier_tst, index(X_lin_tst))
mu_frontier_tst <- table.AnnualizedReturns(ret_frontier_tst)[1, ]
sd_frontier_tst <- table.AnnualizedReturns(ret_frontier_tst)[2, ]

# plot out-of-sample sd-mu scatter plot
chart.RiskReturnScatter(ret_all_tst,
                        main = "Annualized Return and Risk (out-of-sample)",
                        symbolset = c(rep(21, 3), rep(22, 6)),  
                        colorset = c(rep("darkred", 3), rep("darkblue", 6)),
                        bg = "black",
                        add.sharpe = NA,
                        ylim = c(0, 1))
lines(sd_frontier_tst, mu_frontier_tst)

Again this ex post return-risk scatter plot shows that the winners are EWP, IVP, and MDP.

In the previous performance analysis of the different portfolios, we were implicitly assuming a daily rebalancing and that’s why computing the return of the portfolio was as easy as ret_portf <- X_lin %*% w. Now, we will consider a lower rebalancing frequency (weekly, monthly, quarterly, and yearly) for a given portfolio (in particular, we will just use the EWP portfolio w_EWP). Observe that on the day of the rebalancing, the portfolio held will be the designed one w_EWP but in the subsequent days without rebalancing such portfolio slowly deviates as the prices of the stocks change (recall that the portfolio denotes the dollar or capital allocation).

The package PerformanceAnalytics has the convenient function Return.portfolio() that allows us to compute the portfolio return with different rebalancing schemes. Let’s start with a sanity check for the daily rebalancing:

# choose portfolio for the comparison
w <- w_EWP

# recall the computation for the daily rebalancing (with daily returns)
ret_daily_rebal <- X_lin %*% w

# we can alternativaly use this function:
ret_daily_rebal_ <- Return.portfolio(X_lin, weights = w, rebalance_on = "days")
norm(ret_daily_rebal - ret_daily_rebal_)  # sanity check

R>> [1] 8.420819e-14

Now, let’s use a yearly rebalancing frequency to observe the slow deviation as the prices change:

# let's observe how the portfolio slowly deviates from its original design
tmp <- Return.portfolio(X_lin, weights = w, rebalance_on = "years", verbose = TRUE)
round(head(tmp$BOP.Weight, 15), digits = 4)

R>>              AAPL    AMD    ADI   ABBV   AEZS      A    APD     AA     CF
R>> 2013-01-03 0.1111 0.1111 0.1111 0.1111 0.1111 0.1111 0.1111 0.1111 0.1111
R>> 2013-01-04 0.1103 0.1099 0.1099 0.1107 0.1121 0.1121 0.1113 0.1127 0.1111
R>> 2013-01-07 0.1063 0.1134 0.1071 0.1085 0.1125 0.1134 0.1119 0.1141 0.1127
R>> 2013-01-08 0.1056 0.1168 0.1073 0.1086 0.1133 0.1125 0.1117 0.1120 0.1122
R>> 2013-01-09 0.1067 0.1177 0.1070 0.1070 0.1124 0.1124 0.1127 0.1129 0.1114
R>> 2013-01-10 0.1043 0.1152 0.1060 0.1069 0.1130 0.1147 0.1135 0.1119 0.1145
R>> 2013-01-11 0.1053 0.1144 0.1069 0.1069 0.1122 0.1152 0.1130 0.1102 0.1159
R>> 2013-01-14 0.1037 0.1155 0.1051 0.1068 0.1186 0.1135 0.1127 0.1089 0.1152
R>> 2013-01-15 0.1008 0.1151 0.1059 0.1083 0.1164 0.1147 0.1130 0.1094 0.1164
R>> 2013-01-16 0.0978 0.1170 0.1055 0.1101 0.1174 0.1141 0.1128 0.1096 0.1157
R>> 2013-01-17 0.1015 0.1183 0.1053 0.1127 0.1127 0.1132 0.1119 0.1087 0.1157
R>> 2013-01-18 0.0998 0.1180 0.1067 0.1144 0.1120 0.1142 0.1122 0.1084 0.1143
R>> 2013-01-22 0.1000 0.1068 0.1075 0.1181 0.1124 0.1158 0.1134 0.1099 0.1161
R>> 2013-01-23 0.1004 0.1058 0.1071 0.1146 0.1122 0.1164 0.1137 0.1118 0.1179
R>> 2013-01-24 0.1007 0.1161 0.1060 0.1172 0.1097 0.1138 0.1108 0.1095 0.1162

chart.StackedBar(tmp$BOP.Weight, main = "Evolution of uniform portfolio with yearly rebalancing",
                 ylab = "w", space = 0, border = NA)

We can now compare the different rebalancing frequencies:

# now different rebalancing frequencies
ret_weekly_rebal <- Return.portfolio(X_lin, weights = w, rebalance_on = "weeks")
ret_monthly_rebal <- Return.portfolio(X_lin, weights = w, rebalance_on = "months")
ret_quarterly_rebal <- Return.portfolio(X_lin, weights = w, rebalance_on = "quarters")
ret_yearly_rebal <- Return.portfolio(X_lin, weights = w, rebalance_on = "years")
ret_allfreqs <- cbind(ret_daily_rebal, ret_weekly_rebal, ret_monthly_rebal, ret_quarterly_rebal, ret_yearly_rebal)
colnames(ret_allfreqs) <- c("daily", "weekly", "monthly", "quarterly", "yearly")
round(head(ret_allfreqs, 25), digits = 4)

R>>              daily  weekly monthly quarterly  yearly
R>> 2013-01-03 -0.0050 -0.0050 -0.0050   -0.0050 -0.0050
R>> 2013-01-04  0.0078  0.0079  0.0079    0.0079  0.0079
R>> 2013-01-07  0.0010  0.0010  0.0010    0.0010  0.0010
R>> 2013-01-08 -0.0073 -0.0073 -0.0073   -0.0073 -0.0073
R>> 2013-01-09  0.0064  0.0062  0.0064    0.0064  0.0064
R>> 2013-01-10  0.0032  0.0033  0.0031    0.0031  0.0031
R>> 2013-01-11  0.0087  0.0087  0.0089    0.0089  0.0089
R>> 2013-01-14 -0.0075 -0.0075 -0.0075   -0.0075 -0.0075
R>> 2013-01-15 -0.0019 -0.0018 -0.0016   -0.0016 -0.0016
R>> 2013-01-16  0.0042  0.0041  0.0035    0.0035  0.0035
R>> 2013-01-17  0.0099  0.0101  0.0100    0.0100  0.0100
R>> 2013-01-18 -0.0070 -0.0070 -0.0075   -0.0075 -0.0075
R>> 2013-01-22  0.0055  0.0055  0.0054    0.0054  0.0054
R>> 2013-01-23  0.0160  0.0158  0.0152    0.0152  0.0152
R>> 2013-01-24 -0.0119 -0.0117 -0.0104   -0.0104 -0.0104
R>> 2013-01-25  0.0055  0.0061  0.0064    0.0064  0.0064
R>> 2013-01-28  0.0034  0.0034  0.0025    0.0025  0.0025
R>> 2013-01-29  0.0198  0.0210  0.0202    0.0202  0.0202
R>> 2013-01-30 -0.0140 -0.0152 -0.0154   -0.0154 -0.0154
R>> 2013-01-31 -0.0106 -0.0114 -0.0116   -0.0116 -0.0116
R>> 2013-02-01  0.0101  0.0101  0.0101    0.0103  0.0103
R>> 2013-02-04 -0.0196 -0.0196 -0.0195   -0.0196 -0.0196
R>> 2013-02-05  0.0102  0.0102  0.0101    0.0092  0.0092
R>> 2013-02-06  0.0032  0.0030  0.0030    0.0033  0.0033
R>> 2013-02-07  0.0113  0.0111  0.0110    0.0105  0.0105

# performance
t(table.AnnualizedReturns(ret_allfreqs))

R>>           Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> daily                0.1224             0.2296                    0.5334
R>> weekly               0.1186             0.2299                    0.5160
R>> monthly              0.1050             0.2268                    0.4629
R>> quarterly            0.0859             0.2066                    0.4156
R>> yearly               0.1310             0.2020                    0.6485

chart.CumReturns(ret_allfreqs, main = "Daily wealth for different rebalancing frequencies",
                 wealth.index = TRUE, legend.loc = "topleft", colorset = rich6equal)

Until now, we have considered static portfolios in the sense that they are first designed based on a training set and then they remain fixed and are used in the test set. In a more realistic setting, however, one wants to implement this procedure in a rolling-window basis. That is, with some frequency the portfolio is reoptimized (since the dynamics of the asset returns may change with time) and rebalanced.

Recall the procedure for the static portfolio where one estimates \(\boldsymbol\mu\) and \(\boldsymbol\Sigma\) from the training set, designs some portfolio (say, Markowitz mean-variance), and then applies it to the test set with some rebalancing frequency (say, monthly):

# recall the procedure for the static portfolio (Markowitz, for example):
T_trn <- round(0.5*T)
X_log_trn <- X_log[1:T_trn, ]
X_log_tst <- X_log[(T_trn+1):T, ]
X_lin_trn <- X_lin[1:T_trn, ]
X_lin_tst <- X_lin[(T_trn+1):T, ]
mu <- colMeans(X_log_trn)
Sigma <- cov(X_log_trn)
w_MVP_static <- MVP(mu, Sigma, lmd = 2)
ret_static <- Return.portfolio(X_lin_tst, weights = w_MVP_static, rebalance_on = "months")

Now we want to reoptimize the portfolio on a rolling-window basis. R provides some options:

• The package PerformanceAnalytics provides two functions to handle the two main rolling-window strategies: apply.rolling() is for a rolling window of fixed length, whereas apply.fromstart() is for an expanding rolling window. However, those only work independently on each column of the multivariate time series, which is not what we want.
• The package xts has the function apply.rolling() but it only works with nonoverlapping periods, which is not what we want. What we want is to choose a window of, say, 6 months (about 120 days) and shift it, say, every month.

We will then do it by ourselves with a simple loop:

# create empty portfolio matrix
w_MVP_rolling <- X_log
w_MVP_rolling[] <- NA
head(w_MVP_rolling)

R>>            AAPL AMD ADI ABBV AEZS  A APD AA CF
R>> 2013-01-03   NA  NA  NA   NA   NA NA  NA NA NA
R>> 2013-01-04   NA  NA  NA   NA   NA NA  NA NA NA
R>> 2013-01-07   NA  NA  NA   NA   NA NA  NA NA NA
R>> 2013-01-08   NA  NA  NA   NA   NA NA  NA NA NA
R>> 2013-01-09   NA  NA  NA   NA   NA NA  NA NA NA
R>> 2013-01-10   NA  NA  NA   NA   NA NA  NA NA NA

# find rebalancing indices
rebal_indices <- T_trn + endpoints(X_log_tst, on = "months")
rebal_indices

R>>  [1]  504  523  542  564  585  605  627  649  670  691  713  733  755  774  794  816  837  858  880
R>> [20]  900  923  944  965  986 1007

index(X_log)[rebal_indices]

R>>  [1] "2015-01-02" "2015-01-30" "2015-02-27" "2015-03-31" "2015-04-30" "2015-05-29" "2015-06-30"
R>>  [8] "2015-07-31" "2015-08-31" "2015-09-30" "2015-10-30" "2015-11-30" "2015-12-31" "2016-01-29"
R>> [15] "2016-02-29" "2016-03-31" "2016-04-29" "2016-05-31" "2016-06-30" "2016-07-29" "2016-08-31"
R>> [22] "2016-09-30" "2016-10-31" "2016-11-30" "2016-12-30"

# run the rolling window loop
lookback <- 10*20  # maximum value is: floor(T_trn/20)*20
for (i in 1:length(rebal_indices)) {
  # estimate moments
  X_ <- X_log[(rebal_indices[i] - lookback + 1):rebal_indices[i], ]
  mu <- colMeans(X_)
  Sigma <- cov(X_)
  # design portfolio
  w_MVP_rolling[rebal_indices[i], ] <- MVP(mu, Sigma, lmd = 2)
}
w_MVP_rolling <- na.omit(w_MVP_rolling)
w_MVP_rolling

R>>                    AAPL          AMD          ADI         ABBV         AEZS            A
R>> 2015-01-02 1.000000e+00 1.580691e-10 5.050330e-10 1.202175e-09 1.105367e-10 3.573050e-10
R>> 2015-01-30 1.000000e+00 3.539053e-10 7.050078e-10 2.239440e-09 2.554891e-10 6.893228e-10
R>> 2015-02-27 9.999998e-01 1.590879e-09 5.307164e-09 5.757445e-09 1.121228e-09 3.135692e-09
R>> 2015-03-31 7.476072e-01 2.305288e-09 1.534448e-08 1.164569e-08 1.782426e-09 6.584706e-09
R>> 2015-04-30 6.529623e-01 1.929058e-09 2.915103e-08 1.511416e-02 1.810234e-09 9.038248e-09
R>> 2015-05-29 2.315103e-01 6.080842e-10 2.774413e-01 4.580175e-08 3.180671e-10 2.338672e-09
R>> 2015-06-30 4.207176e-08 1.175655e-09 3.356819e-01 1.311964e-08 4.809086e-10 2.857770e-09
R>> 2015-07-31 8.392660e-08 3.044102e-09 5.745496e-01 2.863960e-01 8.534075e-10 1.274881e-08
R>> 2015-08-31 1.592957e-08 1.396224e-09 4.926272e-01 1.703659e-08 3.254074e-10 3.791135e-09
R>> 2015-09-30 5.014320e-01 1.780290e-09 3.394764e-01 6.258564e-09 3.285026e-10 6.633887e-09
R>> 2015-10-30 5.959166e-01 4.504492e-09 4.040833e-01 2.368998e-08 3.681050e-10 2.538993e-08
R>> 2015-11-30 1.573484e-08 3.875573e-09 2.973573e-01 6.790072e-02 9.198252e-10 6.347419e-01
R>> 2015-12-31 7.729590e-10 7.191776e-02 3.307433e-09 2.178741e-01 4.398168e-11 7.102081e-01
R>> 2016-01-29 1.447055e-08 1.663432e-08 4.935063e-02 2.961420e-01 4.236551e-10 4.289580e-01
R>> 2016-02-29 5.134762e-08 1.863314e-07 1.944425e-07 3.292062e-06 2.568049e-09 2.192275e-01
R>> 2016-03-31 2.747988e-08 1.267846e-01 3.120421e-08 3.602898e-08 1.447580e-09 2.777190e-01
R>> 2016-04-29 6.306332e-10 2.687206e-01 1.237174e-09 1.758966e-09 1.033592e-10 1.059529e-08
R>> 2016-05-31 4.182577e-09 4.299935e-01 7.054894e-09 5.718436e-09 5.535112e-10 5.700065e-01
R>> 2016-06-30 9.988696e-10 4.381333e-01 1.585390e-09 4.372585e-09 3.829997e-10 5.618667e-01
R>> 2016-07-29 2.242834e-09 4.987464e-01 4.403256e-09 4.133897e-02 5.280437e-10 4.599146e-01
R>> 2016-08-31 8.031004e-10 5.430246e-01 1.242330e-09 2.333549e-09 1.985657e-10 4.569754e-01
R>> 2016-09-30 1.911442e-08 4.377899e-01 4.943627e-08 1.596992e-01 2.734715e-09 6.312447e-08
R>> 2016-10-31 2.248867e-07 5.312681e-01 2.115053e-01 2.539852e-08 1.146040e-06 5.974307e-08
R>> 2016-11-30 1.266082e-08 6.462092e-01 3.537907e-01 1.425341e-08 5.633842e-09 1.241059e-08
R>> 2016-12-30 1.215596e-08 6.751752e-01 2.282577e-01 9.656699e-02 3.525486e-09 3.534123e-08
R>>                     APD           AA           CF
R>> 2015-01-02 1.037021e-09 1.582697e-09 8.306670e-10
R>> 2015-01-30 2.004292e-09 1.293498e-09 5.901921e-09
R>> 2015-02-27 1.693633e-08 3.335767e-09 1.733711e-07
R>> 2015-03-31 3.507708e-07 4.284927e-09 2.523924e-01
R>> 2015-04-30 1.669919e-08 4.136112e-09 3.319235e-01
R>> 2015-05-29 3.046181e-09 9.222353e-10 4.910484e-01
R>> 2015-06-30 5.425461e-09 1.349488e-09 6.643180e-01
R>> 2015-07-31 3.138748e-08 2.726832e-09 1.390543e-01
R>> 2015-08-31 6.386453e-02 9.593828e-10 4.435083e-01
R>> 2015-09-30 1.020273e-01 1.448415e-09 5.706429e-02
R>> 2015-10-30 5.743667e-08 1.693634e-09 7.838458e-09
R>> 2015-11-30 1.291114e-08 2.487516e-09 5.077051e-09
R>> 2015-12-31 1.205238e-09 4.518472e-10 4.028364e-10
R>> 2016-01-29 2.255493e-01 2.856509e-09 2.805078e-09
R>> 2016-02-29 7.807688e-01 1.891941e-08 1.650022e-08
R>> 2016-03-31 5.954963e-01 1.086565e-08 4.094872e-09
R>> 2016-04-29 7.312794e-01 2.159469e-09 2.712549e-10
R>> 2016-05-31 1.034081e-08 4.375698e-09 1.199165e-09
R>> 2016-06-30 4.289797e-09 9.780217e-10 3.301711e-10
R>> 2016-07-29 4.645468e-09 1.922269e-09 6.623378e-10
R>> 2016-08-31 1.264529e-08 1.402590e-09 2.218040e-10
R>> 2016-09-30 4.025108e-01 9.409569e-09 2.342977e-09
R>> 2016-10-31 2.572251e-01 4.424621e-08 9.895388e-09
R>> 2016-11-30 1.963114e-08 2.841839e-08 5.565432e-09
R>> 2016-12-30 4.998978e-08 1.334315e-08 3.895068e-09

We can now plot the evolution over time of the portfolio:

# compute portfolio returns
tmp <- Return.portfolio(X_lin_tst, weights = w_MVP_rolling, verbose = TRUE)
ret_rolling <- tmp$returns
chart.StackedBar(tmp$BOP.Weight, main = "Evolution of rolling-window Markowitz's mean-variance portfolio",
                 ylab = "w", space = 0, border = NA)

We can now compare the static and rolling window versions:

# performance
ret_MVP <- cbind(ret_static, ret_rolling)
colnames(ret_MVP) <- c("Markowitz MVP - static", "Markowitz MVP - rolling window")

t(table.AnnualizedReturns(ret_MVP))

R>>                                Annualized Return Annualized Std Dev Annualized Sharpe (Rf=0%)
R>> Markowitz MVP - static                    0.0350             0.1978                    0.1770
R>> Markowitz MVP - rolling window            0.3586             0.2957                    1.2126

chart.CumReturns(ret_MVP, main = "Daily wealth for static vs rolling window",
                 wealth.index = TRUE, legend.loc = "topleft", colorset = rich6equal)

chart.Drawdown(ret_MVP, main = "Drawdown for static vs rolling window", 
               legend.loc = "bottomleft", colorset = rich6equal)

• We have backtested a number of different portfolio designs.

• However, we have only used one set of market data and we should not draw any conclusion from that.

• Also, we have used a training data window to design the portfolio and evaluated it in a subsequent test data window.

• A serious backtesting requires:

  • multiple datasets
  • evaluation on a rolling-window basis.

• We will use the R package portfolioBacktest which makes this very simple.

First, let’s load the market data using the package portfolioBacktest:

library(portfolioBacktest)

# load the SP500 assets
data("SP500_symbols")
SP500_YAHOO <- stockDataDownload(stock_symbols = SP500_symbols, from = "2010-12-01", to = "2018-12-01")

# generate 100 random samples each containing 50 random assets over a random window of two years
N_datasets <- 100
mydataset <- stockDataResample(SP500_YAHOO, N = 50, T = 252*2, num_datasets = N_datasets)

Second, let’s define the portfolios that we want to backtest:

library(CVXR)

QuintP <- function(dataset) {
  prices <- dataset$adjusted
  N <- ncol(prices)
  X <- diff(log(prices))[-1]  # returns
  mu <- colMeans(X)
  Sigma <- cov(X)
  
  # portfolio design
  idx <- sort(mu/diag(Sigma), decreasing = TRUE, index.return = TRUE)$ix
  w <- rep(0, N)
  w[idx[1:round(N/5)]] <- 1/round(N/5)
  return(w)
}

IVP <- function(dataset) {
  prices <- dataset$adjusted
  N <- ncol(prices)
  X <- diff(log(prices))[-1]  # returns
  Sigma <- cov(X)

  # portfolio design
  sigma <- sqrt(diag(Sigma))
  w <- 1/sigma
  w <- w/sum(w)
  return(w)
}

GMVP <- function(dataset) {
  prices <- dataset$adjusted
  N <- ncol(prices)
  X <- diff(log(prices))[-1]  # returns
  Sigma <- cov(X)
  
  # portfolio design
  w <- Variable(N)
  prob <- Problem(Minimize(quad_form(w, Sigma)),
                  constraints = list(sum(w) == 1, w >= 0))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

MVP <- function(dataset) {
  prices <- dataset$adjusted
  N <- ncol(prices)
  X <- diff(log(prices))[-1]  # returns
  mu <- colMeans(X)
  Sigma <- cov(X)
  
  # portfolio design
  lmd = 0.5
  w <- Variable(N)
  prob <- Problem(Maximize(t(mu) %*% w - lmd*quad_form(w, Sigma)),
                  constraints = list(sum(w) == 1, w >= 0))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

portfolio_list <- list("QuintP" = QuintP,
                       "IVP"    = IVP,
                       "GMVP"   = GMVP,
                       "MVP"    = MVP)

Then, let’s proceed with the backtesting based on 100 datasets randomly chosen from market data:

bt_all_port <- portfolioBacktest(portfolio_funs = portfolio_list,
                                 dataset = mydataset,
                                 benchmark = c("uniform", "index"),
                                 T_rolling_window = 252*2/3, optimize_every = 20, rebalance_every = 1, 
                                 show_progress_bar = FALSE,
                                 paral_datasets = 5)
res_summary_median <- backtestSummary(bt_all_port)

Now we are ready to compare the portfolios via tables and plots.

A leaderboard is a table with a comparison:

summaryTable(res_summary_median, type = "DT", order_col = 2, order_dir = "desc")

A barplot provides the same information as the leaderboard but in a visual way:

summaryBarPlot(res_summary_median, measures = c("Sharpe ratio", "max drawdown"))

A boxplot is probably the best way to properly compare the performance of different portfolios with a single performance measure:

backtestBoxPlot(bt_all_port, "Sharpe Ratio")

We can confirm the good performance of the EWP and the bad performance of the Markowitz MVP.