run_all_tables.py

import os
import pickle
import numpy as np
import pandas as pd
from dataclasses import replace
from cppi import (
    Params,
    one_path_bs,
    bs_closed_form,
    cppi_bs_discrete,
    cppi_heston,
    cppi_kou,
    cppi_kou_var_ewma,
    cppi_kou_hedged,
    p_lsf_bs,
    p_bf_kou,
    metrics,
)

os.makedirs("results", exist_ok=True)
p = Params()
ARR = {}


# Table 1: impact of rebalancing frequency under BS

print("=" * 60)
print("TABLE 1: discrete BS, impact of rebalancing frequency")
print("=" * 60)
rows = []
for label, rebal in [
    ("Daily (N=1260)", 1),
    ("Weekly (N=252)", 5),
    ("Monthly (N=60)", 21),
]:
    V = cppi_bs_discrete(p, rebal_every=rebal)
    plsf, pbf_t = p_lsf_bs(p, rebal)
    rows.append(
        {
            "Frequency": label,
            "E[V_T]": f"{V.mean():.2f}",
            "Std(V_T)": f"{V.std():.2f}",
            "PBF empirical": f"{np.mean(V < p.G) * 100:.3f}%",
            "PLSF theoretical": f"{plsf:.2e}",
            "PBF theoretical": f"{pbf_t:.2e}",
        }
    )
df1 = pd.DataFrame(rows)
print(df1.to_string(index=False))
ARR["bs_rb1"] = cppi_bs_discrete(p, rebal_every=1)


# Table 2: Heston standard vs stress

print("\n" + "=" * 60)
print("TABLE 2: Heston")
print("=" * 60)
V_bs = cppi_bs_discrete(p, rebal_every=1)
V_hst, rv_hst = cppi_heston(p, v0=0.04, kappa=2.0, theta=0.04, xi=0.4, rho=-0.7)
V_hst_s, rv_hst_s = cppi_heston(p, v0=0.04, kappa=0.5, theta=0.09, xi=1.0, rho=-0.9)
for label, V, rv in [
    ("Daily BS (sigma=20%)", V_bs, 0.20),
    ("Heston standard", V_hst, rv_hst.mean()),
    ("Heston stress", V_hst_s, rv_hst_s.mean()),
]:
    pbf = np.mean(V < p.G) * 100
    print(
        f"{label:35s} | E={V.mean():.2f} | std={V.std():.2f} | "
        f"PBF={pbf:.2f}% | RV={rv:.2f}"
    )
ARR["heston_std"] = V_hst
ARR["heston_stress"] = V_hst_s


# Table 3: Kou theoretical vs empirical

print("\n" + "=" * 60)
print("TABLE 3: Kou theoretical vs empirical")
print("=" * 60)
p_table = replace(p, N_mc=10_000)
m_list = [3, 4, 5, 6, 8, 10]
rows_kou = []
for m_val in m_list:
    p_m = replace(p_table, m=float(m_val))
    pbf_t = p_bf_kou(p_m, lam=1.0, p_neg=0.6, eta_minus=10.0) * 100
    V_kou = cppi_kou(p_m, rebal_every=1)
    pbf_e = np.mean(V_kou < p.G) * 100
    rows_kou.append(
        {
            "m": m_val,
            "PBF theoretical": f"{pbf_t:.2f}%",
            "PBF empirical": f"{pbf_e:.2f}%",
            "Gap (pts)": f"+{pbf_t - pbf_e:.2f}",
        }
    )
df3 = pd.DataFrame(rows_kou)
print(df3.to_string(index=False))
df3.to_csv("results/table_kou_theo_emp.csv", index=False)

# Base simulation m=4 with N=20000 for fig 5
ARR["kou_base"] = cppi_kou(p, rebal_every=1)


# Table 4: summary of regimes

