"""
Realized GARCH (Hansen, Huang & Shek 2012), log-linear specification, from scratch.

Standard GARCH updates the latent variance h_t from the noisy squared return r_{t-1}^2.
Realized GARCH instead lets the informative realized measure RV drive it, and ties the latent h_t
to the observed RV_t through a MEASUREMENT equation with a leverage term:

    return:       r_t = sqrt(h_t) z_t,                          z_t ~ N(0,1)
    GARCH eq:     log h_t = omega + beta log h_{t-1} + gamma log RV_{t-1}
    measurement:  log RV_t = xi + phi log h_t + tau1 z_t + tau2 (z_t^2 - 1) + u_t,   u_t ~ N(0, sigma_u^2)

theta = (omega, beta, gamma, xi, phi, tau1, tau2, log sigma_u).  Effective persistence = beta + phi*gamma.
The likelihood is the JOINT density of returns and realized measures.
"""
import numpy as np
from scipy.optimize import minimize


def _unpack(theta):
    omega, beta, gamma, xi, phi, tau1, tau2, lsu = theta
    return omega, beta, gamma, xi, phi, tau1, tau2, np.exp(lsu)


def filter_logh(theta, r, logRV, logh0):
    omega, beta, gamma, xi, phi, tau1, tau2, su = _unpack(theta)
    T = len(r); logh = np.empty(T); logh[0] = logh0
    for t in range(T - 1):
        logh[t + 1] = omega + beta * logh[t] + gamma * logRV[t]      # RV drives the variance recursion
    return logh


def loglik(theta, r, logRV, logh0):
    omega, beta, gamma, xi, phi, tau1, tau2, su = _unpack(theta)
    if su <= 1e-6 or beta < 0 or beta >= 1 or abs(beta + phi * gamma) >= 1:
        return -np.inf
    logh = filter_logh(theta, r, logRV, logh0)
    h = np.exp(logh); z = r / np.sqrt(h)
    u = logRV - xi - phi * logh - tau1 * z - tau2 * (z ** 2 - 1)
    ll_r = -0.5 * (np.log(2 * np.pi) + logh + r ** 2 / h)           # return density
    ll_x = -0.5 * (np.log(2 * np.pi) + 2 * np.log(su) + u ** 2 / su ** 2)   # measurement density
    ll = ll_r + ll_x
    return np.sum(ll) if np.all(np.isfinite(ll)) else -np.inf


def mle(r, RV):
    """Fit Realized GARCH by joint ML. r = returns, RV = realized variance (same units as r^2)."""
    r = np.asarray(r, float); logRV = np.log(np.asarray(RV, float)); logh0 = float(np.log(np.var(r)))
    x0 = np.array([0.1, 0.55, 0.40, 0.0, 1.0, -0.05, 0.05, np.log(0.4)])
    nll = lambda th: -loglik(th, r, logRV, logh0) if np.isfinite(loglik(th, r, logRV, logh0)) else 1e12
    res = minimize(nll, x0, method="Nelder-Mead", options=dict(xatol=1e-7, fatol=1e-7, maxiter=40000))
    return res.x, logh0


def names():
    return ["omega", "beta", "gamma", "xi", "phi", "tau1", "tau2", "sigma_u"]


def params_readable(theta):
    o, b, g, xi, phi, t1, t2, su = _unpack(theta)
    return dict(omega=o, beta=b, gamma=g, xi=xi, phi=phi, tau1=t1, tau2=t2, sigma_u=su,
                persistence=b + phi * g)