HW 2 Part 2¶

By: Chengyi (Jeff) Chen

%load_ext autotime
%load_ext nb_black
%matplotlib inline

import matplotlib.pyplot as plt

plt.rcParams["figure.dpi"] = 300
plt.rcParams["figure.figsize"] = (16, 12)

import pandas as pd
import numpy as np
import cvxpy as cp
import scipy as sp
from collections import namedtuple
from statsmodels.graphics.gofplots import qqplot
from datetime import datetime, timedelta
from tqdm import tqdm

/opt/anaconda3/envs/ml/lib/python3.7/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm

SPY and US Dollar Index Monthly Returns¶

spy_monthly_returns = pd.read_excel(
    "./data/HW2.xlsx", sheet_name=2, usecols=[0, 1], header=[0], index_col=[0]
)
spy_monthly_returns.index = pd.to_datetime(spy_monthly_returns.index)
spy_monthly_returns.head()

	S&P - 500 Index
2001-07-31	-0.984416
2001-08-31	-6.260163
2001-09-30	-8.075231
2001-10-31	1.906875
2001-11-30	7.670629

time: 872 ms

us_dollar_monthly_returns = pd.read_excel(
    "./data/HW2.xlsx", sheet_name=2, usecols=[0, 2], header=[0], index_col=[0]
)
us_dollar_monthly_returns.index = pd.to_datetime(us_dollar_monthly_returns.index)
us_dollar_monthly_returns.head()

	US Dollar Index
2001-07-31	-1.916799
2001-08-31	-3.208739
2001-09-30	-0.008817
2001-10-31	1.278547
2001-11-30	1.105694

time: 855 ms

S&P 500 Equities Monthly Returns¶

equity_monthly_returns = (
    pd.read_excel("./data/HW2.xlsx", sheet_name=1, header=[1]).iloc[:, 4:].T
)
equity_monthly_returns.columns = pd.MultiIndex.from_frame(
    pd.read_excel("./data/HW2.xlsx", sheet_name=1, header=[1], usecols=[0, 1, 2, 3])
)
equity_monthly_returns.index = pd.to_datetime(equity_monthly_returns.index)
equity_monthly_returns.head()

NAME	3m Co	Abbott Labs	Abbvie Inc	Abiomed Inc	Accenture Plc Ireland	Activision Blizzard Inc	Adobe Sys Inc	Advance Auto Parts	Advanced Micro Devic	Aes Corp	...	Willis Towers Watson Pu	Wynn Resorts Ltd	Xcel Energy Inc	Xerox Corp	Xilinx Inc	Xylem Inc	Yum Brands Inc	Zimmer Hldgs Inc	Zions Bancorp	Zoetis Inc
TICKER	MMM	ABT	ABBV	ABMD	ACN	ATVI	ADBE	AAP	AMD	AES	...	WLTW	WYNN	XEL	XRX	XLNX	XYL	YUM	ZBH	ZION	ZTS
SECTOR	INDU	HLTH	HLTH	HLTH	INFT	TCOM	INFT	DSCR	INFT	UTIL	...	FINA	DSCR	UTIL	INFT	INFT	INDU	DSCR	HLTH	FINA	HLTH
Market Cap ($ Mil)	110949	127036	138674	14640	89887	35535	110435	11478	18449	9577	...	19733	10755	25326	4708	21552	11991	28707	21156	7830	41098
2001-07-31	-1.945100	12.112976	NaN	-19.423240	NaN	-13.577718	-20.234087	NaN	-36.816609	-11.033682	...	-2.816904	NaN	-4.077114	-16.614480	-3.006790	NaN	4.214608	NaN	-0.932129	NaN
2001-08-31	-6.446533	-7.258951	NaN	-1.315789	-0.401070	9.226306	-10.349165	NaN	-25.794085	-13.524804	...	9.043619	NaN	1.707327	15.288163	-2.400000	NaN	-6.842003	-4.895105	-1.698511	NaN
2001-09-30	-5.475902	4.326061	NaN	-6.826667	-14.429530	-26.532794	-28.616782	NaN	-39.852399	-61.292271	...	24.348706	NaN	4.106134	-15.760939	-39.728484	NaN	-7.977550	2.022059	-6.287271	NaN
2001-10-31	6.077282	2.599966	NaN	24.441900	37.803922	32.808711	10.091781	NaN	20.736196	8.034321	...	-0.427533	NaN	0.461629	-9.677419	29.281768	NaN	28.990256	11.387387	-10.696854	NaN
2001-11-30	10.339559	3.812730	NaN	-12.741490	28.628344	3.113693	21.515301	NaN	37.804878	19.277978	...	1.674522	NaN	-3.429815	20.000215	18.704799	NaN	-6.206153	4.367519	1.164570	NaN

