A simple design for sparse index tracking is to decompose the task into two steps: asset selection (choosing the active weights) and capital allocation (optimization of weights). We will assume we know the portfolio implemented by the actual index, which will be a dense portfolio (as opposed to sparse) \(\mathbf{b}\) such that \(\mathbf{X}\mathbf{b} = \mathbf{r}^b\). In practice, it is expensive to get the information on \(\mathbf{b}\), so we will just regress it.

• Step 1: Asset selection We can define the selection of \(K\) assets through the vector \(\mathbf{s}\), defined as \[s_i=\begin{cases} 1\quad & \text{if the }j\text{-th asset is selected}\\ 0 & \text{otherwise}, \end{cases}\] with \(\mathbf{1}^T\mathbf{s}=K\). For example, we can simply choose the \(K\) largest assets of \(\mathbf{b}\).

• Step 2: Capital allocation Once the active assets are known, then the portfolio design is simple. One option is to define the portfolio in the same proportion as the index composition (but only using the active assets): \[\mathbf{w}=\frac{\mathbf{b}\odot\mathbf{s}}{\mathbf{1}^T(\mathbf{b}\odot\mathbf{s})}\] where \(\odot\) denotes Hadamard (elementwise) product. Another option is to fit the index returns subject to the sparsity pattern given by the asset selection: \[\begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \frac{1}{T}\big\|\mathbf{r}^b - \mathbf{X}(\mathbf{w}\odot\mathbf{s})\big\|_2^2\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{w}\ge\mathbf{0}. \end{array}\]

Loading data
We start by loading the market data of the index S&P 500 and its underlying assets during 2010 (from the package sparseIndexTracking):

library(xts)
library(sparseIndexTracking)  #install.packages("sparseIndexTracking")
data(INDEX_2010)
T <- nrow(INDEX_2010$X)
N <- ncol(INDEX_2010$X)
X <- INDEX_2010$X
r <- INDEX_2010$SP500

# training and test sets
T_trn <- round(0.7*T)  # define the training set
T_tst <- T - T_trn
X_trn <- X[1:T_trn, ]
r_trn <- r[1:T_trn, ]
X_tst <- X[-(1:T_trn), ]
r_tst <- r[-(1:T_trn), ]

The data INDEX_2010 contains a list with two xts objects:
1. X: A \(T\times N\) xts with the daily linear returns of the \(N\) assets that were in the index during the year 2010 (total \(T\) trading days)
2. SP500: A \(T\times 1\) xts with the daily linear returns of the index S&P 500 during the same period.

Regressing the dense index portfolio
Since we don’t know the index composition \(\mathbf{b}\), we will get an estimate by solving the constrained regression \[\begin{array}{ll} \underset{\mathbf{b}}{\text{minimize}} & \frac{1}{T}\big\|\mathbf{r}^b - \mathbf{X}\mathbf{b}\big\|_2^2\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{b}=1\\ & \mathbf{b}\ge\mathbf{0}. \end{array}\]

library(CVXR)
b <- Variable(N)
prob <- Problem(Minimize(sum((as.matrix(r) - as.matrix(X) %*% b)^2)), 
                constraints = list(b >= 0, sum(b)==1))
result <- solve(prob)
b_SP500 <- as.vector(result$getValue(b))
barplot(b_SP500, main = "SP500 composition", xlab = "asset index")

Two-stage approach: weights proportional to \(\mathbf{b}\)
To obtain \(\mathbf{s}\), we can simply choose the largest \(K\) elements of \(\mathbf{b}\):

K <- 20
idx_Klargest <- sort(b_SP500, decreasing = TRUE, index.return = TRUE)$ix[1:K]
s <- rep(0, N)
s[idx_Klargest] <- 1

Then we can proceed to define the portfolio in the same proportion as the index composition:

