// ===============================================================
// Heston Vol-Arb (lite-C, Zorro)
// - Calibrate Heston (?v, ?) to ~30D slice via MC RMSE grid
// - Buy underpriced calls, Sell overpriced calls by model-vs-mid
// - Uses Zorro option functions (contractXXX) per manual
// ===============================================================

#include <contract.c>   // required for all contract* functions 
#define UNDERLYING     "SPY"
#define DAYS_TARGET    30
#define DELTA_ATM      0.50
#define HISTVOL_DFLT   0.20
#define LOTS           1
#define USE_EUROPEAN   1

// ---- Mispricing thresholds (tune) ----
var TAU_UNDER  = 0.02;   // buy if model >= mid*(1+2%)
var TAU_OVER   = 0.02;   // sell if model <= mid*(1-2%)
int MAX_BUYS   = 1;      // cap # of long calls this bar
int MAX_WRITES = 1;      // cap # of short calls this bar

// ===== Heston params (kappa, theta, v0 anchored; fit sigmaV, rho) =====
var kappa = 2.0;
var theta = 0.04;   // will be reset from ATM IV^2
var v0    = 0.04;   // will be reset from ATM IV^2
var sigmaV = 0.30;  // fitted
var rho    = -0.70; // fitted

// ===== MC controls =====
int  NPaths = 15000;  // raise in live/research
int  NSteps = 32;     // ~1M maturity
int  ANT    = 1;

// ===== globals refreshed per bar =====
var S0 = 0, r = 0, q = 0;  // market inputs

// ===== Math / RNG =====
var max0(var x){ return ifelse(x>0, x, 0); }

int _bm_hasSpare=0; var _bm_spare=0;
var gauss01(){
  if(_bm_hasSpare){ _bm_hasSpare=0; return _bm_spare; }
  var u,v,s;
  do{ u=2.0*random()-1.0; v=2.0*random()-1.0; s=u*u+v*v; }while(s<=1e-12||s>=1.0);
  var m=sqrt(-2.0*log(s)/s);
  _bm_spare=v*m; _bm_hasSpare=1; return u*m;
}

void corrGauss(var* outZ1,var* outZ2){
  var z1=gauss01(), zP=gauss01();
  var z2 = rho*z1 + sqrt(max0(1.0 - rho*rho))*zP;
  *outZ1=z1; *outZ2=z2;
}

// ===== Heston MC (full truncation) =====
var heston_call_MC(var S, var K, var r_, var q_, var T){
  var dt = T/(var)NSteps;
  var sum=0; int n;
  _bm_hasSpare=0; _bm_spare=0;

  int p;
  for(p=0;p<NPaths;p++){
    var S1=S, v1=v0, S2=S, v2=v0; // antithetic
    int i;
    for(i=0;i<NSteps;i++){
      var z1=gauss01(), zperp=gauss01();
      var z2=rho*z1 + sqrt(max0(1.0-rho*rho))*zperp;
      var z1a=-z1, z2a=-z2;

      // path 1
      { var vPos=max0(v1);
        v1 += kappa*(theta - vPos)*dt + sigmaV*sqrt(vPos)*sqrt(dt)*z2;
        if(v1<0) v1=0;
        S1 += (r_ - q_)*S1*dt + sqrt(max0(vPos))*S1*sqrt(dt)*z1;
        if(S1<1e-6) S1=1e-6;
      }
      // path 2 (antithetic)
      { var vPos2=max0(v2);
        v2 += kappa*(theta - vPos2)*dt + sigmaV*sqrt(vPos2)*sqrt(dt)*z2a;
        if(v2<0) v2=0;
        S2 += (r_ - q_)*S2*dt + sqrt(max0(vPos2))*S2*sqrt(dt)*z1a;
        if(S2<1e-6) S2=1e-6;
      }
    }
    var p1=max0(S1-K), p2=max0(S2-K);
    sum += 0.5*(p1+p2);
  }
  return exp(-r_*T)*(sum/(var)NPaths);
}

// ===== Helpers for options chain =====
int select_30d_call_near(var strike){
  int type = CALL; if(USE_EUROPEAN) type |= EUROPEAN;
  return (contract(type, DAYS_TARGET, strike) != 0);  // ?Days and nearest strike  // :contentReference[oaicite:5]{index=5}
}

int step_to_next_strike_same_expiry(){ // walk strikes at same expiry
  CONTRACT* c = contractNext(CALL | (USE_EUROPEAN?EUROPEAN:0), ContractExpiry, ContractStrike);
  if(!c) return 0;
  return 1;
}

// ===== Collect a small calibration set around ATM =====
int build_calib_slice(var S, var frac, int maxPts, var* outK, var* outMid, int* outN){
  // start from ATM strike by delta target
  var kATM = contractStrike(CALL, DAYS_TARGET, S, HISTVOL_DFLT, r, DELTA_ATM);  // :contentReference[oaicite:6]{index=6}
  if(!select_30d_call_near(kATM)) return 0;

  int cnt=0;
  // walk down strikes
  CONTRACT* start = contract(0); start=start; // keep compiler happy
  var kMin = S*(1.0-frac), kMax = S*(1.0+frac);

  // go to lower strikes within range
  while(ContractStrike > kMin){
    CONTRACT* c = contractNext(CALL | (USE_EUROPEAN?EUROPEAN:0), ContractExpiry, -ContractStrike);
    if(!c) break;
  }
  // move up, collect mids
  do{
    if(ContractStrike >= kMin && ContractStrike <= kMax){
      var mid = contractPrice(0);  // fetch bid/ask/underlying; returns mid  // :contentReference[oaicite:7]{index=7}
      if(mid > 0){
        outK[cnt]=ContractStrike; outMid[cnt]=mid; cnt++;
        if(cnt>=maxPts) break;
      }
    }
  } while(step_to_next_strike_same_expiry());

  *outN = cnt;
  return (cnt>=5); // need at least 5 points
}

