//+------------------------------------------------------------------+
//|                                               EP_WalkForward.mq5 |
//|                                      MMQ — Muhammad Minhas Qamar |
//|                                   www.mql5.com/en/articles/24757 |
//+------------------------------------------------------------------+
#property copyright "MMQ — Muhammad Minhas Qamar"
#property link      "https://www.mql5.com/en/articles/24757"
#property version   "1.00"
#property strict
#property script_show_inputs
#property description "Walks forward through a basket's history and forecasts the book's next-bar"
#property description "VaR and ES from each probability vector, using only the bars before it."
#property description "Scores coverage, the joint VaR-ES loss, and the variance forecast."

#include <EntropyPooling\Scenarios.mqh>
#include <EntropyPooling\Probabilities.mqh>
#include <EntropyPooling\Views.mqh>
#include <EntropyPooling\Posterior.mqh>
#include <Math\Stat\ChiSquare.mqh>
#include <Math\Stat\Normal.mqh>

input group "Data"
input string          InpSymbols     = "EURUSD,GBPUSD,AUDUSD,USDJPY,USDCHF"; // basket
input string          InpWeights     = "0.2,0.2,0.2,-0.2,-0.2"; // book notional per symbol, same order
input int             InpBars        = 5000;               // D1 bars requested per symbol
input group "Test"
input int             InpWarmup      = 1000;               // rows before the first forecast
input double          InpLevel       = 0.025;              // tail level for VaR and ES
input group "Methods"
input int             InpWindow      = 250;                // long rolling window
input int             InpShortWindow = 60;                 // short rolling window
input double          InpHalfLife    = 250.0;              // decay half-life, bars
input int             InpStateWindow = 20;                 // trailing volatility of the book, bars
input double          InpKernelBand  = 0.5;                // kernel bandwidth, share of the state's std
input double          InpMinEns      = 250;                // scenario budget for conditioning

//--- the forecasters, in table order
enum ENUM_METHOD
  {
   M_ROLLING,                          // rolling window, long
   M_ROLLING_SHORT,                    // rolling window, short
   M_DECAY,                            // exponential decay
   M_EWMA,                             // RiskMetrics variance with a normal tail
   M_KERNEL,                           // decay times a Gaussian kernel on the state
   M_STATE_FREE,                       // entropy pooling on the state, no budget
   M_STATE,                            // entropy pooling on the state, ENS budget
   M_COUNT
  };

//--- running scores of one forecaster
struct SScore
  {
   int               hits;             // bars at or below the VaR
   int               degenerate;       // bars whose tail held no loss to average
   double            ens_sum;
   double            ens_min;
   double            fz[];             // joint VaR-ES loss per bar
   double            ql[];             // QLIKE loss of the variance per bar
  };

//+------------------------------------------------------------------+
//| Split a comma list, trimming spaces.                             |
//+------------------------------------------------------------------+
int SplitList(const string text,string &out[])
  {
   const int n=StringSplit(text,',',out);
   for(int i=0;i<n;i++)
     {
      StringTrimLeft(out[i]);
      StringTrimRight(out[i]);
     }
   return n;
  }

//+------------------------------------------------------------------+
//| FZ0 loss of Patton, Ziegel and Chen, joint for VaR and ES.       |
//+------------------------------------------------------------------+
double Fz0(const double y,const double var,const double es,const double level)
  {
   const double hit=(y<=var ? 1.0 : 0.0);
   return -hit*(var-y)/(level*es)+var/es+MathLog(-es)-1.0;
  }

//+------------------------------------------------------------------+
//| Kupiec's p-value for the VaR hit rate against the level.         |
//+------------------------------------------------------------------+
double KupiecP(const int hits,const int n,const double level)
  {
   const double pi=(double)hits/n;
   double lr=-2.0*((n-hits)*MathLog(1.0-level)+hits*MathLog(level));
   if(hits>0 && hits<n)
      lr+=2.0*((n-hits)*MathLog(1.0-pi)+hits*MathLog(pi));
   int err=0;
   return 1.0-MathCumulativeDistributionChiSquare(lr,1.0,err);
  }