w_2step_a <- rep(0, N)
w_2step_a[idx_Klargest] <- b_SP500[idx_Klargest]
w_2step_a <- w_2step_a/sum(w_2step_a)
barplot(w_2step_a, main = "Two-step design with proportional weights", xlab = "asset index")

Two-stage approach: fitted weights
We can also obtain the least squares design given the previous selection:

w <- Variable(N)
prob <- Problem(Minimize(sum((as.matrix(r_trn) - as.matrix(X_trn) %*% (w*s))^2)), 
                constraints = list(w >= 0, sum(w)==1))
result <- solve(prob)
w_2step_b <- as.vector(result$getValue(w))
barplot(w_2step_b, main = "Two-step design with fitted weights", xlab = "asset index")

Comparison We can now compute the tracking errors of the two designs compared with the index value:

# tracking error
ETE <- c("dense b"                = (1/T)*norm(r - X %*% b_SP500, "2")^2,
         "two-step sparse prop b" = (1/T)*norm(r - X %*% w_2step_a, "2")^2,
         "two-step sparse LS"     = (1/T)*norm(r - X %*% w_2step_b, "2")^2)
ETE <- as.matrix(ETE)
colnames(ETE) <- "ETE"
print(ETE)

R>>                                 ETE
R>> dense b                8.753776e-09
R>> two-step sparse prop b 5.949458e-06
R>> two-step sparse LS     4.335993e-06

We can get a better graph of the tracking errors visually:

barplot(c(ETE), main = "Tracking error", ylab = "ETE", col = rainbow(4))

We can also plot the Cum-PnL:

library(PerformanceAnalytics)

# plot of wealth
ret_all <- cbind("SP500"                  = r,
                 "dense b"                = X %*% b_SP500,
                 "two-step sparse prop b" = X %*% w_2step_a,
                 "two-step sparse LS"     = X %*% w_2step_b)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

Recall our nonconvex and nondifferentiable optimization problem: \[\begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \frac{1}{T}\big\|\mathbf{r}^b - \mathbf{X}\mathbf{w}\big\|_2^2 + \lambda\|\mathbf{w}\|_0\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}. \end{array}\]

The idea now is to employ the Majorization-Minimization (MM) algorithm to iteratively solve a sequence of simpler surrogate problems.

The first step is to get rid of the nondifferentiability by approximating the cardinality operator \(\|\mathbf{w}\|_0\) by \(\mathbf{1}^{\top}\mathbf{\boldsymbol{\rho}}{}_{p,\gamma}(\mathbf{w})\), where we have approximated the \(\ell_{0}\)-norm (indicator function) with \[\rho_{p,\gamma}(w)=\frac{\log(1+|w|/p)}{\log(1+\gamma/p)}.\]

Now we can derive a majorizer of the nonconvex term \(\mathbf{1}^{\top}\mathbf{\boldsymbol{\rho}}{}_{p,\gamma}(\mathbf{w})\) (i.e., derive a tangent upper bound) easily by simply linearizing the log function. At iteration \(k\), the resulting problem is the convex QP: \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\big\|\mathbf{X}\mathbf{w}-\mathbf{r}^{b}\big\|_{2}^{2}+\lambda{\mathbf{d}_{p,u}^{(k)}}^\top\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\] where \(d_{p,\gamma}(w^{(k)})=\frac{1}{\log(1+\gamma/p)(p+w^{(k)})}\).

We are ready to implement the iterative MM approach:

