spIndexTrack(): Sparse portfolio construction

Assume that an index is composed of \(N\) assets. We denote by \(\mathbf{r}^b=[r_1^b,\dots,r_T^b]^\top\in\mathbb{R}^T\) and \(\mathbf{X}=[\mathbf{r}_1,\dots,\mathbf{r}_T]^\top\in\mathbb{R}^{T\times N}\) the (arithmetic) net returns of the index and the \(N\) assets in the past \(T\) days, respectively, with \(\mathbf{r}_t\in\mathbb{R}^N\) denoting the net returns of the \(N\) assets at the \(t\)-th day.

The goal of spIndexTrack() is the design of a (sparse) portfolio \(\mathbf{w}\in\mathbb{R}_+^N\), with \(\mathbf{w}^\top\mathbf{1} = 1\), that tracks closely the index, i.e., \(\mathbf{X}\mathbf{w} \approx \mathbf{r}^b\), based on [1]. The underlying optimization problem that is solved is

\[\begin{equation} \begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \text{TE}(\mathbf{w}) + \lambda\|\mathbf{w}\|_0\\ \textsf{subject to} & \mathbf{w}^\top\mathbf{1}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array} \tag{1} \end{equation}\]

where \(\text{TE}(\mathbf{w})\) is a general tracking error (we will see specific tracking errors shortly), \(\lambda\) is a regularization parameter that controls the sparsity of the portfolio, and \(u\) is an upper bound on the weights of the portfolio.

The \(\ell_0\)-“norm” is approximated by the continuous and differentiable (for \(w \geq 0\)) function

\[\begin{equation} \rho_{p,u}(w) = \frac{\log(1 + w/p)}{\log(1 + u/p)}, \end{equation}\]

where \(p>0\) is a parameter that controls the approximation. This leads to the following approximate problem:

\[\begin{equation} \begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \text{TE}(\mathbf{w}) + \lambda\mathbf{1}^\top\boldsymbol{\rho}_{p,u}(\mathbf{w})\\ \textsf{subject to} & \mathbf{w}^\top\mathbf{1}=1\\ & \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}, \end{array} \tag{2} \end{equation}\]

where \(\boldsymbol{\rho}_{p,u}(\mathbf{w})=[\mathbf{\rho}_{p,u}(w_1),\dots,\rho_{p,u}(w_N)]^\top\).

There are four available tracking errors \(\text{TE}(\mathbf{w}\)) in spIndexTrack():

• Empirical tracking error (ETE):

\[ \text{ETE}(\mathbf{w}) = \frac{1}{T}\big\|\mathbf{r}^b - \mathbf{X}\mathbf{w}\big\|_2^2 \]

• Downside risk (DR):

\[ \text{DR}(\mathbf{w}) = \frac{1}{T}\big\|(\mathbf{r}^b-\mathbf{X}\mathbf{w})^+\big\|_2^2 \]

• Huber empirical tracking error (HETE):

\[ \text{HETE}(\mathbf{w}) = \frac{1}{T}\mathbf{1}^\top\boldsymbol{\phi}\left(\mathbf{r}^b - \mathbf{X}\mathbf{w}\right) \]

• Huber downside risk (HDR): \[ \text{HDR}(\mathbf{w}) = \frac{1}{T}\mathbf{1}^\top\boldsymbol{\phi}\left((\mathbf{r}^b-\mathbf{X}\mathbf{w})^+\right) \] where \(\boldsymbol{\phi}(\mathbf{x}) = [\phi(x_1), \dots, \phi(x_T)]^\top\) and \[ \phi(x) = \begin{cases} x^2 &\quad |x| \leq M\\ M(2|x| - M) &\quad |x| > M, \end{cases} \] with \(M>0\) being the Huber parameter.

Regardless of the selected tracking error measure, problem (2) can be solved via Majorization-Minimization (MM) [3] with an iterative closed-form update algorithm (with iterations denoted by \(k\)). It can be shown that all of the above variations boil down to the iterative optimization of the following convex problem:

\[\begin{equation} \begin{array}{ll} \underset{\mathbf{w}}{\text{minimize}} & \mathbf{w}^\top\mathbf{w} + {\mathbf{q}^{(k)}}^\top\mathbf{w}\\ \textsf{subject to} & \mathbf{w}\in\mathcal{W}_{u}, \end{array} \tag{3} \end{equation}\]

where \[ \mathcal{W}_{u} = \big\{\mathbf{w} \big| \mathbf{w}^\top\mathbf{1} = 1, \mathbf{0}\leq\mathbf{w}\leq u\mathbf{1}\big\}, \] and \(\mathbf{q}^{(k)}\in\mathbb{R}^N\).

