//+------------------------------------------------------------------+
//|                                                 distribution.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"base.mqh"
#include<Math\Stat\T.mqh>
#include<Math\Stat\Gamma.mqh>
#include"..\Utility\igami.mqh"
#include"..\Utility\igam.mqh"
//+------------------------------------------------------------------+
//| constraints container                                            |
//+------------------------------------------------------------------+
struct Constraints
  {
   matrix            _one;
   vector            _two;

                     Constraints(void)
     {
      _one = matrix::Zeros(0,0);
      _two = vector::Zeros(0);
     }
                    ~Constraints(void)
     {
     }
                     Constraints(Constraints &other)
     {
      _one = other._one;
      _two = other._two;
     }
   void              operator=(Constraints &other)
     {
      _one = other._one;
      _two = other._two;
     }
  };

//+------------------------------------------------------------------+
//| base class for distributions                                     |
//+------------------------------------------------------------------+
class CDistribution
  {
protected:
   ulong             m_num_params;
   vector            m_parameters;
   vector            m_vector;
   string            m_name;
   uint              m_seed;
   CHighQualityRandStateShell m_generator;
   ENUM_DISTRIBUTION_MODEL m_dist_model;
   bool              m_initialized;
   vector            _check_constraints(vector &params_)
     {
      matrix bounds_ = bounds(vector::Zeros(0));
      ulong nparams = params_.Size();
      if(nparams!=bounds_.Rows())
        {
         Print(__FUNCTION__, " error ", "parameters must have ",bounds_.Rows()," elements");
         return vector::Zeros(0);
        }
      if(!bounds_.Rows())
         return vector::Zeros(0);
      for(ulong i = 0; i<bounds_.Rows(); ++i)
        {
         if(!(bounds_[i,0]<=params_[i]<=bounds_[i,1]))
           {
            Print(__FUNCTION__, " error ", ". Bounds violated ");
            return vector::Zeros(0);
           }
        }
      return params_;
     }

   bool              _initialize_generator(int rstate, bool use_alglib=true)
     {
      if(use_alglib)
        {
         m_seed = uint(MathAbs(rstate));

         if(rstate)
            CHighQualityRand::HQRndSeed((int)m_seed,(int)m_seed+1,m_generator.GetInnerObj());
         else
            CHighQualityRand::HQRndRandomize(m_generator.GetInnerObj());
         return true;

        }
      else
        {
         m_seed = uint(MathAbs(rstate));
         MathSrand(m_seed);
        }
      return true;
     }

   virtual bool      _initialize(vector& distribution_params, int seed=0)
     {
      return false;
     }
public:
                     CDistribution(void):m_name(NULL),
                     m_dist_model(WRONG_VALUE),
                     m_num_params(0),
                     m_initialized(false),
                     m_parameters(vector::Zeros(0))
     {
     }

                     CDistribution(vector& params, int seed=0)
     {
      initialize(params,seed);
     }

                     CDistribution(CDistribution& other)
     {
      m_name = other.name();
      m_dist_model = other.distribution_type();
      m_num_params = other.numParams();
      m_parameters = other.get_params();
      m_seed = other.randomState();
      m_initialized = other.is_initialized();
     }

   void              operator=(CDistribution& other)

     {
      m_name = other.name();
      m_dist_model = other.distribution_type();
      m_num_params = other.numParams();
      m_parameters = other.get_params();
      m_seed = other.randomState();
      m_initialized = other.is_initialized();
     }

                    ~CDistribution(void)
     {

     }
   virtual bool      initialize(vector& distribution_params, int seed=0)
     {
      if(is_initialized())
        {
         m_initialized = false;
        }

      return _initialize(distribution_params,seed);
     }
   virtual ENUM_DISTRIBUTION_MODEL distribution_type(void)
     {
      return m_dist_model;
     }
   virtual bool      is_initialized(void) { return m_initialized; }
   virtual CHighQualityRandStateShell* generator(void) { return GetPointer(m_generator); }
   virtual uint      randomState(void) { return m_seed; }
   virtual ulong     numParams(void)   { return m_num_params; }
   virtual string    name(void)        { return m_name; }
   virtual vector    get_params(void)  { return m_parameters; }
   void              set_params(vector &parameters)
     {
      m_parameters = parameters;
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    rng(ulong nsize)
     {
      matrix out = simulator(nsize);
      return out.Col(0);
     }
   virtual matrix    rng(ulong nrows,ulong ncols)
     {
      return simulator(nrows,ncols);
     }
   virtual Constraints constraints(void)
     {
      return Constraints();
     }
   virtual matrix    bounds(vector& resids)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individial = false)
     {
      return vector::Zeros(0);
     }
   virtual vector    startingValues(vector &resids)
     {
      return vector::Zeros(0);
     }
   virtual vector    ppf(vector &pits, vector &parameters)
     {
      return vector::Zeros(0);
     }
   virtual vector    cdf(vector &resids, vector &parameters)
     {
      return vector::Zeros(0);
     }
   virtual double    moment(int n, vector& parameters)
     {
      return EMPTY_VALUE;
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      return EMPTY_VALUE;
     }
  };
