//+------------------------------------------------------------------+
//|                                                      RunPBO.mq5   |
//|              Probability of Backtest Overfitting - Astralys LLC   |
//|                                                                   |
//| Builds the returns matrix of a full parameter sweep, then runs    |
//| the CSCV engine on it.                                            |
//|                                                                   |
//| Strategy under test, following the full market exposure design:   |
//| one threshold on MACD-on-price switches the position between long |
//| and short. The account is always in the market, so strategy       |
//| volatility equals market volatility and only the sign changes.    |
//| That keeps the returns matrix trivial to build and the comparison |
//| against buy and hold honest.                                      |
//|                                                                   |
//| The signal is read on the previous closed bar and applied to the  |
//| current bar's return. No look ahead.                              |
//|                                                                   |
//| The threshold range is not hard coded. It is derived from the     |
//| observed distribution of the indicator, with the rare extremes    |
//| trimmed, so the grid covers levels the market actually reaches.   |
//+------------------------------------------------------------------+
#property copyright "Astralys LLC"
#property link      "https://pulsar-terminal.com"
#property version   "1.00"
#property script_show_inputs

#include <PBO/MACDp.mqh>
#include <PBO/BearsPower.mqh>
#include <PBO/CSCVEngine.mqh>

// Prefixed because MQL5 already defines IND_BEARS and friends in its
// own built-in ENUM_INDICATOR.
enum ENUM_PBO_INDICATOR
  {
   PBO_IND_MACDP,    // MACD on price
   PBO_IND_BEARS     // Bears Power
  };

input ENUM_PBO_INDICATOR InpIndicator = PBO_IND_MACDP;  // Indicator under test
input string InpSymbol      = "";     // Symbol (empty = chart symbol)
input int    InpBars        = 100000; // Max bars to load
input int    InpPeriodMin   = 2;      // Indicator period, from
input int    InpPeriodMax   = 14;     // Indicator period, to
input int    InpLevels      = 100;    // Number of thresholds tested
input bool   InpNormalise   = true;   // Percentage deviation instead of raw points
input double InpTrimPct     = 1.0;    // Percent trimmed at each tail of the range
input int    InpPartitions  = 16;     // CSCV partitions (S)
input bool   InpExportCsv   = true;   // Write logits and degradation to MQL5\Files

//+------------------------------------------------------------------+
//| Pearson correlation of two equally sized arrays.                  |
//+------------------------------------------------------------------+
double Correlation(const double &a[], const double &b[])
  {
   const int n = ArraySize(a);
   if(n < 2 || ArraySize(b) != n)
      return(0.0);

   double ma = 0.0, mb = 0.0;
   for(int i = 0; i < n; i++) { ma += a[i]; mb += b[i]; }
   ma /= n; mb /= n;

   double num = 0.0, da = 0.0, db = 0.0;
   for(int i = 0; i < n; i++)
     {
      const double u = a[i] - ma, v = b[i] - mb;
      num += u * v; da += u * u; db += v * v;
     }
   if(da <= 0.0 || db <= 0.0)
      return(0.0);

   return(num / MathSqrt(da * db));
  }

//+------------------------------------------------------------------+
//| Value at a given percentile of an already sorted array.           |
//+------------------------------------------------------------------+
double PercentileSorted(const double &sorted[], const double pct)
  {
   const int n = ArraySize(sorted);
   if(n == 0) return(0.0);

   int idx = (int)MathRound(pct / 100.0 * (n - 1));
   if(idx < 0)      idx = 0;
   if(idx > n - 1)  idx = n - 1;
   return(sorted[idx]);
  }

