//+------------------------------------------------------------------+
//|                                                      EP_Demo.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 "Builds one basket's joint history, compares the probability vectors"
#property description "a risk number can stand on, conditions on today's volatility, adds"
#property description "three views, and reports what each step does to the book's risk."

#include <EntropyPooling\Scenarios.mqh>
#include <EntropyPooling\Probabilities.mqh>
#include <EntropyPooling\MinRelEntropy.mqh>
#include <EntropyPooling\Views.mqh>
#include <EntropyPooling\Posterior.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 ENUM_TIMEFRAMES InpTimeframe   = PERIOD_D1;          // bar size of one scenario
input int             InpBars        = 5000;               // bars requested per symbol
input datetime        InpAsOf        = 0;                  // replay as of this date, 0 for today
input group "Prior"
input int             InpWindow      = 250;                // rolling window it is compared with
input double          InpHalfLife    = 250.0;              // exponential decay half-life, bars
input group "State"
input int             InpStateWindow = 20;                 // trailing volatility of the book, bars
input double          InpMinEns      = 250;                // fewest effective scenarios conditioning may leave
input double          InpKernelBand  = 0.5;                // kernel bandwidth, share of the state's std
input group "Views"
input string          InpVolSymbol   = "USDJPY";           // volatility view on
input double          InpVolScale    = 1.5;                // times its conditioned volatility
input string          InpTailSymbol  = "USDCHF";           // tail view on
input double          InpTailMove    = -0.02;              // daily log return at or below
input double          InpTailProb    = 0.01;               // has at least this probability
input string          InpCorrA       = "EURUSD";           // correlation view between
input string          InpCorrB       = "GBPUSD";           // and
input double          InpCorr        = 0.40;               // at most this correlation
input bool            InpCorrPinVol  = true;               // hold both vols, so only the correlation moves
input double          InpConfidence  = 1.0;                // weight on the views, 0..1
input group "Risk"
input double          InpLevel       = 0.025;              // tail level for VaR and ES

//--- one probability vector and its label
struct SStage
  {
   string            name;
   string            code;             // short label for wide tables
   vector            p;
  };

SStage g_stage[];

//+------------------------------------------------------------------+
void AddStage(const string name,const string code,const vector &p)
  {
   const int n=ArraySize(g_stage);
   ArrayResize(g_stage,n+1);
   g_stage[n].name=name;
   g_stage[n].code=code;
   g_stage[n].p   =p;
  }

//+------------------------------------------------------------------+
//| 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;
  }

//+------------------------------------------------------------------+
//| History may still be downloading when a script starts.           |
//+------------------------------------------------------------------+
bool LoadScenarios(CEpScenarios &scen,const string &symbols[])
  {
   for(int attempt=0;attempt<40;attempt++)
     {
      bool ready=true;
      for(int i=0;i<ArraySize(symbols);i++)
        {
         SymbolSelect(symbols[i],true);
         if(!SeriesInfoInteger(symbols[i],InpTimeframe,SERIES_SYNCHRONIZED))
            ready=false;
        }
      if(ready && scen.LoadReturns(symbols,InpTimeframe,InpBars))
         return true;
      Sleep(500);
     }
   return scen.LoadReturns(symbols,InpTimeframe,InpBars);
  }

//+------------------------------------------------------------------+
//| Per stage: ENS, KL to the decay prior, and the book's risk.      |
//+------------------------------------------------------------------+
void ReportStages(const vector &book,const int decay)
  {
   const double ann=MathSqrt(252.0);
   Print("--- the book under each probability vector ---------------------");
   PrintFormat("  %-24s %7s %8s %9s %9s %9s %8s","vector","ENS","KL","vol ann","VaR","ES","size");
   double es_ref=0.0;
   for(int s=0;s<ArraySize(g_stage);s++)
     {
      SEpMoments m;
      EpDescribe(book,g_stage[s].p,InpLevel,m);
      if(s==1)
         es_ref=m.es;
      const double kl=EpRelativeEntropy(g_stage[s].p,g_stage[decay].p);
      PrintFormat("  %-24s %7.0f %8s %8.2f%% %8.3f%% %8.3f%% %8s",g_stage[s].name,
                  EpEffectiveScenarios(g_stage[s].p),
                  (kl==EMPTY_VALUE ? "inf" : StringFormat("%.4f",kl)),100.0*m.vol*ann,
                  100.0*m.var,100.0*m.es,(s>=1 && m.es<0.0 ? StringFormat("%.2fx",es_ref/m.es) : ""));
     }
   PrintFormat("  size: position multiple that keeps the %.1f%% ES of the rolling window",
               100.0*InpLevel);
  }

