//+------------------------------------------------------------------+
//|                                                    EdgeStats.mqh |
//|                                                     AnimateDread |
//|                                                                  |
//| HOW SURE ARE WE, in Lopez de Prado's terms (Advances in Financial |
//| Machine Learning, ch. 8 and 14). Stateless math, no trading.      |
//|                                                                   |
//|   Phi / PhiInv     the normal CDF and its inverse                 |
//|   PSR(SR*)         Probabilistic Sharpe Ratio: the probability    |
//|                    that the true Sharpe exceeds SR*, given the    |
//|                    sample's length, skew and kurtosis. A Sharpe   |
//|                    of 1 on 40 trades with fat tails is far less   |
//|                    certain than the number suggests.              |
//|   DSR(N)           Deflated Sharpe Ratio: PSR against the Sharpe   |
//|                    the BEST of N independent trials would show by |
//|                    luck. N is how many configurations were tried   |
//|                    to arrive at this one - honestly, including the |
//|                    ones that were discarded.                       |
//|                                                                    |
//| The Sharpe here is PER TRADE (mean R / sd R), not annualised: PSR  |
//| is defined on the sampling frequency of the returns it is given.   |
//| The trial-Sharpe variance uses its large-sample form 1/(T-1) - the |
//| caller has one run, not N, to estimate it from. That makes DSR an  |
//| approximation; it is reported as one.                              |
//+------------------------------------------------------------------+
#ifndef WARRIOR_EDGESTATS_MQH
#define WARRIOR_EDGESTATS_MQH

class CEdgeStats
  {
public:
   //--- Abramowitz & Stegun 7.1.26, |error| < 1.5e-7.
   static double     Phi(const double x)
     {
      const double t = 1.0 / (1.0 + 0.2316419 * MathAbs(x));
      const double d = 0.3989422804 * MathExp(-0.5 * x * x);
      const double p = d * t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 +
                       t * (-1.821255978 + t * 1.330274429))));
      return (x >= 0.0) ? 1.0 - p : p;
     }
   //--- Acklam's rational approximation, relative error < 1.2e-9.
   static double     PhiInv(const double p)
     {
      if(p <= 0.0 || p >= 1.0)
         return 0.0;
      static const double a[6] = {-3.969683028665376e+01, 2.209460984245205e+02,
                                  -2.759285104469687e+02, 1.383577518672690e+02,
                                  -3.066479806614716e+01, 2.506628277459239e+00};
      static const double b[5] = {-5.447609879822406e+01, 1.615858368580409e+02,
                                  -1.556989798598866e+02, 6.680131188771972e+01,
                                  -1.328068155288572e+01};
      static const double c[6] = {-7.784894002430293e-03, -3.223964580411365e-01,
                                  -2.400758277161838e+00, -2.549732539343734e+00,
                                  4.374664141464968e+00, 2.938163982698783e+00};
      static const double d[4] = {7.784695709041462e-03, 3.224671290700398e-01,
                                  2.445134137142996e+00, 3.754408661907416e+00};
      const double lo = 0.02425;
      if(p < lo)
        {
         const double q = MathSqrt(-2.0 * MathLog(p));
         return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
                ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);
        }
      if(p > 1.0 - lo)
        {
         const double q = MathSqrt(-2.0 * MathLog(1.0 - p));
         return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
                 ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);
        }
      const double q = p - 0.5, r = q * q;
      return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q /
             (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0);
     }
   //--- Moments of a sample. Kurtosis is the plain (non-excess) fourth standardised moment.
   static bool       Moments(const double &x[], double &mean, double &sd, double &skew, double &kurt)
     {
      const int n = ArraySize(x);
      if(n < 4)
         return false;
      double s = 0.0;
      for(int i = 0; i < n; i++)
         s += x[i];
      mean = s / n;
      double m2 = 0.0, m3 = 0.0, m4 = 0.0;
      for(int i = 0; i < n; i++)
        {
         const double d = x[i] - mean;
         m2 += d * d;
         m3 += d * d * d;
         m4 += d * d * d * d;
        }
      m2 /= n; m3 /= n; m4 /= n;
      if(m2 <= 0.0)
         return false;
      sd = MathSqrt(m2 * n / (n - 1));
      skew = m3 / MathPow(m2, 1.5);
      kurt = m4 / (m2 * m2);
      return true;
     }
   //--- PSR(SR*) for per-trade returns. -1 when the sample is too small to say anything.
   static double     PSR(const double &x[], const double srStar = 0.0)
     {
      double mean, sd, skew, kurt;
      if(!Moments(x, mean, sd, skew, kurt) || sd <= 0.0)
         return -1.0;
      const int n = ArraySize(x);
      const double sr = mean / sd;
      const double den = 1.0 - skew * sr + (kurt - 1.0) / 4.0 * sr * sr;
      if(den <= 0.0)
         return -1.0;
      return Phi((sr - srStar) * MathSqrt(n - 1.0) / MathSqrt(den));
     }
   //--- The Sharpe the best of `trials` unskilled strategies would show (AFML 14.1, eq. 14.1).
   static double     ExpectedMaxSharpe(const int trials, const int n)
     {
      if(trials < 2 || n < 3)
         return 0.0;
      const double g = 0.5772156649;
      const double sigma = 1.0 / MathSqrt(n - 1.0);
      return sigma * ((1.0 - g) * PhiInv(1.0 - 1.0 / trials) +
                      g * PhiInv(1.0 - 1.0 / (trials * M_E)));
     }
   static double     DSR(const double &x[], const int trials)
     {
      return PSR(x, ExpectedMaxSharpe(trials, ArraySize(x)));
     }
  };
#endif // WARRIOR_EDGESTATS_MQH