# Mini Example Source Code

_Previous: [Chapter 18: Appendix: Reference Material](18-appendix.md)_

---

This computes every number quoted across the chapters.
Save it as `mini_example.py` and run `python mini_example.py`.

```python
"""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, an EWMA weighting demo, risk decomposition, performance
attribution with multi-period linking, the hedging example, the optimization
example, a bias-statistic demo, the CROWD factor case study, and the
month-over-month risk change attribution.

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)]
print("\nmonthly returns (%): portfolio", np.round(rp, 3), " benchmark", np.round(rb, 3))
print("Carino kappa_t per month:", np.round(km, 4), " cumulative kappa:", round(kk, 4))
print("linking scale kappa_t/kappa:", np.round(np.array(km) / kk, 4))
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}%")

# ---------------------------------------------------------------------------
# 6b. EWMA weighting demo (chapter 8)
# ---------------------------------------------------------------------------
section("6b. EWMA WEIGHTING DEMO (chapter 8)")
print("half-life -> lambda -> effective observations (1+l)/(1-l):")
for h_l in [20, 60, 90, 250]:
    lam = 2 ** (-1 / h_l)
    print(f"  h = {h_l:3d}: lambda = {lam:.4f}, "
          f"effective obs = {(1 + lam) / (1 - lam):6.1f}  (2.89h = {2.89 * h_l:6.1f})")

# The MOM factor's three monthly returns, weighted two ways (zero-mean form).
mom_f = np.array([f1[5], f2[5], f3[5]])  # MOM factor returns, months 1..3 (%)
lam = 2 ** (-1 / 1.0)                    # half-life = 1 month
wts = lam ** np.arange(len(mom_f) - 1, -1, -1)
wts = wts / wts.sum()
var_eq = (mom_f ** 2).mean()             # zero-mean convention
var_ew = wts @ (mom_f ** 2)
print(f"\nMOM factor returns, months 1-3 (%): {np.round(mom_f, 3)}")
print(f"EWMA weights (h = 1 month, oldest first): {np.round(wts, 3)}")
print(f"monthly variance (%^2): equal-weighted {var_eq:.3f}  EWMA {var_ew:.3f}")
print(f"annualized vol: equal-weighted {np.sqrt(12 * var_eq):.2f}%  "
      f"EWMA {np.sqrt(12 * var_ew):.2f}%  (stipulated MOM vol: 6%)")

# ---------------------------------------------------------------------------
# 7. Hedging example (chapter 13)
# ---------------------------------------------------------------------------
section("7. HEDGING (chapter 13)")
# 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 12): neutralize MOM, keep VALUE, min TE
# ---------------------------------------------------------------------------
section("8. OPTIMIZATION (chapter 12)")
# 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}%")

# ---------------------------------------------------------------------------
# 10. Alpha research: is a candidate signal spanned by the model? (chapter 14)
# ---------------------------------------------------------------------------
section("10. ALPHA RESEARCH (chapter 14)")
# A research team proposes a raw "value-quality composite" score as alpha (higher
# = more attractive). Built here to look novel but mostly track known factors:
# cheap names, with a tilt away from recent winners. This is the common failure
# the factor model is meant to catch before the signal reaches production.
alpha_raw = btp - 0.30 * mom_raw + np.array(
    [0.18, -0.22, 0.10, 0.30, -0.15, 0.05, -0.30, 0.22, -0.12, 0.16])
a = standardize(alpha_raw)
print("candidate standardized exposure a:", np.round(a, 3))
for k in [4, 5, 6]:
    print(f"  corr(a, {factors[k]:5s}) = {np.corrcoef(a, X[:, k])[0, 1]:+.3f}")

# Project a onto the factor space (column space of X): how much is spanned?
coef, *_ = np.linalg.lstsq(X, a, rcond=None)
a_fit = X @ coef
a_perp = a - a_fit
spanned = 1 - a_perp.var(ddof=0) / a.var(ddof=0)
print("signal factor loadings (lstsq):", dict(zip(factors, np.round(coef, 3))))
print(f"spanned fraction (R^2 of a on X): {spanned:.3f}")
print(f"residual (true-alpha) share:      {1 - spanned:.3f}")

# Information coefficient against month-1 returns: raw vs residualized signal
a_perp_z = a_perp / a_perp.std(ddof=0)
ic_raw = np.corrcoef(a, r1)[0, 1]
ic_perp = np.corrcoef(a_perp_z, r1)[0, 1]
print(f"IC raw signal vs r1:   {ic_raw:+.3f}")
print(f"IC residualized vs r1: {ic_perp:+.3f}")

# Decompose the raw signal portfolio's realized return: factor vs specific.
g_sig = a / np.abs(a).sum()        # dollar-neutral long-short, gross = 1
xg = X.T @ g_sig
sig_ret = g_sig @ r1
fac_part = xg @ f1
spec_part = g_sig @ resid
print("signal portfolio exposures:", dict(zip(factors, np.round(xg, 3))))
print(f"raw signal return {sig_ret:+.3f}% = factor {fac_part:+.3f}% "
      f"+ specific {spec_part:+.3f}%")
g_perp = a_perp_z / np.abs(a_perp_z).sum()
print(f"residualized signal return: {g_perp @ r1:+.3f}% "
      f"(factor part {(X.T @ g_perp) @ f1:+.3f}%)")

# ---------------------------------------------------------------------------
# 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}%")

# ---------------------------------------------------------------------------
# 11. Bias statistic on the quarter (chapter 15)
# ---------------------------------------------------------------------------
section("11. BIAS STATISTIC DEMO (chapter 15)")
te_m = sig_a * 100 / np.sqrt(12)      # monthly TE forecast, %
z_bias = np.array(ra_m) / te_m
b_stat = z_bias.std(ddof=1)
print(f"monthly TE forecast: {te_m:.3f}%")
print("standardized active returns z_t:", np.round(z_bias, 3))
print(f"bias statistic b = {b_stat:.2f}   "
      f"(T = 3 acceptance band: [{1 - np.sqrt(2/3):.2f}, {1 + np.sqrt(2/3):.2f}])")

# ---------------------------------------------------------------------------
# 12. CROWD factor case study (chapter 16)
# ---------------------------------------------------------------------------
section("12. CROWD FACTOR CASE STUDY (chapter 16)")
# Stipulated crowding data: hedge-fund ownership share and short interest/float.
hf_own = np.array([0.193, 0.124, 0.125, 0.217, 0.148,
                   0.102, 0.058, 0.090, 0.215, 0.034])
short_int = np.array([0.0483, 0.0311, 0.0311, 0.0541, 0.0369,
                      0.0254, 0.0145, 0.0224, 0.0538, 0.0084])
crowd = standardize(0.6 * standardize(hf_own) + 0.4 * standardize(short_int))
print("CROWD exposures:", dict(zip(names, np.round(crowd, 2))))
coef_c, *_ = np.linalg.lstsq(X, crowd, rcond=None)
r2_aux = 1 - (crowd - X @ coef_c).var(ddof=0) / crowd.var(ddof=0)
print(f"R^2 of CROWD on existing X = {r2_aux:.3f}  ->  VIF = {1/(1-r2_aux):.2f}")
print(f"corr with MOM exposures: {np.corrcoef(crowd, X[:, 5])[0, 1]:+.2f}")
print(f"active CROWD exposure: {crowd @ w_a:+.3f}")

# MiniModel-v2: 8th factor with stipulated vol and correlations; the three
# crowded names' specific vols fall as CROWD explains part of their "specific".
crowd_vol = 0.05
crowd_corr = np.array([0.10, 0.25, -0.10, -0.05, -0.20, 0.30, 0.10])
C8 = np.eye(8)
C8[:7, :7] = C
C8[7, :7] = crowd_corr
C8[:7, 7] = crowd_corr
print("C8 eigenvalues (must be > 0):", np.round(np.linalg.eigvalsh(C8), 3))
vols8 = np.append(vols, crowd_vol)
F8 = np.outer(vols8, vols8) * C8
X8 = np.column_stack([X, crowd])
spec2 = spec_vol.copy()
spec2[[0, 3, 8]] = [16.0, 35.0, 24.0]        # AXIOM, DIGIT, INDIGO
Delta2 = np.diag((spec2 / 100.0) ** 2)
x_a8 = X8.T @ w_a
te_v2 = np.sqrt(x_a8 @ F8 @ x_a8 + w_a @ Delta2 @ w_a)
print(f"TE: v1 {sig_a*100:.2f}%  ->  v2 {te_v2*100:.2f}%")

# Month-1 regression re-run with 8 factors (constraint unchanged, CROWD free).
R8 = np.zeros((8, 7))
R8[:7, :6] = R
R8[7, 6] = 1
XR8 = X8 @ R8
f1_v2 = R8 @ np.linalg.solve(XR8.T @ W @ XR8, XR8.T @ W @ r1)
print("month-1 factor returns, v2 (%):")
for k in range(8):
    print(f"  {(factors + ['CROWD'])[k]:6s} {f1_v2[k]:8.4f}")

# ---------------------------------------------------------------------------
# 13. Risk change attribution: month-1 open -> month-2 open (chapter 11)
# ---------------------------------------------------------------------------
section("13. RISK CHANGE ATTRIBUTION (chapter 11)")
# Between the two reports: (a) month-1 returns reprice everything (drift);
# (b) the manager buys back half the AXIOM active underweight, funded
# pro-rata from the other nine names (the chapter-9 lever); (c) descriptors
# update -- B/P reprices, SIZE = ln(new cap), and the 12-1 momentum window
# rolls one month forward (r1 becomes the *excluded* recent month and does
# not enter; a stipulated month falls off the far end and the previously
# excluded month rolls in); (d) factor vols get one EWMA update from f1
# (half-life 12 months). Correlations and specific risk held fixed.

g1 = 1 + r1 / 100.0
cap2 = cap * g1
cap_w2 = cap2 / cap2.sum()
w_b2 = cap_w2                            # benchmark drifts; no index rebalance
w_drift = w_p * g1 / (1 + (w_p @ r1) / 100.0)

# The trade: halve the drifted AXIOM active underweight, funded pro-rata.
gap_ax = w_drift[0] - w_b2[0]
buy_ax = -0.5 * gap_ax
trade = -buy_ax * w_drift / (1 - w_drift[0])
trade[0] = buy_ax
w_p2 = w_drift + trade
w_a2 = w_p2 - w_b2
print(f"AXIOM active: month-1 {w_a[0]:+.4f} -> drifted {gap_ax:+.4f} -> "
      f"post-trade {w_a2[0]:+.4f}")
print(f"buy {buy_ax:.4f} AXIOM, one-way turnover {np.abs(trade).sum()*50:.1f}%")

# Descriptor updates. Stipulated window-roll months (%): r_drop fell off the
# far end of every stock's 12-1 window, r_add is the previously excluded
# near month that rolled in. 12-1 updated additively (log-return approx.).
r_drop = np.array([-9.0, 1.0, 1.5, -1.0, 3.0, 3.5, 2.5, 0.5, 1.0, 1.5])
r_add = np.array([7.0, 3.0, 1.5, 2.0, -1.0, -1.5, -1.5, 1.5, 1.0, -0.5])
mom_raw2 = mom_raw + (r_add - r_drop) / 100.0
btp2 = btp / g1
size_raw2 = np.log(cap2)

def standardize_w(d, wts):
    z = d - wts @ d
    return z / z.std(ddof=0)

X2 = X.copy()
for k, d2 in [(4, btp2), (5, mom_raw2), (6, size_raw2)]:
    X2[:, k] = standardize_w(d2, cap_w2)

# EWMA vol update from f1; correlations fixed.
lam_v = 2 ** (-1 / 12)
vols2 = np.sqrt(lam_v * vols ** 2 + (1 - lam_v) * 12 * (f1 / 100.0) ** 2)
F2 = np.outer(vols2, vols2) * C
print("factor vols old (% ann):", vols * 100)
print("factor vols new (% ann):", np.round(vols2 * 100, 2))

# --- Exposure change decomposition (exact, bilinear) -----------------------
dw_drift = (w_drift - w_p) - (w_b2 - w_b)   # active-weight change: drift
dw_trade = trade                            # active-weight change: the trade
dw = dw_drift + dw_trade
assert np.allclose(w_a2, w_a + dw)
dX = X2 - X

x_a2 = X2.T @ w_a2
t_drift = X.T @ dw_drift
t_trade = X.T @ dw_trade
t_data = dX.T @ w_a
t_inter = dX.T @ dw
assert np.allclose(x_a2 - x_a, t_drift + t_trade + t_data + t_inter)

print(f"\n{'factor':7s} {'x_a old':>8s} {'drift':>8s} {'trade':>8s} "
      f"{'data':>8s} {'interact':>9s} {'x_a new':>8s}")
for k in range(K):
    print(f"{factors[k]:7s} {x_a[k]:8.3f} {t_drift[k]:8.3f} {t_trade[k]:8.3f} "
          f"{t_data[k]:8.3f} {t_inter[k]:9.3f} {x_a2[k]:8.3f}")

# Split the data term: descriptors re-scored through the OLD cap weights vs
# the frame shift from re-standardizing with NEW cap weights. The frame
# shift is uniform within each column, so it cancels exactly against
# zero-sum active weights (and does NOT cancel for total exposures).
X_mid = X.copy()
for k, d2 in [(4, btp2), (5, mom_raw2), (6, size_raw2)]:
    X_mid[:, k] = standardize_w(d2, cap_w)
assert np.allclose((X2 - X_mid).std(axis=0), 0)   # uniform per column
t_desc = (X_mid - X).T @ w_a
t_frame = (X2 - X_mid).T @ w_a
print("\ndata term split: descriptor", np.round(t_desc[4:], 4),
      " frame", np.round(t_frame[4:], 6), " (styles)")
print("frame shift per style column:", np.round((X2 - X_mid)[0, 4:], 4))
print("frame shift seen by TOTAL portfolio exposures:",
      dict(zip(factors[4:], np.round(((X2 - X_mid).T @ w_p2)[4:], 4))))

# The MOM line, told as the chapter tells it.
kM = 5
x_pre = x_a[kM] + t_drift[kM] + t_trade[kM]
print(f"\nMOM active exposure: month-1 report {x_a[kM]:+.3f}; pre-trade "
      f"calculator promised {x_pre:+.3f}; month-2 report prints {x_a2[kM]:+.3f}")
print(f"the gap: data {t_data[kM]:+.4f}, interaction {t_inter[kM]:+.4f}")
print(f"AXIOM 12-1 momentum {mom_raw[0]:+.2f} -> {mom_raw2[0]:+.2f} "
      f"(dropped {r_drop[0]:+.1f}%, added {r_add[0]:+.1f}%)")
print("MOM z-scores old -> new:")
for i in range(N):
    print(f"  {names[i]:10s} {X[i, kM]:+.3f} -> {X2[i, kM]:+.3f}  "
          f"(delta {X2[i, kM] - X[i, kM]:+.3f})")

# --- TE waterfall (sequential substitution, pipeline order) ----------------
def te_of(wa_s, X_s, F_s):
    xx = X_s.T @ wa_s
    return np.sqrt(xx @ F_s @ xx + wa_s @ Delta @ wa_s)

fwd = [("start (month-1 report)", w_a, X, F),
       ("+ drift", w_a + dw_drift, X, F),
       ("+ trade", w_a2, X, F),
       ("+ exposure updates", w_a2, X2, F),
       ("+ covariance update", w_a2, X2, F2)]
rev = [("start (month-1 report)", w_a, X, F),
       ("+ covariance update", w_a, X, F2),
       ("+ exposure updates", w_a, X2, F2),
       ("+ trade", w_a + dw_trade, X2, F2),
       ("+ drift", w_a2, X2, F2)]
for tag, chain in [("forward (pipeline order)", fwd), ("reverse order", rev)]:
    print(f"\nTE waterfall, {tag}:")
    prev = None
    for lbl, wa_s, X_s, F_s in chain:
        v = te_of(wa_s, X_s, F_s)
        step = "" if prev is None else f"   step {100 * (v - prev):+.2f}"
        print(f"  {lbl:24s} TE {v * 100:6.2f}%{step}")
        prev = v

# Linear (marginal) estimates per forward step vs exact.
sig0 = te_of(w_a + dw_drift, X, F)
mcr0 = Sigma @ (w_a + dw_drift) / sig0
exact_trade = te_of(w_a2, X, F) - sig0
print(f"\ntrade step:      linear {100 * (mcr0 @ dw_trade):+.2f}   "
      f"exact {100 * exact_trade:+.2f}   "
      f"(ch.9 one-name extrapolation MCR_AXIOM*buy = {100 * mcr[0] * buy_ax:+.2f})")
sig1 = te_of(w_a2, X, F)
x1 = X.T @ w_a2
exact_data = te_of(w_a2, X2, F) - sig1
print(f"exposure step:   linear {100 * ((F @ x1) @ (X2.T @ w_a2 - x1) / sig1):+.2f}   "
      f"exact {100 * exact_data:+.2f}")
sig2 = te_of(w_a2, X2, F)
x2v = X2.T @ w_a2
exact_F = te_of(w_a2, X2, F2) - sig2
print(f"covariance step: linear {100 * (x2v @ (F2 - F) @ x2v) / (2 * sig2):+.2f}   "
      f"exact {100 * exact_F:+.2f}")

# The 2x2 first cut: {book} x {model}.
print("\n2x2 first cut (TE %):")
print(f"  old book, old model {te_of(w_a, X, F) * 100:6.2f}    "
      f"old book, new model {te_of(w_a, X2, F2) * 100:6.2f}")
print(f"  new book, old model {te_of(w_a2, X, F) * 100:6.2f}    "
      f"new book, new model {te_of(w_a2, X2, F2) * 100:6.2f}")

# Month-2 attribution shares, for comparison with the month-1 report.
var2 = x_a2 @ F2 @ x_a2 + w_a2 @ Delta @ w_a2
Fx2 = F2 @ x_a2
print(f"\nmonth-2 report: TE {np.sqrt(var2) * 100:.2f}%")
print("factor shares of active variance (%):",
      dict(zip(factors, np.round(100 * x_a2 * Fx2 / var2, 1))))
print(f"specific share: {100 * (w_a2 @ Delta @ w_a2) / var2:.1f}%")
print(f"MOM share of active variance: month-1 "
      f"{100 * x_a[kM] * (F @ x_a)[kM] / sig_a ** 2:.1f}% -> month-2 "
      f"{100 * x_a2[kM] * Fx2[kM] / var2:.1f}%")

# ---------------------------------------------------------------------------
# 14. Hero scatter: month-1 return vs predicted beta (site landing page)
# ---------------------------------------------------------------------------
# Feeds site/src/components/BetaScatter.astro: each stock's predicted beta
# against the cap-weighted benchmark, its month-1 return, and a simple OLS
# line through the ten points. The stock farthest above the line gets the
# "alpha?" mark.
section("14. HERO SCATTER (return vs predicted beta)")
var_b = w_b @ Sigma @ w_b
beta_pred = Sigma @ w_b / var_b
beta_p = w_p @ beta_pred
b_slope, b_int = np.polyfit(beta_pred, r1, 1)
dev = r1 - (b_int + b_slope * beta_pred)
i_alpha = int(np.argmax(dev))
print(f"{'name':10s} {'beta':>7s} {'r1%':>6s} {'dev%':>7s}")
for i in range(N):
    print(f"{names[i]:10s} {beta_pred[i]:7.3f} {r1[i]:6.1f} {dev[i]:7.2f}")
print(f"\nfit through the ten points: r = {b_int:+.3f} + {b_slope:.3f} * beta (%)")
print(f"portfolio predicted beta vs benchmark: {beta_p:.3f}")
print(f"farthest above the line: {names[i_alpha]} ({dev[i_alpha]:+.2f}%)")
```