#!/usr/bin/env python3
"""
Generatore dei dataset usati dai grafici pgfplots della monografia.

Tutti i numeri del volume sono riproducibili eseguendo questo script:

    python3 figure/genera-dati.py

Parametri di riferimento (dichiarati anche nel testo):
    S      = 6500      punti indice
    sigma  = 12%       volatilita' implicita annualizzata
    r = q  = 0
    M      = 100       moltiplicatore del contratto
    anno   = 252 x 6.5 = 1638 ore di borsa
"""

import math
import os

S0 = 6500.0
SIG = 0.12
M = 100
HRS_YR = 252 * 6.5          # 1638 ore di borsa in un anno
HERE = os.path.dirname(os.path.abspath(__file__))

SQRT2PI = math.sqrt(2.0 * math.pi)


def phi(x):
    return math.exp(-0.5 * x * x) / SQRT2PI


def Ncdf(x):
    return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))


def d1d2(S, K, T, sig):
    a = (math.log(S / K) + 0.5 * sig * sig * T) / (sig * math.sqrt(T))
    return a, a - sig * math.sqrt(T)


def gamma(S, K, T, sig):
    a, _ = d1d2(S, K, T, sig)
    return phi(a) / (S * sig * math.sqrt(T))


def charm(S, K, T, sig):          # r = q = 0
    a, b = d1d2(S, K, T, sig)
    return phi(a) * b / (2.0 * T)


def vanna(S, K, T, sig):          # r = q = 0
    a, b = d1d2(S, K, T, sig)
    return -phi(a) * b / sig


def delta_call(S, K, T, sig):
    a, _ = d1d2(S, K, T, sig)
    return Ncdf(a)


def dollar_gamma_1pct(g, S):
    """Gamma in dollari per movimento dell'1%, per singolo contratto."""
    return g * M * S * S * 0.01


def write(name, header, rows):
    path = os.path.join(HERE, name)
    with open(path, "w") as f:
        f.write("# " + header + "\n")
        for r in rows:
            f.write(" ".join(f"{v:.8g}" for v in r) + "\n")
    print(f"scritto {path}  ({len(rows)} righe)")


# --------------------------------------------------------------------------
# 1. Profilo del gamma in dollari per 1%, in milioni, al variare dello strike
#    per quattro valori di tempo residuo.
# --------------------------------------------------------------------------
ORE = [6.5, 2.0, 0.5, 1.0 / 12.0]      # apertura, meta' seduta, 30', 5'
rows = []
K = 6250.0
while K <= 6750.0001:
    vals = []
    for h in ORE:
        T = h / HRS_YR
        vals.append(dollar_gamma_1pct(gamma(S0, K, T, SIG), S0) / 1e6)
    rows.append([K] + vals)
    K += 1.0
write("gamma_strike.dat",
      "K  gamma$1pct_Mdollari  (6.5h  2h  0.5h  5min)", rows)

# --------------------------------------------------------------------------
# 2. Delta della call al variare dello spot, strike fisso 6500
# --------------------------------------------------------------------------
rows = []
S = 6400.0
while S <= 6600.0001:
    vals = []
    for h in ORE:
        T = h / HRS_YR
        vals.append(delta_call(S, 6500.0, T, SIG))
    rows.append([S] + vals)
    S += 0.5
write("delta_spot.dat", "S  delta  (6.5h  2h  0.5h  5min)", rows)

# --------------------------------------------------------------------------
# 3. Profili di charm e vanna lungo la catena, a 30 minuti dalla scadenza.
#    Normalizzati al proprio massimo assoluto per rendere leggibile il
#    confronto di forma (l'ampiezza relativa e' discussa nel testo).
# --------------------------------------------------------------------------
T30 = 0.5 / HRS_YR
ks, chs, vas, gas = [], [], [], []
K = 6400.0
while K <= 6600.0001:
    ks.append(K)
    chs.append(charm(S0, K, T30, SIG))
    vas.append(vanna(S0, K, T30, SIG))
    gas.append(gamma(S0, K, T30, SIG))
    K += 1.0
