LaTeX macros (hidden cell) $ \newcommand{\Q}{\mathcal{Q}} \newcommand{\ECov}{\boldsymbol{\Sigma}} \newcommand{\EMean}{\boldsymbol{\mu}} \newcommand{\EAlpha}{\boldsymbol{\alpha}} \newcommand{\EBeta}{\boldsymbol{\beta}} $
import sys
import os
import re
import datetime as dt
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from mosek.fusion import *
import mosek.fusion.pythonic # Requires MOSEK >= 10.2
from notebook.services.config import ConfigManager
from portfolio_tools import data_download, DataReader, compute_inputs
# Version checks
print(sys.version)
print('matplotlib: {}'.format(matplotlib.__version__))
# Jupyter configuration
c = ConfigManager()
c.update('notebook', {"CodeCell": {"cm_config": {"autoCloseBrackets": False}}})
# Numpy options
np.set_printoptions(precision=5, linewidth=120, suppress=True)
# Pandas options
pd.set_option('display.max_rows', None)
# Matplotlib options
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['figure.dpi'] = 200
3.9.7 (default, Sep 16 2021, 13:09:58) [GCC 7.5.0] matplotlib: 3.7.2
Here we load the raw data that will be used to compute the yearly return observation series used for the optimization. The data consists of daily stock prices of $8$ stocks from the US market, and SPY as the benchmark.
# Data downloading:
# If the user has an API key for alphavantage.co, then this code part will download the data.
# The code can be modified to download from other sources. To be able to run the examples,
# and reproduce results in the cookbook, the files have to have the following format and content:
# - File name pattern: "daily_adjusted_[TICKER].csv", where TICKER is the symbol of a stock.
# - The file contains at least columns "timestamp", "adjusted_close", and "volume".
# - The data is daily price/volume, covering at least the period from 2016-03-18 until 2021-03-18,
# - Files are for the stocks PM, LMT, MCD, MMM, AAPL, MSFT, TXN, CSCO.
list_stocks = ["PM", "LMT", "MCD", "MMM", "AAPL", "MSFT", "TXN", "CSCO"]
list_factors = ["SPY"]
alphaToken = None
list_tickers = list_stocks + list_factors
if alphaToken is not None:
data_download(list_tickers, alphaToken)
We load the daily stock price data from the downloaded CSV files. The data is adjusted for splits and dividends. Then a selected time period is taken from the data.
investment_start = "2016-03-18"
investment_end = "2021-03-18"
# The files are in "stock_data" folder, named as "daily_adjusted_[TICKER].csv"
dr = DataReader(folder_path="stock_data", symbol_list=list_tickers)
dr.read_data()
df_prices, _ = dr.get_period(start_date=investment_start, end_date=investment_end)
Found data files: stock_data/daily_adjusted_AAPL.csv stock_data/daily_adjusted_PM.csv stock_data/daily_adjusted_CSCO.csv stock_data/daily_adjusted_TXN.csv stock_data/daily_adjusted_MMM.csv stock_data/daily_adjusted_IWM.csv stock_data/daily_adjusted_MCD.csv stock_data/daily_adjusted_SPY.csv stock_data/daily_adjusted_MSFT.csv stock_data/daily_adjusted_LMT.csv Using data files: stock_data/daily_adjusted_PM.csv stock_data/daily_adjusted_LMT.csv stock_data/daily_adjusted_MCD.csv stock_data/daily_adjusted_MMM.csv stock_data/daily_adjusted_AAPL.csv stock_data/daily_adjusted_MSFT.csv stock_data/daily_adjusted_TXN.csv stock_data/daily_adjusted_CSCO.csv stock_data/daily_adjusted_SPY.csv
Below we implement the optimization model in Fusion API. We create it inside a function so we can call it later.
def absval(M, x, z):
M.constraint(z >= x)
M.constraint(z >= -x)
def norm1(M, x, t):
z = M.variable(x.getSize(), Domain.greaterThan(0.0))
absval(M, x, z)
M.constraint(Expr.sum(z) == t)
def MinTrackingError(N, R, r_bm, x0, lambda_1, lambda_2, beta=1.5):
with Model("Case study") as M:
# Settings
M.setLogHandler(sys.stdout)
# Variables
# The variable x is the fraction of holdings in each security.
# It is restricted to be positive, which imposes the constraint of no short-selling.
x = M.variable("x", N, Domain.greaterThan(0.0))
xt = x - x0
# The variable t models the OLS objective function term (tracking error).
t = M.variable("t", 1, Domain.unbounded())
# The variables u and v model the regularization terms (transaction cost penalties).
u = M.variable("u", 1, Domain.unbounded())
v = M.variable("v", N, Domain.unbounded())
# Budget constraint
M.constraint('budget', Expr.sum(x) == 1.0)
# Objective
penalty_lin = lambda_1 * u
penalty_32 = lambda_2 * Expr.sum(v)
M.objective('obj', ObjectiveSense.Minimize, t + penalty_lin + penalty_32)
# Constraints for the penalties
norm1(M, xt, u)
M.constraint('market_impact', Expr.hstack(v, Expr.constTerm(N, 1.0), xt), Domain.inPPowerCone(1.0 / beta))
# Constraint for the tracking error
residual = R.T @ x - r_bm
M.constraint('tracking_error', Expr.vstack(t, 0.5, residual), Domain.inRotatedQCone())
# Create DataFrame to store the results. Last security name (the SPY) is removed.
columns = ["track_err", "lin_tcost", "mkt_tcost"] + df_prices.columns[:N].tolist()
df_result = pd.DataFrame(columns=columns)
# Solve optimization
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception("Unexpected solution status!")
# Save results
tracking_error = t.level()[0]
linear_tcost = u.level()[0]
market_impact_tcost = np.sum(v.level())
row = pd.Series([tracking_error, linear_tcost, market_impact_tcost] + list(x.level()), index=columns)
df_result = pd.concat([df_result, pd.DataFrame([row])], ignore_index=True)
return df_result
Here we use the loaded daily price data to compute the corresponding yearly mean return and covariance matrix for logarithmic returns, and compute linear return observations from that. The benchmark will be the last data series, corresponding to SPY.
# Number of securities and observations
N = df_prices.shape[1] - 1
T = 1000
# Mean and covariance of historical yearly log-returns.
m_log, S_log = compute_inputs(df_prices, return_log=True)
# Generate logarithmic return observations assuming normal distribution
scenarios_log = np.random.default_rng().multivariate_normal(m_log, S_log, T)
# Convert logarithmic return observations to linear return observations
scenarios_lin = np.exp(scenarios_log) - 1
We center and normalize the data matrices.
# Center the return data
centered_return = scenarios_lin - scenarios_lin.mean(axis=0)
# Security return scenarios
security_return = scenarios_lin[:, :N] / np.sqrt(T - 1)
# Benchmark return scenarios
benchmark_return = scenarios_lin[:, -1] / np.sqrt(T - 1)
We run the optimization for the given penalty coefficients, and initial portfolio.
lambda_1 = 0.0001
lambda_2 = 0.0001
x0 = np.ones(N) / N
df_result = MinTrackingError(N, security_return.T, benchmark_return, x0, lambda_1, lambda_2)
Problem Name : Case study Objective sense : minimize Type : CONIC (conic optimization problem) Constraints : 18 Affine conic cons. : 9 (1026 rows) Disjunctive cons. : 0 Cones : 0 Scalar variables : 27 Matrix variables : 0 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 1 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - primal attempts : 1 successes : 1 Lin. dep. - dual attempts : 0 successes : 0 Lin. dep. - primal deps. : 0 dual deps. : 0 Presolve terminated. Time: 0.00 Optimizer - threads : 64 Optimizer - solved problem : the dual Optimizer - Constraints : 16 Optimizer - Cones : 10 Optimizer - Scalar variables : 1060 conic : 1028 Optimizer - Semi-definite variables: 0 scalarized : 0 Factor - setup time : 0.00 Factor - dense det. time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 52 after factor : 52 Factor - dense dim. : 0 flops : 7.25e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 1.0e+00 1.3e+00 5.1e+00 0.00e+00 7.086559745e-01 -3.432754826e+00 1.0e+00 0.01 1 4.8e-01 6.2e-01 8.7e-01 2.71e+00 4.120265400e-01 -5.852995592e-01 4.8e-01 0.01 2 9.4e-02 1.2e-01 5.8e-02 1.65e+00 6.234989703e-02 -8.267879019e-02 9.4e-02 0.01 3 1.2e-02 1.6e-02 2.6e-03 1.16e+00 1.277717889e-02 -4.930056889e-03 1.2e-02 0.01 4 4.6e-03 6.0e-03 5.9e-04 1.01e+00 7.690849598e-03 1.100286126e-03 4.6e-03 0.01 5 1.5e-03 1.9e-03 8.8e-05 1.17e+00 7.162884540e-03 5.228693323e-03 1.5e-03 0.01 6 3.2e-04 4.2e-04 4.8e-06 1.43e+00 6.614174541e-03 6.288874887e-03 3.2e-04 0.01 7 5.7e-05 7.4e-05 3.2e-07 1.93e+00 6.352699725e-03 6.316106202e-03 5.7e-05 0.01 8 1.3e-05 1.7e-05 3.6e-08 1.11e+00 6.338763707e-03 6.330495870e-03 1.3e-05 0.01 9 8.8e-07 1.1e-06 5.7e-10 1.02e+00 6.335649771e-03 6.335109751e-03 8.8e-07 0.01 10 2.4e-08 2.4e-08 1.8e-12 1.00e+00 6.335321684e-03 6.335310105e-03 1.9e-08 0.01 11 7.2e-10 9.4e-10 1.3e-14 1.00e+00 6.335315628e-03 6.335315186e-03 7.3e-10 0.01 Optimizer terminated. Time: 0.02 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: 6.3353156286e-03 nrm: 1e+00 Viol. con: 1e-10 var: 0e+00 acc: 0e+00 Dual. obj: 6.3353151856e-03 nrm: 1e+00 Viol. con: 2e-11 var: 2e-11 acc: 0e+00
df_result
track_err | lin_tcost | mkt_tcost | PM | LMT | MCD | MMM | AAPL | MSFT | TXN | CSCO | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.006301 | 0.286418 | 0.060916 | 0.090031 | 0.166661 | 0.168084 | 0.171201 | 0.048057 | 0.101428 | 0.117276 | 0.137262 |