//+------------------------------------------------------------------+
//| Diebold-Mariano of a - b, Newey-West lag 5; negative favours a.  |
//+------------------------------------------------------------------+
double DieboldMariano(const double &a[],const double &b[])
  {
   const int n=ArraySize(a);
   double mu=0.0;
   for(int i=0;i<n;i++)
      mu+=a[i]-b[i];
   mu/=n;
   double lr=0.0;
   for(int i=0;i<n;i++)
      lr+=MathPow(a[i]-b[i]-mu,2);
   lr/=n;
   for(int lag=1;lag<=5;lag++)
     {
      double g=0.0;
      for(int i=lag;i<n;i++)
         g+=(a[i]-b[i]-mu)*(a[i-lag]-b[i-lag]-mu);
      lr+=2.0*(1.0-lag/6.0)*g/n;
     }
   return (lr>0.0 ? mu/MathSqrt(lr/n) : 0.0);
  }

//+------------------------------------------------------------------+
//| Probability vectors over the rows before r only.                 |
//+------------------------------------------------------------------+
void RollingBefore(const int T,const int r,const int window,vector &p)
  {
   p.Init(T);
   p.Fill(0.0);
   const int a=MathMax(0,r-window);
   for(int i=a;i<r;i++)
      p[i]=1.0/(r-a);
  }

//+------------------------------------------------------------------+
void DecayBefore(const int T,const int r,const double half_life,vector &p)
  {
   p.Init(T);
   p.Fill(0.0);
   const double k=MathLog(2.0)/half_life;
   for(int i=0;i<r;i++)
      p[i]=MathExp(-k*(r-1-i));
   EpNormalise(p);
  }

