//+------------------------------------------------------------------+
//|                                                 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 structure for optimization                 |
//| Holds linear inequality constraint matrices (Matrix * X >= Vector)|
//+------------------------------------------------------------------+
struct Constraints
  {
   matrix            _one; // Constraint coefficient mapping matrix (A)
   vector            _two; // Boundary value threshold vector (b)

   // Default Constructor
   Constraints(void)
     {
      _one = matrix::Zeros(0,0);
      _two = vector::Zeros(0);
     }
   // Destructor  
   ~Constraints(void)
     {
     }
   // Copy Constructor  
   Constraints(Constraints &other)
     {
      _one = other._one;
      _two = other._two;
     }
   // Overloaded Assignment Operator  
   void              operator=(Constraints &other)
     {
      _one = other._one;
      _two = other._two;
     }
  };

//+------------------------------------------------------------------+
//| Abstract base class defining the standard interface for continuous |
//| probability distributions used in conditional volatility models  |
//+------------------------------------------------------------------+
class CDistribution
  {
protected:
   ulong                      m_num_params;   // Total estimated parameters in the model
   vector                     m_parameters;   // Vector array housing parameter values
   vector                     m_vector;       // Multi-use workspace vector
   string                     m_name;         // Identifiable distribution title literal string
   uint                       m_seed;         // Active seed initialization flag state
   CHighQualityRandStateShell m_generator;    // Alglib high-quality random engine wrapper shell
   ENUM_DISTRIBUTION_MODEL    m_dist_model;   // Target model classification identifier
   bool                       m_initialized;  // Status flag verifying functional configuration

   // Internal sanity validator checking that parameters fall within valid logical bounds
   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);
         
      // Validate lower/upper boundaries sequentially: Matrix columns 0 and 1
      for(ulong i = 0; i < bounds_.Rows(); ++i)
        {
         if(!(bounds_[i,0] <= params_[i] && params_[i] <= bounds_[i,1]))
           {
            Print(__FUNCTION__, " error ", ". Bounds violated ");
            return vector::Zeros(0);
           }
        }
      return params_;
     }

   // Set seed states across standard RTL or high-quality native math environments
   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;
     }

   // Subclass hook targeting parameterized setup sequences
   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)
     {
     }

   // Prepares operational parameters and safety tracks before model evaluation loops
   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; }

   // Returns simulation matrices filled with generated random numbers
   virtual matrix    simulator(ulong nsize, ulong ncols = 1) { return matrix::Zeros(0,0); }
   
   // Random Number Generation (RNG) outputs formatted to vector or matrix shapes
   virtual vector    rng(ulong nsize)
     {
      matrix out = simulator(nsize);
      return out.Col(0);
     }
   virtual matrix    rng(ulong nrows, ulong ncols)
     {
      return simulator(nrows, ncols);
     }

   // Default stubs for bounds, likelihood evaluation and distribution moments
   virtual Constraints constraints(void) { return Constraints(); }
   virtual matrix    bounds(vector& resids) { return matrix::Zeros(0,0); }
   
   // Computes log-likelihood vectors over residuals vectors given their local variance sigmas
   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); } // Percent Point Function (Inverse CDF)
   virtual vector    cdf(vector &resids, vector &parameters) { return vector::Zeros(0); } // Cumulative Distribution Function
   virtual double    moment(int n, vector& parameters) { return EMPTY_VALUE; }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0) { return EMPTY_VALUE; } // Integrated partial density bounds
  };

