Star | Watch | Fork |
Home |
---|
User Guide |
Examples |
About |
This robust portfolio management model is proposed by Zhu and Fukushima (2009). The portfolio allocation is determined via minimizing the worst-case conditional value-at-risk (CVaR) under ambiguous distribution information. The generic formulation is given as
\[\begin{align} \min~&\max\limits_{\pmb{\pi}\in \Pi} \alpha + \frac{1}{1-\beta}\pmb{\pi}^{\top}\pmb{u} &\\ \text{s.t.}~& u_k \geq -\pmb{y}_k^{\top}\pmb{x} - \alpha, &\forall k = 1, 2, ..., s \\ &u_k \geq 0, &\forall k=1, 2, ..., s \\ &\sum\limits_{k=1}^s\pi_k\pmb{y}_k^{\top}\pmb{x} \geq \mu, &\forall \pmb{\pi} \in \Pi \\ &\underline{\pmb{x}} \leq \pmb{x} \leq \overline{\pmb{x}} \\ &\pmb{1}^{\top}\pmb{x} = w_0 \end{align}\]with investment decisions \(\pmb{x}\in\mathbb{R}^n\) and auxiliary variables \(\pmb{u}\in\mathbb{R}^s\) and \(\alpha\in\mathbb{R}\), where \(n\) is the number of stocks, and \(s\) is the number of samples. The array \(\pmb{\pi}\) denotes the probabilities of samples, and \(\Pi\) is the uncertainty set that captures the distributional ambiguity of probabilities. The constant array \(\pmb{y}_k\in\mathbb{R}^n\) indicates the \(k\)th sample of stock return rates, and \(\bar{x}\) and \(\underline{x}\) are the lower and upper bounds of \(x\). The worst-case minimum expected overall return rate is set to be \(\mu=0.001\), the confidence level is \(\beta=0.95\), and the budget of investment is set to be \(w_0=1\). In this case study, we consider the sample data of eight stocks “JPM”, “AMZN”, “TSLA”, “AAPL”, “GOOG”, “ZM”, “META”, and “MCD”, in the year of 2021, and the other parameters are specified by the following code segment.
import pandas as pd
import yfinance as yf
stocks = ['JPM', 'AMZN', 'TSLA', 'AAPL', 'GOOG', 'ZM', 'META', 'MCD']
start = '2021-1-2' # starting date of historical data
end='2021-12-31' # end date of historical data
data = pd.DataFrame([])
for stock in stocks:
each = yf.Ticker(stock).history(start=start, end=end)
close = each['Close'].values
returns = (close[1:] - close[:-1]) / close[:-1]
data[stock] = returns
data
JPM | AMZN | TSLA | AAPL | GOOG | ZM | META | MCD | |
---|---|---|---|---|---|---|---|---|
0 | 0.005441 | 0.010004 | 0.007317 | 0.012364 | 0.007337 | 0.002361 | 0.007548 | 0.005994 |
1 | 0.046956 | -0.024897 | 0.028390 | -0.033662 | -0.003234 | -0.045506 | -0.028269 | -0.002270 |
2 | 0.032839 | 0.007577 | 0.079447 | 0.034123 | 0.029943 | -0.005546 | 0.020622 | 0.004645 |
3 | 0.001104 | 0.006496 | 0.078403 | 0.008631 | 0.011168 | 0.020759 | -0.004354 | 0.018351 |
4 | 0.014924 | -0.021519 | -0.078214 | -0.023249 | -0.022405 | -0.034038 | -0.040102 | -0.007597 |
... | ... | ... | ... | ... | ... | ... | ... | ... |
245 | 0.003574 | 0.000184 | 0.057619 | 0.003644 | 0.001317 | -0.007663 | 0.014495 | 0.003812 |
246 | 0.005723 | -0.008178 | 0.025248 | 0.022975 | 0.006263 | -0.021967 | 0.032633 | 0.008610 |
247 | 0.003035 | 0.005844 | -0.005000 | -0.005767 | -0.010914 | -0.019580 | 0.000116 | -0.001342 |
248 | -0.000504 | -0.008555 | -0.002095 | 0.000502 | 0.000386 | -0.010666 | -0.009474 | 0.002277 |
249 | -0.000505 | -0.003289 | -0.014592 | -0.006578 | -0.003427 | 0.047907 | 0.004141 | -0.004767 |
250 rows × 8 columns
import numpy as np
y = data.values # stock data as an array
s, n = y.shape # s: sample size; n: number of stocks
x_lb = np.zeros(n) # lower bounds of investment decisions
x_ub = np.ones(n) # upper bounds of investment decisions
beta =0.95 # confidence interval
w0 = 1 # investment budget
mu = 0.001 # target minimum expected return rate
In the nominal model, the CVaR and expected returns are evaluated assuming the exact distribution of stock returns is accurately represented by the historical samples without any distributional ambiguity. In other words, \(\Pi\) is written as a singleton uncertainty \(\Pi = \{\pmb{\pi}^0 \}\), where \(\pi_k^0=1/s\), with \(k=1, 2, …, s\). The Python code for implementing the nominal model is given below.
from rsome import ro
from rsome import msk_solver as msk
model = ro.Model()
pi = np.ones(s) / s
x = model.dvar(n)
u = model.dvar(s)
alpha = model.dvar()
model.min(alpha + 1/(1-beta) * (pi@u))
model.st(u >= -y@x - alpha)
model.st(u >= 0)
model.st(pi@y@x >= mu)
model.st(x >= x_lb, x <= x_ub, x.sum() == w0)
model.solve(msk)
Being solved by Mosek...
Solution status: Optimal
Running time: 0.0208s
The portfolio decision for the nominal model is retrieved by the following code.
x_nom = x.get()
model.get()
0.018195323668177367
Now we consider a box uncertainty set
\[\Pi = \left\{\pmb{\pi}: \pmb{\pi} = \pmb{\pi}^0 + \pmb{\eta}, \pmb{1}^{\top}\pmb{\eta}=0, \underline{\pmb{\eta}}\leq \pmb{\eta} \leq \bar{\pmb{\eta}} \right\}.\]In this case study, we assume that \(-\underline{\pmb{\eta}}=\bar{\pmb{\eta}}=0.0001\), and the Python code for implementation is provided below.
from rsome import ro
from rsome import msk_solver as msk
model = ro.Model()
eta_ub = 0.0001 # upper bound of eta
eta_lb = -0.0001 # lower bound of eta
eta = model.rvar(s) # eta as random variables
uset = (eta.sum() == 0,
eta >= eta_lb,
eta <= eta_ub)
pi = 1/s + eta # pi as inexact probabilities
x = model.dvar(n)
u = model.dvar(s)
alpha = model.dvar()
model.minmax(alpha + 1/(1-beta) * (pi@u), uset)
model.st(u >= -y@x - alpha)
model.st(u >= 0)
model.st(pi@y@x >= mu)
model.st(x >= x_lb, x <= x_ub, x.sum() == w0)
model.solve(msk)
Being solved by Mosek...
Solution status: Optimal
Running time: 0.0208s
x_box = x.get()
model.get()
0.018541181700867923
In cases that \(\Pi\) is an ellipsoidal uncertainty set
\[\Pi = \left\{\pmb{\pi}: \pmb{\pi} = \pmb{\pi}^0 + \rho\pmb{\eta}, \pmb{1}^{\top}\pmb{\eta}=0, \pmb{\pi}^0 + \rho\pmb{\eta} \geq \pmb{0}, \|\pmb{\eta}\| \leq 1 \right\},\]where the nominal probability \(\pmb{\pi}^0\) is the center of the ellipsoid, and the constant \(\rho=0.001\), then the model can be implemented by the code below.
from rsome import ro
from rsome import msk_solver as msk
import rsome as rso
model = ro.Model()
rho = 0.001
eta = model.rvar(s)
uset = (eta.sum() == 0, 1/s + rho*eta >= 0,
rso.norm(eta) <= 1)
pi = 1/s + rho*eta
x = model.dvar(n)
u = model.dvar(s)
alpha = model.dvar()
model.minmax(alpha + 1/(1-beta) * (pi@u), uset)
model.st(u >= -y@x - alpha)
model.st(u >= 0)
model.st(pi@y@x >= mu)
model.st(x >= x_lb, x <= x_ub, x.sum() == w0)
model.solve(grb)
Being solved by Mosek...
Solution status: Optimal
Running time: 0.0284s
x_ellip = x.get()
model.get()
0.018270742978340224
Here, we consider the KL divergence-constrained ambiguity of probabilities
\[\Pi = \left\{\boldsymbol{\pi}: \boldsymbol{\pi} \geq 0, \boldsymbol{1}^{\top}\boldsymbol{\pi} = 1, \sum_{k=1}^s\pi_k\log(\pi_k/\hat{\pi}_k) \leq \epsilon \right\},\]where \(\hat{\pi}_k = 1/s\) is the empirical probability of each sample. Assume that the constant \(\epsilon=0.001\), the code for implementing such a robust model is given below.
from rsome import ro
from rsome import msk_solver as msk
import rsome as rso
model = ro.Model()
epsilon = 0.001
pi = model.rvar(s)
uset = (pi.sum() ==1, pi >= 0,
rso.kldiv(pi, 1/s, epsilon)) # uncertainty set of pi
x = model.dvar(n)
u = model.dvar(s)
alpha = model.dvar()
model.minmax(alpha + 1/(1-beta) * (pi@u), uset)
model.st(u >= -y@x - alpha)
model.st(u >= 0)
model.st(pi@y@x >= mu)
model.st(x >= x_lb, x <= x_ub, x.sum() == w0)
model.solve(msk)
Being solved by Mosek...
Solution status: Optimal
Running time: 0.0936s
x_kld = x.get()
model.get()
0.02128303758055805
Decisions in terms of the allocations of capital in each stock are shown by the figure below.
import matplotlib.pyplot as plt
xdata = np.arange(n)
width = 0.15
plt.figure(figsize=(8, 3.5))
plt.bar(xdata - 1.5*width, x_nom,
width=width, color='b', alpha=0.6, label='Nominal Distribution')
plt.bar(xdata - 0.5*width, x_box,
width=width, color='r', alpha=0.6, label='Box Uncertainty')
plt.bar(xdata + 0.5*width, x_ellip,
width=width, color='m', alpha=0.6, label='Ellipsoidal Uncertainty')
plt.bar(xdata + 1.5*width, x_kld,
width=width, color='g', alpha=0.6, label='KL Divergence')
plt.legend()
plt.ylabel('Fraction of Investment', fontsize=12)
plt.xticks(xdata, data.columns)
plt.show()
In this example, we show that data acquisition tools provided in the Python ecosystem (e.g., pandas-datareader
) can be readily used to collect and feed real data into RSOME models. Apart from acquiring data, rich machine learning tools in the Python ecosystem can also be used to develop data-driven optimization models. More such examples will be provided in introducing the dro
module for modeling distributionally robust optimization problems.
Zhu, Shushang, and Masao Fukushima. 2009. Worst-case conditional value-at-risk with application to robust portfolio management. Operations Research 57(5) 1155-1168.