make_figures.py

import os
import pickle
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy.stats import norm
from dataclasses import replace
from cppi import (
    Params,
    p_lsf_bs,
    cppi_kou,
)

os.makedirs("figures", exist_ok=True)

mpl.rcParams.update(
    {
        "font.family": "serif",
        "font.size": 10,
        "axes.grid": True,
        "grid.alpha": 0.25,
        "grid.linestyle": "--",
        "axes.spines.top": False,
        "axes.spines.right": False,
        "figure.dpi": 130,
        "savefig.bbox": "tight",
        "savefig.format": "pdf",
        "text.usetex": False,
    }
)

C_PORT = "#1f4e79"
C_FLOOR = "#c0392b"
C_CUSH = "#5b9bd5"
C_EXP = "#2e7d32"
C_KOU = "#d95f02"

p = Params()
with open("results/arrays.pkl", "rb") as f:
    ARR = pickle.load(f)


# CPPI mechanism

print("[fig 1]: CPPI mechanism")
t = ARR["path_t"]
V = ARR["path_V"]
F = ARR["path_F"]
E = ARR["path_E"]

fig, ax = plt.subplots(1, 2, figsize=(11, 3.6))
ax[0].plot(t, V, label=r"$V_t$ (portfolio value)", lw=2, color=C_PORT, zorder=3)
ax[0].plot(t, F, label=r"$F_t$ (floor)", lw=1.8, color=C_FLOOR, ls="--", zorder=2)
ax[0].fill_between(
    t, F, V, where=V >= F, alpha=0.20, color=C_CUSH, label=r"cushion $C_t$", zorder=1
)
ax[0].set_xlabel("time (years)")
ax[0].set_ylabel("value")
ax[0].set_title(r"(a) CPPI dynamics, $m = 4$", fontsize=11)
ax[0].legend(loc="upper left", fontsize=9, framealpha=0.9)

alloc = E / np.maximum(V, 1e-9)
ax[1].fill_between(t, 0, alloc, alpha=0.4, color=C_EXP, label="risky asset")
ax[1].fill_between(t, alloc, 1.0, alpha=0.3, color="#bdbdbd", label="risk-free asset")
ax[1].plot(t, alloc, lw=1.5, color=C_EXP)
ax[1].set_xlabel("time (years)")
ax[1].set_ylabel(r"$e_t / V_t$")
ax[1].set_title("(b) Dynamic allocation", fontsize=11)
ax[1].set_ylim(0, 1.05)
ax[1].legend(loc="upper right", fontsize=9)
plt.tight_layout()
plt.savefig("figures/fig1_mechanism.pdf")
plt.close()


# CPPI vs Buy and Hold vs 100% risk-free

print("[fig 2]: CPPI vs BH vs 100% risky asset")
rng = np.random.default_rng(11)
n_paths = 40
M = p.M
dt = p.dt
t_grid = np.linspace(0, p.T, M + 1)
Z = rng.standard_normal((n_paths, M))
logret = (p.mu - 0.5 * p.sigma**2) * dt + p.sigma * np.sqrt(dt) * Z
S = p.V0 * np.exp(np.cumsum(logret, axis=1))
S = np.concatenate([np.full((n_paths, 1), p.V0), S], axis=1)

V_cppi = np.full((n_paths, M + 1), p.V0)
Vc = np.full(n_paths, p.V0)
for k in range(M):
    t_k = k * dt
    F_ = p.G * np.exp(-p.r * (p.T - t_k))
    C = np.maximum(Vc - F_, 0.0)
    e = np.minimum(p.m * C, Vc)
    cash = Vc - e
    Sret = np.exp(logret[:, k])
    Vc = e * Sret + cash * np.exp(p.r * dt)
    V_cppi[:, k + 1] = Vc

V_bh = 0.5 * S + 0.5 * p.V0 * np.exp(p.r * t_grid)

fig, ax = plt.subplots(1, 3, figsize=(13, 3.8), sharey=True)

# CPPI

for pth in V_cppi:
    ax[0].plot(t_grid, pth, color=C_PORT, alpha=0.35, lw=0.7)
ax[0].axhline(p.G, color="black", ls="--", lw=1.2, label=f"floor G = {p.G:.0f}")
ax[0].plot(t_grid, V_cppi.mean(axis=0), color=C_PORT, lw=2.2, label="mean")
ax[0].set_title(r"CPPI ($m = 4$)", fontsize=11)
ax[0].set_xlabel("years")
ax[0].set_ylabel("portfolio value")
ax[0].legend(fontsize=9)

