//+------------------------------------------------------------------+
//|                                                DeflatedSharpe.mqh  |
//|                                   @temonba - binaryforexea.com    |
//|  Tells you whether a strategy's Sharpe ratio is real or luck.     |
//|  A raw Sharpe from a backtest is optimistic for three reasons: a  |
//|  short track record, fat tails, and the fact that you kept the    |
//|  best of many variants you tried. This class answers all three.   |
//|  The Probabilistic Sharpe Ratio (PSR) is the probability the true |
//|  Sharpe is above a benchmark, given the sample size and the       |
//|  higher moments. The Deflated Sharpe Ratio (DSR) is the same test |
//|  against a benchmark raised for the number of trials, so the best |
//|  of N random strategies is no longer mistaken for a real edge.    |
//|  Everything is native: the moments, the normal CDF, and its       |
//|  inverse are written from scratch. No Python, no DLL, no library. |
//+------------------------------------------------------------------+
#property copyright "@temonba"
#property link      "https://www.mql5.com/en/users/temonba"
#property strict

//+------------------------------------------------------------------+
//| Probabilistic and Deflated Sharpe Ratio, from scratch            |
//+------------------------------------------------------------------+
class CDeflatedSharpe
  {
public:
   //--- sample moments of a return series
   static double     Mean(const double &x[]);           // arithmetic mean
   static double     Std(const double &x[]);            // sample standard deviation (divides by n - 1)
   static double     Skew(const double &x[]);           // skewness, third standardized moment
   static double     Kurtosis(const double &x[]);       // Pearson kurtosis, normal = 3

   //--- native normal distribution
   static double     NormalCDF(const double z);         // standard normal cumulative probability
   static double     NormalInv(const double p);         // inverse of the standard normal CDF

   //--- Sharpe and its significance
   static double     Sharpe(const double &r[]);         // observed Sharpe, per observation, not annualized
   static double     PSR(const double &r[], const double srStar = 0.0);                 // probabilistic Sharpe vs a benchmark
   static double     ExpectedMaxSharpe(const double varSR, const int nTrials);          // deflation hurdle for nTrials
   static double     DSR(const double &rBest[], const double varSR, const int nTrials); // deflated Sharpe ratio
  };
//+------------------------------------------------------------------+
//| Arithmetic mean                                                  |
//+------------------------------------------------------------------+
double CDeflatedSharpe::Mean(const double &x[])
  {
   int n = ArraySize(x);
   if(n == 0)
      return(0.0);
   double s = 0.0;
   for(int i = 0; i < n; i++)
      s += x[i];
   return(s / n);
  }
//+------------------------------------------------------------------+
//| Sample standard deviation (divides by n - 1)                     |
//+------------------------------------------------------------------+
double CDeflatedSharpe::Std(const double &x[])
  {
   int n = ArraySize(x);
   if(n < 2)
      return(0.0);
   double m = Mean(x);
   double s = 0.0;
   for(int i = 0; i < n; i++)
      s += (x[i] - m) * (x[i] - m);
   return(MathSqrt(s / (n - 1)));
  }
//+------------------------------------------------------------------+
//| Skewness (third standardized moment)                             |
//+------------------------------------------------------------------+
double CDeflatedSharpe::Skew(const double &x[])
  {
   int n = ArraySize(x);
   if(n < 3)
      return(0.0);
   double m = Mean(x);
   double sd = Std(x);
   if(sd == 0.0)
      return(0.0);
   double s = 0.0;
   for(int i = 0; i < n; i++)
     {
      double z = (x[i] - m) / sd;
      s += z * z * z;
     }
   return(s / n);
  }
//+------------------------------------------------------------------+
//| Kurtosis (fourth standardized moment, normal distribution = 3)   |
//+------------------------------------------------------------------+
double CDeflatedSharpe::Kurtosis(const double &x[])
  {
   int n = ArraySize(x);
   if(n < 4)
      return(3.0);
   double m = Mean(x);
   double sd = Std(x);
   if(sd == 0.0)
      return(3.0);
   double s = 0.0;
   for(int i = 0; i < n; i++)
     {
      double z = (x[i] - m) / sd;
      s += z * z * z * z;
     }
   return(s / n);
  }