//+------------------------------------------------------------------+
//|Standard normal distribution for use with ARCH models             |
//+------------------------------------------------------------------+
class CNormal: public CDistribution
  {
private:

   int               factorial2(int k)
     {
      int out = k;
      for(int i = k-2; i>1; out*=i, i-=2);
      return out;
     }

   double            munp(ulong n)
     {
      if(!n)
         return 1.;
      if(MathMod(double(n),2.)==0)
         return (double)factorial2(int(n)-1);
      else
         return 0.;
     }
   vector            moment_from_stats(ulong n, double mu,double mu2, double g1, double g2)
     {
      vector out(1);
      if(n==0)
        {
         out[0] = 1.0;
         return out;
        }
      else
         if(n==1)
           {
            if(!mu)
               out[0] = munp(1);
            else
               out[0] = mu;
           }
         else
            if(n==2)
              {
               if(!mu2 || !mu)
                  out[0] = munp(2);
               else
                  out[0] = mu2+mu*mu;
              }
            else
               if(n == 3)
                 {
                  if(!g1 || !mu2)
                     out[0] = munp(3);
                  else
                    {
                     double mu3 = g1*pow(mu2,1.5);
                     out[1] = mu3+3.0*mu*mu2+mu*mu*mu;
                    }
                 }
               else
                  if(n==4)
                    {
                     if(!g1 || !g2)
                        out[0] = munp(4);
                     else
                       {
                        double mu3, mu4;
                        mu4 = (g2+3.0)*pow(mu2,2.0);
                        mu3 = g1*pow(mu2, 1.5);
                        out[1] = mu4+4*mu*mu3+6*mu*mu*mu2+mu*mu*mu*mu;
                       }
                    }
                  else
                     out[1] = munp(n);

      return out;

     }
   vector            norm_moment(ulong order, double loc = 0.0, double scale = 1.)
     {
      ulong n = order;
      bool i0 = (scale>0.);
      bool i1 = bool(!loc);
      bool i2 = bool(loc);

      double mu,mu1,g1, g2;
      mu = mu1 = g1 = g2 = 0.0;
      if(n>0 && n<5)
        {
         mu = 0.0;
         mu1 = 1.0;
         g1 = 0.0;
         g2 = 0.0;
        }
      vector val(1);
      val = moment_from_stats(n,mu,mu1,g1,g2);
      val[0] = pow(scale, order)*val[0];
      return val;
     }
   virtual bool      _initialize(vector& distribution_params, int seed=0) override
     {
      m_name = "Normal";
      m_parameters = distribution_params;
      m_num_params = 0;
      m_initialized =  _initialize_generator(seed,true);
      return m_initialized;
     }
public:
                     CNormal(void)
     {
      m_dist_model = DIST_NORMAL;
     }
                     CNormal(vector& params, int seed=0)
     {
      initialize(params,seed);
     }

                    ~CNormal(void)
     {
     }
   virtual Constraints  constraints(void)  override
     {
      Constraints out;
      return out;
     }
   virtual matrix    bounds(vector &resids)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individual = false) override
     {
      vector lls = -0.5 * (log(2 * M_PI) + log(sigmas) + pow(resids, 2.0) / sigmas) ;
      if(individual)
         return lls;
      else
        {
         vector out(1);
         out[0] = lls.Sum();
         return out;
        }
     }
   virtual vector    startingValues(vector& resids) override
     {
      return vector::Zeros(0);
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1) override
     {
      CMatrixDouble out;
      out.Resize(nsize,ncols);
      CHighQualityRand::HQRndNormalM(m_generator.GetInnerObj(),int(nsize),int(ncols),out);
      return out.ToMatrix();
     }
   virtual vector    ppf(vector &pits, vector &parameters) override
     {
      _check_constraints(parameters);
      vector out(pits.Size());
      for(ulong i = 0; i<out.Size(); CNormalDistr::InvNormalCDF(pits[i]), ++i);
      return out;
     }
   virtual vector    cdf(vector &resids, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);
      vector out(resids.Size());
      for(ulong i = 0; i<out.Size(); out[i] = CNormalDistr::NormalCDF(resids[i]), ++i);
      return out;
     }
   virtual double    moment(int n, vector& parameters)
     {
      set_params(parameters);
      if(n<0)
         return double("nan");
      vector out =  norm_moment(ulong(n));
      return out[0];
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      set_params(parameters);
      if(n<0)
         return double("nan");
      else
         if(n == 0)
            return CNormalDistr::NormalCDF(z);
         else
            if(n == 1)
               return -1.*CNormalDistr::NormalPDF(z);
            else
               return -pow(z,double(n-1))*CNormalDistr::NormalPDF(z) + double(n-1)* partialMoment(n-2,parameters,z);
     }
  };