What differentiates the various tracking errors is the exact form of \(\mathbf{q}^{(k)}\) that we need to compute at each iteration \(k\) of the algorithm: \[ \begin{aligned} \mathbf{q}_{\text{ETE}}^{(k)} & = \frac{1}{\lambda_{\text{max}}^{(\mathbf{L}_1)}}(2(\mathbf{L}_1 - \lambda_{\text{max}}^{(\mathbf{L}_1)}\mathbf{I})\mathbf{w}^{(k)} + \lambda{\mathbf{d}_{p,u}^{(k)}} -\frac{2}{T}\mathbf{X}^\top\mathbf{r}^b),\\ \mathbf{q}_{\text{DR}}^{(k)} & = \frac{1}{\lambda_{\text{max}}^{(\mathbf{L}_1)}} (\frac{2}{T} 2(\mathbf{L}_1 - \lambda_{\text{max}}^{(\mathbf{L}_1)}\mathbf{I})\mathbf{w}^{(k)} + \lambda\mathbf{d}_{p,u}^{(k)} + \frac{2}{T}\mathbf{X}^\top(\mathbf{y}^{(k)} - \mathbf{r}^b)),\\ \mathbf{q}_{\text{HETE}}^{(k)} & = \frac{1}{\lambda_{\text{max}}^{(\mathbf{L}_2)}}(2(\mathbf{L}_2 - \lambda_{\text{max}}^{(\mathbf{L}_2)}\mathbf{I})\mathbf{w}^{(k)} + \lambda{\mathbf{d}_{p,u}^{(k)}} -\frac{2}{T}\mathbf{X}^\top\text{Diag}(\mathbf{a}^{(k)})\mathbf{r}^b),\\ \mathbf{q}_{\text{HDR}}^{(k)} & = \frac{1}{\lambda_{\text{max}}^{(\mathbf{L}_3)}}(2(\mathbf{L}_3 - \lambda_{\text{max}}^{(\mathbf{L}_3)}\mathbf{I})\mathbf{w}^{(k)} + \lambda{\mathbf{d}_{p,u}^{(k)}} +\frac{2}{T}\mathbf{X}^\top\text{Diag}(\mathbf{b}^{(k)})(\mathbf{c}^{(k)} - \mathbf{r}^b)), \end{aligned} \] where \(\lambda_{\text{max}}^{(\mathbf{A})}\) denotes the maximum eigenvalue of a matrix \(\mathbf{A}\), \(\mathbf{I}\) denotes the identity matrix, \(\text{Diag}(\mathbf{x})\) is a diagonal matrix with the vector \(\mathbf{x}\) at its principal diagonal, and
\[\begin{align} \mathbf{d}_{p,u}^{(k)} & = \left[d_{p,u}(w_1^{(k)}),\dots,d_{p,u}(w_N^{(k)})\right]^\top,\\ d_{p,u}(w^{(k)}) & = \frac{1}{\log(1 + u/p)(p+w^{(k)})},\\ \mathbf{y}^{(k)} & = -(\mathbf{X}\mathbf{w}^{(k)} - \mathbf{r}^b)^+,\\ \mathbf{a}^{(k)} & = [a([\mathbf{r}^b - \mathbf{X}\mathbf{w}^{(k)}]_1),\dots,a([\mathbf{r}^b - \mathbf{X}\mathbf{w}^{(k)}]_T)]^\top,\\ a(x) & = \begin{cases} 1 &\quad |x|\leq M\\ \frac{M}{|x|} &\quad |x| > M,\end{cases}\\ \mathbf{b}^{(k)} & = [b([\mathbf{r}^b - \mathbf{X}\mathbf{w}^{(k)}]_1),\dots,b([\mathbf{r}^b - \mathbf{X}\mathbf{w}^{(k)}]_T)]^\top,\\ b(x) & = \begin{cases} \frac{M}{M - 2x} &\quad x < 0\\ 1 &\quad0\leq x\leq M\\ \frac{M}{x} &\quad x > M,\end{cases}\\ \mathbf{c}^{(k)} & = [c([\mathbf{r}^b - \mathbf{X}\mathbf{w}^{(k)}]_1),\dots,c([\mathbf{r}^b - \mathbf{X}\mathbf{w}^{(k)}]_T)]^\top,\\ c(x) & = \begin{cases} x &\quad x < 0\\ 0 &\quad x \geq 0,\end{cases}\\ \mathbf{L}_1 & = \frac{1}{T}\mathbf{X}^\top\mathbf{X},\\ \mathbf{L}_2 & = \frac{1}{T}\mathbf{X}^\top\text{Diag}(\mathbf{a}^{(k)})\mathbf{X},\\ \mathbf{L}_3 & = \frac{1}{T}\mathbf{X}^\top\text{Diag}(\mathbf{b}^{(k)})\mathbf{X}. \end{align}\]

The following propositions provide a waterfilling structured solution of problem (3), considering two special cases, namely, \(u=1\) and \(u<1\).

We refer to the iterative procedure of Proposition 1 as AS\(_{1}(\mathbf{q})\) (Active-Set for \(u=1\)) and of Proposition 2 as AS\(_u(\mathbf{q})\) (Active-Set for general \(u<1\)). The iterative closed-form update algorithm is given in Algorithm 1 (where AS\(_{1|u}(\mathbf{q})\) means AS\(_1(\mathbf{q})\) or AS\(_u(\mathbf{q})\)).

Algorithm 1
1. Set \(k=0\) and choose an initial point \(\mathbf{w}^{(0)}\) (by default set to \(\mathbf{w}^{(0)} = \frac{1}{N}\mathbf{1}\))
2. Compute \(\mathbf{q}\) according to the selected tracking error
3. Find the optimal solution \(\mathbf{w}^\star\) with AS\(_{1|u}(\mathbf{q})\) and set it equal to \(\mathbf{w}^{(k+1)}\)
4. \(k \gets k+1\)
5. Repeat steps 2-4 until convergence
6. Return \(\mathbf{w}^{(k)}\)

Finally, note that the approximate problem is controlled by the parameter \(p\), and in particular, as \(p\rightarrow0\) we get \(\rho_{p,u}\rightarrow\ell_0\). However, by setting small values to \(p\), it is likely that the algorithm will get stuck to a local minimum. To solve this issue we start with large values for \(p\), i.e., a “loose” approximation, and solve the corresponding optimization problem. Then, we sequentially decrease \(p\), i.e., we “tighten” the approximation, and solve the problem again using the previous solution as an initial point. In practice we are interested only in the last, “tightest” problem. For each problem that is solved (i.e., for fixed \(p\)) we utilize an acceleration scheme that increases the convergence speed of the MM algorithm. For details, please refer to [4].