//+------------------------------------------------------------------+
void ReportSymbols(const CEpScenarios &scen,const int n)
  {
   const double ann=MathSqrt(252.0);
   Print("--- annualised volatility per symbol ----------------------------");
   string head=StringFormat("  %-24s",""),line;
   for(int k=0;k<n;k++)
      head+=StringFormat(" %8s",scen.Name(k));
   Print(head);
   vector x;
   for(int s=0;s<ArraySize(g_stage);s++)
     {
      line=StringFormat("  %-24s",g_stage[s].name);
      for(int k=0;k<n;k++)
        {
         scen.Column(k,x);
         line+=StringFormat(" %7.2f%%",100.0*ann*EpVolatility(x,g_stage[s].p));
        }
      Print(line);
     }
  }

//+------------------------------------------------------------------+
//| The book's worst days and each vector's weight on them, times T. |
//+------------------------------------------------------------------+
void ReportWorstDays(const CEpScenarios &scen,const vector &book,const int count)
  {
   const int T=scen.Rows();
   double a[][2];
   ArrayResize(a,T);
   for(int t=0;t<T;t++)
     {
      a[t][0]=book[t];
      a[t][1]=t;
     }
   ArraySort(a);
   Print("--- the book's worst days: weight as a multiple of 1/T ----------");
   string head=StringFormat("  %-10s %8s","date","return");
   for(int s=0;s<ArraySize(g_stage);s++)
      head+=StringFormat(" %8s",g_stage[s].code);
   Print(head);
   for(int i=0;i<count && i<T;i++)
     {
      const int t=(int)a[i][1];
      string line=StringFormat("  %-10s %7.2f%%",TimeToString(scen.Time(t),TIME_DATE),100.0*a[i][0]);
      for(int s=0;s<ArraySize(g_stage);s++)
         line+=StringFormat(" %8.3f",g_stage[s].p[t]*T);
      Print(line);
     }
  }

//+------------------------------------------------------------------+
void ReportSolve(const string what,const SEpSolve &r,const CEpViews &v)
  {
   PrintFormat("  %-24s %s: %d pass(es), %d Newton steps, worst miss %.1e, KL %.4f, ENS %.0f, %u ms",
               what,EnumToString(r.status),v.Passes(),r.iterations,r.violation,r.entropy,r.ens,r.ms);
  }