//+------------------------------------------------------------------+
//| Standardized Student's distribution for use with ARCH models     |
//+------------------------------------------------------------------+
class CStudentsT: public CDistribution
  {
private:
   double            moment_from_stats(int n, double _mu,double _mu2, double _g1, double _g2, double df, double loc=0., double scale = 1.)
     {
      double out;
      if(n==0)
         return 1.;
      else
         if(n==1)
           {
            if(_mu == EMPTY_VALUE)
               out = double("nan");
            else
               out = _mu;
           }
         else
            if(n==2)
              {
               if(_mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE)
                  out = double("nan");
               else
                  out = _mu2 + _mu*_mu;
              }
            else
               if(n==3)
                 {
                  if(_g1  == EMPTY_VALUE|| _mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE)
                     out = double("nan");
                  else
                    {
                     double mu3 = _g1*pow(_mu2,1.5);
                     out = mu3+3.*_mu*_mu2+_mu*_mu*_mu;
                    }
                 }
               else
                  if(n==4)
                    {
                     if(_g1  == EMPTY_VALUE||_g2  == EMPTY_VALUE||_mu2  == EMPTY_VALUE||_mu == EMPTY_VALUE)
                        out = double("nan");
                     else
                       {
                        double mu4 =(_g2+3.)*pow(_mu2,2.);
                        double mu3 = _g1*pow(_mu2,1.5);
                        out = mu4+4.*_mu*mu3+6.*_mu*_mu*_mu2+_mu*pow(_mu,3);
                       }
                    }
                  else
                     out = double("nan");

      return out;
     }
   vector            _moment(int order, double df,double loc = 0.0,double scale = 1.)
     {
      vector out(0);

      ulong n = ulong(order);
      bool i0,i1,i2;
      i0 = i1 = i2 = false;
      i0 = (scale>0);
      i1 = (i0 & loc == 0.0);
      i2 = (i0 & loc != 0.0);

      if(order<0)
        {
         Print(__FUNCTION__, " Moment must be positive ");
         return vector::Zeros(0);
        }
      double mu,mu2,g1,g2;
      mu = mu2=g1 = g2 = EMPTY_VALUE;
      if(0<n && n<5)
        {
         mu = df>1.?0.0:double("inf");
         if(df>1. && df<=2.0)
            mu2 = double("inf");
         else
            if(df>2. && MathClassify(df) == FP_NORMAL)
               mu2 = df/(df-2.);
            else
               if(MathClassify(df) == FP_INFINITE)
                  mu2 = 1.;

         g1 = (df>3.)?0.0:double("nan");
         if(df>2. && df<=4.)
            g2 = double("inf");
         if(df>4. && MathClassify(df) == FP_NORMAL)
            g2 = 6./(df-4.);
         else
            if(MathClassify(df) == FP_INFINITE)
               g2 = 0.;

        }
      vector val(1);
      val[0] = moment_from_stats(int(n),mu,mu2,g1,g2,df,loc,scale);
      if(!i0)
         val[0] = double("nan");
      if(i1 && loc == 0.0)
         val[0] = pow(scale,n)*val[0];
      if(i2)
        {
         vector mom(4);
         mom[0] = mu;
         mom[1] = mu2;
         mom[2] = g1;
         mom[3] = g2;
         double res = 0.;
         double fac = scale/loc;
         double valk;
         for(int i = 0; i<int(n); ++i)
           {
            valk = moment_from_stats(i,mu,mu2,g1,g2,loc,scale);
            res += double(np::comb(n,i,true))*pow(fac,i)*valk;
           }
         res+=pow(fac,n)*val[0];
         res*=pow(loc,n);
         val[0] = res;
        }
      return val;
     }
   double            _ord_t_partial_moment(int n, double z, double nu)
     {
      double mom = double("nan");
      int ecode = 0;
      if(n<0 || double(n)>nu)
        {
         Print(__FUNCTION__, " n out of bounds ");
         return mom;
        }
      else
         if(n==0)
           {
            mom = MathCumulativeDistributionT(z,nu,ecode);
            if(ecode)
              {
               Print(__FUNCTION__, " error with cdf ", ecode);
               return double("nan");
              }
           }
         else
            if(n==1)
              {
               double c = CGammaFunc::GammaFunc(0.5+(n+1.))/(sqrt(nu*M_PI)*CGammaFunc::GammaFunc(0.5*nu));
               double e = 0.5*(nu+1.);
               mom = (0.5*(c*nu)/(1.-e))*(pow(1.+pow(z,2.)/nu,1.-e));
              }
            else
              {
               ecode = 0;
               double t1 = pow(z,n-1.)*(nu+pow(z,2.))*MathProbabilityDensityT(z,nu,ecode);
               if(ecode)
                 {
                  Print(__FUNCTION__, " error with cdf ", ecode);
                  return double("nan");
                 }
               double t2 = (n-1)*nu*_ord_t_partial_moment(n-2,z,nu);
               mom = (1./(double(n) - nu))*(t1-t2);
              }
      return mom;
     }
   virtual bool      _initialize(vector& distribution_params, int seed=0) override
     {
      m_name = "Student's T";
      m_num_params = 1;
      m_parameters = distribution_params;
      m_initialized= _initialize_generator(seed,false);
      return m_initialized;
     }
public:
                     CStudentsT(void)
     {
      m_dist_model = DIST_STUDENT;
     }
                     CStudentsT(vector& params, int seed=0)
     {
      initialize(params,seed);
     }
                    ~CStudentsT(void)
     {
     }
   virtual Constraints      constraints(void)  override
     {
      Constraints out;
      matrix x = {{1.0},{-1.}};
      vector y = {2.05,-500.};
      out._one = x;
      out._two = y;
      return out;
     }
   virtual matrix    bounds(vector &resids) override
     {
      matrix out = {{2.05},{500.}};
      return out;
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individual = false) override
     {
      set_params(parameters);
      double nu = m_parameters[0];
      double ll = log(CGammaFunc::GammaFunc((nu+1.0)/2.0)) - log(CGammaFunc::GammaFunc((nu)/2.0)) - log(M_PI*(nu-2.0))/2.;
      vector lls = ll - (0.5*log(sigmas));
      lls-=((nu + 1.) / 2.) * (log(1. + pow(resids, 2.0) / (sigmas * (nu - 2.))));
      if(individual)
         return lls;
      else
        {
         vector out(1);
         out[0] = lls.Sum();
         return out;
        }
     }
   virtual vector    startingValues(vector& resids) override
     {
      double k = np::kurtosis(resids,false,true);
      double temp = k>3.75?(4.0 * k - 6.0) / (k - 3.0):12;
      double sv = MathMax(temp,4.0);
      vector out(1);
      out[0] = sv;
      return out;
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1) override
     {
      matrix out(nsize*ncols,1);
      vector temp;
      double vals[];
      double std_dev = sqrt(m_parameters[0]/(m_parameters[0] - 2.));
      if(!MathRandomT(m_parameters[0],int(nsize*ncols),vals) || !temp.Assign(vals) ||
         !out.Col(temp,0) || !out.Reshape(nsize,ncols))
        {
         Print(__FUNCTION__, " error sampling from student distribution ", GetLastError());
         return matrix::Zeros(0,0);
        }

      return out/(std_dev);
     }

   virtual vector    ppf(vector &pits, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);
      vector out(pits.Size());
      int ecode = 0;
      double var = m_parameters[0]/(m_parameters[0] -2.0);
      double scale = 1.0/sqrt(var);
      for(ulong i = 0; i<out.Size();  ++i)
        {
         out[i] = MathQuantileT(pits[i],m_parameters[0],ecode)*scale;
         if(ecode)
           {
            Print(__FUNCTION__, " ppf() error ", ecode);
            return vector::Zeros(0);
           }
        }
      return out;
     }
   virtual vector    cdf(vector &resids, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);
      vector out(resids.Size());
      int ecode = 0;
      double var = m_parameters[0]/(m_parameters[0] -2.0);
      double scale = 1./sqrt(var);
      vector y = resids/scale;
      for(ulong i = 0; i<out.Size();  ++i)
        {
         out[i] = (MathCumulativeDistributionT(y[i],m_parameters[0],ecode));
         if(ecode)
           {
            Print(__FUNCTION__, " cdf() error ", ecode);
            return vector::Zeros(0);
           }
        }
      return out;
     }
   virtual double    moment(int n, vector& parameters)
     {
      set_params(parameters);
      if(n<0)
         return double("nan");
      _check_constraints(parameters);
      double nu = parameters[0];
      double var =nu/(nu-2.);
      double scale = 1.0/sqrt(var);
      vector out = _moment(n,nu,0.,scale);
      return out[0];
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      set_params(parameters);
      double nu=parameters[0];
      double var = nu/(nu-2.);
      double scale = 1./sqrt(var);
      return pow(scale,n)*_ord_t_partial_moment(n,z/scale,nu);
     }
  };