# Buy and Hold 50/50

for pth in V_bh:
    ax[1].plot(t_grid, pth, color=C_EXP, alpha=0.35, lw=0.7)
ax[1].axhline(p.G, color="black", ls="--", lw=1.2)
ax[1].plot(t_grid, V_bh.mean(axis=0), color=C_EXP, lw=2.2)
ax[1].set_title("Buy and Hold 50/50", fontsize=11)
ax[1].set_xlabel("years")

# 100% risky asset

for pth in S:
    ax[2].plot(t_grid, pth, color=C_KOU, alpha=0.35, lw=0.7)
ax[2].axhline(p.G, color="black", ls="--", lw=1.2)
ax[2].plot(t_grid, S.mean(axis=0), color=C_KOU, lw=2.2)
ax[2].set_title("100% risky asset", fontsize=11)
ax[2].set_xlabel("years")

plt.suptitle(
    "Comparison of 40 simulated paths under the same Gaussian shocks",
    y=1.02,
    fontsize=11,
)
plt.tight_layout()
plt.savefig("figures/fig2_paths.pdf")
plt.close()


# P_LSF and P_BF under BS

print("[fig 3]: P_LSF + P_BF under BS")
rebal_days = np.arange(1, 60)
p_lsf_arr = np.array([p_lsf_bs(p, rb)[0] for rb in rebal_days])
p_bf_arr = np.array([p_lsf_bs(p, rb)[1] for rb in rebal_days])

fig, ax = plt.subplots(1, 2, figsize=(11, 3.6))
ax[0].semilogy(rebal_days, p_lsf_arr, color=C_PORT, lw=1.8)
ax[0].set_xlabel("rebalancing interval (days)")
ax[0].set_ylabel(r"$P^{\mathrm{LSF}} = \mathcal{N}(-d_2)$")
ax[0].set_title("(a) Local breach probability", fontsize=11)
ax[1].semilogy(rebal_days, p_bf_arr, color=C_FLOOR, lw=1.8)
ax[1].set_xlabel("rebalancing interval (days)")
ax[1].set_ylabel(r"$P^{\mathrm{BF}}$")
ax[1].set_title(r"(b) Global probability over $[0, T]$", fontsize=11)
plt.suptitle(
    r"Black-Scholes analytical formula: $m = 4$, $\sigma = 20\%$, $T = 5$ years",
    y=1.02,
    fontsize=11,
)
plt.tight_layout()
plt.savefig("figures/fig3_plsf_bs.pdf")
plt.close()


# Heston stress

print("[fig 4]: Heston stress")
rng = np.random.default_rng(5)
v0_h, kappa_h, theta_h, xi_h, rho_h = 0.04, 0.5, 0.09, 1.0, -0.9
n_paths = 4
S_h = np.full(n_paths, p.V0)
v_h = np.full(n_paths, v0_h)
Ss_h = np.zeros((n_paths, M + 1))
Vs_h = np.zeros((n_paths, M + 1))
Ss_h[:, 0] = S_h
Vs_h[:, 0] = v_h

for k in range(M):
    Z1 = rng.standard_normal(n_paths)
    Z2 = rng.standard_normal(n_paths)
    dW1 = np.sqrt(dt) * Z1
    dW2 = np.sqrt(dt) * (rho_h * Z1 + np.sqrt(1 - rho_h**2) * Z2)
    v_pos = np.maximum(v_h, 0)
    v_h = v_h + kappa_h * (theta_h - v_pos) * dt + xi_h * np.sqrt(v_pos) * dW2
    v_h = np.maximum(v_h, 0)
    S_h = S_h * np.exp((p.mu - 0.5 * v_pos) * dt + np.sqrt(v_pos) * dW1)
    Ss_h[:, k + 1] = S_h
    Vs_h[:, k + 1] = v_h

cols4 = ["#1f4e79", "#d95f02", "#2e7d32", "#7b3294"]
fig, ax = plt.subplots(1, 2, figsize=(11, 3.6))
for i in range(n_paths):
    ax[0].plot(
        t_grid, Ss_h[i], lw=1.3, color=cols4[i], alpha=0.9, label=f"path {i + 1}"
    )
ax[0].set_title(r"(a) Risky asset $S_t$ under Heston stress", fontsize=11)
ax[0].set_xlabel("years")
ax[0].set_ylabel(r"$S_t$")
ax[0].legend(fontsize=8, loc="upper left")
for i in range(n_paths):
    ax[1].plot(t_grid, np.sqrt(Vs_h[i]) * 100, lw=1.3, color=cols4[i], alpha=0.9)
