//+------------------------------------------------------------------+
//|                                                VolEstimators.mqh |
//|        Garman-Klass / Parkinson / Rogers-Satchell / Yang-Zhang   |
//+------------------------------------------------------------------+
//| sigma2_GK = 0.5*(ln(H/L))^2 - (2ln2 - 1)*(ln(C/O))^2              |
//|                                                                  |
//| NON-NEGATIVE BY CONSTRUCTION for any valid bar: H >= max(O,C) and |
//| L <= min(O,C) give ln(H/L) >= |ln(C/O)|, so                       |
//|     GK >= (0.5 - 0.3863)*ln(C/O)^2 = 0.1137*ln(C/O)^2 >= 0.       |
//| Do NOT add a clamp and call it a safety net: on a CORRUPT bar     |
//| (H < L, or an EMPTY_VALUE reading as DBL_MAX) GK does not go      |
//| negative - it returns a plausible POSITIVE variance. The only     |
//| defence is validating O/H/L/C, which is what GKValid() does and   |
//| what every entry point here calls first.                          |
//|                                                                  |
//| GK IS BLIND TO GAPS. It assumes the bar opened where the previous |
//| one closed. Index CFDs gap across the daily maintenance break and |
//| the weekend, so on D1 GK UNDERSTATES the variance a multi-day     |
//| hold is actually exposed to. For sizing an intra-week hold use    |
//| YangZhangSigma(), which prices the gap explicitly.                |
//+------------------------------------------------------------------+
#ifndef VOLESTIMATORS_MQH
#define VOLESTIMATORS_MQH

#define VOL_LN2 0.69314718055994531

//+------------------------------------------------------------------+
//| Bar validity. Rejects the EMPTY_VALUE/DBL_MAX case explicitly -    |
//| it passes MathIsValidNumber and would poison every estimator.     |
//+------------------------------------------------------------------+
bool GKValid(const double o, const double h, const double l, const double c)
  {
   if(!MathIsValidNumber(o) || !MathIsValidNumber(h) ||
      !MathIsValidNumber(l) || !MathIsValidNumber(c))
      return false;
   if(o == EMPTY_VALUE || h == EMPTY_VALUE || l == EMPTY_VALUE || c == EMPTY_VALUE)
      return false;
   if(o <= 0.0 || h <= 0.0 || l <= 0.0 || c <= 0.0)
      return false;
   //--- the inversion that produces a plausible positive variance
   if(h < l)
      return false;
   if(h < o || h < c || l > o || l > c)
      return false;
   return true;
  }

//--- Per-bar GK variance. Returns false (not 0.0) on a bad bar, so a
//--- caller can drop the bar instead of averaging a fabricated zero in.
bool GarmanKlassVar(const double o, const double h, const double l,
                    const double c, double &v)
  {
   if(!GKValid(o, h, l, c))
      return false;
   double hl = MathLog(h / l);
   double co = MathLog(c / o);
   v = 0.5 * hl * hl - (2.0 * VOL_LN2 - 1.0) * co * co;
   return true;
  }

//+------------------------------------------------------------------+
//| Rolling GK sigma over `window` bars ending at `shift`.            |
//|                                                                  |
//| Reads via Copy* directly rather than the stdlib series wrappers:  |
//| a 48-72h window on H4 is 12-18 bars so the 1024 ceiling is not in |
//| play here, but keeping one convention across the vol/FFD layer    |
//| removes the chance of a later window quietly crossing it.         |
//+------------------------------------------------------------------+
bool GarmanKlassSigma(const string symbol, const ENUM_TIMEFRAMES period,
                      const int shift, const int window, double &sigma)
  {
   if(window <= 0 || shift < 0)
      return false;
   double o[], h[], l[], c[];
   if(CopyOpen(symbol, period, shift, window, o)  != window ||
      CopyHigh(symbol, period, shift, window, h)  != window ||
      CopyLow(symbol, period, shift, window, l)   != window ||
      CopyClose(symbol, period, shift, window, c) != window)
      return false;
   double acc = 0.0;
   int used = 0;
   for(int i = 0; i < window; i++)
     {
      double v;
      if(!GarmanKlassVar(o[i], h[i], l[i], c[i], v))
         continue;                        // drop the bar, do not zero it
      acc += v;
      used++;
     }
   //--- Refuse a reading built from a thin remnant: a "volatility" averaged
   //--- over 3 of 18 bars is not the same statistic and must not be gated on.
   if(used < (window + 1) / 2)
      return false;
   sigma = MathSqrt(MathMax(acc / used, 0.0));
   return true;
  }

