"""
har_iv_check.py

Independent Python implementation of the HAR-IV fit, driven by the same
deterministic generator as the self test in VolatilityForecastCheck.mq5, so
the two can be compared coefficient by coefficient without exchanging a data
file.

The point is not that Python can do arithmetic. It is that the two
implementations reach the answer by different routes: this one uses numpy's
least squares, while the MQL5 engine applies Householder reflections written
out by hand, because MQL5 has no linear algebra library. If that code were
wrong it would not announce itself. It would return plausible coefficients.
This is how that gets caught.

    python har_iv_check.py

Then run VolatilityForecastCheck.mq5 with InpSelfTest set to true and compare
the two beta lines. Both should report 573 observations.

Reads only. Never trades.
"""

import math

import numpy as np

TRADING_DAYS = 252
NPARAM = 5
H = 5


class Lcg:
    """
    Mirrors the generator in VolatilityForecastCheck.mq5 exactly, including
    the Box-Muller call order. Two uniforms are drawn per normal and only the
    first variate is kept, which wastes half of them and keeps the two
    sequences aligned.
    """

    def __init__(self, seed):
        """Start the sequence. Assumes a non-negative integer seed."""
        self.state = seed & 0x7fffffff

    def next(self):
        """Next uniform variate in [0, 1)."""
        self.state = (1103515245 * self.state + 12345) & 0x7fffffff
        return self.state / 2147483648.0

    def normal(self):
        """Next normal variate, Box-Muller, first value only."""
        u1 = self.next()
        u2 = self.next()
        if u1 < 1e-12:
            u1 = 1e-12
        return math.sqrt(-2.0 * math.log(u1)) * math.cos(6.283185307179586 * u2)


def daily_variance(prev_close, o, h, l, c):
    """Overnight squared plus Garman-Klass. Returns -1 when a price is missing."""
    if min(prev_close, o, h, l, c) <= 0.0:
        return -1.0
    overnight = math.log(o / prev_close)
    ln_hl = math.log(h / l)
    ln_co = math.log(c / o)
    garman = 0.5 * ln_hl ** 2 - (2.0 * math.log(2.0) - 1.0) * ln_co ** 2
    return overnight ** 2 + garman


def trailing(v, i, n):
    """Annualized mean of the n daily variances ending at i, or -1 if incomplete."""
    if i + 1 < n:
        return -1.0
    window = v[i - n + 1:i + 1]
    if any(x < 0.0 for x in window):
        return -1.0
    return sum(window) / n * TRADING_DAYS


def build_series(n=600, seed=11):
    """
    The shared synthetic market. Assumes 600 finished bars, which is what the
    MQL5 self test holds after it discards its newest, still-forming bar.
    """
    rng = Lcg(seed)
    bars, iv = [], []
    price, sigma = 100.0, 0.011
    for _ in range(n):
        sigma = math.exp(0.97 * math.log(sigma) + 0.03 * math.log(0.011)
                         + 0.10 * rng.normal())
        gap = sigma * 0.6 * rng.normal()
        path = sigma * rng.normal()
        opn = price * math.exp(gap)
        cls = opn * math.exp(path)
        wick = abs(sigma * rng.normal())
        bars.append((opn, max(opn, cls) * math.exp(wick),
                     min(opn, cls) * math.exp(-wick), cls))
        iv.append(sigma * sigma * TRADING_DAYS * 1.35)
        price = cls
    return bars, iv


def main():
    """Fit the shared series and print the coefficients for comparison."""
    print()
    print("Python checker for the HAR-IV self test")
    print("=======================================")
    print()

    bars, iv = build_series()
    n = len(bars)

    dv = [-1.0] * n
    for i in range(1, n):
        o, h, l, c = bars[i]
        dv[i] = daily_variance(bars[i - 1][3], o, h, l, c)

    rows, targets = [], []
    for i in range(22, n - H):
        r1 = trailing(dv, i, 1)
        r5 = trailing(dv, i, 5)
        r22 = trailing(dv, i, 22)
        if min(r1, r5, r22, iv[i]) <= 0.0:
            continue
        fwd = sum(dv[i + 1:i + 1 + H]) / H * TRADING_DAYS
        if fwd <= 0.0:
            continue
        rows.append([1.0, math.log(r1), math.log(r5), math.log(r22),
                     math.log(iv[i])])
        targets.append(math.log(fwd))

    X = np.array(rows)
    y = np.array(targets)
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    forecast = float(X[-1] @ beta)

    print("  observations %d" % len(rows))
    print("  beta " + " ".join("%.10f" % b for b in beta))
    print("  fitted log variance, last design row %.10f" % forecast)
    print("  the same, as annualized volatility    %.6f%%"
          % (100.0 * math.sqrt(math.exp(forecast))))
    print()
    print("  Compare against the beta line printed by the MQL5 self test.")
    print("  Agreement to ten decimal places means the hand-written QR in the")
    print("  engine matches a library least squares.")
    print()

    #--- a second check: the same coefficients through the normal equations
    A = X.T @ X
    b = X.T @ y
    beta2 = np.linalg.solve(A, b)
    worst = float(np.max(np.abs(beta - beta2)))
    print("  numpy lstsq against numpy normal equations: %.3e" % worst)
    print("  %s" % ("consistent" if worst < 1e-6 else
                    "INCONSISTENT, the system is badly conditioned"))
    print()
    return 0


if __name__ == "__main__":
    raise SystemExit(main())