//+------------------------------------------------------------------+
//| Standard normal CDF via a rational approximation of erf          |
//+------------------------------------------------------------------+
double CDeflatedSharpe::NormalCDF(const double z)
  {
   double x = z / MathSqrt(2.0);
   int    sign = (x < 0.0) ? -1 : 1;
   x = MathAbs(x);
   double t = 1.0 / (1.0 + 0.3275911 * x);
   double y = 1.0 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t
                      - 0.284496736) * t + 0.254829592) * t * MathExp(-x * x);
   double erf = sign * y;
   return(0.5 * (1.0 + erf));
  }
//+------------------------------------------------------------------+
//| Inverse standard normal CDF (Acklam's rational approximation)    |
//+------------------------------------------------------------------+
double CDeflatedSharpe::NormalInv(const double p)
  {
   if(p <= 0.0)
      return(-1.0e10);
   if(p >= 1.0)
      return(1.0e10);
   double a1 = -3.969683028665376e+01, a2 =  2.209460984245205e+02;
   double a3 = -2.759285104469687e+02, a4 =  1.383577518672690e+02;
   double a5 = -3.066479806614716e+01, a6 =  2.506628277459239e+00;
   double b1 = -5.447609879822406e+01, b2 =  1.615858368580409e+02;
   double b3 = -1.556989798598866e+02, b4 =  6.680131188771972e+01;
   double b5 = -1.328068155288572e+01;
   double c1 = -7.784894002430293e-03, c2 = -3.223964580411365e-01;
   double c3 = -2.400758277161838e+00, c4 = -2.549732539343734e+00;
   double c5 =  4.374664141464968e+00, c6 =  2.938163982698783e+00;
   double d1 =  7.784695709041462e-03, d2 =  3.224671290700398e-01;
   double d3 =  2.445134137142996e+00, d4 =  3.754408661907416e+00;
   double plow = 0.02425, phigh = 1.0 - 0.02425;
   double q, r;
   if(p < plow)
     {
      q = MathSqrt(-2.0 * MathLog(p));
      return((((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) /
             ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0));
     }
   if(p > phigh)
     {
      q = MathSqrt(-2.0 * MathLog(1.0 - p));
      return(-(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) /
             ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0));
     }
   q = p - 0.5;
   r = q * q;
   return((((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q /
          (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1.0));
  }
//+------------------------------------------------------------------+
//| Observed Sharpe of a return series (per observation)             |
//+------------------------------------------------------------------+
double CDeflatedSharpe::Sharpe(const double &r[])
  {
   double sd = Std(r);
   if(sd == 0.0)
      return(0.0);
   return(Mean(r) / sd);
  }
//+------------------------------------------------------------------+
//| Probabilistic Sharpe Ratio: P(true Sharpe > srStar)              |
//| Corrects for sample length, skewness and kurtosis.               |
//+------------------------------------------------------------------+
double CDeflatedSharpe::PSR(const double &r[], const double srStar = 0.0)
  {
   int n = ArraySize(r);
   if(n < 4)
      return(0.0);
   double sr = Sharpe(r);
   double g3 = Skew(r);
   double g4 = Kurtosis(r);
   double denom = 1.0 - g3 * sr + ((g4 - 1.0) / 4.0) * sr * sr;
   if(denom <= 0.0)
      return(0.0);
   double z = (sr - srStar) * MathSqrt((double)(n - 1)) / MathSqrt(denom);
   return(NormalCDF(z));
  }
//+------------------------------------------------------------------+
//| Expected maximum Sharpe of nTrials independent random strategies |
//| given the variance of their Sharpe estimates. This is the bar a  |
//| real edge must clear once you admit how many variants you tried. |
//+------------------------------------------------------------------+
double CDeflatedSharpe::ExpectedMaxSharpe(const double varSR, const int nTrials)
  {
   if(nTrials <= 1 || varSR <= 0.0)
      return(0.0);
   double emc = 0.5772156649015329;      // Euler-Mascheroni constant
   double e   = 2.718281828459045;
   double q1  = NormalInv(1.0 - 1.0 / nTrials);
   double q2  = NormalInv(1.0 - 1.0 / (nTrials * e));
   return(MathSqrt(varSR) * ((1.0 - emc) * q1 + emc * q2));
  }
//+------------------------------------------------------------------+
//| Deflated Sharpe Ratio: PSR of the best series against the        |
//| expected-maximum benchmark for nTrials. varSR is the variance of |
//| the Sharpe ratios across all trials that were tried.             |
//+------------------------------------------------------------------+
double CDeflatedSharpe::DSR(const double &rBest[], const double varSR, const int nTrials)
  {
   double srStar = ExpectedMaxSharpe(varSR, nTrials);
   return(PSR(rBest, srStar));
  }
//+------------------------------------------------------------------+