# parameters
gamma <- 0.5; p <- 1e-3  # parameters for the rho approximation of log
lambda <- 6e-7  # this achieves a sparsity level of K=20
# MM loop
w <- rep(1/N, N)
obj_value <- (1/T_trn)*norm(r_trn - X_trn %*% w, "2")^2 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p))
while(1) {
  d <- 1/(log(1+gamma/p)*(p+w))
  w_prev <- w
  w <- Variable(N)
  prob <- Problem(Minimize((1/T_trn)*sum((as.matrix(r_trn) - as.matrix(X_trn) %*% w)^2) 
                           + lambda * t(d) %*% w), 
                  constraints = list(w >= 0, w <= 1, sum(w)==1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  obj_value <- c(obj_value, 
                 (1/T_trn)*norm(r_trn - X_trn %*% w, "2")^2 
                 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p)))
  if(norm(w-w_prev, "2")/norm(w_prev, "2") < 1e-3) break
}
w_MM_single <- w
plot(obj_value, type = "b", pch = 19, cex = .8, 
     main = "Convergence", xlab = "iteration", ylab = "objective value")

cat("Sparsity level K =", sum(w_MM_single > 1e-2))

R>> Sparsity level K = 20

barplot(w_MM_single, xlab = "asset index", 
        main = sprintf("Sparse portfolio via MM (sparsity level K=%d)", sum(w_MM_single > 1e-2)))

We can now compute the tracking error:

# tracking error
ETE <- rbind(ETE, 
             "MM-single" = (1/T)*norm(r - X %*% w_MM_single, "2")^2)
print(ETE)

R>>                                 ETE
R>> dense b                8.753776e-09
R>> two-step sparse prop b 5.949458e-06
R>> two-step sparse LS     4.335993e-06
R>> MM-single              1.895406e-06

barplot(c(ETE), main = "Tracking error", ylab = "ETE", col = rainbow(4))

As expected, the tracking error is much lower than the two-step approaches.

We can also plot the Cum-PnL:

# plot of wealth
ret_all <- cbind(ret_all, 
                 "MM-single" = X %*% w_MM_single)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

Recall that after the application of the first majorization, we obtained the following convex optimization problem: \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\big\|\mathbf{X}\mathbf{w}-\mathbf{r}^{b}\big\|_{2}^{2}+\lambda{\mathbf{d}_{p,u}^{(k)}}^\top\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\] which required the use of a solver at each iteration.

However, it is desirable to be able to solve the majorized problem in a simpler way. In our case, we can obtain a closed-form solution (we will focus on the simpler case \(u=1\)).

We first expand the objective: \[\frac{1}{T}\big\|\mathbf{X}\mathbf{w}-\mathbf{r}^{b}\big\|_{2}^{2}+\lambda{\mathbf{d}_{p,u}^{(k)}}^{\top}\mathbf{w}=\frac{1}{T}\mathbf{w}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{w}+\left(\lambda\mathbf{d}_{p,u}^{(k)}-\frac{2}{T}\mathbf{X}^{\top}\mathbf{r}^{b}\right)^{\!\!\top}\!\!\mathbf{w}+\mathsf{\mathit{const.}}\] and then we majorize a second time: \[\begin{aligned} \frac{1}{T}\big\|\mathbf{X}\mathbf{w}-\mathbf{r}^{b}\big\|_{2}^{2} &+ \lambda{\mathbf{d}_{p,u}^{(k)}}^{\top}\mathbf{w} \\ & \le \lambda_{\text{max}}^{(\frac{1}{T}\mathbf{X}^{\top}\mathbf{X})}\mathbf{w}^{\top}\mathbf{w}+\!\left(\!2\left(\frac{1}{T}\mathbf{X}^{\top}\mathbf{X}-\lambda_{\text{max}}^{(\frac{1}{T}\mathbf{X}^{\top}\mathbf{X})}\mathbf{I}\right)\mathbf{w}^{(k)}+\lambda{\mathbf{d}_{p,u}^{(k)}}-\frac{2}{T}\mathbf{X}^{\top}\mathbf{r}^{b}\!\right)^{\!\!\top}\mathbf{\!\!w}+const. \end{aligned}\]