//+------------------------------------------------------------------+
void OnStart(void)
  {
   string symbols[],wtext[];
   const int n=SplitList(InpSymbols,symbols);
   if(n<2 || SplitList(InpWeights,wtext)!=n)
     {
      Print("EP_WalkForward: give at least two symbols and one weight per symbol");
      return;
     }
   int    cols[];
   double weights[];
   ArrayResize(cols,n);
   ArrayResize(weights,n);
   for(int i=0;i<n;i++)
     {
      cols[i]   =i;
      weights[i]=StringToDouble(wtext[i]);
     }

   CEpScenarios scen;
   bool loaded=false;
   for(int attempt=0;attempt<40 && !loaded;attempt++)
     {
      loaded=scen.LoadReturns(symbols,PERIOD_D1,InpBars);
      if(!loaded)
         Sleep(500);
     }
   if(!loaded)
      return;

   vector book;
   EpPortfolio(scen,cols,weights,book);
   double values[];
   ArrayResize(values,scen.Rows());
   for(int t=0;t<scen.Rows();t++)
      values[t]=book[t];
   const int book_col =scen.AddColumn("Book",EP_COLUMN_RETURN,values,EMPTY_VALUE);
   const int state_col=scen.AddTrailingVolatility(book_col,InpStateWindow);
   const int T=scen.Rows();
   if(state_col<0 || T<=InpWarmup+100)
     {
      Print("EP_WalkForward: too little joint history for this warm-up");
      return;
     }
   scen.Column(book_col,book);
   vector z;
   scen.Column(state_col,z);
   int order[];
   EpSortOrder(book,order);

   const int N=T-InpWarmup;
   SScore score[M_COUNT];
   for(int m=0;m<M_COUNT;m++)
     {
      score[m].hits      =0;
      score[m].degenerate=0;
      score[m].ens_sum   =0.0;
      score[m].ens_min   =DBL_MAX;
      ArrayResize(score[m].fz,N);
      ArrayResize(score[m].ql,N);
     }

//--- EWMA variance seeded over the warm-up, as a desk would run it
   double ewma=book[0]*book[0];
   for(int i=1;i<InpWarmup;i++)
      ewma=0.94*ewma+0.06*book[i-1]*book[i-1];
   int err=0;
   const double zq =MathQuantileNormal(InpLevel,0.0,1.0,err);
   const double phi=MathExp(-0.5*zq*zq)/MathSqrt(2.0*M_PI);

   PrintFormat("EP_WalkForward: %d rows %s to %s, forecasting the last %d at the %.1f%% tail",
               T,TimeToString(scen.Time(0),TIME_DATE),TimeToString(scen.Time(T-1),TIME_DATE),N,100.0*InpLevel);
   const uint t0=GetTickCount();
   vector p[M_COUNT];
   int bound=0;                        // forecasts where the budget pulled the target back
   for(int s=0;s<N;s++)
     {
      const int    r=InpWarmup+s;
      const double y=book[r];
      RollingBefore(T,r,InpWindow,p[M_ROLLING]);
      RollingBefore(T,r,InpShortWindow,p[M_ROLLING_SHORT]);
      DecayBefore(T,r,InpHalfLife,p[M_DECAY]);
      ewma=0.94*ewma+0.06*book[r-1]*book[r-1];

      //--- the state for row r was known before its bar opened
      const double sd=EpVolatility(z,p[M_DECAY]);
      EpKernel(z,z[r],InpKernelBand*sd,p[M_DECAY],p[M_KERNEL]);
      SEpSolve res;
      EpConditionState(scen,p[M_DECAY],state_col,z[r],0.0,p[M_STATE_FREE],res);
      if(EpConditionState(scen,p[M_DECAY],state_col,z[r],InpMinEns,p[M_STATE],res)!=z[r])
         bound++;

      for(int m=0;m<M_COUNT;m++)
        {
         double var,es,vol;
         if(m==M_EWMA)
           {
            vol=MathSqrt(ewma);
            var=zq*vol;
            es =-vol*phi/InpLevel;
           }
         else
           {
            EpTailOrdered(book,order,p[m],InpLevel,var,es);
            vol=EpVolatility(book,p[m]);
            const double ens=EpEffectiveScenarios(p[m]);
            score[m].ens_sum+=ens;
            score[m].ens_min =MathMin(score[m].ens_min,ens);
           }
         if(y<=var)
            score[m].hits++;
         if(es>=0.0 || vol<=0.0)
           {
            score[m].degenerate++;
            score[m].fz[s]=0.0;
            score[m].ql[s]=0.0;
            continue;
           }
         score[m].fz[s]=Fz0(y,var,es,InpLevel);
         score[m].ql[s]=MathLog(vol*vol)+y*y/(vol*vol);
        }
     }

   const string names[M_COUNT]={StringFormat("rolling %d",InpWindow),StringFormat("rolling %d",InpShortWindow),
                                StringFormat("decay %.0f",InpHalfLife),"EWMA 0.94 normal","decay x kernel",
                                "state, no budget",StringFormat("state, ENS >= %.0f",InpMinEns)};
   PrintFormat("EP_WalkForward: %d forecasts in %u ms",N,GetTickCount()-t0);
   PrintFormat("%-18s %6s %6s %7s %8s %9s %8s %8s %8s %9s","method","ENS","minENS","hits",
               "Kupiec p","FZ0","DM roll","DM decay","DM kern","QLIKE");
   for(int m=0;m<M_COUNT;m++)
     {
      double fz=0.0,ql=0.0;
      for(int s=0;s<N;s++)
        {
         fz+=score[m].fz[s];
         ql+=score[m].ql[s];
        }
      const bool   usable=(score[m].degenerate==0);
      const string ens   =(m==M_EWMA ? "-" : StringFormat("%.0f",score[m].ens_sum/N));
      const string ens_lo=(m==M_EWMA ? "-" : StringFormat("%.0f",score[m].ens_min));
      PrintFormat("%-18s %6s %6s %6.2f%% %8.3f %9s %8s %8s %8s %9s",names[m],ens,ens_lo,100.0*score[m].hits/N,
                  KupiecP(score[m].hits,N,InpLevel),
                  (usable ? StringFormat("%.4f",fz/N) : "n/a"),
                  (usable ? StringFormat("%.2f",DieboldMariano(score[m].fz,score[M_ROLLING].fz)) : "n/a"),
                  (usable ? StringFormat("%.2f",DieboldMariano(score[m].fz,score[M_DECAY].fz)) : "n/a"),
                  (usable ? StringFormat("%.2f",DieboldMariano(score[m].fz,score[M_KERNEL].fz)) : "n/a"),
                  (usable ? StringFormat("%.4f",ql/N) : "n/a"));
      if(!usable)
         PrintFormat("%-18s %d of %d forecasts had a non-negative ES, where FZ0 is undefined",
                     "",score[m].degenerate,N);
     }
   PrintFormat("EP_WalkForward: the budget pulled the state target back on %d of %d forecasts",bound,N);
  }
//+------------------------------------------------------------------+