mc = max(abs(x) for x in chs)
mv = max(abs(x) for x in vas)
mg = max(abs(x) for x in gas)
write("charm_vanna.dat", "K  charm_norm  vanna_norm  gamma_norm",
      [[ks[i], chs[i] / mc, vas[i] / mv, gas[i] / mg] for i in range(len(ks))])

# --------------------------------------------------------------------------
# 4. Gamma ATM in funzione delle ore residue (scala 1/sqrt(T))
# --------------------------------------------------------------------------
rows = []
h = 6.5
while h >= 0.02:
    T = h / HRS_YR
    rows.append([h, dollar_gamma_1pct(gamma(S0, S0, T, SIG), S0) / 1e6])
    h -= 0.02
write("gamma_atm_tempo.dat", "ore_residue  gamma$1pct_Mdollari", rows)

# --------------------------------------------------------------------------
# 5. Profilo GEX(S) di una catena stilizzata, per illustrare lo zero gamma.
#    Catena: call vendute dai dealer sopra, put comprate dai dealer sotto,
#    con open interest concentrato sugli strike tondi.
# --------------------------------------------------------------------------
# (strike, OI call, OI put) --- catena illustrativa aggregata su piu' scadenze.
# Calibrata per produrre i tre riferimenti usati nel testo:
#   spot 6508 ; zero gamma 6492 ; call wall 6550 ; put wall 6475 ; NetGEX +3,1 mld
CATENA = [
    (6400,  2000, 12000), (6425,  2500,  9000), (6450,  4000, 13000),
    (6475,  5000, 58000), (6500, 16000, 15000), (6525, 22000,  6000),
    (6550, 40000,  3000), (6575, 15000,  1800), (6600, 26000,  1400),
    (6625,  9000,   700), (6650, 14000,   500),
]
# tempo residuo effettivo della catena aggregata: ~3 giorni di borsa
T_PROF = 3.0 / 252.0
SPOT = 6508.0


def gex_profilo(S):
    """GEX netto (miliardi di dollari per 1%) valutato allo spot ipotetico S."""
    tot = 0.0
    for (K, oic, oip) in CATENA:
        g = gamma(S, K, T_PROF, SIG)
        tot += g * (oic - oip) * M * S * S * 0.01
    return tot / 1e9


rows = []
S = 6380.0
while S <= 6670.0001:
    rows.append([S, gex_profilo(S)])
    S += 1.0
write("gex_profilo.dat", "S  GEX_miliardi_per_1pct", rows)

# zero gamma per bisezione
lo, hi = 6380.0, 6670.0
for _ in range(80):
    mid = 0.5 * (lo + hi)
    if gex_profilo(lo) * gex_profilo(mid) <= 0.0:
        hi = mid
    else:
        lo = mid
print(f"  -> zero gamma = {0.5*(lo+hi):.2f} ; NetGEX(spot) = {gex_profilo(SPOT):.3f} mld")

# GEX per strike valutato allo spot corrente (istogramma)
rows = []
for (K, oic, oip) in CATENA:
    g = gamma(SPOT, K, T_PROF, SIG)
    c = g * oic * M * SPOT ** 2 * 0.01 / 1e9
    p = -g * oip * M * SPOT ** 2 * 0.01 / 1e9
    rows.append([K, c, p, c + p])
write("gex_per_strike.dat", "K  GEX_call  GEX_put  GEX_netto  (miliardi)", rows)

# 6. Moltiplicatore di volatilita' del feedback: m = 1/(1 + lambda*Gamma_D)
# --------------------------------------------------------------------------
rows = []
x = -0.95
while x <= 2.0001:
    rows.append([x, 1.0 / (1.0 + x)])
    x += 0.005