// ===== ATM IV anchor (uses RQuantLib via manual function) =====
var estimate_ATM_IV(var S){
  var kATM = contractStrike(CALL, DAYS_TARGET, S, HISTVOL_DFLT, r, DELTA_ATM);   // :contentReference[oaicite:8]{index=8}
  if(!select_30d_call_near(kATM)) return 0.0;
  var mid = contractPrice(0); if(mid<=0) return 0.0;                            // :contentReference[oaicite:9]{index=9}
  var T = (var)contractDays(0)/365.0;                                           // :contentReference[oaicite:10]{index=10}
  // contractVol: IV from market price (needs R + RQuantLib)                    // :contentReference[oaicite:11]{index=11}
  return contractVol(0, S, HISTVOL_DFLT, mid, q, r);
}

// ===== Calibrate sigmaV & rho on small grid (RMSE) =====
void calibrate_sigma_rho(var S, var T, var* Ks, var* Mids, int N){
  // anchor theta = v0 = ATM^2
  var ivATM = estimate_ATM_IV(S);
  if(ivATM > 0){
    theta = ivATM*ivATM;
    v0    = theta;
  }

  var bestErr = 1e100, bestSig = sigmaV, bestRho = rho;

  var sMin=0.10, sMax=1.00, sStep=0.10;
  var rMin=-0.90, rMax=-0.05, rStep=0.05;

  var s, rr;
  for(s=sMin; s<=sMax+1e-12; s+=sStep){
    sigmaV = s;
    for(rr=rMin; rr<=rMax+1e-12; rr+=rStep){
      rho = rr;
      // RMSE over the set
      var se=0; int i;
      for(i=0;i<N;i++){
        var mdl = heston_call_MC(S, Ks[i], r, q, T);
        var diff = mdl - Mids[i];
        se += diff*diff;
      }
      var rmse = sqrt(se/ifelse(N>0,N,1));
      if(rmse < bestErr){ bestErr=rmse; bestSig=s; bestRho=rr; }
    }
  }
  sigmaV = bestSig; rho = bestRho;
  printf("\n[Calib] ATM IV=%.3f  theta=v0=%.5f  sigmaV=%.3f  rho=%.2f  RMSE=%.4f",
         (double)ivATM,(double)theta,(double)sigmaV,(double)rho,(double)bestErr);
}

// ===== Trading: scan ~30D calls, buy underpriced, sell overpriced =====
void scan_and_trade_30d_calls(var S){
  // start near ATM
  var kStart = contractStrike(CALL, DAYS_TARGET, S, HISTVOL_DFLT, r, DELTA_ATM); // :contentReference[oaicite:12]{index=12}
  if(!select_30d_call_near(kStart)) return;

  int nBuys=0, nWrites=0;
  int sameExp = ContractExpiry;
  do{
    if(ContractExpiry != sameExp) break; // stay in same maturity bucket
    var mid = contractPrice(0);
    if(mid <= 0) continue;
    var T = (var)contractDays(0)/365.0;

    var model = heston_call_MC(S, ContractStrike, r, q, T);

    // underpriced? buy
    if(nBuys < MAX_BUYS && model >= mid*(1.0 + TAU_UNDER) && !NumOpenLong){
      Lots = LOTS; enterLong();  // buys the selected CALL  // :contentReference[oaicite:13]{index=13}
      nBuys++;
      printf("\nBUY CALL K=%.2f mid=%.3f mdl=%.3f", (double)ContractStrike,(double)mid,(double)model);
    }
    // overpriced? write
    if(nWrites < MAX_WRITES && model <= mid*(1.0 - TAU_OVER) && !NumOpenShort){
      Lots = LOTS; enterShort(); // writes (sells) the selected CALL  // :contentReference[oaicite:14]{index=14}
      nWrites++;
      printf("\nSELL CALL K=%.2f mid=%.3f mdl=%.3f", (double)ContractStrike,(double)mid,(double)model);
    }

  } while(step_to_next_strike_same_expiry());
}

// ===== Zorro hooks =====
function init(){
  set(LOGFILE);
  asset(UNDERLYING);
  BarPeriod = 60;     // 1h
  LookBack  = 60;     // ensure chain loads after warmup (manual)  // :contentReference[oaicite:15]{index=15}
  Multiplier = 100;   // typical US equity options; adjust per market (manual) // :contentReference[oaicite:16]{index=16}
}

function run(){
  if(Bar <= LookBack) return;

  // 1) Load options chain (CALL|PUT) once per day ideally (manual)  // :contentReference[oaicite:17]{index=17}
  if(!contractUpdate(Asset,0,CALL|PUT)) return;

  // 2) Market inputs
  S0 = contractUnderlying(); if(S0<=0) S0 = priceClose();              // :contentReference[oaicite:18]{index=18}
  r  = yield()/100;   // annual risk-free
  q  = 0.00;          // set dividend if needed

  // 3) Build a small calibration slice around ATM (±15% strikes)
  var Ks[32], Mids[32]; int N=0;
  if(!build_calib_slice(S0, 0.15, 32, Ks, Mids, &N)) return;

  var T_bar = (var)DAYS_TARGET/365.0;
  calibrate_sigma_rho(S0, T_bar, Ks, Mids, N);

  // 4) Scan ~30D calls and trade mispricings
  scan_and_trade_30d_calls(S0);
}