//+------------------------------------------------------------------+
//|                                                    EP_Checks.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 "Checks the entropy pooling library against results known in closed form,"
#property description "on synthetic Gaussian scenarios: identities, the exponential tilt, bounds,"
#property description "infeasibility, marginal cost, second-moment views and the scenario budget."

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

input int             InpRows        = 200000;             // synthetic scenarios
input int             InpSeed        = 7;                  // random seed

int g_failed=0;

//+------------------------------------------------------------------+
//| One line per check: the measured miss against its tolerance.     |
//+------------------------------------------------------------------+
void Check(const string what,const double miss,const double tol)
  {
   const bool ok=(MathIsValidNumber(miss) && MathAbs(miss)<=tol);
   if(!ok)
      g_failed++;
   PrintFormat("  %-52s %10.2e  (tol %.0e)  %s",what,miss,tol,(ok ? "ok" : "FAILED"));
  }

//+------------------------------------------------------------------+
//| Three correlated Gaussian returns, T rows.                       |
//+------------------------------------------------------------------+
void Gaussian(const int T,const int seed,CEpScenarios &scen)
  {
   MathSrand(seed);
   const double L[3][3]={{0.010,0.0,0.0},{0.006,0.008,0.0},{-0.002,0.003,0.012}};
   matrix x;
   x.Init(T,3);
   int err=0;
   for(int t=0;t<T;t++)
     {
      double e[3];
      for(int i=0;i<3;i++)
         e[i]=MathRandomNormal(0.0,1.0,err);
      for(int i=0;i<3;i++)
         x[t][i]=L[i][0]*e[0]+L[i][1]*e[1]+L[i][2]*e[2];
     }
   string names[]={"A","B","C"};
   datetime none[];
   scen.SetData(x,names,none);
  }

//+------------------------------------------------------------------+
double MaxDiff(const vector &a,const vector &b)
  {
   double d=0.0;
   for(ulong i=0;i<a.Size();i++)
      d=MathMax(d,MathAbs(a[i]-b[i]));
   return d;
  }

//+------------------------------------------------------------------+
//| Posterior entropy for a mean view on A at a given target.        |
//+------------------------------------------------------------------+
double EntropyAt(const CEpScenarios &scen,const vector &p,CEpSolver &solver,const double target,
                 double &cost)
  {
   CEpViews v;
   v.AddMean(0,EP_EQUAL,target);
   vector q;
   SEpSolve r;
   v.Pool(scen,p,solver,q,r);
   cost=solver.MarginalCost(0);
   return r.entropy;
  }

//+------------------------------------------------------------------+
void OnStart(void)
  {
   CEpScenarios scen;
   Gaussian(InpRows,InpSeed,scen);
   const int T=scen.Rows();
   vector p,a,b,q;
   EpUniform(T,p);
   scen.Column(0,a);
   scen.Column(1,b);
   const double ma=EpMean(a,p),sa=EpVolatility(a,p),sb=EpVolatility(b,p);
   CEpSolver solver;
   SEpSolve  r;
   PrintFormat("EP_Checks: %d Gaussian scenarios",T);

//--- 1. with nothing to say, the posterior is the prior
   CEpViews slack;
   slack.AddMean(0,EP_AT_MOST,ma+5.0*sa);
   slack.Pool(scen,p,solver,q,r);
   Check("slack view leaves the prior, max |q - p| * T",MaxDiff(q,p)*T,1e-9);
   Check("slack view multiplier",solver.Multiplier(0),1e-12);

//--- 2. a half-sigma mean view on Gaussian rows is an exponential tilt
   CEpViews tilt;
   tilt.AddMean(0,EP_EQUAL,ma+0.5*sa);
   tilt.Pool(scen,p,solver,q,r);
   const double shift=0.5*sa;
   const double mb_bl=EpMean(b,p)+EpCorrelation(a,b,p)*sb/sa*shift;
   Check("mean view met",EpMean(a,q)-(ma+shift),1e-12);
   Check("relative entropy against z^2/2",r.entropy-0.125,5e-3);
   Check("mean of B against the Black-Litterman mean, / shift",(EpMean(b,q)-mb_bl)/shift,5e-3);
   Check("vol of A unchanged by a mean view, relative",EpVolatility(a,q)/sa-1.0,5e-3);
   Check("vol of B unchanged by a mean view, relative",EpVolatility(b,q)/sb-1.0,5e-3);

//--- 3. a binding >= view gives the equality posterior; a slack one the prior
   vector qe=q;
   CEpViews bind;
   bind.AddMean(0,EP_AT_LEAST,ma+0.5*sa);
   bind.Pool(scen,p,solver,q,r);
   Check("binding >= view against the equality posterior, * T",MaxDiff(q,qe)*T,1e-9);

//--- 4. a view no reweighting can reach is reported
   CEpViews reach;
   reach.AddMean(0,EP_EQUAL,a.Max()+0.001);
   const bool pooled=reach.Pool(scen,p,solver,q,r);
   Check("unreachable mean view reported as infeasible",(!pooled && r.status==EP_INFEASIBLE ? 0.0 : 1.0),0.0);

//--- 5. the multiplier is the price of moving the target
   const double c0=ma+0.3*sa,h=1e-6*sa;
   double cost=0.0,unused=0.0;
   const double up=EntropyAt(scen,p,solver,c0+h,unused);
   const double dn=EntropyAt(scen,p,solver,c0-h,unused);
   EntropyAt(scen,p,solver,c0,cost);
   Check("marginal cost against a finite difference, relative",cost/((up-dn)/(2.0*h))-1.0,1e-6);

//--- 6. volatility and correlation views land after re-referencing passes
   CEpViews moments;
   moments.AddVolatility(0,EP_EQUAL,1.3*sa);
   moments.AddCorrelation(0,1,EP_EQUAL,0.80);
   moments.Pool(scen,p,solver,q,r);
   Check("vol view met, relative",EpVolatility(a,q)/(1.3*sa)-1.0,1e-9);
   Check("correlation view met",EpCorrelation(a,b,q)-0.80,1e-9);
   PrintFormat("  second-moment views took %d passes and %u ms",moments.Passes(),r.ms);

//--- 7. a tail probability view is linear and exact
   CEpViews tail;
   tail.AddTail(0,ma-2.0*sa,EP_EQUAL,0.05);
   tail.Pool(scen,p,solver,q,r);
   Check("tail probability view met",EpTailProbability(a,q,ma-2.0*sa)-0.05,1e-12);

//--- 8. conditioning never leaves fewer scenarios than the budget
   double absa[];
   ArrayResize(absa,T);
   for(int t=0;t<T;t++)
      absa[t]=MathAbs(a[t]);
   CEpScenarios state;
   Gaussian(T,InpSeed,state);
   const int zc=state.AddColumn("|A|",EP_COLUMN_STATE,absa,0.0);
   vector z;
   state.Column(zc,z);
   const double used=EpConditionState(state,p,zc,z.Max()+0.01,500.0,q,r);
   Check("state beyond history keeps ENS at the budget, relative",r.ens/500.0-1.0,1e-3);
   Check("pulled-back target met",EpMean(z,q)-used,1e-12);

   PrintFormat("EP_Checks: %s",(g_failed==0 ? "all checks passed" :
               StringFormat("%d check(s) FAILED",g_failed)));
  }
//+------------------------------------------------------------------+