//+------------------------------------------------------------------+
//| Standard Normal Distribution implementation (Gaussian)           |
//+------------------------------------------------------------------+
class CNormal: public CDistribution
  {
private:
   // Double factorial function evaluating product reductions: k!! = k * (k-2) * (k-4)...
   int               factorial2(int k)
     {
      int out = k;
      for(int i = k-2; i > 1; out*=i, i-=2);
      return out;
     }

   // Formulates non-central moments centered at zero for standardized normal terms
   double            munp(ulong n)
     {
      if(!n)
         return 1.;
      if(MathMod(double(n), 2.) == 0)
         return (double)factorial2(int(n)-1); // Even orders match double factorials
      else
         return 0.; // Odd standard symmetric distributions track zero
     }

   // Maps structural moments directly from standard parameters (Mean, Variance, Skewness, Kurtosis)
   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) { out[0] = (!mu) ? munp(1) : mu; }
      else if(n == 2) { out[0] = (!mu2 || !mu) ? munp(2) : 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;
     }

   // Converts standardized normal moment evaluations across given location/scale anchors
   vector            norm_moment(ulong order, double loc = 0.0, double scale = 1.)
     {
      ulong n = order;
      double 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; // Gaussian has zero shape parameters to estimate
      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); }

   // Normal distribution log-likelihood calculation: -0.5 * (ln(2*pi) + ln(sigma) + e^2 / sigma)
   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); }

   // Simulates multiple series vectors via Matrix tracking interfaces
   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();
     }

   // Inverse Normal CDF calculations map probabilities directly onto standardized values
   virtual vector    ppf(vector &pits, vector &parameters) override
     {
      _check_constraints(parameters);
      vector out(pits.Size());
      for(ulong i = 0; i < out.Size(); ++i)
         out[i] = CNormalDistr::InvNormalCDF(pits[i]);
      return out;
     }

   // Computes cumulative density vectors using traditional mathematical error integrations
   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];
     }

   // Evaluates conditional distributions inside partial integration bounds using recursive derivations
   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 t-Distribution for fat-tailed residuals   |
//+------------------------------------------------------------------+
class CStudentsT: public CDistribution
  {
private:
   // Maps empirical parameters to generalized moment profiles specifically tailored to Student's t domains
   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) { out = (_mu == EMPTY_VALUE) ? double("nan") : _mu; }
      else if(n == 2) { out = (_mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE) ? double("nan") : _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;
     }

   // Tracks internal moment values depending directly on active degrees of freedom (df) thresholds
   vector            _moment(int order, double df, double loc = 0.0, double scale = 1.)
     {
      vector out(0);
      ulong n = ulong(order);
      bool 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 = EMPTY_VALUE, mu2 = EMPTY_VALUE, g1 = EMPTY_VALUE, 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)
        {
         double res = 0., fac = scale / loc, 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;
     }

   // Integrates a localized snippet of the lower partial t-distribution curve
   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; // 1 shape parameter: degrees of freedom
      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) {}

   // Optimization bounds requiring df > 2.05 to maintain valid variance definitions
   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;
     }

   // Student's t Log-Likelihood accounting for fat-tailed kurtosis parameterizations
   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;
        }
     }

   // Estimates robust starting parameters using the sample kurtosis metrics
   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;
     }

   // Samples from the distribution and standardizes variance down to unity scale limits
   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 t-Distribution (Hansen 1994 setup)    |
