"""
horizon.py -- The invariance quest and horizon projection.

Backs the notebooks in  "The Invariance Quest & Horizon Projection".
Follows Meucci, "Risk and Asset Allocation" (Springer 2005), chapter 3-A/B.

THE TWO STEPS OF MARKET MODELLING
---------------------------------
1. INVARIANCE QUEST.  Before you can estimate anything you must find the quantity
   whose distribution repeats identically through time -- the "invariant".  You
   cannot model a stock PRICE (it trends and is perfectly persistent) but you can
   model its log-RETURN; you cannot model a bond YIELD but you can model its daily
   CHANGE.  An invariant is identically and independently distributed (i.i.d.):
   its autocorrelation is ~0 and its distribution is stable across sub-samples.

2. HORIZON PROJECTION.  You estimate the invariant at a short base horizon (say
   daily), but you invest over a long one (say a year).  Because the invariant is
   i.i.d., the horizon quantity is the SUM of T independent copies, whose
   distribution is the T-fold self-convolution.  Its moments scale in a precise
   way and, by the CLT, it drifts toward normal -- but at short-to-medium horizons
   it keeps the skew and fat tails that a naive "square-root-of-time" rule ignores.

CONVENTIONS
-----------
x : (n,) a candidate invariant (e.g. daily log-returns).  T : projection horizon
(number of base periods).  All in the invariant's own units.
"""

import numpy as np
from scipy.stats import skew, kurtosis, norm


# --------------------------------------------------------------------------- #
#  Step 1 -- the invariance quest                                               #
# --------------------------------------------------------------------------- #

def autocorr(x, lags=20):
    """Sample autocorrelation function up to `lags`. An invariant has acf ~ 0."""
    x = np.asarray(x, float) - np.mean(x)
    n = len(x); denom = np.dot(x, x)
    return np.array([1.0 if k == 0 else np.dot(x[:-k], x[k:]) / denom for k in range(lags + 1)])


def ljung_box(x, lags=20):
    """Ljung-Box Q statistic testing joint zero autocorrelation (independence)."""
    n = len(x); ac = autocorr(x, lags)[1:]
    q = n * (n + 2) * np.sum(ac ** 2 / (n - np.arange(1, lags + 1)))
    return q          # compare to chi-square(lags); large Q -> reject i.i.d.


def split_sample_moments(x, k=2):
    """Mean/std of x over k equal consecutive sub-samples. A stable (invariant)
    series has near-identical moments across blocks."""
    x = np.asarray(x, float); parts = np.array_split(x, k)
    return np.array([[p.mean(), p.std()] for p in parts])


# --------------------------------------------------------------------------- #
#  Moments and their horizon scaling                                            #
# --------------------------------------------------------------------------- #

def sample_moments(x):
    """(mean, std, skewness, excess-kurtosis) of the invariant."""
    x = np.asarray(x, float)
    return dict(mean=x.mean(), std=x.std(ddof=1),
                skew=float(skew(x)), exkurt=float(kurtosis(x)))   # kurtosis() is excess by default


def project_moments(m, T):
    """Moments of the sum of T i.i.d. copies of the invariant.
       mean_T = T*mean;  var_T = T*var;  skew_T = skew/sqrt(T);  exkurt_T = exkurt/T.
    The last two shrink toward 0 -- the central-limit drift to normality."""
    return dict(mean=T * m["mean"], std=np.sqrt(T) * m["std"],
                skew=m["skew"] / np.sqrt(T), exkurt=m["exkurt"] / T)


# --------------------------------------------------------------------------- #
#  Step 2 -- horizon projection of the whole distribution                       #
# --------------------------------------------------------------------------- #

def project_density_fft(x, T, bins=2000):
    """EXACT i.i.d. horizon distribution via the characteristic function.
    Bin the invariant into a density p, then the T-fold self-convolution is
    IFFT( FFT(p)^T ) -- the horizon quantity's density on a grid. Returns
    (grid, density)."""
    x = np.asarray(x, float)
    lo, hi = x.min(), x.max()
    dx = (hi - lo) / (bins - 1)
    p, _ = np.histogram(x, bins=bins, range=(lo - dx / 2, hi + dx / 2))
    p = p.astype(float); p /= p.sum()
    L = T * (bins - 1) + 1
    ph = np.fft.rfft(p, n=L)
    dens = np.fft.irfft(ph ** T, n=L)
    dens = np.maximum(dens, 0.0)
    grid = T * lo + np.arange(L) * dx
    dens /= (dens.sum() * dx)
    return grid, dens


def simulate_horizon(x, T, n, rng):
    """Bootstrap: n samples of the horizon sum = sum of T i.i.d. draws of x."""
    idx = rng.integers(0, len(x), size=(n, T))
    return np.asarray(x)[idx].sum(axis=1)


# --------------------------------------------------------------------------- #
#  Quantiles / VaR at the horizon                                               #
# --------------------------------------------------------------------------- #

def cornish_fisher_z(alpha, skew, exkurt):
    """Cornish-Fisher adjusted standard normal quantile for skew & excess
    kurtosis -- a moment-based correction to the Gaussian quantile."""
    z = norm.ppf(alpha)
    return (z + (z ** 2 - 1) * skew / 6
            + (z ** 3 - 3 * z) * exkurt / 24
            - (2 * z ** 3 - 5 * z) * skew ** 2 / 36)


def var_cornish_fisher(m, T, alpha=0.01):
    """Horizon VaR from projected moments + Cornish-Fisher (captures skew/kurtosis)."""
    mT = project_moments(m, T)
    zcf = cornish_fisher_z(alpha, mT["skew"], mT["exkurt"])
    return mT["mean"] + zcf * mT["std"]


def var_sqrt_rule(m, T, alpha=0.01):
    """The naive Gaussian 'square-root-of-time' VaR: normal, moments scaled."""
    return T * m["mean"] + norm.ppf(alpha) * np.sqrt(T) * m["std"]


def var_from_density(grid, dens, alpha=0.01):
    """Quantile of a density given on a grid (for the exact FFT projection)."""
    cdf = np.cumsum(dens) * (grid[1] - grid[0])
    return float(np.interp(alpha, cdf, grid))