## Portfolio Optimization Example

In portfolio optimization [1], [2], we are interested in finding an optimal portfolio $w\in \mathbb{R}^n$ of $n$ assets. The return $r$ has expected value $\alpha$ and covariance $\Sigma \in \mathbb{S}_{++}^n$. The risk aversion factor is denoted by $\gamma \geq 0$, and transaction (short-selling) cost is denoted by $\kappa_\mathrm{tc(sh)} \geq 0$. We solve the optimization

\begin{equation}
\begin{array}{ll}
\text{maximize} \quad &\alpha^T w - \gamma w^T \Sigma w - \kappa_\mathrm{tc}^T |w-w^\mathrm{prev}| + \kappa_\mathrm{sh}^T (w)_- \\
\text{subject to} \quad &\mathbb{1}^T w = 1 \\
&\Vert w \Vert_1 \leq L,
\end{array}
\end{equation}

where $w \in \mathbb{R}^n$ is the variable. The previous portfolio is $w^\mathrm{prev}$, $(\cdot)_-$ represents the argument's negative part, and $L \geq 1$ is the maximum value of total long positions plus the total magnitude of short positions.

With the factor model $\Sigma = F \Sigma^f F^T + D$ [3], we can rewrite an equivalent and [DPP-compliant](https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming) problem, i.e.,

\begin{equation}
\begin{array}{ll}
\text{maximize} \quad &a^T w - \Vert \left(\Sigma^f\right)^{1/2} f \Vert_2^2 - \Vert D^{1/2} w\Vert_2^2 - k_\mathrm{tc}^T |\Delta w| + k_\mathrm{sh}^T (w)_- \\
\text{subject to} &f = F^T w \\
&\mathbb{1}^T w = 1 \\
&\Vert w \Vert_1 \leq L \\
&\Delta w = w-w^\mathrm{prev}, \\
\end{array}
\end{equation}

where $w \in \mathbb{R}^n$ and $f \in \mathbb{R}^m$ are variables. To summarize, the parameters are:

\begin{equation}
\begin{array}{ll}
a = \frac{\alpha}{\gamma} \\
F \\
\left(\Sigma^f\right)^{1/2} \\
D^{1/2} \\
k_\mathrm{tc} = \frac{\kappa_\mathrm{tc}}{\gamma} \\
k_\mathrm{sh} = \frac{\kappa_\mathrm{sh}}{\gamma} \\
w^\mathrm{prev} \\
L.
\end{array}
\end{equation}

Note that we divided the objective function by the risk aversion factor $\gamma$. This way, updating the value of $\gamma$ only affects the linear part of the objective function, avoiding to compute a matrix factorization when solving the problem repeatedly.

Let's define the corresponding CVXPY problem. Note that we define the diagonal of $D^{1/2}$ as a vector.


```python
import cvxpy as cp
import numpy as np

# define dimensions
n, m = 100, 10

# define variables
w = cp.Variable(n, name='w')
delta_w = cp.Variable(n, name='delta_w')
f = cp.Variable(m, name='f')

# define parameters
a = cp.Parameter(n, name='a')
F = cp.Parameter((n, m), name='F')
Sig_f_sqrt = cp.Parameter((m, m), name='Sig_f_sqrt')
d_sqrt = cp.Parameter(n, name='d_sqrt')
k_tc = cp.Parameter(n, nonneg=True, name='k_tc')
k_sh = cp.Parameter(n, nonneg=True, name='k_sh')
w_prev = cp.Parameter(n, name='w_prev')
L = cp.Parameter(nonneg=True, name='L')

# define objective
objective = cp.Maximize(a@w
                        -cp.sum_squares(Sig_f_sqrt@f)
                        -cp.sum_squares(cp.multiply(d_sqrt, w))
                        -k_tc@cp.abs(delta_w)
                        +k_sh@cp.minimum(0, w))

# define constraints
constraints = [f == F.T@w,
               np.ones(n)@w == 1, 
               cp.norm(w, 1) <= L, 
               delta_w == w-w_prev]

# define problem
problem = cp.Problem(objective, constraints)
```

Assign parameter values and solve the problem.


```python
np.random.seed(0)
gamma = 1
alpha = np.random.randn(n)
kappa_tc = 0.01*np.ones(n)
kappa_sh = 0.05*np.ones(n)

a.value = alpha/gamma
F.value = np.random.randn(n, m)
Sig_f_sqrt.value = np.random.rand(m, m)
d_sqrt.value = np.random.rand(n)
k_tc.value = kappa_tc/gamma
k_sh.value = kappa_sh/gamma
w_prev_unnormalized = np.random.rand(n)
w_prev.value = w_prev_unnormalized/np.linalg.norm(w_prev_unnormalized)
L.value = 1.6

val = problem.solve()
```

Generating C source for the problem is as easy as:


```python
from cvxpygen import cpg

cpg.generate_code(problem, code_dir='portfolio_code')
```

Now, you can use a python wrapper around the generated code as a custom CVXPY solve method.


```python
from portfolio_code.cpg_solver import cpg_solve
import numpy as np
import pickle
import time

# load the serialized problem formulation
with open('portfolio_code/problem.pickle', 'rb') as f:
    prob = pickle.load(f)

# assign parameter values
n, m = 100, 10
np.random.seed(1)
gamma = 1
alpha = np.random.randn(n)
kappa_tc = 0.01*np.ones(n)
kappa_sh = 0.05*np.ones(n)
prob.param_dict['a'].value = alpha/gamma
prob.param_dict['F'].value = np.round(np.random.randn(n, m))
prob.param_dict['Sig_f_sqrt'].value = np.diag(np.random.rand(m))
prob.param_dict['d_sqrt'].value = np.random.rand(n)
prob.param_dict['k_tc'].value = kappa_tc/gamma
prob.param_dict['k_sh'].value = kappa_sh/gamma
prob.param_dict['w_prev'].value = np.zeros(n)
prob.param_dict['L'].value = 1.6

# solve problem conventionally
t0 = time.time()
# CVXPY chooses eps_abs=eps_rel=1e-5, max_iter=10000, polish=True by default,
# however, we choose the OSQP default values here, as they are used for code generation as well
val = prob.solve(eps_abs=1e-3, eps_rel=1e-3, max_iter=4000, polish=False)
t1 = time.time()
print('\nCVXPY\nSolve time: %.3f ms' % (1000 * (t1 - t0)))
print('Objective function value: %.6f\n' % val)

# solve problem with C code via python wrapper
prob.register_solve('CPG', cpg_solve)
t0 = time.time()
val = prob.solve(method='CPG')
t1 = time.time()
print('\nCVXPYgen\nSolve time: %.3f ms' % (1000 * (t1 - t0)))
print('Objective function value: %.6f\n' % val)
```

### References

[1] Lobo, Miguel Sousa, Maryam Fazel, and Stephen Boyd. "Portfolio optimization with linear and fixed transaction costs." Annals of Operations Research 152.1 (2007): 341-365.

[2] Moehle, N., Kochenderfer, M. J., Boyd, S., and Ang, A. "Tax-Aware Portfolio Construction via Convex Optimization." Journal of Optimization Theory and Applications 189.2 (2021): 364-383.

[3] Ng, Victor, Robert F. Engle, and Michael Rothschild. "A multi-dynamic-factor model for stock returns." Journal of Econometrics 52.1-2 (1992): 245-266.
