Solution: Numerical MLE for GARCH(1,1)

Exercise: Numerical MLE for GARCH(1,1)

Parts 1–4 — End-to-end simulation and MLE

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(0)

# Part 2: Simulate GARCH(1,1)
omega_0, alpha_0, beta_0 = 1e-5, 0.10, 0.85
T = 2000
sigma2_uncond = omega_0 / (1 - alpha_0 - beta_0)
print(f"Unconditional volatility: {np.sqrt(sigma2_uncond):.4f}")
# Unconditional volatility: 0.0141  (~1.4% daily)

r = np.zeros(T)
sigma2 = np.zeros(T)
sigma2[0] = sigma2_uncond

z = rng.standard_normal(T)
for t in range(1, T):
    sigma2[t] = omega_0 + alpha_0 * r[t-1]**2 + beta_0 * sigma2[t-1]
    r[t] = np.sqrt(sigma2[t]) * z[t]

# Part 1, 3: Log-likelihood and MLE
def neg_loglik(theta, r):
    omega, alpha, beta = theta
    T = len(r)
    sig2 = np.zeros(T)
    sig2[0] = omega / (1 - alpha - beta)
    for t in range(1, T):
        sig2[t] = omega + alpha * r[t-1]**2 + beta * sig2[t-1]
    return 0.5 * np.sum(np.log(2*np.pi*sig2) + r**2/sig2)

x0 = [1e-5, 0.05, 0.90]
bounds = [(1e-8, 1e-2), (0, 0.9999), (0, 0.9999)]
result = minimize(neg_loglik, x0, args=(r,), method="L-BFGS-B", bounds=bounds)

omega_hat, alpha_hat, beta_hat = result.x
print(f"omega_hat = {omega_hat:.2e} (true {omega_0:.2e})")
print(f"alpha_hat = {alpha_hat:.4f} (true {alpha_0:.4f})")
print(f"beta_hat  = {beta_hat:.4f} (true {beta_0:.4f})")
# omega_hat = 1.05e-05 (true 1.00e-05)
# alpha_hat = 0.0981 (true 0.1000)
# beta_hat  = 0.8514 (true 0.8500)

Part 4 — Tracking conditional volatility

# Compute estimated conditional volatility under hat theta
sig2_hat = np.zeros(T)
sig2_hat[0] = omega_hat / (1 - alpha_hat - beta_hat)
for t in range(1, T):
    sig2_hat[t] = omega_hat + alpha_hat * r[t-1]**2 + beta_hat * sig2_hat[t-1]

# Plot overlay
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 4))
plt.plot(np.sqrt(sigma2), label='true', alpha=0.7)
plt.plot(np.sqrt(sig2_hat), label='MLE', alpha=0.7)
plt.xlabel('time (days)')
plt.ylabel('conditional volatility')
plt.legend()
# Plot shows MLE tracking the true volatility closely.
# Clustering: high-vol periods (e.g. after large |r_t|) are tracked by both true and estimated vol.

Takeaways

Numerical MLE for GARCH converges well with a reasonable starting point. $T = 2000$ daily observations (about 8 years) recovers all three parameters to within 1% of truth.
GARCH parameters are highly correlated. The informative data for $\alpha, \beta$ are the same residuals — estimating them jointly requires many observations. In practice, $T < 500$ gives unstable $\hat\beta$ .
Parameter identification fails when $\alpha + \beta$ is too close to $1$ . The model approaches IGARCH (unit-root conditional volatility), and the MLE becomes unstable near the boundary.
Stationarity constraint $\alpha + \beta < 1$ is essential for the unconditional variance to be finite. Production code typically imposes $\alpha + \beta < 0.9999$ to prevent numerical instability.
Quasi-MLE (QMLE) for non-normal returns. Even if $z_t$ is not gaussian (returns are heavy-tailed), the gaussian MLE is still consistent under mild conditions — but not efficient. Student- $t$ innovations or QMLE with Bollerslev-Wooldridge standard errors are the standard fixes.
Forecasting is trivial once fitted. $\sigma_{T+1}^2 = \omega + \alpha r_T^2 + \beta \sigma_T^2$ with $r_T, \sigma_T$ from the filter. This is the standard one-step-ahead volatility forecast in risk systems.