"""
FIGARCH (Fractionally Integrated GARCH; Baillie, Bollerslev & Mikkelsen 1996) from scratch.

Long-memory conditional variance: shocks to volatility decay HYPERBOLICALLY (~ j^{d-1}), not
geometrically as in GARCH. FIGARCH(1,d,1):

    phi(L)(1-L)^d eps_t^2 = omega + (1 - beta L) nu_t,   nu_t = eps_t^2 - sigma_t^2

which gives the ARCH(inf) form  sigma_t^2 = omega/(1-beta) + sum_{j>=1} lambda_j eps_{t-j}^2,
with the lambda_j computed from the fractional-differencing weights of (1-L)^d.
d in (0,1) interpolates between GARCH (d=0, short memory) and IGARCH (d=1, infinite persistence).
"""
import numpy as np
from scipy.optimize import minimize
from scipy.special import gammaln


def weights(d, phi, beta, J=1000):
    """ARCH(inf) coefficients lambda_j of FIGARCH(1,d,1) (lambda[0]=0)."""
    pi = np.empty(J + 1); pi[0] = 1.0
    for j in range(1, J + 1):
        pi[j] = pi[j - 1] * (j - 1 - d) / j                 # (1-L)^d fractional-diff weights
    c = pi.copy(); c[1:] -= phi * pi[:-1]                    # (1-phi L)(1-L)^d
    g = np.empty(J + 1); g[0] = c[0]
    for j in range(1, J + 1):
        g[j] = c[j] + beta * g[j - 1]                        # divide by (1-beta L)
    lam = -g; lam[0] = 0.0                                   # lambda(L) = 1 - g(L)
    return lam


def var(r2, omega, d, phi, beta, J=1000):
    lam = weights(d, phi, beta, J)
    mbar = r2.mean()
    padded = np.concatenate([np.full(J, mbar), r2])         # pad pre-sample with unconditional mean
    conv = np.convolve(padded, lam)[J:J + len(r2)]
    return omega / (1.0 - beta) + conv


def _ok(d, phi, beta, lam):
    return (0 < d < 1) and (0 <= phi < 1) and (0 <= beta < 1) and np.all(lam[1:200] >= -1e-8)


def loglik(theta, r2, student=False, J=1000):
    omega, d, phi, beta = theta[0], theta[1], theta[2], theta[3]
    nu = theta[4] if student else np.inf
    if omega <= 0 or nu <= 2:
        return -np.inf
    lam = weights(d, phi, beta, J)
    if not _ok(d, phi, beta, lam):
        return -np.inf
    s2 = omega / (1 - beta) + np.convolve(np.concatenate([np.full(J, r2.mean()), r2]), lam)[J:J + len(r2)]
    if np.any(s2 <= 0) or not np.all(np.isfinite(s2)):
        return -np.inf
    if np.isfinite(nu):
        sc2 = s2 * (nu - 2) / nu
        return np.sum(gammaln((nu + 1) / 2) - gammaln(nu / 2) - 0.5 * np.log(nu * np.pi * sc2)
                      - (nu + 1) / 2 * np.log1p(r2 / (nu * sc2)))
    return -0.5 * np.sum(np.log(2 * np.pi * s2) + r2 / s2)


def mle(r, student=False, J=1000):
    r = np.asarray(r, float); r = r - r.mean(); r2 = r ** 2
    x0 = [0.05 * np.var(r), 0.4, 0.2, 0.6] + ([8.0] if student else [])
    def nll(th):
        ll = loglik(th, r2, student, J)
        return 1e12 if not np.isfinite(ll) else -ll
    res = minimize(nll, np.array(x0), method="Nelder-Mead",
                   options=dict(xatol=1e-7, fatol=1e-7, maxiter=20000))
    return res.x, r2


def names(student=False):
    return ["omega", "d", "phi", "beta"] + (["nu"] if student else [])