Appendix · source
Mini Example Source Code
This computes every number quoted across the chapters.
Save it as mini_example.py and run python mini_example.py.
"""Mini example for the equity factor model primer.
Defines the 10-stock "MiniModel" used in every chapter's worked example and
computes all quantities quoted in the text: standardized exposures, the
constrained WLS cross-sectional regression, pure factor portfolios, the risk
model assembly, risk decomposition, performance attribution with multi-period
linking, the hedging example, and the optimization example.
Run: python mini_example.py
"""
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=160)
def section(title):
print("\n" + "=" * 78)
print(title)
print("=" * 78)
# ---------------------------------------------------------------------------
# 1. The universe
# ---------------------------------------------------------------------------
names = ["AXIOM", "BINARY", "CIPHER", "DIGIT", "EVERGREEN",
"FIDELIS", "GUARDIAN", "HARVEST", "INDIGO", "JUNIPER"]
industry = ["Tech", "Tech", "Tech", "Tech", "Fin",
"Fin", "Fin", "Cons", "Cons", "Cons"]
cap = np.array([150.0, 80, 40, 10, 120, 60, 20, 90, 30, 15]) # $bn
N = len(names)
cap_w = cap / cap.sum() # benchmark (cap) weights
# Raw style descriptors
btp = np.array([0.15, 0.25, 0.45, 0.60, 0.85, 0.95, 1.10, 0.40, 0.55, 0.70]) # book/price
mom_raw = np.array([0.32, 0.18, -0.05, 0.40, 0.06, -0.02, -0.12, 0.10, 0.02, -0.08]) # 12-1m ret
size_raw = np.log(cap)
# One-month realized returns (%, month 1)
r1 = np.array([4.2, 2.8, 0.5, 6.0, 0.8, -0.6, -1.8, 1.2, 2.0, -0.5])
# Annualized specific vols (%)
spec_vol = np.array([18.0, 22, 30, 38, 16, 20, 28, 17, 26, 32])
# Manager portfolio (long-only, value-tilted)
w_p = np.array([0.10, 0.08, 0.10, 0.03, 0.22, 0.14, 0.06, 0.15, 0.08, 0.04])
assert abs(w_p.sum() - 1) < 1e-12
w_b = cap_w.copy()
w_a = w_p - w_b
section("1. UNIVERSE")
print(f"{'name':10s} {'ind':5s} {'cap':>6s} {'capw%':>7s} {'B/P':>5s} {'mom':>6s} "
f"{'lncap':>6s} {'r1%':>5s} {'specvol':>7s} {'w_p':>6s} {'w_b':>7s} {'w_a':>8s}")
for i in range(N):
print(f"{names[i]:10s} {industry[i]:5s} {cap[i]:6.0f} {100*cap_w[i]:7.2f} "
f"{btp[i]:5.2f} {mom_raw[i]:6.2f} {size_raw[i]:6.2f} {r1[i]:5.1f} "
f"{spec_vol[i]:7.0f} {w_p[i]:6.2f} {w_b[i]:7.4f} {w_a[i]:8.4f}")
print("industry cap weights:",
{j: round(cap_w[[i for i in range(N) if industry[i] == j]].sum(), 4)
for j in ["Tech", "Fin", "Cons"]})
# ---------------------------------------------------------------------------
# 2. Standardized exposures: cap-weighted mean 0, equal-weighted std 1
# ---------------------------------------------------------------------------
def standardize(d):
z = d - cap_w @ d
return z / z.std(ddof=0)
val = standardize(btp)
mom = standardize(mom_raw)
size = standardize(size_raw)
factors = ["MKT", "TECH", "FIN", "CONS", "VALUE", "MOM", "SIZE"]
K = len(factors)
X = np.zeros((N, K))
X[:, 0] = 1.0
for i in range(N):
X[i, 1 + ["Tech", "Fin", "Cons"].index(industry[i])] = 1.0
X[:, 4], X[:, 5], X[:, 6] = val, mom, size
section("2. STANDARDIZED EXPOSURES X")
print(f"{'name':10s}" + "".join(f"{f:>8s}" for f in factors))
for i in range(N):
print(f"{names[i]:10s}" + "".join(f"{X[i, k]:8.3f}" for k in range(K)))
print("cap-weighted style means:", cap_w @ X[:, 4:])
print("equal-weighted style stds:", X[:, 4:].std(axis=0, ddof=0))
# ---------------------------------------------------------------------------
# 3. Constrained WLS cross-sectional regression (month 1)
# ---------------------------------------------------------------------------
# Constraint: cap-weighted industry factor returns sum to zero.
ind_w = np.array([cap_w[[i for i in range(N) if industry[i] == j]].sum()
for j in ["Tech", "Fin", "Cons"]])
# Restriction matrix R (K x K-1): free factors MKT,TECH,FIN,VALUE,MOM,SIZE;
# CONS = -(wT*TECH + wF*FIN)/wC
R = np.zeros((K, K - 1))
R[0, 0] = 1 # MKT
R[1, 1] = 1 # TECH
R[2, 2] = 1 # FIN
R[3, 1] = -ind_w[0] / ind_w[2] # CONS from TECH
R[3, 2] = -ind_w[1] / ind_w[2] # CONS from FIN
R[4, 3] = 1 # VALUE
R[5, 4] = 1 # MOM
R[6, 5] = 1 # SIZE
regw = np.sqrt(cap)
regw = regw / regw.sum()
W = np.diag(regw)
XR = X @ R
g = np.linalg.solve(XR.T @ W @ XR, XR.T @ W @ r1)
f1 = R @ g
P = R @ np.linalg.solve(XR.T @ W @ XR, XR.T @ W) # K x N: rows = pure factor ptfs
resid = r1 - X @ f1
section("3. MONTH-1 CONSTRAINED WLS")
print("regression weights (sqrt-cap, normalized):", regw)
print("\nfactor returns f1 (%):")
for k in range(K):
print(f" {factors[k]:6s} {f1[k]:8.4f}")
print("constraint check (cap-w industry sum):", ind_w @ f1[1:4])
print("\nresiduals e (%):")
for i in range(N):
print(f" {names[i]:10s} {resid[i]:8.4f}")
wss_tot = regw @ (r1 ** 2)
wss_res = regw @ (resid ** 2)
rbar = regw @ r1
wss_tot_c = regw @ ((r1 - rbar) ** 2)
print(f"\nweighted R^2 (uncentered): {1 - wss_res / wss_tot:.4f}")
print(f"weighted R^2 (centered): {1 - wss_res / wss_tot_c:.4f}")
section("3b. PURE FACTOR PORTFOLIOS P (rows: factors, cols: stocks)")
print(f"{'':8s}" + "".join(f"{n:>10s}" for n in names))
for k in range(K):
print(f"{factors[k]:8s}" + "".join(f"{P[k, i]:10.4f}" for i in range(N)))
print("\nP X (should be I on the constrained subspace):")
print(P @ X)
print("\nVerify f1 = P r1:", P @ r1)
print("VALUE pure portfolio: sum of weights =", P[4].sum(),
" gross leverage =", np.abs(P[4]).sum())
# ---------------------------------------------------------------------------
# 4. Risk model: factor covariance F (annualized) and specific risk
# ---------------------------------------------------------------------------
vols = np.array([16.0, 9, 7, 5, 4, 6, 4]) / 100.0 # annualized factor vols
C = np.array([
# MKT TECH FIN CONS VALUE MOM SIZE
[1.00, 0.10, -0.05, -0.10, -0.20, 0.05, 0.15],
[0.10, 1.00, -0.40, -0.30, -0.35, 0.30, 0.05],
[-0.05, -0.40, 1.00, -0.10, 0.40, -0.15, 0.00],
[-0.10, -0.30, -0.10, 1.00, 0.05, -0.05, -0.05],
[-0.20, -0.35, 0.40, 0.05, 1.00, -0.45, 0.10],
[0.05, 0.30, -0.15, -0.05, -0.45, 1.00, 0.05],
[0.15, 0.05, 0.00, -0.05, 0.10, 0.05, 1.00],
])
assert np.allclose(C, C.T)
F = np.outer(vols, vols) * C # annualized covariance (decimal^2)
eig = np.linalg.eigvalsh(C)
Delta = np.diag((spec_vol / 100.0) ** 2)
Sigma = X @ F @ X.T + Delta
section("4. RISK MODEL")
print("factor vols (% ann):", vols * 100)
print("correlation eigenvalues (must be > 0):", eig)
def risk_report(w, label):
x = X.T @ w
var_f = x @ F @ x
var_s = w @ Delta @ w
var = var_f + var_s
sig = np.sqrt(var)
print(f"\n--- {label} ---")
print("exposures x:", dict(zip(factors, np.round(x, 4))))
print(f"factor var {var_f:.6f} specific var {var_s:.6f} total var {var:.6f}")
print(f"vol: factor {np.sqrt(var_f)*100:.2f}% specific {np.sqrt(var_s)*100:.2f}% "
f"total {sig*100:.2f}% (ann)")
contrib = x * (F @ x) / var # share of total variance from each factor
print("factor contributions to TOTAL variance (%):",
dict(zip(factors, np.round(100 * contrib, 2))))
print(f"specific share of variance: {100*var_s/var:.2f}%")
return x, sig
x_p, sig_p = risk_report(w_p, "PORTFOLIO (total risk)")
x_b, sig_b = risk_report(w_b, "BENCHMARK (total risk)")
x_a, sig_a = risk_report(w_a, "ACTIVE (tracking error)")
# Stock-level contributions to tracking error
mcr = Sigma @ w_a / sig_a
ctr = w_a * mcr
section("4b. STOCK-LEVEL TE CONTRIBUTIONS (annualized)")
print(f"{'name':10s} {'w_a':>8s} {'MCR':>9s} {'CTR':>9s} {'CTR%':>7s}")
for i in np.argsort(-np.abs(ctr)):
print(f"{names[i]:10s} {w_a[i]:8.4f} {mcr[i]:9.4f} {ctr[i]:9.5f} {100*ctr[i]/sig_a:7.2f}")
print("sum of CTR:", ctr.sum(), " = TE:", sig_a)
# Factor-block decomposition of TE
blocks = {"Market": [0], "Industries": [1, 2, 3], "Styles": [4, 5, 6]}
var_a = sig_a ** 2
print("\nTE variance decomposition by block (x_k * (F x)_k summed in block):")
Fx = F @ x_a
for b, idx in blocks.items():
v = sum(x_a[k] * Fx[k] for k in idx)
print(f" {b:11s} {v:.6f} ({100*v/var_a:.2f}% of active variance)")
print(f" {'Specific':11s} {w_a @ Delta @ w_a:.6f} "
f"({100*(w_a @ Delta @ w_a)/var_a:.2f}% of active variance)")
# ---------------------------------------------------------------------------
# 5. Stress test: VALUE -2 sigma with correlated propagation
# ---------------------------------------------------------------------------
section("5. STRESS TEST (chapter 9)")
shock_size = -2 * vols[4] # -2 sigma on VALUE, annualized terms
kv = 4
f_cond = F[:, kv] / F[kv, kv] * shock_size # E[f | f_VALUE = shock]
print(f"VALUE shock: {shock_size*100:.1f}%")
print("propagated factor moves (%):", dict(zip(factors, np.round(100 * f_cond, 2))))
print(f"naive portfolio impact (only VALUE moves): {x_a[kv]*shock_size*100:+.2f}% active")
print(f"correlated impact (all factors move): {x_a @ f_cond*100:+.2f}% active")
print(f"portfolio (total) correlated impact: {x_p @ f_cond*100:+.2f}%")
# ---------------------------------------------------------------------------
# 6. Performance attribution, 3 months with Carino linking (chapter 10)
# ---------------------------------------------------------------------------
section("6. ATTRIBUTION (3 months, constant exposures assumed)")
def adjust_constraint(f):
f = f.copy()
f[3] = -(ind_w[0] * f[1] + ind_w[1] * f[2]) / ind_w[2]
return f
f2 = adjust_constraint(np.array([-2.0, -1.5, 1.0, 0.0, 1.2, -0.8, 0.5]))
f3 = adjust_constraint(np.array([3.0, 0.8, -0.5, 0.0, -0.6, 1.0, -0.3]))
print("f2 (%):", dict(zip(factors, np.round(f2, 3))))
print("f3 (%):", dict(zip(factors, np.round(f3, 3))))
spec_active = [w_a @ resid, 0.30, -0.10] # month-1 computed, months 2-3 stipulated
fs = [f1, f2, f3]
ra_m, fac_contrib_m = [], []
for m in range(3):
fc = x_a * fs[m] # per-factor active contribution (%)
ra = fc.sum() + spec_active[m]
fac_contrib_m.append(fc)
ra_m.append(ra)
print(f"\nmonth {m+1}: active return {ra:+.3f}% "
f"(factor {fc.sum():+.3f}%, specific {spec_active[m]:+.3f}%)")
print(" per-factor (%):", dict(zip(factors, np.round(fc, 3))))
# Portfolio and benchmark total returns per month (factor part + specific part)
# month 1 actuals:
rp1, rb1 = w_p @ r1, w_b @ r1
print(f"\nmonth-1 portfolio return {rp1:.3f}%, benchmark {rb1:.3f}%, "
f"active {rp1-rb1:+.3f}% (check vs {ra_m[0]:+.3f}%)")
# For months 2-3 stipulate benchmark returns to do Carino linking
rb = [rb1, x_b @ f2 + 0.0, x_b @ f3 + 0.0] # benchmark specific ~ 0
rp = [rb[m] + ra_m[m] for m in range(3)]
Rp = np.prod([1 + v / 100 for v in rp]) - 1
Rb = np.prod([1 + v / 100 for v in rb]) - 1
print(f"\ncumulative: portfolio {Rp*100:.3f}% benchmark {Rb*100:.3f}% "
f"active {100*(Rp-Rb):+.3f}%")
print(f"sum of monthly active returns (arithmetic): {sum(ra_m):+.3f}% <- does not match")
# Carino linking coefficients
kk = (np.log(1 + Rp) - np.log(1 + Rb)) / (Rp - Rb)
km = [(np.log(1 + rp[m] / 100) - np.log(1 + rb[m] / 100)) / ((rp[m] - rb[m]) / 100)
for m in range(3)]
linked = np.zeros(K)
linked_spec = 0.0
for m in range(3):
scale = km[m] / kk
linked += fac_contrib_m[m] * scale
linked_spec += spec_active[m] * scale
print("\nCarino-linked factor contributions (%):", dict(zip(factors, np.round(linked, 3))))
print(f"Carino-linked specific: {linked_spec:+.3f}%")
print(f"linked total: {linked.sum()+linked_spec:+.3f}% vs true active {100*(Rp-Rb):+.3f}%")
# ---------------------------------------------------------------------------
# 7. Hedging example (chapter 12)
# ---------------------------------------------------------------------------
section("7. HEDGING (chapter 12)")
# Instruments: broad index future (exposures = benchmark) and small-cap future
x_h1 = x_b.copy()
x_h2 = np.array([1.05, 0.35, 0.30, 0.35, 0.10, -0.05, -1.20])
print("index future exposures:", dict(zip(factors, np.round(x_h1, 4))))
print("small-cap future exposures:", dict(zip(factors, np.round(x_h2, 4))))
# Solve for h to zero portfolio MKT and SIZE exposures
A = np.array([[x_h1[0], x_h2[0]], [x_h1[6], x_h2[6]]])
b = -np.array([x_p[0], x_p[6]])
h = np.linalg.solve(A, b)
print(f"\nhedge notionals (per $1 of portfolio): index {h[0]:+.4f}, small-cap {h[1]:+.4f}")
x_hedged = x_p + h[0] * x_h1 + h[1] * x_h2
print("hedged exposures:", dict(zip(factors, np.round(x_hedged, 4))))
var_h = x_hedged @ F @ x_hedged + w_p @ Delta @ w_p
print(f"pre-hedge total vol {sig_p*100:.2f}% -> post-hedge {np.sqrt(var_h)*100:.2f}% "
f"(factor {np.sqrt(x_hedged @ F @ x_hedged)*100:.2f}%, "
f"specific {np.sqrt(w_p @ Delta @ w_p)*100:.2f}%)")
# Minimum-variance single-instrument hedge with the index future
cov_ph = x_p @ F @ x_h1
var_hh = x_h1 @ F @ x_h1
h_mv = -cov_ph / var_hh
x_mv = x_p + h_mv * x_h1
var_mv = x_mv @ F @ x_mv + w_p @ Delta @ w_p
print(f"\nmin-variance index-only hedge: h = {h_mv:+.4f}; "
f"post-hedge vol {np.sqrt(var_mv)*100:.2f}%")
print(f"(portfolio 'beta' to index future = {cov_ph/var_hh:.4f})")
# ---------------------------------------------------------------------------
# 8. Optimization example (chapter 11): neutralize MOM, keep VALUE, min TE
# ---------------------------------------------------------------------------
section("8. OPTIMIZATION (chapter 11)")
# Decision: active weights w. Minimize w' Sigma w s.t. 1'w = 0,
# (X'w)_MOM = 0, (X'w)_VALUE = current active value exposure.
Aeq = np.vstack([np.ones(N), X[:, 5], X[:, 4]])
beq = np.array([0.0, 0.0, x_a[4]])
KKT = np.block([[2 * Sigma, Aeq.T], [Aeq, np.zeros((3, 3))]])
sol = np.linalg.solve(KKT, np.concatenate([np.zeros(N), beq]))
w_opt = sol[:N]
x_opt = X.T @ w_opt
te_opt = np.sqrt(w_opt @ Sigma @ w_opt)
print("current active weights:", np.round(w_a, 4))
print("optimal active weights:", np.round(w_opt, 4))
print("implied portfolio weights:", np.round(w_b + w_opt, 4),
" min:", (w_b + w_opt).min())
print("optimal active exposures:", dict(zip(factors, np.round(x_opt, 4))))
print(f"TE: current {sig_a*100:.2f}% -> optimized {te_opt*100:.2f}%")
print(f"turnover from current active to optimal (one-way): "
f"{0.5*np.abs(w_opt - w_a).sum()*100:.1f}%")
# ---------------------------------------------------------------------------
# 9. Chapter 2 mini example: 5 stocks, 2 factors, hand-checkable
# ---------------------------------------------------------------------------
section("9. CHAPTER-2 MINI EXAMPLE (5 stocks, 2 factors)")
X5 = np.array([[1.0, 1.2], [1.0, 0.5], [1.0, -0.3], [1.0, -1.0], [1.0, -0.4]])
F5 = np.array([[0.0256, -0.00128], [-0.00128, 0.0016]]) # MKT 16%, VAL 4%, rho -0.2
d5 = np.diag(np.array([0.20, 0.25, 0.18, 0.30, 0.22]) ** 2)
w5 = np.array([0.30, 0.25, 0.20, 0.15, 0.10])
x5 = X5.T @ w5
vf = x5 @ F5 @ x5
vs = w5 @ d5 @ w5
print("exposures:", x5)
print(f"factor var {vf:.6f} specific var {vs:.6f} total vol {np.sqrt(vf+vs)*100:.2f}%")
print(f"factor vol {np.sqrt(vf)*100:.2f}% specific vol {np.sqrt(vs)*100:.2f}%")
print(f"factor share {100*vf/(vf+vs):.1f}%")