//+---------------------------------------------------------------------+
//|Standardized Skewed Student's distribution for use with ARCH models  |
//+---------------------------------------------------------------------+
class CSkewStudent: public CDistribution
  {
private:
   double            _const_c(vector& parameters)
     {
      double eta = parameters[0];
      return log(CGammaFunc::GammaFunc((eta+1.)/2.)) - log(CGammaFunc::GammaFunc(eta/2.)) - log(M_PI*(eta-2.))/2.;
     }
   double            _const_b(vector& parameters)
     {
      double lam = parameters[1];
      double a = _const_a(parameters);
      return pow(1.+3.* pow(lam,2.)-pow(a,2.),0.5);
     }
   double            _const_a(vector& parameters)
     {
      double eta = parameters[0];
      double lam = parameters[1];
      double c = _const_c(parameters);
      return (4.*lam*exp(c)*(eta-2.)/(eta-1.));
     }

   double            _ord_t_partial_moment(int n, double z, double nu)
     {
      double mom = double("nan");
      int ecode = 0;
      if(n<0 || double(n)>nu)
        {
         Print(__FUNCTION__, " n out of bounds ");
         return mom;
        }
      else
         if(n==0)
           {
            mom = MathCumulativeDistributionT(z,nu,ecode);
            if(ecode)
              {
               Print(__FUNCTION__, " error with cdf ", ecode);
               return double("nan");
              }
           }
         else
            if(n==1)
              {
               double c = CGammaFunc::GammaFunc(0.5+(n+1.))/(sqrt(nu*M_PI)*CGammaFunc::GammaFunc(0.5*nu));
               double e = 0.5*(nu+1.);
               mom = (0.5*(c*nu)/(1.-e))*(pow(1.+pow(z,2.)/nu,1.-e));
              }
            else
              {
               ecode = 0;
               double t1 = pow(z,n-1.)*(nu+pow(z,2.))*MathProbabilityDensityT(z,nu,ecode);
               if(ecode)
                 {
                  Print(__FUNCTION__, " error with cdf ", ecode);
                  return double("nan");
                 }
               double t2 = (n-1)*nu*_ord_t_partial_moment(n-2,z,nu);
               mom = (1./(double(n) - nu))*(t1-t2);
              }
      return mom;
     }
   virtual bool      _initialize(vector& distribution_params, int seed=0) override
     {
      m_name = "Standardized Skew Student's T";
      m_num_params = 2;
      m_parameters = distribution_params;
      m_initialized= _initialize_generator(seed,true);
      return m_initialized;
     }
public:
                     CSkewStudent(void)
     {
      m_dist_model = DIST_SKEW_STUDENT;
     }
                     CSkewStudent(vector& params, int seed=0)
     {
      initialize(params,seed);
     }
                    ~CSkewStudent(void)
     {

     }
   virtual Constraints      constraints(void)  override
     {
      Constraints out;
      matrix x = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}};
      vector y = {2.05, -300.0, -1, -1};
      out._one = x;
      out._two = y;
      return out;
     }
   virtual matrix    bounds(vector &resids) override
     {
      matrix out = {{-1.,2.05},{1.,300.}};
      return out;
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individual = false) override
     {
      set_params(parameters);
      double eta = parameters[0];
      double lam = parameters[1];
      double const_c = _const_c(parameters);
      double const_a = _const_a(parameters);
      double const_b = _const_b(parameters);

      resids=resids/pow(sigmas,0.5);
      vector lls = log(const_b)+const_c - log(sigmas)/2.;
      if(MathAbs(lam)>=1.)
         lam = np::sign(lam)*(1.-1e-6);
      vector signs(resids.Size());
      for(ulong i = 0; i<signs.Size(); signs[i] = np::sign(resids[i]+const_a/const_b), ++i);
      vector llf_resid=pow((const_b*resids+const_a)/(1.+signs*lam),2.);
      lls -= (eta+1.)/2.*log(1.+llf_resid/(eta-2.));
      if(individual)
         return lls;
      else
        {
         vector out(1);
         out[0] = lls.Sum();
         return out;
        }
     }
   virtual vector    startingValues(vector& resids) override
     {
      double k = np::kurtosis(resids,false,true);
      double temp = k>3.75?(4.0 * k - 6.0) / (k - 3.0):12;
      double sv = (k>3.75)?MathMax((4.*k -6.)/(k-3.),4.):MathMax(12.,4.);
      vector out(2);
      out[0] = sv;
      out[1] = 0.0;
      return out;
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1) override
     {
      matrix out(nsize*ncols,1);
      vector uniforms(out.Rows());
      CHighQualityRand::HQRndContinuous(m_generator.GetInnerObj(),int(out.Rows()),uniforms);
      vector _ppf = ppf(uniforms,m_parameters);
      out.Col(_ppf,0);
      out.Reshape(nsize,ncols);
      return out;
     }

   virtual vector    ppf(vector &pits, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);
      vector out(pits.Size());
      double eta = m_parameters[0];
      double lam = m_parameters[1];

      //double c = _const_c(parameters);
      double a = _const_a(parameters);
      double b = _const_b(parameters);

      vector icdf = vector::Ones(pits.Size()) * -999.99;
      vector signs = vector::Zeros(pits.Size());
      int ecode = 0;
      for(ulong i  = 0; i<pits.Size(); ++i)
        {
         if(pits[i]<(1.-lam)/2.)
            icdf[i] = MathQuantileT((pits[i]/(1.-lam)),eta,ecode);
         else
            icdf[i] = MathQuantileT(((pits[i] - (1.-lam)/2)/(1.+lam)),eta,ecode);
         icdf[i] = icdf[i]*(1.+np::sign(pits[i] - (1.-lam)/2.)*lam)*pow(1.-2./eta,0.5-a);
        }

      icdf = icdf/b;
      return icdf;
     }
   virtual vector    cdf(vector &resids, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);

      double eta = m_parameters[0];
      double lam = m_parameters[1];

      double a = _const_a(parameters);
      double b = _const_b(parameters);
      double var = eta/(eta-2.);

      vector y1 = (b*resids+a)/(1.-lam)*sqrt(var);
      vector y2 = (b*resids+a)/(1.+lam)*sqrt(var);
      vector p = vector::Zeros(resids.Size());
      int ecode = 0;
      for(ulong i = 0; i<p.Size(); ++i)
        {
         if(resids[i]<(-a/b))
            p[i] = (1.-lam)*MathCumulativeDistributionT(y1[i],eta,ecode);
         else
            p[i] += ((1. - lam) / 2. + (1 + lam) * (MathCumulativeDistributionT(y2[i],eta,ecode) - 0.5));
        }
      return p;
     }
   virtual double    moment(int n, vector& parameters)
     {
      set_params(parameters);
      if(n<0)
         return double("nan");
      _check_constraints(parameters);

      double eta = m_parameters[0];
      double lam = m_parameters[1];

      double a = _const_a(parameters);
      double b = _const_b(parameters);

      double loc = -a/b;
      double lscale = sqrt(1.-2./eta)*(1.-lam)/b;
      double rscale = sqrt(1.-2./eta)*(1.+lam)/b;
      double rpmom,lpmom,lhs,rhs,mom = 0.;
      for(int k = 0; k<(n+1); ++k)
        {
         rpmom = (0.5 * (CGammaFunc::GammaFunc(0.5 * (k + 1)) * CGammaFunc::GammaFunc(0.5 * (eta - k)) * pow(eta,(0.5 * k)))/ (sqrt(M_PI) * CGammaFunc::GammaFunc(0.5 * eta)));
         lpmom = pow(-1, k) * rpmom;
         lhs = (1 - lam) * pow(lscale,k) * pow(loc, (n - k)) * lpmom;
         rhs = (1 + lam) * pow(rscale, k) * pow(loc, (n - k)) * rpmom;
         mom += np::comb(ulong(n),double(k)) * (lhs + rhs);
        }
      return mom;
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      set_params(parameters);
      if(n<0 || double(n)>=m_parameters[0])
         return double("nan");
      _check_constraints(parameters);

      double eta = m_parameters[0];
      double lam = m_parameters[1];

      double a = _const_a(parameters);
      double b = _const_b(parameters);

      double loc = -a/b;
      double lscale = sqrt(1.-2./eta)*(1.-lam)/b;
      double rscale = sqrt(1.-2./eta)*(1.+lam)/b;
      double lbound,lhs,rhs,mom = 0.;
      for(int k = 0; k<(n+1); ++k)
        {
         lbound = MathMin(z,loc);
         lhs = ((1 - lam)* pow(loc, (n - k))* pow(lscale,k)*_ord_t_partial_moment(k,(lbound - loc) / lscale, eta));
         if(z>loc)
            rhs = ((1 + lam)* pow(loc, (n - k))* pow(rscale, k)* (_ord_t_partial_moment(k, (z - loc) / rscale, eta) - _ord_t_partial_moment(k,0.0,eta)));
         else
            rhs = 0.0;
         mom +=np::comb(ulong(n),double(k))*(lhs+rhs);
        }
      return mom;
     }

  };