5 rows × 505 columns

time: 1.9 s

intersect = lambda a, b: sorted(set(a) & set(b))


def β_hat(X, y, dropna=True):
    if dropna:
        X, y = X.dropna(), y.dropna()
        common_index = intersect(X.index, y.index)
        X, y = X.loc[common_index], y.loc[common_index]
    else:
        X, y = X.fillna(0.0), y.fillna(0.0)
    return (
        pd.DataFrame(np.linalg.inv(X.T @ X), columns=X.columns, index=X.columns)
        @ X.T
        @ y.values
    )


# Factor Exposures to SPY and US Dollar
factor_exposures = equity_monthly_returns.apply(
    lambda equity: β_hat(
        X=pd.concat(
            [
                spy_monthly_returns,
                us_dollar_monthly_returns,
                pd.Series(1, name="α", index=spy_monthly_returns.index),
            ],
            axis=1,
        ),
        y=equity,
        dropna=True,
    ),
    axis=0,
).T
factor_exposures.reset_index(level=[1, 2, 3], drop=True, inplace=True)
factor_exposures.rename(
    {
        "S&P - 500 Index": "β (equity market factor exposure)",
        "US Dollar Index": "β (us dollar factor exposure)",
    },
    axis=1,
    inplace=True,
)

# Sector Exposure
sector_factor = pd.concat(
    [
        pd.Series(equity_monthly_returns.columns.get_level_values("NAME")),
        pd.Series(equity_monthly_returns.columns.get_level_values("SECTOR")),
    ],
    axis=1,
)
sector_factor.set_index("NAME", inplace=True)
sector_factor = pd.get_dummies(sector_factor).drop(["SECTOR_UTIL"], axis=1)

# Size Exposure
size_factor = (
    pd.concat(
        [
            pd.Series(equity_monthly_returns.columns.get_level_values("NAME")),
            pd.Series(equity_monthly_returns.columns.get_level_values("SECTOR")),
            pd.Series(
                equity_monthly_returns.columns.get_level_values("Market Cap ($ Mil)")
            ),
        ],
        axis=1,
    )
    .groupby("SECTOR")["Market Cap ($ Mil)"]
    .transform(lambda x: (np.log(x) - np.log(x).mean()) / np.log(x).std())
    .to_frame()
)
size_factor.index = equity_monthly_returns.columns.get_level_values("NAME")
size_factor.rename(
    {"Market Cap ($ Mil)": "β (size factor exposure)"}, axis=1, inplace=True
)

# Pure Monthly Returns of each factor
γ_pure_factor_monthly_returns = equity_monthly_returns.apply(
    lambda returns: β_hat(
        X=pd.concat(
            [
                factor_exposures.drop(["α"], axis=1),
                sector_factor,
                size_factor,
                pd.Series(1, index=size_factor.index, name="γ_residual"),
            ],
            axis=1,
        ),
        y=returns,
        dropna=False,
    ),
    axis=1,
)
γ_pure_factor_monthly_returns.head()

	β (equity market factor exposure)	β (us dollar factor exposure)	SECTOR_DSCR	SECTOR_ENER	SECTOR_FINA	SECTOR_HLTH	SECTOR_INDU	SECTOR_INFT	SECTOR_MATS	SECTOR_REAL	SECTOR_STPL	SECTOR_TCOM	β (size factor exposure)	γ_residual
2001-07-31	-5.081906	0.546521	9.689164	4.497985	6.003830	3.868261	7.997561	5.971431	7.795002	3.298414	5.834931	4.522445	-0.666501	-0.172946
2001-08-31	-5.351213	-2.530621	-0.865683	-1.974351	-0.172676	0.152704	0.283117	-0.786473	2.069103	2.237107	0.399686	-0.210663	-0.887846	3.563947
2001-09-30	-16.448408	-4.242000	6.115425	5.129572	15.719172	7.981975	9.199334	8.180248	9.807025	11.983056	5.616929	8.494224	-0.404314	0.634008
2001-10-31	7.730491	4.661245	-0.237472	5.161065	-8.600545	2.401637	-3.667805	5.580222	-2.020707	-6.257782	-0.408720	0.168053	-0.612423	-2.110206
2001-11-30	9.140003	4.696506	3.084856	-7.999442	-0.934668	0.446165	1.300240	2.209795	3.643307	-0.218926	2.664973	-1.297612	1.092705	-3.667960

time: 4.07 s

Problems¶

Black Litterman Model¶

Preliminaries¶

From HW2 Part 1, you have located 10 sector factor returns, i.e., 10 monthly return time series. Let’s practise Black-Litterman model with the following assumptions:

For the 10 sector factors, equlibrium expected monthly returns $\Pi$ are assumed to equal historical average monthly returns.
Using historical monthly return covariance matrix $\Sigma$
Uncertainty metric $\tau$ for the prior distribution of the factor expected returns equals 0.05.
Assuming we have the views as follows:
- View 1: Sector ENER will have return of 0.20% in next month.
- View 2: Sector FINA will outperform Sector INDU by 0.05% in next month.
- View 3: Sector INFT and Sector HLTH will outperform Sector REAL and Sector STPL by 0.30% in next month. (Assuming equal weighting on both long and short legs)
Assuming view uncertainty for each view portfolio equals the historical variance of each view portfolio. ($\Omega=\tau P \Sigma P^\top$).

δ_eq = lambda R_mkt, R_f, σ_mkt: (R_mkt - R_f) / (σ_mkt ** 2)  # Risk Aversion measure
Π = lambda δ_eq, Σ, w_eq: δ_eq * (Σ @ w_eq)  # Equilibrium Excess Return vector
C = lambda τ, Σ: τ * Σ  # Covariance matrix for Prior Normal
H = lambda P, Ω, C: (P.T @ np.linalg.inv(Ω)) @ P + np.linalg.inv(
    C
)  # Precision Matrix of Posterior Normal
ν = lambda H, P, Ω, Q, C, Π: np.linalg.inv(H) @ (
    (P.T @ np.linalg.inv(Ω)) @ Q + np.linalg.inv(C) @ Π
)  # Mean of Posterior Normal

μ = γ_pure_factor_monthly_returns[sector_factor.columns].mean()
Π = μ  # 1. Arithmetic Historical average of monthly returns (N x 1)
Σ = γ_pure_factor_monthly_returns[
    sector_factor.columns
].cov()  # 2. Historical monthly covariance matrix (N x N)
τ = 0.05  # 3. Uncertainty metric for prior distribution of θ's covariance matrix

P = pd.DataFrame(
    np.zeros((3, len(sector_factor.columns))),
    index=["view 1", "view 2", "view 3"],
    columns=sector_factor.columns,
)  # View Portfolio Matrix (k x N)
Q = pd.DataFrame(
    np.zeros((3, 1)), index=["view 1", "view 2", "view 3"], columns=["expected return"]
)  # View Expected Return Vector (k x 1)

# Active View 1
P.loc["view 1"]["SECTOR_ENER"] = 1
Q.loc["view 1"]["expected return"] = 0.2

# Active View 2
P.loc["view 2"]["SECTOR_FINA"], P.loc["view 2"]["SECTOR_INDU"] = 1, -1
Q.loc["view 2"]["expected return"] = 0.05

# Active View 3
(
    P.loc["view 3"]["SECTOR_INFT"],
    P.loc["view 3"]["SECTOR_HLTH"],
    P.loc["view 3"]["SECTOR_REAL"],
    P.loc["view 3"]["SECTOR_STPL"],
) = (0.5, 0.5, -0.5, -0.5)
Q.loc["view 3"]["expected return"] = 0.3

Ω = lambda τ, P, Σ: np.diag(
    np.diag(τ * ((P @ Σ) @ P.T))
)  # 5. View uncertainty for each view portfolio equals the historical variance of each view portfolio

time: 9.5 ms

A.¶

Without active views, given $\Pi$ (from 1) and $\Sigma$ (from 2), locate the optimal portfolio of the 10 sectors which has the highest expected return and at the same time its annualized volatility no greater than 10%. Note: this is an optimization problem with only one constraint on the total variance of the portfolio. Also notice that you need to transform the annualized vol to monthly vol.

w_eq = cp.Variable((len(sector_factor.columns),), integer=False)
constraints = [
    cp.quad_form(w_eq, Σ)
    <= (10 ** 2)
    / 12,  # (annualized volatility no greater than 10%)^2 / 12 => monthly variance
]
obj = cp.Maximize(w_eq.T @ Π)
prob = cp.Problem(obj, constraints)
prob.solve()

print(f"Status: {prob.status}")
print(f"The optimal value is: {np.round(prob.value, 2)}")
print(f"The optimal solution is: w = {[np.round(w_i, 2) for w_i in w_eq.value]}")

Status: optimal
The optimal value is: 0.68
The optimal solution is: w = [0.23, 0.06, -0.64, 0.77, -0.28, 0.1, 0.44, 0.22, 0.02, -0.62]
time: 15.3 ms

B.¶

Now let’s incorporate the active views on 4 using Black-Litterman method,

show the updated factor expected return distribution and compare it with the original one, i.e., you need to calculate $\nu$, which is the update return vector; and $(H^{-1}+\Sigma)$, which is the updated covariance matrix.

(64)¶\[\begin{align} R\vert Q \sim \mathcal{N}(\nu, H^{-1} + \Sigma) \end{align}\]

