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                                   # reuse the from-scratch variance recursion for plotting
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)); print("%d daily returns" % len(r))

def fit_garch(student, draws=800):
    with pm.Model() as m:
        omega   = pm.HalfNormal("omega", 0.5)
        persist = pm.Beta("persist", 20, 1.5)               # alpha+beta = p < 1  (stationarity)
        wshare  = pm.Beta("wshare", 1, 8)                   # alpha's share of p
        alpha = pm.Deterministic("alpha", persist * wshare)
        beta  = pm.Deterministic("beta",  persist * (1 - wshare))
        def step(r2m, hprev, om, al, be):
            return om + al * r2m + be * hprev
        h, _ = pytensor.scan(step, sequences=[pt.as_tensor(r2[:-1])],
                             outputs_info=[pt.as_tensor(s0)], non_sequences=[omega, alpha, 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.9,
                          random_seed=1, progressbar=False)
    return idata

idata_n = fit_garch(False)
idata_t = fit_garch(True)
print("done")

5523 daily returns

C:\Users\user\AppData\Local\Temp\ipykernel_63600\1739718354.py:19: 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])],

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [omega, persist, wshare]

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 97 seconds.

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

C:\Users\user\AppData\Local\Temp\ipykernel_63600\1739718354.py:19: 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])],

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)
Initializing NUTS using jitter+adapt_diag...

def summ(idata, extra=()):
    return az.summary(idata, var_names=["omega", "alpha", "beta", *extra])[["mean", "sd", "r_hat"]]
print("PyMC GARCH(1,1) — Normal innovations");    print(summ(idata_n).to_string())
print("\nPyMC GARCH(1,1) — Student-t innovations"); print(summ(idata_t, extra=["nu"]).to_string())
pn = idata_n.posterior; pt_ = idata_t.posterior
persist_n = float((pn["alpha"] + pn["beta"]).mean()); persist_t = float((pt_["alpha"] + pt_["beta"]).mean())
print("\npersistence  Normal %.4f   Student-t %.4f   nu %.2f" % (persist_n, persist_t, float(pt_["nu"].mean())))
print("cross-check (garch_python from-scratch):  Normal a+b~0.993 ;  t a+b~0.997, nu~6.2")

PyMC GARCH(1,1) — Normal innovations
         mean      sd r_hat
omega  0.0147  0.0026  1.00
alpha  0.0901   0.008  1.00
beta   0.9013  0.0088  1.00

PyMC GARCH(1,1) — Student-t innovations
         mean       sd r_hat
omega  0.0071  0.00175  1.00
alpha  0.0626   0.0068  1.00
beta   0.9326    0.007  1.01
nu        6.4     0.52  1.00

persistence  Normal 0.9914   Student-t 0.9952   nu 6.40
cross-check (garch_python from-scratch):  Normal a+b~0.993 ;  t a+b~0.997, nu~6.2

# conditional volatility (posterior-mean params, reuse from-scratch recursion) + 1% VaR + param posteriors
from scipy.stats import t as tdist
mt = np.array([float(pt_[k].mean()) for k in ["omega", "alpha", "beta"]])
h_t = G.garch_var(mt, r2, s0); vol_t = np.sqrt(h_t); nu_hat = float(pt_["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="steelblue", lw=.7, label="PyMC GARCH(1,1)-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 GARCH(1,1)-t - persistence %.3f, nu %.1f, next-day 1%% VaR %.2f%%" % (mt[1]+mt[2], nu_hat, var99))
for j, p in enumerate(["omega", "alpha", "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")
    a.set_title(p); a.set_yticks([])
    if j == 0: a.legend(fontsize=8)
plt.tight_layout(); plt.savefig("garch_pymc.png", dpi=120, bbox_inches="tight"); plt.show()

GARCH(1,1) — PyMC¶

The Bayesian-Python engine of the trio (from-scratch · PyMC · R)¶

Results¶