print("\n" + "=" * 60)
print("TABLE 4: summary of regimes")
print("=" * 60)
V_bs_cf = bs_closed_form(p)
V_kou_m4 = ARR["kou_base"]
for label, V in [
    ("Continuous BS (closed form)", V_bs_cf),
    ("Discrete BS daily", V_bs),
    ("Heston standard", V_hst),
    ("Heston stress", V_hst_s),
    ("Kou daily m=4", V_kou_m4),
]:
    m_obj = metrics(V, p.G, label)
    pbf = m_obj["P(V_T<G) % raw"]
    es_c = m_obj["ES_cond"]
    es_t = m_obj["ES_tail5%"]
    print(
        f"{label:30s} | E={m_obj['E[V_T]']:7.2f} | std={m_obj['Std(V_T)']:7.2f} | "
        f"PBF={pbf:6.2f}% | ES_c={es_c:.2f} | ES_t={es_t:.2f}"
    )


# Table 5: dynamic VaR multiplier (residual horizon) vs fixed horizon

print("\n" + "=" * 60)
print("TABLE 5a: dynamic VaR multiplier (tau=T-t), various caps")
print("=" * 60)
p_dyn = replace(p, N_mc=10_000)

# Constant references
V_m3_kou = cppi_kou(replace(p_dyn, m=3.0), rebal_every=1)
V_m4_kou = cppi_kou(p_dyn, rebal_every=1)
ARR["V_m4_kou"] = V_m4_kou
ARR["V_m3_kou"] = V_m3_kou

# VaR with residual horizon (T-t)
V_cap3, mh_cap3 = cppi_kou_var_ewma(p_dyn, m_max=3.0)
V_cap4, mh_cap4 = cppi_kou_var_ewma(p_dyn, m_max=4.0)
V_cap12, mh_cap12 = cppi_kou_var_ewma(p_dyn, m_max=12.0)
ARR["V_cap3"] = V_cap3
ARR["V_cap4"] = V_cap4
ARR["V_cap12"] = V_cap12
ARR["m_hist_cap4"] = mh_cap4
ARR["m_hist_cap12"] = mh_cap12

for label, V, mh in [
    ("m=3 constant", V_m3_kou, np.full(1, 3.0)),
    ("m=4 constant", V_m4_kou, np.full(1, 4.0)),
    ("VaR(T-t), m_max=3", V_cap3, mh_cap3),
    ("VaR(T-t), m_max=4", V_cap4, mh_cap4),
    ("VaR(T-t), m_max=12", V_cap12, mh_cap12),
]:
    m_obj = metrics(V, p.G, label)
    pbf = m_obj["P(V_T<G) % raw"]
    es_c = m_obj["ES_cond"]
    es_t = m_obj["ES_tail5%"]
    print(
        f"{label:25s} | E={m_obj['E[V_T]']:7.2f} | PBF={pbf:6.2f}% | "
        f"ES_c={es_c:.2f} | ES_t={es_t:.2f} | m_avg={np.mean(mh):.3f}"
    )

# Table 5b
print("\nTABLE 5b: residual horizon vs fixed horizon tau=1/12")
# Fixed horizon tau=1/12: passed via tau_fixed parameter only
# Local helper that fixes tau
from scipy.stats import norm as _norm


def _cppi_kou_var_fixed_tau(
    p,
    alpha=0.995,
    tau_fixed=1 / 12,
    lambda_ewma=0.94,
    m_max=None,
    m_min=1.0,
    lam=1.0,
    p_neg=0.6,
    eta_up=25.0,
    eta_dn=10.0,
    seed=42,
):

    from cppi import _kou_compensator

    rng = np.random.default_rng(seed)
    comp = _kou_compensator(lam, p_neg, eta_up, eta_dn)
    z_alpha = _norm.ppf(alpha)
    V = np.full(p.N_mc, p.V0)
    e = np.zeros(p.N_mc)
    cash = np.zeros(p.N_mc)
    sigma_ewma = np.full(p.N_mc, p.sigma**2 * p.dt)
    prev_logret = np.zeros(p.N_mc)
    mh = []
    for k in range(p.M):
        t = k * p.dt
        if k > 0:
            sigma_ewma = lambda_ewma * sigma_ewma + (1 - lambda_ewma) * prev_logret**2
        sigma_t = np.sqrt(sigma_ewma / p.dt)
        tau = tau_fixed  # fixed horizon
        exponent = (p.mu - p.r - 0.5 * sigma_t**2) * tau - sigma_t * z_alpha * np.sqrt(
            tau
        )
        denom = np.where(1 - np.exp(exponent) > 1e-8, 1 - np.exp(exponent), 1e-8)
        m_t = np.clip(1.0 / denom, m_min, m_max if m_max is not None else 1e9)
        mh.append(float(m_t[0]))
        F = p.G * np.exp(-p.r * (p.T - t))
        C = np.maximum(V - F, 0.0)
        e_new = np.minimum(m_t * C, V)
        tc_cost = p.tc * np.abs(e_new - e)
        V = V - tc_cost
        e = e_new
        cash = V - e
        Z = rng.standard_normal(p.N_mc)
        logret = (p.mu - 0.5 * p.sigma**2 - comp) * p.dt + p.sigma * np.sqrt(p.dt) * Z
        Nj = rng.poisson(lam * p.dt, size=p.N_mc)
        mask = Nj > 0
        if mask.any():
            nj = int(mask.sum())
            up = rng.random(nj) < (1 - p_neg)
            J = np.where(
                up, rng.exponential(1 / eta_up, nj), -rng.exponential(1 / eta_dn, nj)
            )
            logret[mask] += J
        S_ret = np.exp(logret)
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        e = e * S_ret
        cash = cash * np.exp(p.r * p.dt)
        prev_logret = logret
    return V, np.array(mh)