Finally, the double-majorized problem can be written as \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \mathbf{w}^{\top}\mathbf{w}+{\mathbf{q}_{1}^{(k)}}^{\top}\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq \mathbf{1}, \end{array}\] where \[\mathbf{q}_{1}^{(k)}=\frac{1}{\lambda_{\text{max}}^{(\frac{1}{T}\mathbf{X}^{\top}\mathbf{X})}}\left(2\left(\frac{1}{T}\mathbf{X}^{\top}\mathbf{X}-\lambda_{\text{max}}^{(\frac{1}{T}\mathbf{X}^{\top}\mathbf{X})}\mathbf{I}\right)\mathbf{w}^{(k)}+\lambda{\mathbf{d}_{p,u}^{(k)}}-\frac{2}{T}\mathbf{X}^{\top}\mathbf{r}^{b}\right).\] This problem can be solved again with a solver but also directly with a simple algorithm in a finite number of iterations.

We can now implement this double iterative MM approach (see https://cvxr.rbind.io/cvxr_functions/ for the list of functions in CVXR):

# parameters
gamma <- 0.5; p <- 1e-3  # parameters for the rho approximation of log
lambda <- 5.7e-7  # this achieves a sparsity level of K=20  (5.5e-7 gives K=22, 5.8e-7 gives K=19)
# MM loop
w <- rep(1/N, N)
obj_value <- (1/T_trn)*norm(r_trn - X_trn %*% w, "2")^2 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p))
while(1) {
  d <- 1/(log(1+gamma/p)*(p+w))
  XX <- t(X_trn) %*% X_trn
  lmd_max <- max(eigen(XX/T_trn)$values)
  q1 <- 1/lmd_max * (2*(XX/T_trn - lmd_max*diag(N)) %*% w + lambda*d - (2/T_trn)*t(X_trn)%*%r_trn)
  w_prev <- w
  w <- Variable(N)
  prob <- Problem(Minimize(sum_squares(w) + t(q1) %*% w), 
                  constraints = list(w >= 0, w <= 1, sum(w)==1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  obj_value <- c(obj_value, 
                 (1/T_trn)*norm(r_trn - X_trn %*% w, "2")^2 
                 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p)))
  if(norm(w-w_prev, "2")/norm(w_prev, "2") < 1e-5) break
}
w_MM_double <- w

plot(obj_value, type = "b", pch = 19, cex = .5, 
     main = "Convergence", xlab = "iteration", ylab = "objective value")

cat("Sparsity level K =", sum(w_MM_double > 1e-2))

R>> Sparsity level K = 21

barplot(w_MM_double, xlab = "asset index", 
        main = sprintf("Sparse portfolio via MM (sparsity level K=%d)", sum(w_MM_double > 1e-2)))

We can now compute the tracking error:

# tracking error
ETE <- rbind(ETE, 
             "MM-double" = (1/T)*norm(r - X %*% w_MM_double, "2")^2)
print(ETE)

R>>                                 ETE
R>> dense b                8.753776e-09
R>> two-step sparse prop b 5.949458e-06
R>> two-step sparse LS     4.335993e-06
R>> MM-single              1.895406e-06
R>> MM-double              2.076981e-06

barplot(c(ETE), main = "Tracking error", ylab = "ETE", col = rainbow(5))

We can also plot the Cum-PnL:

# plot of wealth
ret_all <- cbind(ret_all, 
                 "MM-double" = X %*% w_MM_double)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

Now we will use the package sparseIndexTracking:

library(sparseIndexTracking)
w_pkg <- spIndexTrack(X_trn, r_trn, lambda = 8.3e-7, u = gamma)

Portfolio:

cat("Sparsity level K =", sum(w_pkg > 1e-2))

R>> Sparsity level K = 20

barplot(w_pkg, xlab = "asset index", names.arg = rep("", N),
        main = sprintf("Sparse portfolio via MM (sparsity level K=%d)", sum(w_pkg > 1e-2)))

Tracking error:

# tracking error
ETE <- rbind(ETE, 
             "MM-pkg" = (1/T)*norm(r - X %*% w_pkg, "2")^2)