//| Handles both heavy-tail kurtosis and asymmetric financial profiling  |
//+---------------------------------------------------------------------+
class CSkewStudent: public CDistribution
  {
private:
   // Internal distribution constants configuration (a, b, c) tracking skew adjustments
   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 Honest Student's T";
      m_num_params = 2; // 2 shape parameters: Tail thickness (eta) and Skewness (lambda)
      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) {}

   // Multi-variable optimization parameters restricting asymmetry constraints to [-1, 1] range limits
   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;
     }

   // Complex two-sided conditional likelihood calculations based on sign thresholds
   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; // Start with symmetrical assumption
      return out;
     }

   // Map uniform arrays across functional inverse maps to run simulated sets
   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], lam = m_parameters[1];
      double a = _const_a(parameters), b = _const_b(parameters);

      vector icdf = vector::Ones(pits.Size()) * -999.99;
      int ecode = 0;
      for(ulong i = 0; i < pits.Size(); ++i)
        {
         // Choose calculation fork based on regional probability split bounds
         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], lam = m_parameters[1];
      double a = _const_a(parameters), 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;
     }

   // Blends left-tail and right-tail asymmetric components over combinatoric tracking fields
   virtual double    moment(int n, vector& parameters)
     {
      set_params(parameters);
      if(n < 0) return double("nan");
      _check_constraints(parameters);

      double eta = m_parameters[0], lam = m_parameters[1];
      double a = _const_a(parameters), 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], lam = m_parameters[1];
      double a = _const_a(parameters), 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:
   //--- Calculates the standardized partial moment of the Generalized Normal distribution
   double            _ord_gennorm_partial_moment(int n, int z, double beta)
     {
      if(n<0)
         return double("nan");

      // Math constants and scaling factor using Gamma Function
      double w = 0.5*beta/CGammaFunc::GammaFunc(1./beta);
      
      // Compute left-side integration boundary using the regularized complementary upper incomplete gamma function (igamc)
      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));
      
      // Compute right-side integration boundary
      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;
     }

   //--- Computes statistical properties (mean, variance, skewness, excess kurtosis) using shape parameter beta
   void             stats(double beta, double &a, double &b, double &c, double &d)
     {
      // Use log-gamma variants to avoid numeric overflow
      double c1 = log(CGammaFunc::GammaFunc(1./beta));
      double c3 = log(CGammaFunc::GammaFunc(3./beta));
      double c5 = log(CGammaFunc::GammaFunc(5./beta));

      a = 0.;               // Mean (Always zero due to symmetry)
      b = exp(c3-c1);       // Variance
      c = 0.;               // Skewness (Symmetric distribution)
      d = exp(c5+c1-2.*c3)-3.; // Excess Kurtosis
      return;
     }

   //--- Converts central moments to raw moments given location, scale, and statistical traits
   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;
     }

   //--- Core implementation for raw nth-order central moment calculations
   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];
         
      // Apply binomial expansions if location shifting is required (loc != 0.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;
     }

   //--- Percent Point Function (Inverse CDF) using Cephes inverse incomplete gamma logic
   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);
     }

   //--- Wrapper for scalar PPF calls with out-of-bounds safety checks
   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;
     }

   //--- Vectorized version of the Percent Point Function
   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;
     }

   //--- Private override initialization sequence called via base class interface
   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:
   //+------------------------------------------------------------------+
   //| Constructors & Destructor                                        |
   //+------------------------------------------------------------------+
   CGeneralizedError(void)
     {
      m_dist_model = DIST_GEN_ERROR;
     }

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

   ~CGeneralizedError(void)
     {
     }

   //+------------------------------------------------------------------+
   //| Parameter Constraints & Optimization Bounds                       |
   //+------------------------------------------------------------------+
   virtual Constraints   constraints(void)  override
     {
      Constraints out;
      matrix x = {{1}, {-1}};   // Constraint definition matrix
      vector y = {1.01, -500.0}; // Shape bounds (typically nu > 1.0)
      out._one = x;
      out._two = y;
      return out;
     }

   virtual matrix    bounds(vector &resids) override
     {
      matrix out = {{1.01},{500.}}; // Box constraints for optimizer [Min, Max]
      return out;
     }

   //+------------------------------------------------------------------+
   //| Log-Likelihood: Used during ARCH parameter estimation            |
   //+------------------------------------------------------------------+
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individual = false) override
     {
      set_params(parameters);
      double nu = parameters[0]; // Shape parameter (Degrees of freedom)
      
      // Scaling constants for GED normalization
      double logc = 0.5 * (-2. / nu * log(2.) + log(CGammaFunc::GammaFunc(1. / nu)) - log(CGammaFunc::GammaFunc(3. / nu)));
      double c = exp(logc);
      
      // Constant component of the log-likelihood function
      double lls = log(nu) - logc - log(CGammaFunc::GammaFunc(1./nu)) - (1.+1./nu)*log(2.);
      
      // Apply localized dynamic variance offsets over standard conditional variance (sigmas)
      vector llsv = lls - (0.5 * log(sigmas));
      llsv = llsv - (0.5*pow(MathAbs(resids/(sqrt(sigmas)*c)),nu));
      
      if(individual)
         return llsv; // Return pointwise vector for diagnostics
      else
        {
         vector out(1);
         out[0] = llsv.Sum(); // Return aggregated objective function value
         return out;
        }
     }

   //--- Optimization structural starting value guess
   virtual vector    startingValues(vector& resids) override
     {
      vector out(1);
      out[0] = 1.5; // Standard starting shape value for financial residual optimization
      return out;
     }

   //+------------------------------------------------------------------+
   //| Simulator: Generates random matrix pathways via Gamma conversion |
   //+------------------------------------------------------------------+
   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;
      
      // Generate standard gamma distributions to transform into GED random paths
      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);

      // Randomize direction signs uniformly to guarantee distribution symmetry
      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; // Rescale back to unit variance base
     }

   //+------------------------------------------------------------------+
   //| Public Distribution Interfaces: PPF, CDF, Moments                |
   //+------------------------------------------------------------------+
   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;

      // Classify outputs and safely compute localized probability values via igamc
      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);

      // Compute partial localized moment translated against current bounds scale
      double mom = pow(scale,n)*_ord_gennorm_partial_moment(n,int(z/scale),nu);
      return mom;
     }
  };