V_fix4, mh_fix4 = _cppi_kou_var_fixed_tau(p_dyn, m_max=4.0)
V_fix12, mh_fix12 = _cppi_kou_var_fixed_tau(p_dyn, m_max=12.0)
ARR["V_fix4"] = V_fix4
ARR["V_fix12"] = V_fix12
ARR["m_hist_fix4"] = mh_fix4
ARR["m_hist_nocap"] = mh_cap12  # same m_hist as cap12, renamed for table 5b

for label, V, mh in [
    ("VaR(T-t), m_max=4", V_cap4, mh_cap4),
    ("VaR(1/12), m_max=4", V_fix4, mh_fix4),
    ("VaR(T-t), m_max=12", V_cap12, mh_cap12),
    ("VaR(1/12), m_max=12", V_fix12, mh_fix12),
]:
    pbf = np.mean(V < p.G) * 100
    print(
        f"{label:25s} | E={V.mean():7.2f} | PBF={pbf:6.2f}% | m_avg={np.mean(mh):.3f}"
    )


# Table 6: OTM put hedging under Kou

print("\n" + "=" * 60)
print("TABLE 6: OTM put hedging under Kou")
print("=" * 60)
p_hedge = replace(p, N_mc=10_000)
V_naked_base = cppi_kou(p_hedge, rebal_every=1)
V_hedge_base = cppi_kou_hedged(p_hedge, rebal_every=1)
V_naked_fat = cppi_kou(p_hedge, rebal_every=1, eta_dn=5.0)
V_hedge_fat = cppi_kou_hedged(p_hedge, rebal_every=1, eta_dn=5.0)
ARR["V_naked_base"] = V_naked_base
ARR["V_hedge_base"] = V_hedge_base
ARR["V_naked_fat"] = V_naked_fat
ARR["V_hedge_fat"] = V_hedge_fat

MATERIALITY = 1.0
for label, V in [
    ("Kou base, no hedge", V_naked_base),
    ("Kou base, with puts", V_hedge_base),
    ("Kou fat tails, no hedge", V_naked_fat),
    ("Kou fat tails, with puts", V_hedge_fat),
]:
    m_obj = metrics(V, p.G, label, materiality=MATERIALITY)
    pbf_b = m_obj["P(V_T<G) % raw"]
    pbf_m = m_obj["P(V_T<G) % material"]
    es_c = m_obj["ES_cond"]
    es_t = m_obj["ES_tail5%"]
    print(
        f"{label:30s} | E={m_obj['E[V_T]']:7.2f} | "
        f"PBF_raw={pbf_b:6.2f}% | PBF_mat={pbf_m:6.2f}% | "
        f"ES_c={es_c:.2f} | ES_t={es_t:.2f}"
    )


# Single path

t, Vs, Fs, Es, Ss = one_path_bs(p, seed=2)
ARR["path_t"] = t
ARR["path_V"] = Vs
ARR["path_F"] = Fs
ARR["path_E"] = Es


# Save

with open("results/arrays.pkl", "wb") as f:
    pickle.dump(ARR, f)

print("Saved to results/arrays.pkl -> OK")