//+------------------------------------------------------------------+
void OnStart(void)
  {
   string symbols[],wtext[];
   const int n=SplitList(InpSymbols,symbols);
   if(n<2 || SplitList(InpWeights,wtext)!=n)
     {
      Print("EP_Demo: 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]);
     }

   Print("================================================================");
   PrintFormat("Entropy pooling on %s, %s%s",InpSymbols,EnumToString(InpTimeframe),
               (InpAsOf>0 ? ", as of "+TimeToString(InpAsOf,TIME_DATE) : ""));
   Print("================================================================");

   CEpScenarios scen;
   if(!LoadScenarios(scen,symbols))
      return;
   if(InpAsOf>0 && !scen.KeepUntil(InpAsOf))
     {
      Print("EP_Demo: too little history before the replay date");
      return;
     }

//--- the book as a column, then its trailing volatility as the state
   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);
   if(state_col<0)
      return;
   scen.Column(book_col,book);
   const int T=scen.Rows();

   vector z;
   scen.Column(state_col,z);
   const double z_today=scen.Next(state_col);
   PrintFormat("  scenarios        : %d joint rows, %s to %s",T,
               TimeToString(scen.Time(0),TIME_DATE),TimeToString(scen.Time(T-1),TIME_DATE));
   PrintFormat("  book             : %s",InpWeights);
   PrintFormat("  state            : book's %d-bar volatility, today %.3f%% (ann. %.2f%%)",
               InpStateWindow,100.0*z_today,100.0*z_today*MathSqrt(252.0));

//--- priors: what a lookback silently assumes, and two explicit choices
   vector p_full,p_roll,p_decay,p_kern;
   EpUniform(T,p_full);
   EpRollingWindow(T,InpWindow,p_roll);
   EpExponentialDecay(T,InpHalfLife,p_decay);
   AddStage("full history","full",p_full);
   AddStage(StringFormat("rolling %d",InpWindow),"rolling",p_roll);
   AddStage(StringFormat("decay hl %.0f",InpHalfLife),"decay",p_decay);
   const int decay=2;

   const double z_sd=EpVolatility(z,p_decay);
   EpKernel(z,z_today,InpKernelBand*z_sd,p_decay,p_kern);
   AddStage("decay x kernel","kernel",p_kern);

   CEpSolver solver;
   Print("--- solves -----------------------------------------------------");

//--- condition on the state: tilt the decay prior to today's level, within budget
   vector p_state;
   SEpSolve rs;
   const double z_now=EpConditionState(scen,p_decay,state_col,z_today,InpMinEns,p_state,rs);
   PrintFormat("  %-24s %s: %d Newton steps, KL %.4f, ENS %.0f",
               "state conditioning",EnumToString(rs.status),rs.iterations,rs.entropy,rs.ens);
   if(z_now!=z_today)
      PrintFormat("  %-24s today's state leaves under %.0f scenarios; target pulled to %.3f%%",
                  "",InpMinEns,100.0*z_now);
   AddStage("decay | state","state",p_state);

//--- the views, stated against the conditioned distribution
   const int vol_col =scen.Find(InpVolSymbol);
   const int tail_col=scen.Find(InpTailSymbol);
   const int ca=scen.Find(InpCorrA),cb=scen.Find(InpCorrB);
   if(vol_col<0 || tail_col<0 || ca<0 || cb<0)
     {
      Print("EP_Demo: every view symbol must be in the basket");
      return;
     }
   vector xv;
   scen.Column(vol_col,xv);
   const double vol_target=InpVolScale*EpVolatility(xv,p_state);

   CEpViews views;
   views.AddMean(state_col,EP_EQUAL,z_now);
   const int v_vol =views.AddVolatility(vol_col,EP_EQUAL,vol_target);
   const int v_tail=views.AddTail(tail_col,InpTailMove,EP_AT_LEAST,InpTailProb);
   const int v_corr=views.AddCorrelation(ca,cb,EP_AT_MOST,InpCorr);
   if(InpCorrPinVol)
     {
      vector xa,xb;
      scen.Column(ca,xa);
      scen.Column(cb,xb);
      if(ca!=vol_col)
         views.AddVolatility(ca,EP_EQUAL,EpVolatility(xa,p_state));
      if(cb!=vol_col)
         views.AddVolatility(cb,EP_EQUAL,EpVolatility(xb,p_state));
     }
   if(v_vol<0 || v_tail<0 || v_corr<0)
     {
      Print("EP_Demo: a view could not be added, check its inputs");
      return;
     }
   vector p_post;
   SEpSolve rv;
   const bool pooled=views.Pool(scen,p_decay,solver,p_post,rv);
   ReportSolve("state and views",rv,views);
   if(!pooled)
      return;
   AddStage("decay | state + views","views",p_post);
   if(InpConfidence<1.0)
     {
      vector p_mix;
      EpBlend(p_state,p_post,InpConfidence,p_mix);
      AddStage(StringFormat("blend c=%.2f",InpConfidence),"blend",p_mix);
     }

   Print("--- views -------------------------------------------------------");
   PrintFormat("  %-34s %12s %12s %12s %12s","view","before","after","multiplier","nats/unit");
   for(int k=0;k<views.Total();k++)
      PrintFormat("  %-34s %12.6g %12.6g %12.4f %12.4g",views.Describe(k,scen),
                  views.Achieved(k,scen,p_state),views.Achieved(k,scen,p_post),
                  solver.Multiplier(k),solver.MarginalCost(k));
   Print("  before = conditioned on the state only; multiplier 0 = the view was already true");

   ReportStages(book,decay);
   ReportSymbols(scen,n);
   ReportWorstDays(scen,book,8);
   Print("================================================================");
  }
//+------------------------------------------------------------------+