ax[1].axhline(
    np.sqrt(theta_h) * 100,
    color="black",
    ls="--",
    lw=1,
    label=r"long-term vol $\sqrt{\theta} = 30\%$",
)
ax[1].set_title(r"(b) Instantaneous volatility $\sqrt{v_t}$", fontsize=11)
ax[1].set_xlabel("years")
ax[1].set_ylabel("vol (%)")
ax[1].legend(fontsize=9)
plt.tight_layout()
plt.savefig("figures/fig4_heston.pdf")
plt.close()


# BS vs Kou distributions + theoretical/empirical curves under Kou

print("[fig 5]: Kou distributions + theoretical/empirical curves under Kou")
df_te = pd.read_csv("results/table_kou_theo_emp.csv")
V_bs = ARR["bs_rb1"]
V_kou = ARR["kou_base"]

fig, ax = plt.subplots(1, 2, figsize=(11, 3.8))
bins = np.linspace(40, 400, 80)
ax[0].hist(
    np.clip(V_bs, 40, 400),
    bins=bins,
    color=C_PORT,
    alpha=0.55,
    edgecolor="white",
    lw=0.3,
    label="Black-Scholes daily",
)
ax[0].hist(
    np.clip(V_kou, 40, 400),
    bins=bins,
    color=C_KOU,
    alpha=0.65,
    edgecolor="white",
    lw=0.3,
    label="Kou daily",
)
ax[0].axvline(p.G, color="black", ls="--", lw=1.3, label=f"G = {p.G:.0f}")
ax[0].set_xlabel(r"$V_T$")
ax[0].set_ylabel("frequency")
ax[0].set_title(r"(a) Terminal distribution, $m = 4$", fontsize=11)
ax[0].legend(fontsize=9, loc="upper right")
ax[1].plot(
    df_te["m"],
    df_te["PBF theoretical"].str.rstrip("%").astype(float),
    "o-",
    color=C_PORT,
    lw=1.8,
    label="theoretical formula (upper bound)",
    markersize=7,
)
ax[1].plot(
    df_te["m"],
    df_te["PBF empirical"].str.rstrip("%").astype(float),
    "s-",
    color=C_KOU,
    lw=1.8,
    label="empirical simulation",
    markersize=7,
)
ax[1].set_xlabel(r"multiplier $m$")
ax[1].set_ylabel(r"$P^{\mathrm{BF}}$ (%)")
ax[1].set_title("(b) Theoretical vs empirical", fontsize=11)
ax[1].legend(fontsize=9, loc="upper left")
plt.tight_layout()
plt.savefig("figures/fig5_kou.pdf")
plt.close()


# Fig 6: heatmap

print("[fig 6]: heatmap")
ms_grid = [3, 4, 5, 6, 8]
lams_grid = [0.3, 0.5, 1.0, 2.0]
heat = np.zeros((len(ms_grid), len(lams_grid)))
p_small = replace(p, N_mc=3_000)
for i, mval in enumerate(ms_grid):
    for j, lam in enumerate(lams_grid):
        p_m = replace(p_small, m=float(mval))
        Vk = cppi_kou(p_m, rebal_every=1, lam=lam)
        heat[i, j] = 100 * np.mean(Vk < p.G)

fig, ax = plt.subplots(figsize=(7.5, 4.2))
im = ax.imshow(
    heat, aspect="auto", cmap="Reds", origin="lower", vmin=0, vmax=heat.max()
)
ax.set_xticks(range(len(lams_grid)))
ax.set_xticklabels(lams_grid)
ax.set_yticks(range(len(ms_grid)))
ax.set_yticklabels(ms_grid)
ax.set_xlabel(r"jump intensity $\lambda$ (per year)", fontsize=11)
ax.set_ylabel(r"multiplier $m$", fontsize=11)
ax.set_title(r"$P(V_T < G)$ under Kou (%), daily rebalancing", fontsize=11)
for i in range(len(ms_grid)):
    for j in range(len(lams_grid)):
        val = heat[i, j]
        ax.text(
            j,
            i,
            f"{val:.1f}",
            ha="center",
            va="center",
            color="white" if val > heat.max() / 2 else "black",
            fontsize=9.5,
            fontweight="bold",
        )