print(ETE)

R>>                                 ETE
R>> dense b                8.753776e-09
R>> two-step sparse prop b 5.949458e-06
R>> two-step sparse LS     4.335993e-06
R>> MM-single              1.895406e-06
R>> MM-double              2.076981e-06
R>> MM-pkg                 1.993646e-06

barplot(c(ETE), main = "Tracking error", ylab = "ETE", col = rainbow(6))

We can observe that the single MM performs the best and the double MM is a bit worse because it didn’t have enough time to converge (convergence is very slow in this implementation because we are still using a solver instead of the efficient closed-form solution, with more time it would converge though). The package sparseIndexTracking is using the double MM but with the efficient closed-form solution, so convergence time is fast and the final solution is like the single MM.

In the original formulation, we used the empirical tracking error to measure the performance: \[\textsf{ETE}(\mathbf{w}) = \frac{1}{T}\big\|\mathbf{r}^b - \mathbf{X}\mathbf{w}\big\|_2^2.\] However, depending on the goals of the tracking, one could argue that if the performance of the portfolio is better than the index then it should not be penalized. For that purpose, we can use the downside risk (DR): \[\textsf{DR}(\mathbf{w})=\frac{1}{T}\big\|(\mathbf{r}^{b}-\mathbf{X}\mathbf{w})^{+}\big\|_{2}^{2},\] where \((x)^+=\text{max}(0,x)\).

The problem formulation is then \[\begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \frac{1}{T}\big\|(\mathbf{r}^{b}-\mathbf{X}\mathbf{w})^{+}\big\|_{2}^{2} + \lambda\|\mathbf{w}\|_0\\ \textsf{subject to} & \mathbf{w}^\top\mathbf{1}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\] which, after approximating and majorizing the cardinality term \(\|\mathbf{w}\|_0\) as before, becomes the convex problem \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\big\|(\mathbf{r}^{b}-\mathbf{X}\mathbf{w})^{+}\big\|_{2}^{2} + \lambda{\mathbf{d}_{p,u}^{(k)}}^\top\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\]

We can easily solve this convex problem with a solver (note that we could also use a double majorization like before to further simplify the problem):

# parameters
gamma <- 0.5; p <- 1e-3  # parameters for the rho approximation of log
lambda <- 1e-7  # this achieves a sparsity level of K=20
# MM loop
w <- rep(1/N, N)
obj_value <- (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w), "2")^2 +
             lambda * sum(log(1+abs(w)/p)/log(1+gamma/p))
while(1) {
  d <- 1/(log(1+gamma/p)*(p+w))
  w_prev <- w
  w <- Variable(N)
  prob <- Problem(Minimize((1/T_trn)*sum((pos(as.matrix(r_trn) - as.matrix(X_trn) %*% w))^2) 
                           + lambda * t(d) %*% w), 
                  constraints = list(w >= 0, w <= 1, sum(w)==1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  obj_value <- c(obj_value, 
                 (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w), "2")^2 
                 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p)))
  if(norm(w-w_prev, "2")/norm(w_prev, "2") < 1e-3) break
}
w_DR_MM <- w
plot(obj_value, type = "b", pch = 19, cex = .8, 
     main = "Convergence", xlab = "iteration", ylab = "objective value")

cat("Sparsity level K =", sum(w_DR_MM > 1e-2))

R>> Sparsity level K = 20

barplot(w_DR_MM, xlab = "asset index", 
        main = sprintf("Sparse portfolio via MM (sparsity level K=%d)", sum(w_DR_MM > 1e-2)))

We can observe the difference from the Cum-PnL:

# plot of wealth
ret_all <- cbind(ret_all, 
                 "DR-MM-single" = X %*% w_DR_MM)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

Let’s plot the out-of-sample Cum-PnL for a fair comparison:

