# Stratified model Implementation
## Dataset
The experiment utilizes daily adjusted closing prices of stocks in 18 ETFs. The dataset is divided into the following periods:
- **In-sample data**: From January 1st, 2018, to June 30th, 2023, used for model training. The validation set consisting of 20% of them
- **Out-of-sample data:** From July 1st, 2023, to July 31st, 2024, used for evaluating the model's performance on unseen data.
The 18 chosen ETFs span diverse asset classes:
| Ticker | Name | Description | Asset Class |
|--------|------|-------------|-------------|
| AGG | iShares Core U.S. Aggregate Bond ETF | Invests in U.S. investment-grade bonds | Fixed Income |
| DBC | Invesco DB Commodity Index Tracking Fund | Tracks various commodities | Commodities |
| GLD | SPDR Gold Shares | Invests in gold bullion | Commodities |
| IBB | iShares Biotechnology ETF | Tracks biotechnology stocks | Equity |
| ITA | iShares U.S. Aerospace & Defense ETF | Tracks U.S. aerospace and defense companies | Equity |
| PBJ | Invesco Dynamic Food & Beverage ETF | Tracks U.S. food and beverage companies | Equity |
| IEF | iShares 7-10 Year Treasury Bond ETF | Invests in U.S. Treasury bonds (7-10 years) | Fixed Income |
| VNQ | Vanguard Real Estate ETF | Invests in REITs | Real Estate |
| VTI | Vanguard Total Stock Market ETF | Tracks broad U.S. stock market | Equity |
| XLB | Materials Select Sector SPDR Fund | Tracks materials sector stocks | Equity |
| XLE | Energy Select Sector SPDR Fund | Tracks energy sector stocks | Equity |
| XLF | Financial Select Sector SPDR Fund | Tracks financial sector stocks | Equity |
| XLI | Industrial Select Sector SPDR Fund | Tracks industrial sector stocks | Equity |
| XLK | Technology Select Sector SPDR Fund | Tracks technology sector stocks | Equity |
| XLP | Consumer Staples Select Sector SPDR Fund | Tracks consumer staples sector stocks | Equity |
| XLU | Utilities Select Sector SPDR Fund | Tracks utilities sector stocks | Equity |
| XLV | Health Care Select Sector SPDR Fund | Tracks healthcare sector stocks | Equity |
| XLY | Consumer Discretionary Select Sector SPDR Fund | Tracks consumer discretionary sector stocks | Equity |
### Importing Libraries and Loading Data
```python=
import cvxportfolio as cvx
import pandas as pd
import numpy as np
import yfinance as yf
import pandas_datareader as pdr
import matplotlib.pyplot as plt
import networkx as nx
import scipy.sparse as sparse
import scipy.sparse.linalg as splinalg
import multiprocessing as mp
import time
import torch
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
import cvxpy as cp
```
### Data Loading and Preprocessing
```python=
def load_data(start_date='2007-01-01', end_date='2024-07-31'):
tickers = ["AGG", "DBC", "GLD", "IBB", "ITA", "PBJ", "IEF", "VNQ", "VTI", "XLB", "XLE", "XLF", "XLI", "XLK", "XLP", "XLU", "XLV", "XLY"]
returns = pd.DataFrame({ticker: yf.download(ticker, start=start_date, end=end_date)['Adj Close'].pct_change()
for ticker in tickers})
return returns.fillna(method='ffill').iloc[1:]
daily_returns = load_data()
daily_returns = daily_returns.dropna()
```
## Stratified Market Conditions
Each data record includes a market condition variable, denoted as $z$, which is known at the close of the previous trading day. We construct $z$ using several real-valued market indicators:
1. **Market Implied Volatility (VIX):** A 20-day moving average of the CBOE Volatility Index, lagged by one day. This captures market expectations of near-term volatility.
2. **Inflation Rate:** Measured using the 1-month rolling percentage change of the following ETFs, lagged by one day, to proxy for changes in the Consumer Price Index (CPI):
* TIP (iShares TIPS Bond ETF)
* DBA (Invesco DB Agriculture Fund)
* USO (United States Oil Fund)
* REZ (iShares Residential Real Estate ETF)
3. **Treasury 20+ Year Bond Index (TLT):** An 8-week rolling percentage change of the iShares 20+ Year Treasury Bond ETF, lagged by one day. This reflects trends in long-term interest rates.
### ETFs Used for Market Condition Construction
| Ticker | Name | Description | Asset Class |
|--------|------|-------------|-------------|
| VIX | CBOE Volatility Index | Measures the market's expectation of 30-day volatility based on S&P 500 index options prices.s | Index |
| TIP | iShares TIPS Bond ETF | Tracks the performance of U.S. Treasury Inflation-Protected Securities | Bond|
| DBA | Invesco DB Agriculture Fund | Tracks an index of agricultural commodities' futures contracts | Commodity|
| USO | United States Oil Fund | Tracks the daily price movements of West Texas Intermediate (WTI) light, sweet crude oil. | Commodity|
| REZ | iShares Residential Real Estate ETF | Tracks the performance of U.S. real estate investment trusts (REITs) that focus on residential properties| Equity|
| TLT | iShares 20+ Year Treasury Bond ETF | Tracks the performance of U.S. Treasury bonds with remaining maturities greater than 20 years | Bond|
### Data Preparation and Discretization
```python=
def load_data_Z(start_date='2007-01-01', end_date='2024-07-31'):
constituents = ["^VIX", "TIP", "DBA", "USO", "REZ", "TLT"]
Z_price = pd.DataFrame({constituent: yf.download(constituent, start=start_date, end=end_date)['Adj Close']
for constituent in constituents})
Z = pd.DataFrame({
"VIX": Z_price["^VIX"].rolling(window=20).mean().shift(1),
"TIP": Z_price["TIP"].pct_change(20).mul(100).shift(1),
"DBA": Z_price["DBA"].pct_change(20).mul(100).shift(1),
"USO": Z_price["USO"].pct_change(20).mul(100).shift(1),
"REZ": Z_price["REZ"].pct_change(20).mul(100).shift(1),
"TLT": Z_price["TLT"].pct_change(40).mul(100).shift(1)
})
return Z.dropna()
Z_raw = load_data_Z()
```
#### Discretization
Each of the continuous market indicators is discretized into quartiles (four bins), labeled 1 through 4. This results in a total of 4 x 4 x 4 x 4 x 4 x 4 = 4096 possible combinations of market conditions.
```python=
def discretize_data(start_train, end_train, start_backtest, Z_raw):
Z = Z_raw.loc[start_train:end_train, :]
Z_holdout = Z_raw.loc[start_backtest:, :]
for col in Z.columns:
#No holdout data is used in computing the deciles
col_arr, col_bins = pd.qcut(Z.loc[:,col], q=4, labels=range(4), retbins=True)
Z.loc[:, col] = pd.to_numeric(col_arr).astype(int)
def map_to_quantile(val):
return np.argmin(abs(val-col_bins[1:]))
Z_holdout.loc[:,col] = Z_holdout.loc[:,col].apply(map_to_quantile)
Z_train, Z_val = train_test_split(Z, test_size=0.2, random_state=0)
return Z, Z_holdout, Z_train, Z_val
start_train = "2018-01-01"
end_train = "2023-06-30"
start_backtest = "2023-07-01"
Z, Z_holdout, Z_train, Z_val = discretize_data(start_train, end_train, start_backtest, Z_raw)
```
### Correlation Matrix
The correlation matrix below provides insights into the relationships between the discretized market indicators during the training period:
```python=
Z.corr().round(2)
```
|---| VIX| TIP |DBA |USO |REZ| TLT|
| -------- | -------- | -------- | -------- | -------- | -------- |-------- |
VIX| 1.00| -0.07| 0.06| -0.05| -0.11 |-0.25
TIP| -0.07| 1.00| -0.06| -0.01| 0.36| 0.50
DBA| 0.06 |-0.06| 1.00| 0.26 |0.04| -0.20
USO| -0.05 |-0.01 |0.26 |1.00| 0.00 |-0.23
REZ| -0.11| 0.36| 0.04| 0.00| 1.00 |0.19
TLT| -0.25 |0.50 |-0.20 |-0.23 |0.19 |1.00
### Evolution of Market Conditions Over Time
```python=
fig, ax = plt.subplots(1,1,figsize=(15,5))
pd.concat([Z, Z_holdout]).plot(ax=ax, fontsize="large", color=["black", "blue", "red","pink","purple","green"])
ax.set_xlabel("Date", fontsize="x-large")
ax.set_ylabel("Quartile", fontsize="x-large") # Use "Quartile" instead of "Decile"
ax.set_yticks(range(4))
ax.set_yticklabels(range(0,4))
plt.axvline(start_backtest, color="black")
plt.tight_layout()
plt.show()
# Calculate and print average L1 norm of differences between consecutive market conditions
avg_change = Z.diff().dropna().apply(lambda x: np.linalg.norm(x,1), axis=1).mean()
print(f"The average daily change in market conditions is {avg_change:.2f} quartiles.")
```
```
The average daily change in market conditions is 1.27 quartiles.
```
#### Data Scarcity
The total number of possible market conditions, with each indicator having four quartiles, is $4^6 = 4096$. However, within the training and validation datasets, only 724 unique market condition combinations are observed. This means that for roughly 78% of the possible market condition combinations, we have no training data, highlighting the challenge of data scarcity in high-dimensional categorical spaces.
```python=
train_conditions = Z.value_counts().to_dict()
print(f"Number of unique market conditions in training data: {len(train_conditions)}")
```
```
Number of unique market conditions in training data: 724
```
## Modeling Expected Returns and Risk
The optimal parameters for the stratified models (both return and risk) are determined by minimizing an overall objective function:
$$
F(\theta_1, \dots, \theta_k) = \sum_{k=1}^K (l_k(\theta_k) + r(\theta_k)) + \mathcal{L}(\theta_1, \dots, \theta_K)
$$
### Data Preparation
We start by preparing the data for model fitting. This involves splitting the dataset into training, validation, and holdout (out-of-sample) sets.
```python=
start_train = "2018-01-01"
end_train = "2023-06-30"
start_backtest = "2023-07-01"
def Slicedf(start_train, end_train, start_backtest, prices):
df = prices.loc[start_train:end_train,:]
df_holdout = prices.loc[start_backtest:,:]
df_train, df_val = train_test_split(df, test_size=0.2, random_state=0)
return df_holdout, df_train, df_val
df_holdout, df_train, df_val = Slicedf(start_train, end_train, start_backtest, daily_returns)
```
### Stratified Return Model
The stratified return model's core function is to estimate the expected returns for each asset under various market conditions. It achieves this by employing a Huber loss function in conjunction with a Laplacian regularization term derived from the regularization graph.
#### Huber Loss Function
Distinguished by its robustness, the Huber loss function proves less susceptible to the influence of outliers when compared to the squared error loss. Its mathematical definition is as follows:
\begin{split}
\ell_k(\mu_k) &= \sum_{t:z_t=k} \mathbf{1}^T H(\mu_k − y_t) \\
H(z) &= \begin{cases}z^2,&\ |z|\le M\\2M|z| − M^2, &\ |z|\ge M\end{cases} \\
r(\mu_k) &= \gamma_{ret,loc}\|\mu_k\|^2_2
\end{split}
where $M$ is a tuning parameter controlling the transition point between the squared error and linear portions of the loss.
```python=
M = np.percentile(np.abs(df_train.values), 79)
```
```python=
# Import the stratified model, can download at (https://github.com/cvxgrp/strat_models)
import strat_models
torch.manual_seed(0)
def huber_return_prox(Y, nu, theta, t, M):
if Y is None:
return nu
###torch
if np.allclose(nu, 0):
nu_tch = torch.from_numpy(np.random.randn(nu.shape[0])/1000)
else:
nu_tch = torch.from_numpy(nu)
n,nk = Y[0].shape
theta_tch = torch.from_numpy(theta).requires_grad_(True)
loss = torch.nn.SmoothL1Loss(beta=M)
optim = torch.optim.SGD([theta_tch], lr=.01, momentum=0.9)
def closure():
optim.zero_grad()
l = torch.sum( (theta_tch - nu_tch)**2 )/(2*t)
for y in Y[0].T:
l += (1/1)*loss(theta_tch, torch.from_numpy(y))
l.backward()
return l
optim.step(closure)
return theta_tch.data.numpy()
class huber_return_loss(strat_models.Loss):
def __init__(self, M=None):
if M is None:
raise ValueError("M must be a number.")
super().__init__()
self.isDistribution = True
self.M = M
def evaluate(self, theta, data):
assert "Y" in data
return None
def setup(self, data, G):
Y = data["Y"]
Z = data["Z"]
K = len(G.nodes())
shape = (data["n"],)
theta_shape = (K,) + shape
#preprocess data
for y, z in zip(Y, Z):
vertex = G._node[z]
if "Y" in vertex:
vertex["Y"] += [y]
else:
vertex["Y"] = [y]
Y_data = []
for i, node in enumerate(G.nodes()):
vertex = G._node[node]
if 'Y' in vertex:
Y = vertex['Y']
Y_data += [Y]
del vertex['Y']
else:
Y_data += [None]
cache = {"Y": Y_data, "n":data["n"], "theta_shape":theta_shape, "shape":shape, "K":K}
return cache
def prox(self, t, nu, warm_start, pool, cache):
"""
Proximal operator for joint covariance estimation
"""
res = pool.starmap(huber_return_prox, zip(cache["Y"], nu, warm_start, t*np.ones(cache["K"]), self.M*np.ones(cache["K"])))
return np.array(res)
def logprob(self, data, G):
logprobs = []
for y,z in zip(data["Y"], data["Z"]):
n, nk = y.shape
y_bar = np.mean(y, axis=1).flatten()
if (y_bar == np.zeros(n)).all():
continue
mu = G._node[z]["theta"].copy().flatten()
lp = 0
for i in range(nk):
lp += sum((mu - y[:,i])**2) / (2*nk)
logprobs += [-lp]
return logprobs
def sample(self, data, G):
"""
Samples from ~N(mu_z, I)
"""
Z = turn_into_iterable(data["Z"])
mus = [G._node[z]["theta"] for z in Z]
n = mus[0].shape[0]
return [np.random.multivariate_normal(mu, np.eye(n)) for mu in mus]
def joint_cov_prox(Y, nu, theta, t):
"""
Proximal operator for joint covariance estimation
"""
if Y is None:
s, Q = np.linalg.eigh(nu)
s[s <= 0] = 1e-8
return Q @ np.diag(s) @ Q.T
n, nk = Y[0].shape
Yemp = Y[0]@Y[0].T/nk
s, Q = np.linalg.eigh(nu/(t*nk)-Yemp)
w = ((t*nk)*s + np.sqrt(((t*nk)*s)**2 + 4*(t*nk)))/2
return Q @ np.diag(w) @ Q.T
def prox(self, t, nu, warm_start, pool, cache):
"""
Proximal operator for joint covariance estimation
"""
res = pool.starmap(joint_cov_prox, zip(cache["Y"], nu, warm_start, t*np.ones(cache["K"])))
return np.array(res)
def logprob(self, data, G):
logprobs = []
for y,z in zip(data["Y"], data["Z"]):
n, nk = y.shape
Y = (y@y.T)/nk
if (np.zeros((n,n)) == Y).all():
continue
theta = G._node[z]["theta_tilde"]
logprobs += [np.linalg.slogdet(theta)[1] - np.trace(Y@theta)]
return logprobs
def sample(self, data, G):
Z = turn_into_iterable(data["Z"])
sigmas = [np.linalg.inv(G._node[z]["theta"]) for z in Z]
n = sigmas[0].shape[0]
return [np.random.multivariate_normal(np.zeros(n), sigma) for sigma in sigmas]
```
#### Laplacian Regularization
The Laplacian regularization serves to promote smoothness among model parameters across market conditions that exhibit similarities. We establish a regularization graph by taking the Cartesian product of six chain graphs. Each of these chain graphs links consecutive quartiles within each market indicator.
```python=
def make_G(w1, w2, w3, w4, w5, w6):
G_list = [nx.path_graph(4) for _ in range(6)]
weights = [w1, w2, w3, w4, w5, w6]
for G, w in zip(G_list, weights):
for u, v in G.edges():
G[u][v]['weight'] = w
G = nx.cartesian_product(G_list)
return G.copy()
```
#### Hyperparameter Optimization via Grid Search
We leverage grid search, a prevalent hyperparameter optimization method, to ascertain the most suitable combination of hyperparameters for the stratified return model. This entails a systematic evaluation of the model's performance across a predetermined array of hyperparameter values. The hyperparameter combination that yields the highest correlation between the model's estimated returns and the actual returns on the validation set is selected.
The hyperparameters subject to tuning include:
* `local`: This represents the strength of regularization applied to the local loss term.
* `w1`, `w2`, ..., `w6`: These denote the edge weights within the regularization graph, corresponding to each market condition indicator (VIX, TIP, USO, DBA, TLT, REZ).
```python=
def evaluate_return(params, train_return, val_return, Z_val, df_val, M):
local, w1, w2, w3, w4, w5, w6 = params
G = make_G(w1, w2, w3, w4, w5, w6)
loss = HuberReturnLoss(M=M)
reg = strat_models.sum_squares_reg(lambd=local)
bm = strat_models.BaseModel(loss=loss, reg=reg)
sm = strat_models.StratifiedModel(BaseModel=bm, graph=G)
sm.fit(data=train_return, **kwargs)
preds_val = np.vstack([
sm.G._node[tuple(Z_val.loc[date].values)]["theta"] for date in Z_val.index])
corr_val = corr(preds=preds_val, df=df_val)
print(f"STRATIFIED RETURN MODEL: train correlation = {corr_val:.4f}")
return corr_val
# Hyperparameter optimization
from itertools import product
def optimize_hyperparameters(param_grid, train_return, val_return, Z_val, df_val, M):
best_params = None
best_corr = -float('inf')
for params in product(*param_grid.values()):
corr_train = evaluate_return(params, train_return, val_return, Z_val, df_val, M)
if corr_train > best_corr:
best_corr = corr_train
best_params = params
print(f"New best parameters: {best_params}")
print(f"New best correlation: {best_corr:.4f}")
return best_params, best_corr
```
##### First Grid Search
```python=
param_grid = {
'local': [0.05, 0.1],
'w1': [1, 1000, 10000],
'w2': [1, 1000, 10000],
'w3': [1, 1000, 10000],
'w4': [1, 1000, 10000],
'w5': [1, 1000, 10000],
'w6': [1, 1000, 10000]
}
kwargs = dict(verbose=True, abs_tol=1e-6, maxiter=150, rho=1, n_jobs=8)
num_assets = len(daily_returns.columns)
train_return = get_data_dict(df_Y=df_train, df_Z=Z_train, num_assets=num_assets)
val_return = get_data_dict(df_Y=df_val, df_Z=Z_val, num_assets=num_assets)
best_params, best_corr = optimize_hyperparameters(param_grid, train_return, val_return, Z_val, df_val, M)
```
```
Best parameters:(0.05, 1, 1, 1000, 10000, 1000, 10000)
Best corr: 0.02149
```
##### Second Grid Search
```python=
param_grid = {
'local': [0.01, 0.05],
'w1': [1,5,10,20,50],
'w2': [1,5,10,20,50],
'w3': [100, 1000, 5000],
'w4': [2000, 5000, 7500, 10000],
'w5': [100, 1000, 5000],
'w6': [2000, 5000, 7500, 10000]
}
kwargs = dict(verbose=True, abs_tol=1e-6, maxiter=150, rho=1, n_jobs=8)
num_assets = len(daily_returns.columns)
train_return = get_data_dict(df_Y=df_train, df_Z=Z_train, num_assets=num_assets)
val_return = get_data_dict(df_Y=df_val, df_Z=Z_val, num_assets=num_assets)
best_params_2, best_corr = optimize_hyperparameters(param_grid, train_return, val_return, Z_val, df_val, M)
```
```
Best parameters:(0.01, 10, 1, 100, 7500, 5000, 5000)
Best corr: 0.04129
```
##### Explanation of Grid Search Results
The initial grid search encompassed a broader spectrum of hyperparameter values, whereas the second grid search concentrated on a narrower range centered around the promising values revealed by the first search. The optimal hyperparameter combination identified was:
```
local = 0.01, w1 = 10, w2 = 1, w3 = 100, w4 = 7500, w5 = 5000, w6 = 5000
```
This configuration led to a training correlation of 0.4508, signifying a reasonably strong association between the model's predictions and the actual returns within the training data.
#### Tuned Stratified Return Model
Using the optimized hyperparameters, we fit the tuned stratified return model to the training data. This model provides distinct return estimates for each asset under each identified market condition, as captured by the values stored in the `returns` dictionary.
```python=
def fit_stratified_return_model(df_train, Z_train, df_val, Z_val, best_params, **kwargs):
local, w1, w2, w3, w4, w5, w6 = best_params
G = make_G(w1, w2, w3, w4, w5, w6)
loss = huber_return_loss(M=M)
reg = strat_models.sum_squares_reg(lambd=local)
bm = strat_models.BaseModel(loss=loss, reg=reg)
sm = strat_models.StratifiedModel(BaseModel=bm, graph=G)
train_return = get_data_dict(df_Y=df_train, df_Z=Z_train, num_assets=num_assets)
val_return = get_data_dict(df_Y=df_val, df_Z=Z_val, num_assets=num_assets)
sm.fit(data=train_return, **kwargs)
preds_train = np.vstack([sm.G.nodes[tuple(Z_train.loc[date].values)]["theta"] for date in Z_train.index])
preds_val = np.vstack([sm.G.nodes[tuple(Z_val.loc[date].values)]["theta"] for date in Z_val.index])
corr_train = np.corrcoef(preds_train.flatten(), df_train.values.flatten())[0, 1]
corr_val = np.corrcoef(preds_val.flatten(), df_val.values.flatten())[0, 1]
print("STRATIFIED RETURN MODEL:")
print(f"\ttrain correlation = {corr_train:.4f}")
print(f"\tvalidation correlation = {corr_val:.4f}")
returns = {node: sm.G.nodes[node]["theta"].copy() for node in sm.G.nodes()}
return sm, returns, corr_train, corr_val
kwargs = dict(verbose=True, abs_tol=1e-6, maxiter=150, rho=1, n_jobs=4)
# Fit the stratified return model with the best parameters
sm_return, returns_pred, strat_corr_train, strat_corr_val = fit_stratified_return_model(df_train, Z_train, df_val, Z_val, best_params_2, **kwargs)
```
```
STRATIFIED RETURN MODEL:
train correlation = 0.4541
validation correlation = 0.0408
```
#### Common Return Model
To establish a benchmark for comparison, we also fit a common return model. This simpler model estimates a single return value for each asset by calculating the average return across the entire training dataset. Importantly, this baseline model does not factor in any variations in market conditions.
```python=
def fit_common_return_model(df_train, Z_train, df_val, Z_val):
common_return = df_train.mean(axis=0)
preds_train = np.vstack([common_return for _ in Z_train.index])
preds_val = np.vstack([common_return for _ in Z_val.index])
common_corr_train = np.corrcoef(preds_train.flatten(), df_train.values.flatten())[0, 1]
common_corr_val = np.corrcoef(preds_val.flatten(), df_val.values.flatten())[0, 1]
print("COMMON RETURN MODEL:")
print(f"\ttrain correlation = {common_corr_train:.4f}")
print(f"\tvalidation correlation = {common_corr_val:.4f}")
return common_return, common_corr_train, common_corr_val
# Fit the common return model
common_return, common_corr_train, common_corr_val = fit_common_return_model(df_train, Z_train, df_val, Z_val)
```
```
COMMON RETURN MODEL:
train correlation = 0.0188
validation correlation = -0.0057
```
#### Results and Interpretation
```python=
# Compare results
print("\nModel Comparison:")
print(f"{'Model':<20} {'Train Correlation':<20} {'Validation Correlation':<20}")
print(f"{'Stratified':<20} {strat_corr_train:<20.4f} {strat_corr_val:<20.4f}")
print(f"{'Common':<20} {common_corr_train:<20.4f} {common_corr_val:<20.4f}")
# Display return predictions
rets = pd.DataFrame(data=np.vstack([returns[key] for key in returns.keys()]), columns=df_train.columns)
tab = rets.describe().loc[["50%", "min", "max"]].rename(index={"50%":"median"})
tab = tab.T
tab["common"] = common_return
tab = tab[["common", "median", "min", "max"]].drop("VTI")
print("\nReturn Predictions (% daily returns):")
print((tab*100).round(3))
```
| Model | Train correlation |Validation correlation |
| -------- | -------- | -------- |
| Stratified return model |0.45 |0.0408 |
| Common return model| 0.0188 | -0.0057 |
It's readily apparent that the stratified return model demonstrates a substantially higher correlation on the training set. This suggests its proficiency in capturing return patterns that are specific to different market regimes. However, the lower correlation observed on the validation set raises a flag regarding potential overfitting. This is an aspect we will address in subsequent refinements of the model.
|-------- | common |median |min |max|
| -------- | -------- | -------- |-------- |-------- |
AGG | 0.004 |0.006 |-0.009 |0.015
DBC |0.041 |0.054 |-0.008 |0.080
GLD |0.039 |0.021 |-0.031 |0.081
IBB |0.026 |0.014 |-0.010 |0.046
ITA |0.068 |0.040 |-0.002 |0.071
PBJ |0.038 |0.035 |-0.021 |0.062
IEF |-0.001 |0.003 |-0.015 |0.010
VNQ |0.045 |0.031 |-0.025 |0.072
XLB |0.059 |0.047 |-0.026 |0.082
XLE |0.065 |0.045 |-0.032 |0.078
XLF |0.057 |0.043 |-0.004 |0.064
XLI |0.065 |0.052 |0.009 |0.088
XLK |0.114 |0.083 |0.012 |0.127
XLP |0.047 |0.038 |0.005 |0.053
XLU |0.060 |0.049 |-0.013 |0.074
XLV |0.063 |0.051 |0.019 |0.064
XLY |0.073 |0.067 |0.004 |0.110
The table above displays the return predictions (expressed as percentage daily returns) for each asset under both the common and stratified models. The first column presents the estimates from the common return model, while the subsequent columns showcase the median, minimum, and maximum return predictions across the various market conditions for the Laplacian-regularized stratified model. All returns are calculated relative to VTI, which is assigned a return of zero.
### Stratified Risk Model
In addition to modeling expected returns, we also need to estimate the risk associated with each asset under different market conditions. To achieve this, we model the covariance matrix for each market condition using a negative log-likelihood loss function:
$$
\ell_k (\theta_k) = \text{Tr}(S_k \Sigma_k^{−1}) − \log \text{det}(\Sigma_k^{−1})
$$
where
* $\theta_k$ represents the parameters of the covariance matrix $\Sigma_k$ for market condition $k$.
* $S_k = \frac{1}{n_k} \sum_{t:z_t=k}y_t y_t^T$ is the empirical covariance matrix calculated from the observed returns $y_t$ when the market condition is $k$.
* $n_k$ is the number of observations for market condition
#### Model Fitting
To estimate the risk model, we can reasonably assume that the mean return is negligible and treat it as zero. This simplifies the estimation process without significantly impacting the accuracy of the risk model.
```python=
class covariance_max_likelihood_loss(strat_models.Loss):
"""
f(theta) = Trace(theta @ Y) - logdet(theta)
"""
def __init__(self):
super().__init__()
self.isDistribution = True
def evaluate(self, theta, data):
assert "Y" in data
return np.trace(theta @ data["Y"]) - np.linalg.slogdet(theta)[1]
def setup(self, data, G):
Y = data["Y"]
Z = data["Z"]
K = len(G.nodes())
shape = (data["n"], data["n"])
theta_shape = (K,) + shape
#preprocess data
for y, z in zip(Y, Z):
vertex = G._node[z]
if "Y" in vertex:
vertex["Y"] += [y]
else:
vertex["Y"] = [y]
Y_data = []
for i, node in enumerate(G.nodes()):
vertex = G._node[node]
if 'Y' in vertex:
Y = vertex['Y']
Y_data += [Y]
del vertex['Y']
else:
Y_data += [None]
cache = {"Y": Y_data, "n":data["n"], "theta_shape":theta_shape, "shape":shape, "K":K}
return cache
```
##### Hyperparameter Search
Similar to the return model, we employ grid search for hyperparameter optimization. The objective is to find the hyperparameter combination that minimizes the negative log-likelihood (our loss function) on the validation set.
```python=
def evaluate_risk(params, train_cov, val_cov, **kwargs):
w1, w2, w3, w4, w5, w6 = params
G = make_G(w1, w2, w3, w4, w5, w6)
loss = covariance_max_likelihood_loss()
reg = strat_models.trace_reg(lambd=0)
bm = strat_models.BaseModel(loss=loss, reg=reg)
sm = strat_models.StratifiedModel(BaseModel=bm, graph=G)
sm.fit(data=train_cov, **kwargs)
risk_val_log = sm.anll(val_cov)
print(f"STRATIFIED RISK MODEL: train log-likelihood = {risk_val_log:.6f}")
return risk_val_log
# Hyperparameter optimization for risk model
def optimize_risk_hyperparameters(param_grid, train_cov, val_cov, **kwargs):
best_params = None
best_log = float('inf')
for params in product(*param_grid.values()):
risk_train_log = evaluate_risk(params, train_cov, val_cov, **kwargs)
if risk_train_log < best_log:
best_log = risk_train_log
best_params = params
print(f"New best parameters: {best_params}")
print(f"New best log-likelihood: {best_log:.6f}")
return best_params, best_log
```
```python=
param_grid_risk = {
'w1': [0.2, 1, 10, 20, 30, 50],
'w2': [0.2, 1, 10, 20, 30, 50],
'w3': [0.2, 1, 10, 20, 30, 50],
'w4': [0.2, 1, 10, 20, 30, 50],
'w5': [0.2, 1, 10, 20, 30, 50],
'w6': [0.2, 1, 10, 20, 30, 50]
}
kwargs_risk = dict(verbose=True, abs_tol=1e-3, maxiter=50, rho=25, n_jobs=8)
train_cov = get_data_dict(df_Y=df_train, df_Z=Z_train, num_assets=num_assets)
val_cov = get_data_dict(df_Y=df_val, df_Z=Z_val, num_assets=num_assets)
# Scale the data
for data in [train_cov, val_cov]:
for i in range(len(data["Y"])):
if not np.allclose(data["Y"][i], 0):
data["Y"][i] *= 100
best_params_risk, best_log = optimize_risk_hyperparameters(param_grid_risk, train_cov, val_cov, **kwargs_risk)
print(f"Best parameters: {best_params_risk}")
print(f"Best log-likelihood: {best_log:.6f}")
```
```
Best parameters: (0.2, 10, 20, 10, 50, 30)
Best corr: - 0.0905354
```
#### Tuned Stratified risk model
After the hyperparameter tuning process, the optimal values for the stratified risk model are determined to be:
\begin{split}
\gamma^{VIX}= 0.2
,\ \gamma^{TIP} = 10,\ \gamma^{DBA} = 20,\ \gamma^{USO} = 10, \ \gamma^{REZ}= 50,\ \gamma^{TLT}= 30
\end{split}
We then fit the tuned stratified risk model using these optimal hyperparameters. The model estimates a covariance matrix for each identified market condition, and these matrices are stored in the `covs` dictionary.
```python=
def fit_stratified_risk_model(df_train, Z_train, df_val, Z_val, best_params, **kwargs):
w1, w2, w3, w4, w5, w6 = best_params
G = make_G(w1, w2, w3, w4, w5, w6)
loss = covariance_max_likelihood_loss()
reg = strat_models.trace_reg(lambd=0)
bm = strat_models.BaseModel(loss=loss, reg=reg)
sm = strat_models.StratifiedModel(BaseModel=bm, graph=G)
train_cov = get_data_dict(df_Y=df_train, df_Z=Z_train, num_assets=df_train.shape[1])
val_cov = get_data_dict(df_Y=df_val, df_Z=Z_val, num_assets=df_val.shape[1])
# Scale the data
for data in [train_cov, val_cov]:
for i in range(len(data["Y"])):
if not np.allclose(data["Y"][i], 0):
data["Y"][i] *= 100
sm.fit(data=train_cov, **kwargs)
train_log = sm.anll(train_cov)
val_log = sm.anll(val_cov)
print("STRATIFIED RISK MODEL:")
print(f"\t(w1,w2,w3,w4,w5,w6)=({w1:.3f},{w2:.3f},{w3:.3f},{w4:.3f},{w5:.3f},{w6:.3f})")
print(f"\ttrain log-likelihood: {train_log:.6f}")
print(f"\tvalidation log-likelihood: {val_log:.6f}")
covs = {node: np.linalg.inv(sm.G.nodes[node]["theta"].copy()) for node in sm.G.nodes()}
return covs, train_log, val_log
# Fit the stratified risk model with the best parameters
kwargs_fit = dict(verbose=True, abs_tol=1e-3, maxiter=2500, rho=25, n_jobs=16)
covs_strat, train_log_strat, val_log_strat = fit_stratified_risk_model(df_train, Z_train, df_val, Z_val, best_params_risk, **kwargs_fit)
```
```
STRATIFIED RISK MODEL:
(w1,w2,w3,w4,w5,w6)=(0.200,10.000,20.000,10.000,50.000,30.000)
train log-likelihood: -3.356525
validation log-likelihood: -0.090509
```
#### Common Risk Model
As a baseline for comparison, we also fit a common risk model. This model estimates a single covariance matrix by calculating the empirical covariance across all training data, irrespective of market conditions.
```python=
def fit_common_risk_model(df_train, Z_train, df_val, Z_val):
train_cov = get_data_dict(df_Y=df_train, df_Z=Z_train, num_assets=df_train.shape[1])
val_cov = get_data_dict(df_Y=df_val, df_Z=Z_val, num_assets=df_val.shape[1])
# Scale the data
for data in [train_cov, val_cov]:
for i in range(len(data["Y"])):
if not np.allclose(data["Y"][i], 0):
data["Y"][i] *= 100
theta_common = (df_train*100).cov().values
G = make_G(w1=1e10, w2=1e10, w3=1e10,w4=1e10, w5=1e10, w6=1e10)
loss = covariance_max_likelihood_loss()
reg = strat_models.trace_reg(lambd=0)
bm_common = strat_models.BaseModel(loss=loss,reg=reg)
sm_common = strat_models.StratifiedModel(BaseModel=bm_common, graph=G)
for node in G.nodes():
sm_common.G._node[node]["theta"] = np.linalg.inv(theta_common)
sm_common.G._node[node]["theta_tilde"] = np.linalg.inv(theta_common)
sm_common.G._node[node]["theta_hat"] = np.linalg.inv(theta_common)
train_log = sm_common.anll(train_cov)
val_log = sm_common.anll(val_cov)
print("COMMON RISK MODEL:")
print(f"\ttrain log-likelihood: {train_log:.6f}")
print(f"\tvalidation log-likelihood: {val_log:.6f}")
covs_common = {node: np.linalg.inv(sm_common.G._node[node]["theta"].copy()) for node in sm_common.G.nodes()}
return covs_common, train_log, val_log
# Fit the stratified risk model with the best parameters
covs_common, train_log_common, val_log_common = fit_common_risk_model(df_train, Z_train, df_val, Z_val)
```
```
COMMON RISK MODEL:
train log-likelihood: -0.113580
validation log-likelihood: 2.357243
```
#### Results and Comparison
```python=
# Compare results
print("\nModel Comparison:")
print(f"{'Model':<20} {'Train Log-Likelihood':<25} {'Validation Log-Likelihood':<25}")
print(f"{'Stratified':<20} {train_log_strat:<25.6f} {val_log_strat:<25.6f}")
print(f"{'Common':<20} {train_log_common:<25.6f} {val_log_common:<25.6f}")
# Compare volatilities
common_vols = np.sqrt((100*df_train).cov().values.diagonal()/(100*100))
vols = pd.DataFrame(data=np.vstack([np.sqrt(covs_strat[key].diagonal()/(100*100)) for key in covs_strat.keys()]), columns=df_train.columns)
tab = vols.describe().loc[["50%", "min", "max"]].rename(index={"50%":"Median"}).T
tab["Common"] = common_vols
tab = tab[["Common", "Median", "min", "max"]].drop("VTI")
(tab*100).round(3)
```
| Model | Train loss |Validation loss |
| -------- | -------- | -------- |
| Stratified risk model |-3.356 | -0.0905 |
| Common risk model| -0.1135 | 2.3572 |
It's evident from the table that the stratified risk model significantly outperforms the common risk model on both training and validation sets. This suggests its ability to capture the varying risk characteristics associated with different market regimes.
|-------- | common |median |min |max|
| -------- | -------- | -------- |-------- |-------- |
AGG| 0.366 |0.638 |0.500 |2.000
DBC |1.212 |1.009 |0.652 |5.546
GLD |0.893 |0.857 |0.580 |4.167
IBB |1.550 |1.164 |0.724 |4.862
ITA |1.686 |1.054 |0.687 |8.304
PBJ |1.026 |0.831 |0.622 |4.561
IEF |0.420 |0.658 |0.523 |1.340
VNQ |1.418 |0.986 |0.651 |6.749
XLB |1.458 |1.005 |0.708 |6.359
XLE |2.198 |1.346 |0.800 |14.065
XLF |1.592 |1.018 |0.692 |7.981
XLI |1.430 |0.955 |0.598 |6.848
XLK |1.651 |1.132 |0.734 |6.399
XLP |1.004 |0.805 |0.623 |5.336
XLU |1.290 |0.940 |0.665 |6.341
XLV |1.134 |0.894 |0.646 |4.895
XLY |1.509 |1.042 |0.673 |5.888
The table above summarizes key statistics (median, minimum, and maximum) of the asset volatility predictions generated by the stratified risk model across different market conditions. It also includes the volatility estimates from the common risk model for comparison. Notably, the stratified model exhibits substantial variation in its volatility predictions, with some assets showing a nearly seven-fold difference between their minimum and maximum predicted volatilities. This highlights the model's ability to adapt its risk estimates to the prevailing market regime.
## Trading Policy and Backtesting
In this section, we will implement the trading policy outlined in the paper and evaluate its performance through backtesting.
### Data Preparation for Backtesting
We begin by preparing the data for the backtesting phase. This includes:
* Fetching the risk-free rate data using the FRED API.
* Preparing return and covariance estimates from the stratified and common models for the holdout (out-of-sample) period.
```python=
start_date = '2007-01-01'
end_date = '2024-07-31'
# Fetch risk-free rate data
risk_free_rate = pdr.get_data_fred('DFF', start=start_date, end=end_date)['DFF'] / (252 * 100)
BORROW_FEE = risk_free_rate.rolling(window=250).mean().shift(2).dropna()
# Fetch asset returns data
tickers = ["AGG", "DBC", "GLD","IBB","ITA","PBJ","IEF","VNQ","VTI","XLB","XLE","XLF","XLI","XLK","XLP","XLU","XLV","XLY"]
returns = pd.DataFrame(dict([(ticker, yfinance.download(ticker, start=start_date, end=end_date)['Adj Close'].pct_change())for ticker in tickers]))
returns = returns.fillna(method='ffill').iloc[1:]
```
```python=
strat_returns = np.vstack([
sm_return.G.nodes[tuple(Z.loc[date].values)]["theta"]
for date in Z.index])
r_hat_strat = pd.DataFrame(strat_returns, columns=df_holdout.columns, index= Z.index).shift(1).dropna()
# Create a DataFrame directly with the desired structure
r_hat_common = pd.DataFrame(
[common_return] * len(Z), index=Z.index, columns=df_holdout.columns ).shift(1).dropna()
# Prepare covariance estimates
# Create a list to hold individual DataFrames for each date and stock combination
cov_dfs = []
cov_dfc = []
for date in Z.index:
for stock in range(len(df_holdout.columns)):
cov_dfs.append(
pd.Series(covs_strat[tuple(Z.loc[date].values)][stock], index=df_holdout.columns)
)
cov_dfc.append(
pd.Series(covs_common[tuple(Z.loc[date].values)][stock], index=df_holdout.columns)
)
# Concatenate all DataFrames into a single one, transposing it for the desired shape
d_cov_ = pd.concat(cov_dfs, axis=1).T
d_cov_c = pd.concat(cov_dfc, axis=1).T
# Create a MultiIndex for the final output
iterables = [Z.index, df_holdout.columns]
index = pd.MultiIndex.from_product(iterables)
# Set the MultiIndex and assign to out_cov
cov_strat = d_cov_.set_axis(index, axis=0).shift(18)
cov_common = d_cov_c.set_axis(index, axis=0).shift(18)
```
### Trading Policy Implementation
We adopt the trading policy specified in the paper, which aims to minimize a combination of risk and trading costs while targeting a specific level of expected return:
\begin{split}
\text{minimize} &\quad \mu^T w -\gamma_{hold}\phi_{hold}(w) - \gamma_{trade} \phi_{trade}(z) \\
\text{subject to} &\quad w^T\Sigma w \le (\sigma_{tar})^2, \quad 1^Tw = 1,\\
&\quad w_{min} \leq w \leq w_{max}, \quad L \leq L_{tar}
\end{split}
For our implementation, we set $\gamma^{hold} = \gamma^{trade} = 1$ and impose weight and leverage constraints as follows:
* Long position limit $w_{max}$: 40%
* Short position limit $w_{min}$: -25%
* Maximum leverage $L_{tar}$: 2
* Target volatility $\sigma_{tar}$: 10%
These choices align with those in the original paper, allowing for diversification while permitting some degree of short selling.
```python=
# Prepare market data and cost models
market_data = cvx.UserProvidedMarketData(returns=daily_returns, cash_key='USDOLLAR')
tcost_model = cvx.TcostModel(a=0.0005) # Assuming HALF_SPREAD is 0.0005
hcost_model = cvx.HcostModel(short_fees=BORROW_FEE)
# Initial portfolio
init_portfolio = pd.Series(index=market_data.returns.columns, data=1000000/18)
init_portfolio.USDOLLAR = 0
# Trading constraints
max_weight = cvx.MaxWeights(0.4)
min_weight = cvx.MinWeights(-0.25)
l_limit = cvx.LeverageLimit(2)
# Risk targets
risk_target_str = cvx.FullCovariance(cov_stri) <= 0.00004
risk_target_common = cvx.FullCovariance(cov_com) <= 0.00004
# Define initial weights.
w_b = pd.Series(index=market_data.returns.columns, data=1)
w_b.USDOLLAR = 0.
w_b /= sum(w_b)
target_weights = w_b
rebalancing_times = pd.date_range('2018-01-01', '2023-12-31', freq='D')
# Define the market simulator and set teading frequency weekly
market_sim = cvx.MarketSimulator(
market_data = market_data,
costs = [
cvx.TcostModel(a=HALF_SPREAD, b=None),
cvx.HcostModel(short_fees=BORROW_FEE)])
```
We formulate the trading policies for both the stratified and common models using `cvx.SinglePeriodOpt`:
```python=
# Trading policies
gamma_trade, gamma_hold = 1.0, 1.0
spo_policy_stri = cvx.SinglePeriodOpt(
objective=cvx.ReturnsForecast(r_hat_) - gamma_trade * tcost_model - gamma_hold * hcost_model,
constraints=[risk_target_str, l_limit, max_weight, min_weight],
include_cash_return=False
)
spo_policy_common = cvx.SinglePeriodOpt(
objective=cvx.ReturnsForecast(common_return_df) - gamma_trade * tcost_model - gamma_hold * hcost_model,
constraints=[risk_target_common, l_limit, max_weight, min_weight],
include_cash_return=False
)
```
### Backtesting
We now conduct a backtest to compare the performance of the following strategies:
1. **Stratified Model:** Employs return and risk estimates from the stratified models.
2. **Common Model:** Utilizes return and risk estimates from the common models.
3. **Equally Weighted Portfolio:** An equally weighted portfolio with daily rebalancing allocates the same amount of capital to each asset and adjusts these allocations daily to maintain equal weights.
#### In sample
```python=
# Backtesting
in_results = market_sim.run_multiple_backtest(
h=[init_portfolio] * 3,
start_time='2018-01-03', end_time='2023-06-30',
policies=[spo_policy_stri, spo_policy_common, cvx.PeriodicRebalance(target=target_weights, rebalancing_times=rebalancing_times)]
)
# Plotting results
plt.figure(figsize=(12, 8))
for i, label in enumerate(['Stratified model', 'Common model', 'Equal weights']):
results[i].v.plot(label=label)
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.title('In-Sample Performance Comparison')
plt.legend()
plt.show()
```
##### Results
The in sample backtest results are summarized in the table below:
| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight |9.1% | 16.2% | 0.48 | 0.1 | 1 | -29.2% |
| Common model |1.6% |10.6%|0.0| 15.7|1.18 |-28.7%|
| Stratified model | 14.5% |10.2%|1.42| 15.6| 1.95 | -10.7% |
#### Out of sample
```python=
# Backtesting
results = market_sim.run_multiple_backtest(
h=[init_portfolio] * 3,
start_time='2023-07-05', end_time='2024-06-30',
policies=[spo_policy_stri, spo_policy_common, cvx.PeriodicRebalance(target=target_weights, rebalancing_times=rebalancing_times)]
)
# Plotting results
plt.figure(figsize=(12, 8))
for i, label in enumerate(['Stratified model', 'Common model', 'Equal weights']):
results[i].v.plot(label=label)
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.title('Out-of-Sample Performance Comparison')
plt.legend()
plt.show()
```
##### Results and Observations
The backtest results are summarized in the table below:
| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight |11.8% | 8.8% | 0.73 | 0.1 | 1 | -8.7% |
| Common model |11.2% |7.6%| 0.89| 36.2|1.16 |-6.0%|
| Stratified model | 20.3% | 9.3%| 1.62 | 19.7| 1.98 | -4.1% |
### Key Observations
1. **Superiority of the Stratified Model:** The stratified model showcases the best risk-adjusted performance (Sharpe ratio), effective volatility control, and a favorable balance between turnover and active management.
2. **Risk-Adjusted Returns:** The stratified model's highest Sharpe ratio (1.62) signifies its superior risk-adjusted returns compared to the other strategies.
3. **Return and Risk Management:** The stratified model achieves the highest return (20.3%) while maintaining a 9.3% target volatility demonstrates the model's success in limiting risk exposure. .
4. **Turnover and Leverage:** The stratified model's moderate turnover (19.7) reflects a balance between active management and cost efficiency. Its slightly higher leverage suggests it utilizes borrowing to potentially enhance returns.
5. **Overall Assessment:** The stratified model emerges as the most effective strategy, offering the best combination of return, risk management, and drawdown control. The equally weighted portfolio also performs reasonably well, particularly for investors seeking simplicity and minimal turnover. In contrast, the common model exhibits subpar performance across all metrics.
These findings emphasize the value of incorporating market conditions into both return and risk modeling for portfolio construction. The stratified model's ability to adapt to changing market regimes leads to improved investment outcomes compared to strategies that rely on static assumptions about asset behavior.
## Reference
* J. Tuck, S. Barratt, and S. Boyd, Portfolio Construction Using Stratified Models, Chapter 17 in Machine Learning in Financial Markets: A Guide to Contemporary Practice, A. Capponi and C.-A. Lehalle, editors, Cambridge University Press, pages 317–339, 2023.
* J. Tuck, S. Barratt, and S. Boyd, A Distributed Method for Fitting Laplacian Regularized Stratified Models, Journal of Machine Learning Research, (22):1–37, 2021.
* J. Tuck and S. Boyd Fitting Laplacian Regularized Stratified Gaussian Models
* S. Barratt and S. Boyd Fitting Feature-Dependent Markov Chains, Optimization and Engineering, 23:895–915, 2022.
* [lrsm_portfolio](https://github.com/cvxgrp/lrsm_portfolio)
* [strat_models](https://github.com/cvxgrp/strat_models)
Sign in with Wallet
Connect another wallet