//--- Parkinson: high-low only, no drift term. Useful as a cross-check when
//--- GK's close-open term behaves oddly on a gappy series.
bool ParkinsonSigma(const string symbol, const ENUM_TIMEFRAMES period,
                    const int shift, const int window, double &sigma)
  {
   if(window <= 0 || shift < 0)
      return false;
   double h[], l[];
   if(CopyHigh(symbol, period, shift, window, h) != window ||
      CopyLow(symbol, period, shift, window, l)  != window)
      return false;
   double acc = 0.0;
   int used = 0;
   for(int i = 0; i < window; i++)
     {
      if(!MathIsValidNumber(h[i]) || !MathIsValidNumber(l[i]) ||
         h[i] == EMPTY_VALUE || l[i] == EMPTY_VALUE || h[i] <= 0.0 || l[i] <= 0.0 || h[i] < l[i])
         continue;
      double hl = MathLog(h[i] / l[i]);
      acc += hl * hl / (4.0 * VOL_LN2);
      used++;
     }
   if(used < (window + 1) / 2)
      return false;
   sigma = MathSqrt(acc / used);
   return true;
  }

//+------------------------------------------------------------------+
//| Yang-Zhang: overnight gap variance + open-close variance + the    |
//| Rogers-Satchell drift-robust term. THE ONE TO GATE ON for a       |
//| multi-day CFD hold, because the gap it prices is the risk the     |
//| position actually carries overnight.                              |
//+------------------------------------------------------------------+
bool YangZhangSigma(const string symbol, const ENUM_TIMEFRAMES period,
                    const int shift, const int window, double &sigma,
                    const double kAlpha = 1.34)
  {
   if(window < 2 || shift < 0)
      return false;
   //--- one extra OLDER bar: the first overnight term needs the prior close
   int need = window + 1;
   double o[], h[], l[], c[];
   if(CopyOpen(symbol, period, shift, need, o)  != need ||
      CopyHigh(symbol, period, shift, need, h)  != need ||
      CopyLow(symbol, period, shift, need, l)   != need ||
      CopyClose(symbol, period, shift, need, c) != need)
      return false;
   //--- Copy* returns oldest-first; index 0 is the OLDEST of the need bars.
   double sumO = 0.0, sumC = 0.0, sumRS = 0.0;
   double oc[], co[];
   ArrayResize(oc, window);
   ArrayResize(co, window);
   int used = 0;
   for(int i = 1; i < need; i++)
     {
      if(!GKValid(o[i], h[i], l[i], c[i]) || !GKValid(o[i-1], h[i-1], l[i-1], c[i-1]))
         return false;                    // YZ needs a contiguous run, not a subset
      oc[used] = MathLog(o[i] / c[i - 1]);     // overnight
      co[used] = MathLog(c[i] / o[i]);         // open-to-close
      sumO += oc[used];
      sumC += co[used];
      sumRS += MathLog(h[i] / c[i]) * MathLog(h[i] / o[i]) +
               MathLog(l[i] / c[i]) * MathLog(l[i] / o[i]);
      used++;
     }
   if(used < 2)
      return false;
   double mO = sumO / used, mC = sumC / used;
   double vO = 0.0, vC = 0.0;
   for(int i = 0; i < used; i++)
     {
      vO += (oc[i] - mO) * (oc[i] - mO);
      vC += (co[i] - mC) * (co[i] - mC);
     }
   vO /= (used - 1);
   vC /= (used - 1);
   double vRS = sumRS / used;
   double k = kAlpha / (kAlpha + (used + 1.0) / (used - 1.0));
   sigma = MathSqrt(MathMax(vO + k * vC + (1.0 - k) * vRS, 0.0));
   return true;
  }
#endif // VOLESTIMATORS_MQH