//+------------------------------------------------------------------+
//|Generalized Error distribution for use with ARCH models           |
//+------------------------------------------------------------------+
class CGeneralizedError: public CDistribution
  {
private:
   double            _ord_gennorm_partial_moment(int n, int z, double beta)
     {
      if(n<0)
         return double("nan");

      double w = 0.5*beta/CGammaFunc::GammaFunc(1./beta);
      double lz = pow(MathAbs(MathMin(z,0)),beta);
      double lterm = (w* (pow((-1), n))* (1 / beta)* CGammaFunc::GammaFunc((n + 1) / beta)* special::cephes::igamc((n + 1) / beta, lz));
      double rz = pow(MathMax(0,z),beta);
      double rterm = w * (1 / beta) * CGammaFunc::GammaFunc((n + 1) / beta) * special::cephes::igamc((n + 1) / beta, rz);
      return lterm+rterm;
     }
   void              stats(double beta, double &a, double &b, double &c, double &d)
     {
      double c1 = log(CGammaFunc::GammaFunc(1./beta));
      double c3 = log(CGammaFunc::GammaFunc(3./beta));
      double c5 = log(CGammaFunc::GammaFunc(5./beta));

      a = 0.;
      b = exp(c3-c1);
      c = 0.;
      d = exp(c5+c1-2.*c3)-3.;
      return;
     }
   double            moment_from_stats(int n, double _mu,double _mu2, double _g1, double _g2, double df, double loc=0., double scale = 1.)
     {
      double out;
      if(n==0)
         return 1.;
      else
         if(n==1)
           {
            if(_mu == EMPTY_VALUE)
               out = double("nan");
            else
               out = _mu;
           }
         else
            if(n==2)
              {
               if(_mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE)
                  out = double("nan");
               else
                  out = _mu2 + _mu*_mu;
              }
            else
               if(n==3)
                 {
                  if(_g1  == EMPTY_VALUE|| _mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE)
                     out = double("nan");
                  else
                    {
                     double mu3 = _g1*pow(_mu2,1.5);
                     out = mu3+3.*_mu*_mu2+_mu*_mu*_mu;
                    }
                 }
               else
                  if(n==4)
                    {
                     if(_g1  == EMPTY_VALUE||_g2  == EMPTY_VALUE||_mu2  == EMPTY_VALUE||_mu == EMPTY_VALUE)
                        out = double("nan");
                     else
                       {
                        double mu4 =(_g2+3.)*pow(_mu2,2.);
                        double mu3 = _g1*pow(_mu2,1.5);
                        out = mu4+4.*_mu*mu3+6.*_mu*_mu*_mu2+_mu*pow(_mu,3);
                       }
                    }
                  else
                     out = double("nan");

      return out;
     }
   vector            _moment(int order, double df,double loc = 0.0,double scale = 1.)
     {
      vector out(0);

      ulong n = ulong(order);
      bool i0,i1,i2;
      i0 = i1 = i2 = false;
      i0 = (scale>0);
      i1 = (i0 & loc == 0.0);
      i2 = (i0 & loc != 0.0);

      if(order<0)
        {
         Print(__FUNCTION__, " Moment must be positive ");
         return vector::Zeros(0);
        }
      double mu,mu2,g1,g2;
      mu = mu2=g1 = g2 = EMPTY_VALUE;
      if(0<n && n<5)
        {
         stats(df,mu,mu2,g1,g2);
        }
      vector val(1);
      val[0] = moment_from_stats(int(n),mu,mu2,g1,g2,df,loc,scale);
      if(!i0)
         val[0] = double("nan");
      if(i1 && loc == 0.0)
         val[0] = pow(scale,n)*val[0];
      if(i2)
        {
         vector mom(4);
         mom[0] = mu;
         mom[1] = mu2;
         mom[2] = g1;
         mom[3] = g2;
         double res = 0.;
         double fac = scale/loc;
         double valk;
         for(int i = 0; i<int(n); ++i)
           {
            valk = moment_from_stats(i,mu,mu2,g1,g2,loc,scale);
            res += np::comb(n,i,true)*pow(fac,i)*valk;
           }
         res+=pow(fac,n)*val[0];
         res*=pow(loc,n);
         val[0] = res;
        }
      return val;
     }
   double            _ppf(double q, double beta)
     {
      double c = np::sign(q - 0.5);
      return c*pow(special::cephes::igamci(1./beta,(1.0 + c) - 2.0*c*q),1./beta);
     }
   double            gennorm_ppf(double nu,double loc, double scale, double q)
     {
      if(q<=0 || q>=1 || scale<=0)
        {
         double lb = -double("inf");
         double ub = double("inf");
         if(q == 0.)
            return lb;
         else
            if(q==1.)
               return ub;
        }
      return _ppf(q,nu)*scale+loc;
     }
   vector            gennorm_ppfv(double nu,double loc, double scale, vector& q)
     {
      vector out(q.Size());

      for(ulong i = 0; i<out.Size(); ++i)
         out[i] = gennorm_ppf(nu,loc,scale,q[i]);

      return out;
     }
   virtual bool      _initialize(vector& distribution_params, int seed=0) override
     {
      m_name = "Generalized Error Distribution";
      m_num_params = 1;
      m_parameters = distribution_params;
      m_initialized= _initialize_generator(seed,false);
      return m_initialized;
     }
public:
                     CGeneralizedError(void)
     {
      m_dist_model = DIST_GEN_ERROR;
     }
                     CGeneralizedError(vector& params, int seed=0)
     {
      initialize(params,seed);
     }
                    ~CGeneralizedError(void)
     {

     }
   virtual Constraints      constraints(void)  override
     {
      Constraints out;
      matrix x = {{1}, {-1}};
      vector y = {1.01, -500.0};
      out._one = x;
      out._two = y;
      return out;
     }
   virtual matrix    bounds(vector &resids) override
     {
      matrix out = {{1.01},{500.}};
      return out;
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individual = false) override
     {
      set_params(parameters);
      double nu = parameters[0];
      double logc = 0.5 * (-2. / nu * log(2.) + log(CGammaFunc::GammaFunc(1. / nu)) - log(CGammaFunc::GammaFunc(3. / nu)));
      double c = exp(logc);
      double lls = log(nu) - logc - log(CGammaFunc::GammaFunc(1./nu)) - (1.+1./nu)*log(2.);
      vector llsv = lls - (0.5 * log(sigmas));
      llsv = llsv - (0.5*pow(MathAbs(resids/(sqrt(sigmas)*c)),nu));
      if(individual)
         return llsv;
      else
        {
         vector out(1);
         out[0] = llsv.Sum();
         return out;
        }
     }
   virtual vector    startingValues(vector& resids) override
     {
      vector out(1);
      out[0] = 1.5;
      return out;
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1) override
     {
      double nu = m_parameters[0];
      double gammas[];
      double scale = sqrt(CGammaFunc::GammaFunc(3./nu)/CGammaFunc::GammaFunc(1./nu));
      double shape = 1./nu;
      if(!MathRandomGamma(shape,scale,int(nsize*ncols),gammas))
        {
         Print(__FUNCTION__," math gamma error ", GetLastError());
         return matrix::Zeros(0,0);
        }
      vector gammasv(nsize*ncols);
      gammasv.Assign(gammas);
      ArrayFree(gammas);

      for(ulong i = 0; i<(nsize*ncols); gammasv[i]*= (2.*(double)CHighQualityRand::HQRndUniformI(m_generator.GetInnerObj(),2)-1.), ++i);
      matrix out(nsize*ncols,1);

      out.Col(gammasv,0);
      out.Reshape(nsize,ncols);
      return out/scale;//do i have to divide by scale
     }

   virtual vector    ppf(vector &pits, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);
      double nu = m_parameters[0];
      double c1 = log(CGammaFunc::GammaFunc(1./nu));
      double c3 = log(CGammaFunc::GammaFunc(3./nu));

      double var = exp(c3-c1);
      double scale = 1./sqrt(var);

      return gennorm_ppfv(nu,0.0,scale,pits);
     }
   virtual vector    cdf(vector &resids, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);

      double nu = m_parameters[0];
      double c1 = log(CGammaFunc::GammaFunc(1./nu));
      double c3 = log(CGammaFunc::GammaFunc(3./nu));

      double var = exp(c3-c1);
      double scale = 1./sqrt(var);

      vector out = resids/scale;

      for(ulong i = 0; i<out.Size(); ++i)
        {
         switch(MathClassify(out[i]))
           {
            case FP_INFINITE:
               out[i] = 1.;
               break;
            case FP_NAN:
               out[i] = double("nan");
               break;
            default:
               out[i] = (0.5 + (0.5 * np::sign(out[i]))) - (0.5 * np::sign(out[i])) * special::cephes::igamc(1.0/nu, pow(MathAbs(out[i]),nu));
               break;
           }
        }
      return out;
     }
   virtual double    moment(int n, vector& parameters)
     {
      set_params(parameters);
      _check_constraints(parameters);

      double nu = m_parameters[0];
      double c1 = log(CGammaFunc::GammaFunc(1./nu));
      double c3 = log(CGammaFunc::GammaFunc(3./nu));

      double var = exp(c3-c1);
      double scale = 1./sqrt(var);

      vector out = _moment(n,nu,0.0,scale);
      return out[0];
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      set_params(parameters);
      _check_constraints(parameters);

      double nu = m_parameters[0];
      double c1 = log(CGammaFunc::GammaFunc(1./nu));
      double c3 = log(CGammaFunc::GammaFunc(3./nu));

      double var = exp(c3-c1);
      double scale = 1./sqrt(var);

      double mom = pow(scale,n)*_ord_gennorm_partial_moment(n,int(z/scale),nu);
      return mom;
     }
  };