//+------------------------------------------------------------------+
//| RNG                                                              |
//+------------------------------------------------------------------+
class BootstrapRng
  {
protected:
   vector            m_std_resid;     // Vector storing standardized residuals for empirical bootstrapping
   ulong             m_start;         // Initial starting index for empirical sampling
   ulong             m_index;         // Current index position tracked during random generation
   CHighQualityRandStateShell m_rstate;// State wrapper for the High Quality Random Number Generator
   CDistribution* m_distribution;  // Pointer to the chosen parametric distribution model
   bool              m_initialized;   // Flag tracking whether the generator is successfully ready for use

public:
   //+------------------------------------------------------------------+
   //| Default Constructor: Initializes pointers and states to safe defaults|
   //+------------------------------------------------------------------+
   BootstrapRng(void)
     {
      //--- m_rstate = NULL;
      m_distribution = NULL;
      m_initialized = false;
      m_start = m_index = ULONG_MAX;
      m_std_resid = vector::Zeros(0);
     }

   //+------------------------------------------------------------------+
   //| Parametric Constructor: Instantiates an RNG using a specific     |
   //| mathematical distribution model.                                 |
   //+------------------------------------------------------------------+
   BootstrapRng(ENUM_DISTRIBUTION_MODEL type, vector& parameters, int seed)
     {
      m_distribution = NULL;
      
      // Instantiate the specific parametric distribution class based on type
      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;
        }
      
      // Verify successful initialization of the distribution object
      m_initialized = (m_distribution.is_initialized()==true);//((CheckPointer(m_distribution)!=POINTER_INVALID) && (m_distribution.initialize(parameters,seed)));
     }

   //+------------------------------------------------------------------+
   //| Empirical Constructor: Instantiates an RNG for historic          |
   //| bootstrapping using a fixed dataset of standardized residuals.   |
   //+------------------------------------------------------------------+
   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;

      // Seed the High Quality Random Number Generator
      if(seed)
         CHighQualityRand::HQRndSeed(seed,seed+1,m_rstate.GetInnerObj());
      else
         CHighQualityRand::HQRndRandomize(m_rstate.GetInnerObj());
         
      m_initialized = true;
     }

   //+------------------------------------------------------------------+
   //| Destructor: Clean up dynamically allocated distribution memory   |
   //+------------------------------------------------------------------+
   ~BootstrapRng(void)
     {
      if(CheckPointer(m_distribution) == POINTER_DYNAMIC)
         delete m_distribution;
     }

   //+------------------------------------------------------------------+
   //| Check if the instance is ready to generate random data           |
   //+------------------------------------------------------------------+
   bool is_initialized(void)
     {
      return m_initialized;
     }

   //+------------------------------------------------------------------+
   //| Generates a 1D vector of random samples via empirical bootstrap  |
   //+------------------------------------------------------------------+
   vector rng(ulong size)
     {
      // Check if we are using empirical data (residuals vector populated)
      if(m_std_resid.Size())
        {
         // Prevent out-of-bounds tracking if the pointer index exceeds available data
         if(m_index>=m_std_resid.Size())
           {
            Print(__FUNCTION__, " not enough data points ");
            return vector::Zeros(0);
           }
           
         // 1. Generate Uniform random values between [0, 1)
         vector index = vector::Zeros(size);
         for(ulong i = 0; i<size; index[i] = CHighQualityRand::HQRndUniformR(m_rstate.GetInnerObj()), ++i);
         
         // 2. Map uniform samples to empirical indices within current bound constraints
         index = floor(double(m_index)*index);
         m_index+=1; // Increment current tracked horizon point
         
         // 3. Extract the actual residual values from the historic vector matching the generated indices
         vector out = index;
         for(ulong i = 0; i<size; ++i)
            out[i] = m_std_resid[ulong(index[i])];
            
         return out;
        }
      return EMPTY_VECTOR;
     }

   //+------------------------------------------------------------------+
   //| Generates a 2D matrix of random samples using either parametric   |
   //| or empirical methods depending on how the class was instantiated.|
   //+------------------------------------------------------------------+
   matrix rng(ulong size, ulong cols)
     {
      // Route generation to parametric distribution if instantiated
      if(CheckPointer(m_distribution)!=POINTER_INVALID)
         return m_distribution.rng(size,cols);
      else
        {
         // Otherwise, fallback to generating a 1D empirical vector and reshaping it into a matrix
         matrix out(size*cols,1);
         vector col = rng(size*cols);
         out.Col(col,0);
         out.Reshape(size,cols);
         return out;
        }
     }
  };