chart.CumReturns(ret_all[-(1:T_trn), ], 
                 main = sprintf("Out-of-sample cumulative P&L (N=%d, K=%d)", N, K), 
                 wealth.index = TRUE, legend.loc = "topleft")

Recall the original problem formulation is then \[\begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \frac{1}{T}\big\|(\mathbf{r}^{b}-\mathbf{X}\mathbf{w})^{+}\big\|_{2}^{2} + \lambda\|\mathbf{w}\|_0\\ \textsf{subject to} & \mathbf{w}^\top\mathbf{1}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\] and the convex formulation after approximating and majorizing the cardinality term \(\|\mathbf{w}\|_0\): \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\big\|(\mathbf{r}^{b}-\mathbf{X}\mathbf{w})^{+}\big\|_{2}^{2} + \lambda{\mathbf{d}_{p,u}^{(k)}}^\top\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}. \end{array}\] This problem is convex and can be solved with a QP solver.

Nevertheless, we can further simplify the problem by majorizing the DR term as follows: \[\frac{1}{T}\big\|(\mathbf{r}^{b}-\mathbf{X}\mathbf{w})^{+}\big\|_{2}^{2} \le \frac{1}{T}\big\|\mathbf{r}^{b}-\mathbf{X}\mathbf{w}-\mathbf{y}^{(k)}\big\|_{2}^{2}, \] where \[\mathbf{y}^{(k)}=-\left(\mathbf{X}\mathbf{w}^{(k)}-\mathbf{r}^{b}\right)^{+}.\] Effectively, this majorization transforms the DR into a term that looks like the ETE! This means we can use all the methods developed before simply by substituting \(\mathbf{r}^{b}\) with \(\mathbf{r}^{b}-\mathbf{y}^{(k)}\). This is the resulting optimization problem: \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\big\|\mathbf{r}^{b}-\mathbf{X}\mathbf{w}-\mathbf{y}^{(k)}\big\|_{2}^{2} + \lambda{\mathbf{d}_{p,u}^{(k)}}^\top\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\]

# parameters
gamma <- 0.5; p <- 1e-3  # parameters for the rho approximation of log
lambda <- 0.807324e-7  # this achieves a sparsity level close to 20
# MM loop
w <- rep(1/N, N)
obj_value <- (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w), "2")^2 +
             lambda * sum(log(1+abs(w)/p)/log(1+gamma/p))
while(1) {
  d <- 1/(log(1+gamma/p)*(p+w))
  y <- -pmax(0, X_trn %*% w - r_trn)
  w_prev <- w
  w <- Variable(N)
  prob <- Problem(Minimize((1/T_trn)*sum((as.matrix(r_trn) - as.matrix(X_trn) %*% w - y)^2) 
                           + lambda * t(d) %*% w), 
                  constraints = list(w >= 0, w <= 1, sum(w)==1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  obj_value <- c(obj_value, 
                 (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w), "2")^2
                 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p)))
  if(norm(w-w_prev, "2")/norm(w_prev, "2") < 1e-3) break
}
w_DR_MM_bis <- w
plot(obj_value, type = "b", pch = 19, cex = .8, 
     main = "Convergence", xlab = "iteration", ylab = "objective value")

cat("Sparsity level K =", sum(w_DR_MM_bis > 1e-2))

R>> Sparsity level K = 19

There is some difference in the portfolio by majorizing only the cardinality term or both the cardinality and the downside risk tracking error term:

norm(w_DR_MM - w_DR_MM_bis, "2")

R>> [1] 0.2887955

However, both Cum-PnL look very similar:

# plot of wealth
ret_all <- cbind(ret_all, 
                 "DR-MM-single-bis" = X %*% w_DR_MM_bis)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

We can again use the package sparseIndexTracking:

w_DR_pkg <- spIndexTrack(X_trn, r_trn, lambda = 2.78e-7, measure = "dr")
cat("Sparsity level K =", sum(w_DR_pkg > 1e-2))