//+------------------------------------------------------------------+
//|RNG                                                               |
//+------------------------------------------------------------------+
class BootstrapRng
  {
protected:
   vector            m_std_resid;
   ulong             m_start,m_index;
   CHighQualityRandStateShell m_rstate;
   CDistribution*     m_distribution;
   bool              m_initialized;
public:
                     BootstrapRng(void)
     {
      //m_rstate = NULL;
      m_distribution = NULL;
      m_initialized = false;
      m_start = m_index = ULONG_MAX;
      m_std_resid = vector::Zeros(0);
     }
                     BootstrapRng(ENUM_DISTRIBUTION_MODEL type, vector& parameters, int seed)
     {
      m_distribution = NULL;
      switch(type)
        {
         case DIST_NORMAL:
            m_distribution = new CNormal(parameters,seed);
            break;
         case DIST_STUDENT:
            m_distribution = new CStudentsT(parameters,seed);
            break;
         case DIST_SKEW_STUDENT:
            m_distribution = new CSkewStudent(parameters,seed);
            break;
         case DIST_GEN_ERROR:
            m_distribution = new CGeneralizedError(parameters,seed);
            break;
         default:
            Print(__FUNCTION__, " critical error. Failed instantiation of distribution ");
            break;
        }
      m_initialized = (m_distribution.is_initialized()==true);//((CheckPointer(m_distribution)!=POINTER_INVALID) && (m_distribution.initialize(parameters,seed)));
     }
                     BootstrapRng(vector& std_resid, ulong start,int seed = 0)
     {
      m_initialized = false;
      m_std_resid = std_resid;
      m_start = start;
      m_index = m_start;
      m_distribution = NULL;

      if(seed)
         CHighQualityRand::HQRndSeed(seed,seed+1,m_rstate.GetInnerObj());
      else
         CHighQualityRand::HQRndRandomize(m_rstate.GetInnerObj());
      m_initialized = true;
     }
                    ~BootstrapRng(void)
     {
      if(CheckPointer(m_distribution) == POINTER_DYNAMIC)
         delete m_distribution;
     }
   bool              is_initialized(void)
     {
      return m_initialized;
     }
   vector            rng(ulong size)
     {
      if(m_std_resid.Size())
        {
         if(m_index>=m_std_resid.Size())
           {
            Print(__FUNCTION__, " not enough data points ");
            return vector::Zeros(0);
           }
         vector index = vector::Zeros(size);
         for(ulong i = 0; i<size; index[i] = CHighQualityRand::HQRndUniformR(m_rstate.GetInnerObj()), ++i);
         index = floor(double(m_index)*index);
         m_index+=1;
         vector out = index;
         for(ulong i = 0; i<size; ++i)
            out[i] = m_std_resid[ulong(index[i])];
         return out;
        }
      return EMPTY_VECTOR;
     }
   matrix            rng(ulong size, ulong cols)
     {
      if(CheckPointer(m_distribution)!=POINTER_INVALID)
         return m_distribution.rng(size,cols);
      else
        {
         matrix out(size*cols,1);
         vector col = rng(size*cols);
         out.Col(col,0);
         out.Reshape(size,cols);
         return out;
        }
     }

  };
//+------------------------------------------------------------------+