//+------------------------------------------------------------------+
//| Script entry point                                                |
//+------------------------------------------------------------------+
void OnStart(void)
  {
   const string symbol = (InpSymbol == "" ? _Symbol : InpSymbol);

   //--- 1. price history ------------------------------------------------
   MqlRates rates[];
   ArraySetAsSeries(rates, false);              // index 0 = oldest

   const int bars = CopyRates(symbol, PERIOD_D1, 0, InpBars, rates);
   if(bars <= 0)
     {
      PrintFormat("Cannot load history for %s: error %d", symbol, GetLastError());
      return;
     }

   double close[], low[];
   ArrayResize(close, bars);
   ArrayResize(low,   bars);
   for(int i = 0; i < bars; i++)
     {
      close[i] = rates[i].close;
      low[i]   = rates[i].low;
     }

   const string indName = (InpIndicator == PBO_IND_MACDP) ? "MACD on price" : "Bears Power";

   Print("============================================================");
   PrintFormat("PBO run on %s D1   |   %s, %s",
               symbol, indName, InpNormalise ? "normalised" : "raw points");
   PrintFormat("%d bars, %s to %s",
               bars,
               TimeToString(rates[0].time, TIME_DATE),
               TimeToString(rates[bars-1].time, TIME_DATE));

   //--- 2. one indicator series per period ------------------------------
   const int nPeriods = InpPeriodMax - InpPeriodMin + 1;
   const int nLevels  = InpLevels;
   const int N        = nPeriods * nLevels;

   if(nPeriods < 1 || nLevels < 2)
     {
      Print("Empty parameter grid.");
      return;
     }

   // The indicator depends only on the period, so it is computed once
   // per period and reused across every threshold. That is the
   // difference between 13 indicator passes and 1300 of them.
   double ind[];                                 // nPeriods x bars
   ArrayResize(ind, nPeriods * bars);

   for(int p = 0; p < nPeriods; p++)
     {
      double one[];
      const int n  = InpPeriodMin + p;
      const int ok = (InpIndicator == PBO_IND_MACDP)
                     ? MACDPrice(close, n, InpNormalise, one)
                     : BearsPower(low, close, n, InpNormalise, one);

      if(ok < 0)
        {
         PrintFormat("%s failed for period %d", indName, n);
         return;
        }
      for(int i = 0; i < bars; i++)
         ind[p * bars + i] = one[i];
     }

   //--- 3. threshold range, read from the data --------------------------
   // Taken from the longest period, following the thesis method: look at
   // the distribution of the indicator and keep the extremes that recur,
   // not the rare ones.
   double sample[];
   ArrayResize(sample, 0);
   for(int i = InpPeriodMax - 1; i < bars; i++)
     {
      const double v = ind[(nPeriods - 1) * bars + i];
      if(v == EMPTY_VALUE) continue;
      const int k = ArraySize(sample);
      ArrayResize(sample, k + 1);
      sample[k] = v;
     }
   ArraySort(sample);

   const double lo = PercentileSorted(sample, InpTrimPct);
   const double hi = PercentileSorted(sample, 100.0 - InpTrimPct);

   if(hi <= lo)
     {
      Print("Degenerate threshold range.");
      return;
     }

   const double step = (hi - lo) / (nLevels - 1);
   PrintFormat("Threshold range from data: %.4f to %.4f in %d steps of %.4f%s",
               lo, hi, nLevels, step, InpNormalise ? " (percent)" : "");

   //--- 4. returns matrix ------------------------------------------------
   // The first usable bar is the one where every period already has a
   // valid reading on the PREVIOUS bar, so the whole grid shares the
   // same rows and the matrix stays rectangular.
   const int first = InpPeriodMax;
   const int T     = bars - first;

   if(T < InpPartitions * 2)
     {
      PrintFormat("Only %d usable rows, not enough for %d partitions.", T, InpPartitions);
      return;
     }

   double market[];                              // log return of the index
   ArrayResize(market, T);
   for(int t = 0; t < T; t++)
     {
      const int b = first + t;
      market[t] = MathLog(close[b] / close[b-1]);
     }

   PrintFormat("Grid: periods %d..%d x %d levels = %d combinations",
               InpPeriodMin, InpPeriodMax, nLevels, N);
   PrintFormat("Matrix: %d rows x %d columns (%.1f MB)",
               T, N, (double)T * N * 8.0 / 1048576.0);

   // Note: "matrix" is a reserved type name in MQL5, hence retMatrix.
   double retMatrix[];
   if(ArrayResize(retMatrix, T * N) != T * N)
     {
      Print("Cannot allocate the returns matrix.");
      return;
     }

   const uint t0 = GetTickCount();
   int signals[];                                // position changes per column
   ArrayResize(signals, N);
   ArrayInitialize(signals, 0);

   for(int p = 0; p < nPeriods; p++)
     {
      for(int l = 0; l < nLevels; l++)
        {
         const int    n     = p * nLevels + l;   // column index
         const double level = lo + l * step;

         int prevSign = 0;
         for(int t = 0; t < T; t++)
           {
            const int    b   = first + t;
            const double val = ind[p * bars + (b - 1)];   // signal on the closed bar
            const int    s   = (val > level) ? 1 : -1;

            if(prevSign != 0 && s != prevSign)
               signals[n]++;
            prevSign = s;

            retMatrix[t * N + n] = s * market[t];
           }
        }
     }

   PrintFormat("Matrix built in %.1f s", (GetTickCount() - t0) / 1000.0);

   //--- 5. buy and hold reference ---------------------------------------
   double bh = 0.0;
   for(int t = 0; t < T; t++) bh += market[t];
   PrintFormat("Buy and hold over the period: %.1f%% cumulative",
               (MathExp(bh) - 1.0) * 100.0);

   //--- 6. CSCV ----------------------------------------------------------
   CCSCVEngine engine;
   if(!engine.SetPartitions(InpPartitions)) return;
   if(!engine.SetReturns(retMatrix, T, N))  return;
   if(!engine.Run())                        return;

   //--- 7. what the sweep would have selected in sample -----------------
   double is[], oos[];
   engine.GetDegradation(is, oos);

   int selected[];
   engine.GetSelected(selected);

   double bestIS = -DBL_MAX, pairedOOS = 0.0;
   int    bestCol = -1;
   for(int c = 0; c < ArraySize(is); c++)
      if(is[c] > bestIS) { bestIS = is[c]; pairedOOS = oos[c]; bestCol = selected[c]; }

   Print("------------------------------------------------------------");
   PrintFormat("Best in-sample half seen across all splits: %.1f%% cumulative",
               (MathExp(bestIS) - 1.0) * 100.0);
   PrintFormat("Same parameters out of sample on that split: %.1f%%",
               (MathExp(pairedOOS) - 1.0) * 100.0);

   if(bestCol >= 0)
      PrintFormat("That winner is period %d, threshold %+.3f, %d position changes",
                  InpPeriodMin + bestCol / nLevels,
                  lo + (bestCol % nLevels) * step,
                  signals[bestCol]);

   //--- 7b. is the procedure stable, or is it picking at random? ---------
   // How often each column is chosen across the 12,870 splits. A single
   // column dominating means the sweep keeps landing on the same place.
   // A scattered picture means it is chasing noise.
   int picks[];
   ArrayResize(picks, N);
   ArrayInitialize(picks, 0);
   for(int c = 0; c < ArraySize(selected); c++)
      picks[selected[c]]++;

   int distinct = 0;
   for(int n = 0; n < N; n++)
      if(picks[n] > 0) distinct++;

   PrintFormat("Columns ever selected: %d out of %d", distinct, N);

   for(int rank = 0; rank < 5; rank++)
     {
      int top = -1, best = 0;
      for(int n = 0; n < N; n++)
         if(picks[n] > best) { best = picks[n]; top = n; }
      if(top < 0) break;

      PrintFormat("   #%d  period %2d, threshold %+.3f, %5d changes  ->  chosen %.1f%% of splits",
                  rank + 1,
                  InpPeriodMin + top / nLevels,
                  lo + (top % nLevels) * step,
                  signals[top],
                  100.0 * picks[top] / ArraySize(selected));
      picks[top] = 0;
     }

   //--- 7c. how much of the grid is degenerate --------------------------
   // A column with almost no position changes is buy and hold, or its
   // mirror image. In a rising market those score well in sample without
   // being a strategy at all.
   int deg = 0;
   for(int n = 0; n < N; n++)
      if(signals[n] < 10) deg++;
   PrintFormat("Near degenerate columns (fewer than 10 changes): %d of %d (%.1f%%)",
               deg, N, 100.0 * deg / N);

   //--- 8. logit distribution shape --------------------------------------
   double logits[];
   engine.GetLogits(logits);
   const int L = ArraySize(logits);

   double mean = 0.0, sd = 0.0;
   for(int i = 0; i < L; i++) mean += logits[i];
   mean /= L;
   for(int i = 0; i < L; i++) { const double d = logits[i] - mean; sd += d * d; }
   sd = MathSqrt(sd / (L - 1));

   PrintFormat("Logits: mean %+.3f, sd %.3f over %d splits", mean, sd, L);
   PrintFormat("(under no information the null is standard logistic, sd = %.3f)",
               M_PI / MathSqrt(3.0));

   //--- 9. trading activity ----------------------------------------------
   int minSig = INT_MAX, maxSig = 0; double avgSig = 0.0;
   for(int n = 0; n < N; n++)
     {
      if(signals[n] < minSig) minSig = signals[n];
      if(signals[n] > maxSig) maxSig = signals[n];
      avgSig += signals[n];
     }
   PrintFormat("Position changes per combination: min %d, average %.0f, max %d",
               minSig, avgSig / N, maxSig);

   //--- 9b. control: what does a random selection give on the same splits?
   // The identity IS + OOS = total makes any such scatter lean negative
   // on its own. The control measures how much of the slope is just that.
   int poolAll[], poolSel[];
   ArrayResize(poolAll, N);
   for(int n = 0; n < N; n++) poolAll[n] = n;

   ArrayResize(poolSel, 0);
   for(int n = 0; n < N; n++)
     {
      bool used = false;
      for(int c = 0; c < ArraySize(selected) && !used; c++)
         if(selected[c] == n) used = true;
      if(used)
        {
         const int k = ArraySize(poolSel);
         ArrayResize(poolSel, k + 1);
         poolSel[k] = n;
        }
     }

   double ctlAllIS[], ctlAllOOS[], ctlSelIS[], ctlSelOOS[];
   engine.RunControl(poolAll, 20260816, ctlAllIS, ctlAllOOS);
   engine.RunControl(poolSel, 20260816, ctlSelIS, ctlSelOOS);

   Print("------------------------------------------------------------");
   PrintFormat("Correlation of out-of-sample against in-sample, same splits:");
   PrintFormat("   real selection, best in sample : %+.3f", Correlation(is, oos));
   PrintFormat("   control, random over all %4d  : %+.3f", N, Correlation(ctlAllIS, ctlAllOOS));
   PrintFormat("   control, random over the %3d ever selected : %+.3f",
               ArraySize(poolSel), Correlation(ctlSelIS, ctlSelOOS));
   Print("   (a real selection effect shows as a correlation below the control)");

   //--- 10. raw output for plotting --------------------------------------
   // The console gives the summary. The figures need the 12,870 values
   // behind it, so they are written out rather than re-derived by hand.
   if(InpExportCsv)
     {
      const string tag = (InpIndicator == PBO_IND_MACDP) ? "macdp" : "bears";

      const string f1 = StringFormat("PBO_logits_%s.csv", tag);
      int h = FileOpen(f1, FILE_WRITE | FILE_CSV | FILE_ANSI, ',');
      if(h != INVALID_HANDLE)
        {
         FileWrite(h, "logit");
         for(int i = 0; i < L; i++)
            FileWrite(h, DoubleToString(logits[i], 6));
         FileClose(h);
         PrintFormat("Wrote %s (%d rows)", f1, L);
        }
      else
         PrintFormat("Cannot write %s, error %d", f1, GetLastError());

      const string f2 = StringFormat("PBO_degradation_%s.csv", tag);
      h = FileOpen(f2, FILE_WRITE | FILE_CSV | FILE_ANSI, ',');
      if(h != INVALID_HANDLE)
        {
         FileWrite(h, "is_log", "oos_log", "column", "period", "threshold", "changes");
         for(int c = 0; c < ArraySize(is); c++)
           {
            const int col = selected[c];
            FileWrite(h,
                      DoubleToString(is[c],  6),
                      DoubleToString(oos[c], 6),
                      IntegerToString(col),
                      IntegerToString(InpPeriodMin + col / nLevels),
                      DoubleToString(lo + (col % nLevels) * step, 4),
                      IntegerToString(signals[col]));
           }
         FileClose(h);
         PrintFormat("Wrote %s (%d rows)", f2, ArraySize(is));
        }
      else
         PrintFormat("Cannot write %s, error %d", f2, GetLastError());
     }

   Print("============================================================");
  }
//+------------------------------------------------------------------+