plt.colorbar(im, ax=ax, label="%")
plt.tight_layout()
plt.savefig("figures/fig6_heatmap_kou.pdf")
plt.close()


# Dynamic multiplier

print("[fig 7]: dynamic multiplier")
V_m4 = ARR["V_m4_kou"]
V_cap4 = ARR["V_cap4"]
V_cap12 = ARR["V_cap12"]
mh_cap4 = ARR["m_hist_cap4"]
mh_cap12 = ARR["m_hist_nocap"]  # alias
V_fix4 = ARR["V_fix4"]
mh_fix4 = ARR["m_hist_fix4"]

fig, ax = plt.subplots(1, 2, figsize=(11, 3.8))
bins = np.linspace(40, 400, 70)
ax[0].hist(
    np.clip(V_m4, 40, 400),
    bins=bins,
    color=C_PORT,
    alpha=0.55,
    edgecolor="white",
    lw=0.3,
    label=r"$m = 4$ constant",
)
ax[0].hist(
    np.clip(V_cap4, 40, 400),
    bins=bins,
    color=C_EXP,
    alpha=0.55,
    edgecolor="white",
    lw=0.3,
    label=r"VaR$(T\!-\!t)$, cap $= 4$",
)
ax[0].hist(
    np.clip(V_fix4, 40, 400),
    bins=bins,
    color=C_FLOOR,
    alpha=0.45,
    edgecolor="white",
    lw=0.3,
    label=r"VaR$(\tau\!=\!1/12)$, cap $= 4$",
)
ax[0].axvline(p.G, color="black", ls="--", lw=1.3)
ax[0].set_xlabel(r"$V_T$")
ax[0].set_ylabel("frequency")
ax[0].set_title("(a) Distributions under Kou", fontsize=11)
ax[0].legend(fontsize=8, loc="upper right")

t_mh = np.arange(len(mh_cap4)) * dt
ax[1].plot(t_mh, mh_cap4, color=C_EXP, lw=1.3, label=r"VaR$(T\!-\!t)$, cap $= 4$")
ax[1].plot(
    t_mh, mh_fix4, color=C_FLOOR, lw=1.3, label=r"VaR$(\tau\!=\!1/12)$, cap $= 4$"
)
ax[1].axhline(4, color=C_PORT, ls="--", lw=1.2, label=r"$m = 4$ constant")
ax[1].set_xlabel("years")
ax[1].set_ylabel(r"$m_t$")
ax[1].set_title(r"(b) $m_t$ trajectory on a single path", fontsize=11)
ax[1].legend(fontsize=8, loc="upper left")
plt.tight_layout()
plt.savefig("figures/fig7_dyn_m.pdf")
plt.close()


# Put hedging

print("[fig 8]: put hedging")
V_naked = ARR["V_naked_base"]
V_hedge = ARR["V_hedge_base"]
MATERIALITY = 1.0

fig, ax = plt.subplots(1, 2, figsize=(11, 3.8), sharey=True)
bins = np.linspace(40, 400, 70)

pbf_naked = 100 * np.mean(V_naked < p.G - MATERIALITY)
ax[0].hist(
    np.clip(V_naked, 40, 400),
    bins=bins,
    color=C_FLOOR,
    alpha=0.7,
    edgecolor="white",
    lw=0.3,
)
ax[0].axvline(p.G, color="black", ls="--", lw=1.3)
ax[0].set_title(rf"(a) No hedge, $P(V_T < G-1)$ = {pbf_naked:.1f}%", fontsize=11)
ax[0].set_xlabel(r"$V_T$")
ax[0].set_ylabel("frequency")

pbf_hedge = 100 * np.mean(V_hedge < p.G - MATERIALITY)
ax[1].hist(
    np.clip(V_hedge, 40, 400),
    bins=bins,
    color=C_PORT,
    alpha=0.7,
    edgecolor="white",
    lw=0.3,
)
ax[1].axvline(p.G, color="black", ls="--", lw=1.3)
ax[1].set_title(rf"(b) With OTM puts, $P(V_T < G-1)$ = {pbf_hedge:.1f}%", fontsize=11)
ax[1].set_xlabel(r"$V_T$")

plt.suptitle(
    r"Impact of put hedging under Kou "
    r"($m = 4$, $\sigma_{\mathrm{hedge}} = 25\%$)",
    y=1.02,
    fontsize=11,
)
plt.tight_layout()
plt.savefig("figures/fig8_hedging.pdf")
plt.close()

print("\nAll figures generated in figures/")