R>> Sparsity level K = 20

We can compare the DR:

# tracking error
DR_trn <- c("DR MM single"     = (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w_DR_MM), "2")^2,
            "DR MM single-bis" = (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w_DR_MM_bis), "2")^2,
            "DR-MM-pkg"        = (1/T_trn)*norm(pmax(0, r_trn - X_trn %*% w_DR_pkg), "2")^2)
print(DR_trn)

R>>     DR MM single DR MM single-bis        DR-MM-pkg 
R>>     1.247079e-07     1.457982e-07     2.466550e-07

DR <- c("DR MM single"     = (1/T)*norm(pmax(0, r - X %*% w_DR_MM), "2")^2,
        "DR MM single-bis" = (1/T)*norm(pmax(0, r - X %*% w_DR_MM_bis), "2")^2,
        "DR-MM-pkg"        = (1/T)*norm(pmax(0, r - X %*% w_DR_pkg), "2")^2)
print(DR)

R>>     DR MM single DR MM single-bis        DR-MM-pkg 
R>>     8.397384e-07     6.616777e-07     1.241276e-06

barplot(DR, main = "Tracking error", ylab = "DR", col = rainbow(4))

We can observe the difference from the Cum-PnL:

# plot of wealth
ret_all <- cbind(ret_all,
                 "DR-MM-pkg" = X %*% w_DR_pkg)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

In the original formulation, we used the empirical tracking error to measure the performance: \[\textsf{ETE}(\mathbf{w}) = \frac{1}{T}\big\|\mathbf{r}^b - \mathbf{X}\mathbf{w}\big\|_2^2.\] However, the \(\ell_2\)-norm is very sensitive to outliers. A way to make it more robust is to use a more robust penalty function like the Huber penalty function: \[\phi(x)=\begin{cases} x^{2}, & \quad|x|\leq M,\\ M(2|x|-M), & \quad|x|>M. \end{cases}\] The robust tracking error measure becomes \[\textsf{ETE}_\textsf{huber}(\mathbf{w}) = \frac{1}{T}\mathbf{1}^\top\boldsymbol{\phi}(\mathbf{r}^b - \mathbf{X}\mathbf{w}).\]

The problem formulation after majorizing the cardinality term \(\|\mathbf{w}\|_0\) as before is then \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\mathbf{1}^\top\boldsymbol{\phi}(\mathbf{r}^b - \mathbf{X}\mathbf{w}) + \lambda\mathbf{d}_{p,u}^{(k)^\top}\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\]

We can easily solve this convex problem with a solver:

# parameters
gamma <- 0.5; p <- 1e-3  # parameters for the rho approximation of log
lambda <- 5e-7  # this achieves a sparsity level close to K=20
# MM loop
phi <- function(x, M = 0.1) {
  ifelse(abs(x) <= M, x^2, M*(2*abs(x)-M))
}
w <- rep(1/N, N)
obj_value <- (1/T_trn)*sum(phi(r_trn - X_trn %*% w)) + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p))
while(1) {
  d <- 1/(log(1+gamma/p)*(p+w))
  w_prev <- w
  w <- Variable(N)
  prob <- Problem(Minimize((1/T_trn)*sum(huber(as.matrix(r_trn) - as.matrix(X_trn) %*% w)) 
                           + lambda * t(d) %*% w), 
                  constraints = list(w >= 0, w <= 1, sum(w)==1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  obj_value <- c(obj_value, 
                 (1/T_trn)*sum(phi(r_trn - X_trn %*% w)) 
                 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p)))
  if(norm(w-w_prev, "2")/norm(w_prev, "2") < 1e-3) break
}
w_Huber_MM <- w
plot(obj_value, type = "b", pch = 19, cex = .8, 
     main = "Convergence", xlab = "iteration", ylab = "objective value")

cat("Sparsity level K =", sum(w_Huber_MM > 1e-2))

