ภาคผนวก A–B · Python และ Excel

โค้ดและสูตรในหน้านี้คงไว้ตามต้นฉบับเพื่อให้ทำซ้ำได้ คำอธิบายเป็นภาษาไทย ส่วนชื่อฟังก์ชันและคอมเมนต์ในโค้ดเป็นภาษาอังกฤษตามต้นฉบับ

อ้างอิงต้นฉบับ: หน้า 47–48 · เปิด PDF หน้า 47

ภาคผนวก A · Python reference implementation

ต้นฉบับระบุว่าโค้ดด้านล่างรันได้ทันที (numpy, pandas, scipy) และทำซ้ำผลทั้งหมดของ §9.2 (T=998 เพราะการจัดแนวคลาดหนึ่งวัน; จำนวน breach เท่ากัน) รวมถึง reverse stress test ของ §11.3, liquidity coverage ของ §11.4 และสถิติ drawdown ของบทที่ 13 กฎการออกแบบ: การพยากรณ์ของวัน $t$ ใช้ข้อมูลถึง $t-1$ เท่านั้น; ขอบเขตของเขตใช้ “ต่ำกว่า ppf อย่างเคร่งครัด”

ส่วน __main__ ต้องใช้ไฟล์ sp500.csv จาก workbook ของคอร์สต้นทาง ซึ่งไม่ได้แจกจ่ายในเว็บนี้ ฟังก์ชันที่ไม่ต้องใช้ข้อมูลตลาดถูกเรียกใช้ใน Jupyter Notebook ของเว็บนี้ เพื่อตรวจตัวอย่างคำนวณ

import numpy as np, pandas as pd
from scipy.stats import norm, binom, chi2

def log_returns(px: pd.Series) -> pd.Series:
    return np.log(px).diff().dropna()

def ewma_sigma(r: pd.Series, lam=0.94, seed_n=250) -> pd.Series:
    """sigma_{t|t-1}: forecast for day t using returns up to t-1."""
    v = np.empty(len(r)); v[0] = np.mean(r.values[:seed_n] ** 2)
    for t in range(1, len(r)):
        v[t] = lam * v[t-1] + (1 - lam) * r.values[t-1] ** 2
    return pd.Series(np.sqrt(v), index=r.index)

def cf_factor(p, S, K):
    z = norm.ppf(p)
    return z + (z*z-1)*S/6 + (z**3-3*z)*K/24 - (2*z**3-5*z)*S*S/36

def normal_es_factor(c):
    return -norm.pdf(norm.ppf(1-c)) / (1-c)  # -2.6652 at c = 0.99

def var_forecasts(r, c=0.99, window=250, lam=0.94):
    p, z = 1 - c, norm.ppf(1 - c)
    sig = ewma_sigma(r, lam)
    out = pd.DataFrame(index=r.index)
    out["roll21"] = z * r.rolling(21).std().shift(1)
    out["ewma"]   = z * sig
    out["hs"]     = r.rolling(window).quantile(p).shift(1)
    resid = r / sig  # devolatilise
    out["fhs"]    = resid.rolling(window).quantile(p).shift(1) * sig  # revolatilise
    return out

def backtest(r, var, p=0.01):
    m = var.notna(); hit = (r[m] < var[m]).astype(int).values
    T, x = len(hit), hit.sum(); ph = x / T
    ll = lambda q, k0, k1: (k0*np.log(1-q) if k0 else 0) + (k1*np.log(q) if k1 else 0)
    pof = -2 * (ll(p, T-x, x) - ll(ph, T-x, x))
    a, b = hit[:-1], hit[1:]
    n00, n01 = ((a==0)&(b==0)).sum(), ((a==0)&(b==1)).sum()
    n10, n11 = ((a==1)&(b==0)).sum(), ((a==1)&(b==1)).sum()
    pi0, pi1 = n01/(n00+n01), (n11/(n10+n11) if n10+n11 else 0.0)
    pi = (n01+n11)/(n00+n01+n10+n11)
    ind = -2 * (ll(pi, n00+n10, n01+n11) - ll(pi0, n00, n01) - ll(pi1, n10, n11))
    k250 = hit[-250:].sum()
    g, y = binom.ppf(0.95, 250, p), binom.ppf(0.9999, 250, p)  # 5, 10
    zone = "green" if k250 < g else ("yellow" if k250 < y else "red")
    return dict(T=T, x=int(x), rate=ph, LR_pof=pof, p_pof=chi2.sf(pof, 1),
                LR_ind=ind, p_ind=chi2.sf(ind, 1), cc=pof+ind, p_cc=chi2.sf(pof+ind, 2),
                pi1=pi1, last250=int(k250), zone=zone)

# ---------- Part II: stress testing ----------
def reverse_stress(w, cov, loss):
    """Most plausible scenario x* (min Mahalanobis) with w'x = -loss."""
    w, cov = np.asarray(w, float), np.asarray(cov, float)
    var_p = w @ cov @ w
    x = -loss * cov @ w / var_p
    return x, loss / np.sqrt(var_p)  # scenario, Mahalanobis distance

def coverage(buckets, redemption, horizon, stress_mult=1.0):
    """buckets: list of (share_of_nav, days_to_liquidate)."""
    avail = sum(a for a, d in buckets if d * stress_mult <= horizon)
    return avail / redemption

# ---------- Part III: drawdown ----------
def drawdown_stats(nav: pd.Series) -> dict:
    peak = nav.cummax(); dd = nav / peak - 1
    trough = dd.idxmin(); start = nav.loc[:trough].idxmax()
    after = nav.loc[trough:]; rec = after[after >= nav.loc[start]]
    yrs = (nav.index[-1] - nav.index[0]).days / 365.25
    cagr = (nav.iloc[-1] / nav.iloc[0]) ** (1 / yrs) - 1
    return dict(mdd=dd.min(), peak=start, trough=trough,
                recovery=rec.index[0] if len(rec) else None,
                calmar=cagr / -dd.min(), ulcer=np.sqrt((dd ** 2).mean()))

def mc_mdd(sigma=0.16, mu=0.0, years=3, n_paths=10000, steps=252, seed=1):
    rng = np.random.default_rng(seed)
    lr = (mu - 0.5 * sigma**2) / steps + sigma / np.sqrt(steps) * rng.standard_normal((n_paths, years * steps))
    lv = np.concatenate([np.zeros((n_paths, 1)), np.cumsum(lr, axis=1)], axis=1)
    mdd = 1 - np.exp(lv - np.maximum.accumulate(lv, axis=1))
    return mdd.max(axis=1)  # distribution of MDD

if __name__ == "__main__":
    px = pd.read_csv("sp500.csv", index_col=0, parse_dates=True)["Px"]
    r = log_returns(px)
    f = var_forecasts(r).iloc[251:]
    for col in f:
        print(col, {k: (round(v, 4) if isinstance(v, float) else v)
                    for k, v in backtest(r.loc[f.index], f[col]).items()})
    print(drawdown_stats(px))
    m = mc_mdd(years=3, mu=0.06); print("3y MDD median/95%:", np.median(m).round(3), np.quantile(m, .95).round(3))

ภาคผนวก B · สูตร Excel ที่ใช้บ่อย