write("moltiplicatore.dat", "lambda*Gamma_D  moltiplicatore", rows)


# --------------------------------------------------------------------------
# 7. Profilo GEX con radici multiple (controesempio all'unicita' del flip)
# --------------------------------------------------------------------------
CATENA_MULTI = [
    (6420, 30000,  2000), (6450, 35000,  3000), (6480,  3000, 30000),
    (6500,  4000, 45000), (6520,  5000, 20000), (6560, 40000,  2000),
    (6600, 25000,  1500),
]
T_MULTI = 1.0 / 252.0
rows = []
S = 6360.0
while S <= 6660.0001:
    tot = sum(gamma(S, K, T_MULTI, SIG) * (c - p) * M * S * S * 0.01
              for (K, c, p) in CATENA_MULTI)
    rows.append([S, tot / 1e9])
    S += 1.0
write("gex_multiradice.dat", "S  GEX_miliardi", rows)


# --------------------------------------------------------------------------
# 8. Totali della catena illustrativa (riga di totale della tabella del
#    capitolo sull'identificazione) e sensibilita' alla dinamica dello smile
#    (tabella del capitolo sullo zero gamma). Erano gli unici due blocchi di
#    numeri del volume non generati da questo script.
# --------------------------------------------------------------------------
tot_c = sum(gamma(SPOT, K, T_PROF, SIG) * oic for (K, oic, _) in CATENA) \
        * M * SPOT ** 2 * 0.01 / 1e9
tot_p = -sum(gamma(SPOT, K, T_PROF, SIG) * oip for (K, _, oip) in CATENA) \
        * M * SPOT ** 2 * 0.01 / 1e9
print(f"\n  totali catena: call {tot_c:+.2f} ; put {tot_p:+.2f} ; "
      f"netto {tot_c + tot_p:+.2f} mld")

PENDENZA = -1.2          # skew lineare in log-moneyness


def sigma_smile(K, x, modo):
    """Volatilita' di linea sotto le tre convenzioni di dinamica dello smile."""
    if modo == "piatta":
        return SIG
    if modo == "strike":                       # sticky strike: legata a K
        return SIG + PENDENZA * math.log(K / SPOT)
    return SIG + PENDENZA * math.log(K / x)    # sticky moneyness: legata a K/x


def gex_smile(x, modo):
    tot = 0.0
    for (K, oic, oip) in CATENA:
        tot += gamma(x, K, T_PROF, sigma_smile(K, x, modo)) * (oic - oip)
    return tot * M * x * x * 0.01 / 1e9


print("\n  sensibilita' alla dinamica dello smile"
      f" (skew {PENDENZA}, cioe' "
      f"{100*(sigma_smile(6400, SPOT, 'strike') - sigma_smile(6650, SPOT, 'strike')):.1f}"
      " punti fra 6.400 e 6.650)")
print(f"    {'convenzione':22} {'zero gamma':>11} {'GEX 6450':>10}"
      f" {'GEX 6508':>10} {'GEX 6600':>10}")
righe = []
for modo, etichetta in (("piatta", "volatilita' piatta"),
                        ("strike", "smile sticky strike"),
                        ("moneyness", "smile sticky moneyness")):
    lo, hi = 6380.0, 6670.0
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        if gex_smile(lo, modo) * gex_smile(mid, modo) <= 0.0:
            hi = mid
        else:
            lo = mid
    z = 0.5 * (lo + hi)
    a, b, c = gex_smile(6450, modo), gex_smile(6508, modo), gex_smile(6600, modo)
    righe.append([z, a, b, c])
    print(f"    {etichetta:22} {z:>11.1f} {a:>10.2f} {b:>10.2f} {c:>10.2f}")
write("smile_sensibilita.dat",
      "zero_gamma  GEX_6450  GEX_6508  GEX_6600  (piatta, sticky strike, sticky moneyness)",
      righe)

print("\nfatto.")