import os
_lib = r"C:\Users\user\anaconda3\envs\pymc-env\Library\bin"
if _lib not in os.environ.get("PATH", ""): os.environ["PATH"] = _lib + os.pathsep + os.environ.get("PATH", "")
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
import numpy as np, pandas as pd, pymc as pm, pytensor, pytensor.tensor as pt, arviz as az, matplotlib.pyplot as plt
import garch_scratch as G
d = pd.read_csv("sp500ret.csv", parse_dates=["date"]); r = d["ret"].values - d["ret"].values.mean()
r2 = r**2; s0 = float(np.var(r)); neg = (r < 0).astype("float64"); print("%d daily returns" % len(r))

def fit_gjr(student, draws=800):
    with pm.Model() as m:
        omega   = pm.HalfNormal("omega", 0.5)
        persist = pm.Beta("persist", 20, 1.5)                          # alpha + gamma/2 + beta = p < 1
        sh      = pm.Dirichlet("sh", a=np.array([1.0, 1.0, 8.0]))      # shares of (alpha, gamma/2, beta)
        alpha = pm.Deterministic("alpha", persist * sh[0])
        gamma = pm.Deterministic("gamma", 2.0 * persist * sh[1])
        beta  = pm.Deterministic("beta",  persist * sh[2])
        def step(r2m, negm, hprev, om, al, ga, be):
            return om + (al + ga * negm) * r2m + be * hprev
        h, _ = pytensor.scan(step, sequences=[pt.as_tensor(r2[:-1]), pt.as_tensor(neg[:-1])],
                             outputs_info=[pt.as_tensor(s0)], non_sequences=[omega, alpha, gamma, beta])
        h = pt.concatenate([[pt.as_tensor(s0)], h]); sd = pt.sqrt(h)
        if student:
            nu = pm.Uniform("nu", 2.1, 40.0)
            pm.StudentT("obs", nu=nu, mu=0.0, sigma=sd * pt.sqrt((nu - 2) / nu), observed=r)
        else:
            pm.Normal("obs", 0.0, sd, observed=r)
        idata = pm.sample(draws, tune=draws, chains=2, cores=1, target_accept=0.92,
                          random_seed=1, progressbar=False)
    return idata

idata_n = fit_gjr(False); idata_t = fit_gjr(True); print("done")

5523 daily returns

C:\Users\user\AppData\Local\Temp\ipykernel_2336\422053643.py:20: DeprecationWarning: Scan return signature will change. Updates dict will not be returned, only the first argument. Pass `return_updates=False` to conform to the new API and avoid this warning
  h, _ = pytensor.scan(step, sequences=[pt.as_tensor(r2[:-1]), pt.as_tensor(neg[:-1])],

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [omega, persist, sh]

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)

Sampling 2 chains for 800 tune and 800 draw iterations (1_600 + 1_600 draws total) took 354 seconds.

There were 39 divergences after tuning. Increase `target_accept` or reparameterize.

We recommend running at least 4 chains for robust computation of convergence diagnostics

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

def summ(idata, extra=()):
    return az.summary(idata, var_names=["omega", "alpha", "gamma", "beta", *extra])[["mean", "sd", "r_hat"]]
print("PyMC GJR-GARCH — Normal");    print(summ(idata_n).to_string())
print("\nPyMC GJR-GARCH — Student-t"); print(summ(idata_t, extra=["nu"]).to_string())
gt = idata_t.posterior
print("\ngamma (leverage): Normal %.3f   Student-t %.3f   -> P(gamma>0) = %.3f"
      % (float(idata_n.posterior["gamma"].mean()), float(gt["gamma"].mean()), float((gt["gamma"] > 0).mean())))
print("cross-check (garch_asym_python from-scratch GJR-t): omega .012 alpha .008 gamma .115 beta .924 nu 6.8")

PyMC GJR-GARCH — Normal
         mean      sd r_hat
omega  0.0201  0.0029  1.01
alpha  0.0094  0.0052  1.01
gamma   0.136   0.012  1.01
beta    0.905   0.009  1.01

PyMC GJR-GARCH — Student-t
         mean      sd r_hat
omega  0.0129  0.0024  1.00
alpha    0.01  0.0057  1.01
gamma  0.1166  0.0144  1.00
beta   0.9209  0.0082  1.00
nu       6.86    0.59  1.00

gamma (leverage): Normal 0.136   Student-t 0.117   -> P(gamma>0) = 1.000
cross-check (garch_asym_python from-scratch GJR-t): omega .012 alpha .008 gamma .115 beta .924 nu 6.8

# conditional vol at posterior mean (from-scratch GJR recursion) + gamma/alpha/beta posteriors (N vs t)
from scipy.stats import t as tdist
gt = idata_t.posterior
mt = np.array([float(gt[k].mean()) for k in ["omega", "alpha", "beta", "gamma"]])   # garch_var order: o,a,b,g
h_t = G.garch_var(mt, r2, s0, neg=neg, gjr=True); vol_t = np.sqrt(h_t); nu_hat = float(gt["nu"].mean())
var99 = -np.sqrt(h_t[-1]) * np.sqrt((nu_hat - 2)/nu_hat) * tdist.ppf(0.01, nu_hat)
fig = plt.figure(figsize=(13, 7)); gs = fig.add_gridspec(2, 3, height_ratios=[1.3, 1])
axv = fig.add_subplot(gs[0, :])
axv.plot(d["date"], vol_t, color="firebrick", lw=.7, label="PyMC GJR-GARCH-t conditional vol")
axv.plot(d["date"], np.abs(r), color="0.8", lw=.3, zorder=0, label="|return|")
axv.set_ylim(0, 8); axv.set_ylabel("daily vol (%)"); axv.legend(fontsize=8)
axv.set_title("PyMC GJR-GARCH-t - leverage gamma %.3f, persistence %.3f, nu %.1f, next-day 1%% VaR %.2f%%"
              % (mt[3], mt[1] + mt[3]/2 + mt[2], nu_hat, var99))
for j, p in enumerate(["alpha", "gamma", "beta"]):
    a = fig.add_subplot(gs[1, j])
    a.hist(idata_n.posterior[p].values.ravel(), bins=50, density=True, color="steelblue", alpha=.55, label="Normal")
    a.hist(idata_t.posterior[p].values.ravel(), bins=50, density=True, color="firebrick", alpha=.55, label="Student-t")
    if p == "gamma": a.axvline(0, color="k", lw=1)
    a.set_title(p); a.set_yticks([])
    if j == 0: a.legend(fontsize=8)
plt.tight_layout(); plt.savefig("garch_asym_pymc.png", dpi=120, bbox_inches="tight"); plt.show()

GJR-GARCH — PyMC¶

The PyMC corner of the Phase-2 trio (from-scratch · PyMC · R)¶

Results¶