Porfolio Optimization¶
A portfolio manager oversees multiple assets \((i=1,\dots,N)\) and a bank account \((i=N+1)\). For a specified number of stages \(T\), the manager wants to maximize his utility by dynamically reblanacing the portfolio. Let \(\{r_{it}\}\) be the return process of asset \(i\). At the end of each period, the position of \(i^{th}\) asset \(x_{it}\) equals the start position \(x_{i,t-1}\), plus the realized return \(r_{it} x_{i,t-1}\) during the period, plus the newly long positions \(b_{it}\), minus the newly short positions \(s_{it}\). Transaction costs are \(f_b,f_s\) for buying and selling respectively. The capital in the bank account will be adjusted accordingly.
We consider a simple asset pricing model that decomposes the excess return as the return explained by Capital Asset Pricing Model (CAPM), alpha and idiosyncratic risk,
where \(\alpha_i,\beta_i,\sigma_i\) are assumed to be constant. We refer to \(\{r_{ft} + \beta_i (r_{Mt}-r_{ft})\}\) as the market-exposure return process and \(\{\epsilon_{it}\}\) as the idiosyncratic return process. The market return process \(\{r_{Mt}\}\) is modelled as a first-order autoregressive process (AR) with normal generalized autoregressive conditional heteroscedastic GARCH(1,1) innovation due to,
Formulation¶
Data¶
[1]:
import pandas,numpy
rf = 0.0005
fee = 0.001
params = pandas.read_csv("./data/parameters.csv",index_col=0)
coeffs = pandas.read_csv("./data/coefficients.csv",index_col=0)
mu,phi,omega,alpha_1,alpha_2 = params.iloc[:,0]
alpha = numpy.array(coeffs['alpha'])
beta = numpy.array(coeffs['beta'])
sigma = numpy.array(coeffs['epsilon'])
Solution¶
Consider three stage and five assets.
[2]:
T = 3
N = 5
The indiosyncratic return process is straightforward to construct
[3]:
def f(alpha,sigma):
def inner(random_state):
return random_state.normal(alpha+1,sigma)
return inner
The Markovian market-exposure return process \(\{r_{it}\},i=1,\dots,100\) are a bit more involved. One can certainly build directly a 103 dimensional return process and then discretize. But there is a more efficient way. Given the fact that the market-exposure return processes are virtually determined by the market return process, we can first discretize the market return process and then arithmetically compute the discretization of the market-exposure return process. Therefore, we need to manually deal with the discretization a bit. Normally we don’t need to touch the discretize module, but this is a special case.
[4]:
from msppy.discretize import Markovian
[5]:
# Markovian process generator
def generator(random_state, size):
# (r_Mt, epsilon_Mt, sigma^2_Mt)
epsilon = random_state.normal(size=[T,size])
process = numpy.zeros(shape=[size,T,3])
process[:,0,0] = -0.006
process[:,0,2] = omega/(1-alpha_1-alpha_2)
for t in range(1,T):
process[:,t,2] = omega + alpha_1*process[:,t-1,1]**2 + alpha_2*process[:,t-1,2]
process[:,t,1] = numpy.sqrt(process[:,t,2]) * epsilon[t]
process[:,t,0] = mu + phi*process[:,t-1,0] + process[:,t,1]
return process
[6]:
# augmented Markovian process generator
def generator_augmented(random_state, size):
# (r_it, r_Mt, epsilon_Mt, sigma^2_Mt)
process = generator(random_state, size)
market_return = process[:,:,0]
process_aug = numpy.concatenate(
(beta[:N]*(market_return[:,:,numpy.newaxis]-rf) + rf,process),
axis=-1,
)
return process_aug
[7]:
# Markov chain discretization
sample_paths = generator(numpy.random.RandomState(0),size=1000)
return_sample_paths = sample_paths[:,:,0]
var_sample_paths = sample_paths[:,:,2]
price_sample_paths = numpy.cumprod(numpy.exp(return_sample_paths),axis=1)
markovian = Markovian(generator,n_Markov_states=[1]+[100]*(T-1),n_sample_paths=100000)
markovian.SA()
# augment to N+3 dimension
Markov_states = [None for _ in range(T)]
transition_matrix = markovian.transition_matrix
for t in range(T):
market_return = markovian.Markov_states[t][:,0].reshape(-1,1)
asset_return_market_exposure = beta[:N]*(market_return-rf) + rf
Markov_states[t] = numpy.concatenate(
(asset_return_market_exposure,markovian.Markov_states[t]), axis=1)
[8]:
from msppy.msp import MSLP
from msppy.solver import SDDP
import gurobipy
[9]:
AssetMgt = MSLP(T=T, sense=-1, bound=200)
AssetMgt.add_Markovian_uncertainty(generator_augmented)
for t in range(T):
m = AssetMgt[t]
now, past = m.addStateVars(N+1, lb=0, obj=0, name='asset')
if t == 0:
buy = m.addVars(N, name='buy')
sell = m.addVars(N, name='sell')
m.addConstrs(now[j] == buy[j] - sell[j] for j in range(N))
m.addConstr(
now[N] == 100
- (1+fee) * gurobipy.quicksum(buy[j] for j in range(N))
+ (1-fee) * gurobipy.quicksum(sell[j] for j in range(N))
)
elif t != T-1:
sell = m.addVars(N, name='sell')
buy = m.addVars(N, name='buy')
capm = m.addVars(N, lb = -gurobipy.GRB.INFINITY, name='capm')
idio = m.addVars(N, name='idio')
m.addConstr(
now[N] == (
(1+rf) * past[N]
- (1+fee) * gurobipy.quicksum(buy[j] for j in range(N))
+ (1-fee) * gurobipy.quicksum(sell[j] for j in range(N))
)
)
m.addConstrs(
now[j] == capm[j] + idio[j] + buy[j] - sell[j]
for j in range(N)
)
for j in range(N):
m.addConstr(past[j] == capm[j], uncertainty_dependent={past[j]:j})
m.addConstr(past[j] == idio[j], uncertainty={past[j]:f(alpha[j],sigma[j])})
else:
v = m.addVar(obj=1, lb=-gurobipy.GRB.INFINITY, name='wealth')
capm = m.addVars(N, lb = -gurobipy.GRB.INFINITY, name='capm')
idio = m.addVars(N, name='idio')
m.addConstr(v == gurobipy.quicksum(now[j] for j in range(N+1)))
m.addConstrs(
now[j] == capm[j] + idio[j]
for j in range(N)
)
for j in range(N):
m.addConstr(past[j] == capm[j], uncertainty_dependent={past[j]:j})
m.addConstr(past[j] == idio[j], uncertainty={past[j]:f(alpha[j],sigma[j])})
m.addConstr(now[N] == (1+rf) * past[N])
AssetMgt.discretize(
n_samples=100,
method='input',
Markov_states=Markov_states,
transition_matrix=transition_matrix,
random_state=888,
)
AssetMgt.set_AVaR(l=0.5, a=0.25)
AssetMgt_SDDP = SDDP(AssetMgt)
AssetMgt_SDDP.solve(max_iterations=10)
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
----------------------------------------------------------------
SDDP Solver, Lingquan Ding
----------------------------------------------------------------
Iteration Bound Value Time
----------------------------------------------------------------
1 200.000000 100.100025 2.631516
2 200.000000 100.100025 2.325153
3 101.129474 98.520141 2.808576
4 100.592014 100.253018 2.746442
5 100.157229 100.345333 3.842424
6 100.100025 100.262957 2.932778
7 100.100025 100.100025 1.899105
8 100.100025 100.100025 1.759192
9 100.100025 100.100025 1.858058
10 100.100025 100.100025 1.841612
----------------------------------------------------------------
Time: 24.644856214523315 seconds
Algorithm stops since iteration:10 has reached