# Critical slowing down (Scheffer et al. 2009): as a system nears a tipping point it recovers more
# slowly from small shocks, so lag-1 autocorrelation and variance rise. Test it on AR(1) series:
#   x[t+1] = phi[t] x[t] + noise
#   "approaching": phi rises from 0.3 to 0.95 over 1,000 steps (a system drifting towards instability)
#   "control":     phi stays at 0.3 (nothing is changing)
# Indicator: lag-1 autocorrelation in a moving window of 250 steps; its trend is measured by
# Kendall's tau. A "warning" is declared when tau > 0.5. 200 series of each kind.
import random, statistics as st

def series(n, ramp, rnd):
    x, out = 0.0, []
    for t in range(n):
        phi = 0.3 + (0.65 * t / n if ramp else 0.0)
        x = phi * x + rnd.gauss(0, 1); out.append(x)
    return out

def ac1(w):
    m = st.fmean(w)
    return sum((w[i] - m) * (w[i + 1] - m) for i in range(len(w) - 1)) / sum((u - m) ** 2 for u in w)

def kendall(y):
    c = d = 0
    for i in range(len(y)):
        for j in range(i + 1, len(y)):
            if y[j] > y[i]: c += 1
            elif y[j] < y[i]: d += 1
    return (c - d) / (c + d)

def tau_of(x, win=250, step=10):
    return kendall([ac1(x[i - win:i]) for i in range(win, len(x) + 1, step)])

rnd = random.Random(8)
for ramp, lab in ((True, "approaching"), (False, "control")):
    taus = [tau_of(series(1000, ramp, rnd)) for _ in range(200)]
    print(f"{lab:11s}: median tau {st.median(taus):+.2f}; 'warning' (tau > 0.5) in {sum(t > 0.5 for t in taus) / 2:.0f} % of series")
# With a shorter record the same rule loses power and gains false alarms
for ramp, lab in ((True, "approaching"), (False, "control")):
    taus = [tau_of(series(300, ramp, rnd), win=100, step=5) for _ in range(200)]
    print(f"{lab:11s}, only 300 steps: 'warning' in {sum(t > 0.5 for t in taus) / 2:.0f} % of series")
# approaching: median tau +0.91; 'warning' (tau > 0.5) in 100 % of series
# control    : median tau -0.04; 'warning' (tau > 0.5) in 7 % of series
# approaching, only 300 steps: 'warning' in 91 % of series
# control    , only 300 steps: 'warning' in 14 % of series

# Over-fitting a log-periodic power law (LPPL; Johansen, Ledoit & Sornette 2000) to series that
# contain NO bubble: random walks with an upward drift (a steady boom, 250 trading days).
#   ln p(t) = A + B (tc - t)^m + (tc - t)^m [C1 cos(w ln(tc - t)) + C2 sin(w ln(tc - t))]
# For each (tc, m, w) on a grid, A, B, C1, C2 are fitted by least squares; the best grid point wins.
# Grid: tc 5-117 days after the last observation (steps of 8), m 0.1-0.9, w 6-13 (common filter ranges).
import math, random

def solve(M, v):
    n = len(v); M = [row[:] + [v[i]] for i, row in enumerate(M)]
    for c in range(n):
        piv = max(range(c, n), key=lambda i: abs(M[i][c])); M[c], M[piv] = M[piv], M[c]
        if abs(M[c][c]) < 1e-14: return None
        for i in range(n):
            if i != c:
                f = M[i][c] / M[c][c]
                M[i] = [a - f * b for a, b in zip(M[i], M[c])]
    return [M[i][n] / M[i][i] for i in range(n)]

def lstsq(cols, y):
    k = len(cols)
    G = [[sum(a * b for a, b in zip(cols[i], cols[j])) for j in range(k)] for i in range(k)]
    beta = solve(G, [sum(a * b for a, b in zip(cols[i], y)) for i in range(k)])
    if beta is None: return None, float("inf")
    fit = [sum(beta[i] * cols[i][t] for i in range(k)) for t in range(len(y))]
    return beta, sum((a - b) ** 2 for a, b in zip(y, fit))

def fit_lppl(y):
    T = len(y); ts = range(0, T, 2); yy = [y[t] for t in ts]
    best = (float("inf"), None)
    for tc in range(T + 4, T + 120, 8):
        for m in (0.1, 0.3, 0.5, 0.7, 0.9):
            for w in (6, 7.5, 9, 10.5, 12, 13):
                d = [(tc - t) ** m for t in ts]; L = [math.log(tc - t) for t in ts]
                cols = [[1.0] * len(ts), d, [a * math.cos(w * l) for a, l in zip(d, L)], [a * math.sin(w * l) for a, l in zip(d, L)]]
                _, sse = lstsq(cols, yy)
                if sse < best[0]: best = (sse, (tc - (T - 1), m, w))
    tot = sum((u - sum(yy) / len(yy)) ** 2 for u in yy)
    _, lin = lstsq([[1.0] * len(ts), [float(t) for t in ts]], yy)
    return best[1], 1 - best[0] / tot, 1 - lin / tot

rnd = random.Random(21)
print("series   LPPL R^2   straight-line R^2   'critical time' tc   largest 20-day fall in the next 120 days")
for k in range(8):
    x, y = 0.0, []
    for _ in range(370):
        x += 0.0012 + 0.012 * rnd.gauss(0, 1); y.append(x)       # 30 % a year drift, 19 % a year vol
    (tc, m, w), r2, r2lin = fit_lppl(y[:250])
    fut = y[249:]
    fall = max(1 - math.exp(fut[j] - max(fut[max(0, j - 20):j + 1])) for j in range(len(fut)))
    print(f"   {k + 1}       {r2:.3f}          {r2lin:.3f}            {tc:4d} days later           {fall:5.1%}")
print("No series contains a bubble or a crash mechanism, yet every one yields a fitted 'critical time'.")
# series   LPPL R^2   straight-line R^2   'critical time' tc   largest 20-day fall in the next 120 days
#    1       0.819          0.766              77 days later            7.6%
#    2       0.812          0.269             109 days later            6.3%
#    3       0.441          0.076              45 days later           15.5%
#    4       0.951          0.925             117 days later            8.9%
#    5       0.739          0.041              37 days later            7.2%
#    6       0.960          0.542             117 days later           10.3%
#    7       0.964          0.958              37 days later            9.4%
#    8       0.971          0.867               5 days later            7.3%
# No series contains a bubble or a crash mechanism, yet every one yields a fitted 'critical time'.