R>> Sparsity level K = 19

barplot(w_Huber_MM, xlab = "asset index", 
        main = sprintf("Sparse portfolio via MM (sparsity level K=%d)", sum(w_Huber_MM > 1e-2)))

We can observe the difference from the Cum-PnL:

# plot of wealth
ret_all <- cbind(ret_all[, 1:6], 
                 "Huber-MM-single" = X %*% w_Huber_MM)
{ chart.CumReturns(ret_all, main = sprintf("Cumulative P&L (N=%d, K=%d)", N, K), 
                   wealth.index = TRUE, legend.loc = "topleft")
  addEventLines(xts("training", index(X[T_trn])), srt=90, pos=2, lwd = 2, col = "black") }

Recall the problem formulation after majorizing the cardinality term \(\|\mathbf{w}\|_0\): \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\mathbf{1}^\top\boldsymbol{\phi}(\mathbf{r}^b - \mathbf{X}\mathbf{w}) + \lambda\mathbf{d}_{p,u}^{(k)^\top}\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\] This problem is convex and can be solved with a QP solver.

Nevertheless, we can further simplify the problem by majorizing the Huber term as follows: \[\frac{1}{T}\mathbf{1}^\top\boldsymbol{\phi}(\mathbf{r}^b - \mathbf{X}\mathbf{w}) \le \frac{1}{T}\mathbf{a}^{(k)^\top}(\mathbf{r}^b - \mathbf{X}\mathbf{w})^2, \] where \[a^{(k)}_i=\begin{cases} 1, & \quad |x^{(k)}_i|\leq M,\\ \frac{M}{|x^{(k)}_i|}, & \quad|x^{(k)}_i|>M. \end{cases}\]

Effectively, this majorization again transforms the Huber term into a term that looks like the ETE! This means we can use all the methods developed before. We can now solve the problem \[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \frac{1}{T}\mathbf{a}^{(k)^\top}(\mathbf{r}^b - \mathbf{X}\mathbf{w})^2 + \lambda\mathbf{d}_{p,u}^{(k)^\top}\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^\top\mathbf{w}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array}\]

# parameters
gamma <- 0.5; p <- 1e-3  # parameters for the rho approximation of log
lambda <- 5e-7  # this achieves a sparsity level of K=20
# MM loop
ak <- function(x, M = 0.1) {
  ifelse(abs(x) <= M, 1, M/abs(x))
}
w <- rep(1/N, N)
obj_value <- (1/T_trn)*sum(phi(r_trn - X_trn %*% w)) + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p))
while(1) {
  d <- 1/(log(1+gamma/p)*(p+w))
  a <- ak(r_trn - X_trn %*% w)
  w_prev <- w
  w <- Variable(N)
  prob <- Problem(Minimize((1/T_trn)* t(a) %*% (as.matrix(r_trn) - as.matrix(X_trn) %*% w)^2 
                           + lambda * t(d) %*% w), 
                  constraints = list(w >= 0, w <= 1, sum(w)==1))
  result <- solve(prob)
  w <- as.vector(result$getValue(w))
  obj_value <- c(obj_value, 
                 (1/T_trn)*sum(phi(r_trn - X_trn %*% w)) 
                 + lambda * sum(log(1+abs(w)/p)/log(1+gamma/p)))
  if(norm(w-w_prev, "2")/norm(w_prev, "2") < 1e-3) break
}
w_Huber_MM_bis <- w
cat("Sparsity level K =", sum(w_Huber_MM_bis > 1e-2))

R>> Sparsity level K = 19

plot(obj_value, type = "b", pch = 19, cex = .8, 
     main = "Convergence", xlab = "iteration", ylab = "objective value")

Of course we get exactly the same portfolio by majorizing only the cardinality term or both the cardinality and the tracking error term:

norm(w_Huber_MM - w_Huber_MM_bis, "2")

R>> [1] 4.235976e-05