updated_factor_expected_return = ν(
    H(P=P.values, Ω=Ω(τ=τ, P=P.values, Σ=Σ.values), C=C(τ, Σ.values)),
    P=P.values,
    Ω=Ω(τ=τ, P=P.values, Σ=Σ.values),
    Q=Q.values.flatten(),
    C=C(τ, Σ.values),
    Π=Π.values,
)

updated_factor_covariance_matrix = (
    np.linalg.inv(H(P=P.values, Ω=Ω(τ=τ, P=P.values, Σ=Σ.values), C=C(τ, Σ.values))) + Σ
)

time: 4.24 ms

Π.name = "original factor expected return"
pd.concat(
    [
        np.round(Π, 2),
        pd.Series(
            np.round(updated_factor_expected_return, 2),
            index=Π.index,
            name="updated factor expected return",
        ),
    ],
    axis=1,
)

	original factor expected return	updated factor expected return
SECTOR_DSCR	0.18	0.22
SECTOR_ENER	0.07	0.16
SECTOR_FINA	-0.19	-0.01
SECTOR_HLTH	0.37	0.47
SECTOR_INDU	0.09	0.12
SECTOR_INFT	0.16	0.26
SECTOR_MATS	0.19	0.25
SECTOR_REAL	0.37	0.32
SECTOR_STPL	0.05	0.08
SECTOR_TCOM	-0.11	-0.03

time: 8.36 ms

make observations on how active views impact updated returns.

Directly Affected Factors:

View 1: Sector ENER will have return of 0.20% in next month.
- This absolute view has pushed up the original SECTOR_ENER’s monthly expected return from 0.07% to 0.16%.
View 2: Sector FINA will outperform Sector INDU by 0.05% in next month.
- SECTOR_FINA’s monthly expected return has been pushed up from -0.19% to -0.01%.
- However, the relative view pushed up SECTOR_INDU’s monthly expected return from 0.09% to 0.12% even though the relative view is bearish towards it
- Although, SECTOR_FINA’s monthly expected return has increased more than SECTOR_INDU’s
View 3: Sector INFT and Sector HLTH will outperform Sector REAL and Sector STPL by 0.30% in next month. (Assuming equal weighting on both long and short legs)
- This relative view has pushed up SECTOR_HLTH’s monthly expected return from 0.37% to 0.47% and SECTOR_INFT from 0.16% to 0.26%, while pushing down SECTOR_REAL from 0.37% to 0.32%
- However, SECTOR_STPL’s monthly expected return has increased from 0.05% to 0.08% even though the relative view was bearish on it.

Indirectly Affected Factors:

Because Black-Litterman accounts for how changes in one factor’s expected return affects others through the covariance matrix:
- SECTOR_DSCR increased from 0.18% to 0.22%.
- SECTOR_MATS increased from 0.19% to 0.25%.
- SECTOR_TCOM increased from -0.11% to -0.03%.

Re-run the optimization as in A, but with the new return vector and covariance matrix. Please compare the updated optimal portfolio (the weight vector) with the solution from A.

updated_w_eq = cp.Variable((len(sector_factor.columns),), integer=False)
constraints = [
    cp.quad_form(updated_w_eq, updated_factor_covariance_matrix)
    <= (10 ** 2)
    / 12,  # (annualized volatility no greater than 10%)^2 / 12 => monthly variance
]
obj = cp.Maximize(updated_w_eq.T @ updated_factor_expected_return)
prob = cp.Problem(obj, constraints)
prob.solve()

print(f"Status: {prob.status}")
print(f"The optimal value is: {np.round(prob.value, 2)}")
print(
    f"The optimal solution is: w = {[np.round(w_i, 2) for w_i in updated_w_eq.value]}"
)

Status: optimal
The optimal value is: 0.63
The optimal solution is: w = [0.24, 0.07, -0.39, 0.89, -0.55, 0.2, 0.45, 0.12, -0.08, -0.64]
time: 16.7 ms

pd.DataFrame(
    np.round(np.array([w_eq.value, updated_w_eq.value]).T, 4) * 100,
    index=sector_factor.columns,
    columns=[
        "Original Optimal Portfolio Weights - A (%)",
        "Updated Optimal Portfolio Weights - B (%)",
    ],
)

	Original Optimal Portfolio Weights - A (%)	Updated Optimal Portfolio Weights - B (%)
SECTOR_DSCR	23.01	23.53
SECTOR_ENER	5.73	7.48
SECTOR_FINA	-63.93	-38.99
SECTOR_HLTH	77.41	89.05
SECTOR_INDU	-27.68	-54.70
SECTOR_INFT	9.83	19.94
SECTOR_MATS	43.70	44.69
SECTOR_REAL	21.60	12.21
SECTOR_STPL	1.78	-8.06
SECTOR_TCOM	-62.16	-63.57

time: 10 ms