Article-23677-EGARCH-MQL5-Volatility-Modeling/Arch/Univariate/volatility.mqh

//+------------------------------------------------------------------+
//|                                                   volatility.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "recursions.mqh"
#include "distribution.mqh"
//---
//+------------------------------------------------------------------+
//| Container for variance forecasts: Stores point-in-time predictions |
//| and full stochastic simulation paths for GARCH-family models.      |
//+------------------------------------------------------------------+
struct VarianceForecast
  {
   matrix            forecasts;     // Holds point variance estimates [steps x horizons]
   matrix            forecastpaths[]; // Array of simulation paths for each forecast origin
   matrix            shocks[];        // Simulated residual shocks corresponding to paths

   // Default constructor
                     VarianceForecast(void)
     {
      forecasts = matrix::Zeros(0, 0);
     }

   // Parameterized constructor to initialize from existing matrices
                     VarianceForecast(matrix &_forecasts, matrix &_forecastpaths[], matrix &_shocks[])
     {
      forecasts = _forecasts;

      // Copy simulation paths and associated shocks
      ArrayResize(forecastpaths, _forecastpaths.Size());
      ArrayResize(shocks, _shocks.Size());

      for(uint i = 0 ; i < shocks.Size(); shocks[i] = _shocks[i], ++i);
      for(uint i = 0 ; i < forecastpaths.Size(); forecastpaths[i] = _forecastpaths[i], ++i);
     }

   // Copy constructor
                     VarianceForecast(VarianceForecast &other)
     {
      forecasts = other.forecasts;

      // Deep copy array members
      ArrayResize(forecastpaths, other.forecastpaths.Size());
      ArrayResize(shocks, other.shocks.Size());

      for(uint i = 0 ; i < shocks.Size(); shocks[i] = other.shocks[i], ++i);
      for(uint i = 0 ; i < forecastpaths.Size(); forecastpaths[i] = other.forecastpaths[i], ++i);
     }

   // Assignment operator for deep object copying
   void              operator=(VarianceForecast &other)
     {

      forecasts = other.forecasts;

      // Resize and copy buffers from source object
      ArrayResize(forecastpaths, other.forecastpaths.Size());
      ArrayResize(shocks, other.shocks.Size());

      for(uint i = 0 ; i < shocks.Size(); shocks[i] = other.shocks[i], ++i);
      for(uint i = 0 ; i < forecastpaths.Size(); forecastpaths[i] = other.forecastpaths[i], ++i);
     }
  };
//+------------------------------------------------------------------+
//| Base class for all volatility processes (GARCH, EGARCH, etc.)    |
//| Enforces interface consistency for estimation and forecasting.   |
//+------------------------------------------------------------------+
class CVolatilityProcess
  {
protected:
   bool               m_updateable;    // Can model parameters be updated?
   bool               m_initialized;   // Is the model ready for computation?
   ulong              m_num_params;    // Number of model parameters to optimize
   bool               m_closedform;    // Does model have an analytical forecast?
   ulong              m_bootstrap_obs; // Minimum observations required for bootstrap
   long               m_start, m_stop; // Data range indices
   CNormal            m_normal;        // Gaussian distribution reference
   VolatilityUpdater *m_volupdater;    // Polymorphic helper for variance recursion
   string             m_name;          // Model identification string
   int                m_seed;          // Seed for stochastic simulations
   ENUM_VOLATILITY_MODEL m_model;      // Model type identifier

   //--- Compute the variance for a single observation (Overridden by child)
   virtual double     _update(ulong index, vector& parameters, vector &resids, vector &sigma2, vector &backcast, vector &var_bounds)
     {
      return EMPTY_VALUE;
     }

   //--- Verify if a forecasting method (e.g., Analytic vs Simulation) is valid
   virtual bool       _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon)
     {
      return false;
     }

   //--- One-step ahead forecast: Projects variance forward by one period
   virtual bool       _onestepforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, vector&out_sigma, matrix&out_forecasts)
     {
      ulong t = resids.Size();
      vector _resids = resids;
      _resids.Resize(t + 1); // Expand to include future residual
      _resids[t] = 0.;       // Assume zero innovation for point forecast

      matrix _vb = varbounds;
      _vb.Resize(varbounds.Rows() + 1, varbounds.Cols());
      vector temp = {0., double("inf")};
      _vb.Row(temp, varbounds.Rows()); // Extend bounds

      vector sigma2 = vector::Zeros(t + 1);
      computeVariance(parameters, _resids, sigma2, backcast, _vb);

      matrix forecasts = matrix::Zeros(t - _start, horizon);
      temp = np::sliceVector(sigma2, _start + 1);
      forecasts.Col(temp, 0);
      out_sigma = np::sliceVector(sigma2, 0, -1);
      out_forecasts = forecasts;
      return true;
     }

   //--- Computes Gaussian Log-Likelihood, useful for initial parameter estimation
   virtual double     _gaussianloglikelihood(vector &parameters, vector& resids, vector&backcast, matrix& varbounds)
     {
      vector sigma2 = vector::Zeros(resids.Size());
      computeVariance(parameters, resids, sigma2, backcast, varbounds);
      vector empty = vector::Zeros(0);
      vector v = m_normal.loglikelihood(empty, resids, sigma2);
      return v[0];
     }

   //--- Analytic multi-step forecasts: Mathematical projection of variance
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon)
     {
      return VarianceForecast(); // Default: must be implemented by subclasses
     }

   //--- Simulation-based forecasts: Generates Monte Carlo paths of future variance
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      return VarianceForecast();
     }

   //--- Bootstrap forecast: Uses historical residuals to simulate future volatility
   virtual VarianceForecast _bootstrapforecast(vector& parameters, vector& resids, vector& backcast, matrix &varbounds, ulong _start, ulong horizon, ulong simulations, int seed = 0)
     {
      VarianceForecast out;
      vector sigma2 = vector::Zeros(resids.Size());
      computeVariance(parameters, resids, sigma2, backcast, varbounds);

      // Calculate standardized residuals to resample from
      vector std_resid = resids / sqrt(sigma2);

      // Safety check for data sufficiency
      if(_start < m_bootstrap_obs)
        {
         Print(__FUNCTION__, " start must include more than ", m_bootstrap_obs, " observations");
         return out;
        }

      BootstrapRng Rng(std_resid, _start, seed);
      out = _simulationforecast(parameters, resids, backcast, varbounds, _start, horizon, simulations, Rng);
      return out;
     }

   //--- Internal initialization routine to set model hyper-parameters
   virtual bool       _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100)
     {
      m_updateable = updateable;
      m_num_params = nparams;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;
      m_model = volmodel;
      return m_normal.initialize(vector::Zeros(0), seed);
     }

   //+------------------------------------------------------------------+
   //| Public Interface for CVolatilityProcess                          |
   //+------------------------------------------------------------------+
public:
   // Constructor: Initializes defaults for model configuration
                     CVolatilityProcess(void): m_updateable(true),
                     m_num_params(0),
                     m_closedform(false),
                     m_bootstrap_obs(100),
                     m_start(0),
                     m_stop(-1),
                     m_initialized(false),
                     m_name(NULL),
                     m_volupdater(NULL),
                     m_model(WRONG_VALUE)
     {}

   // Copy constructor: Performs shallow copy of configuration properties
                     CVolatilityProcess(CVolatilityProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closedForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
     }

   // Destructor: Safely cleans up the dynamically allocated volupdater
                    ~CVolatilityProcess(void)
     {
      if(CheckPointer(m_volupdater) == POINTER_DYNAMIC)
         delete m_volupdater;
      m_volupdater = NULL;
     }

   // Assignment operator: Updates object state from another process
   void              operator=(CVolatilityProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closedForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
     }

   // --- Accessors ---
   ENUM_VOLATILITY_MODEL volatilityprocess(void) { return m_model; }
   ulong             bootstraps(void) { return m_bootstrap_obs; }

   // Setters with validation
   void              bootstraps(ulong bts) { if(bts > 10) m_bootstrap_obs = bts; }
   string            name(void) { return m_name; }
   bool              is_initialized(void) { return m_initialized; }

   // Index management
   virtual long      start(void) { return m_start; }
   virtual long      stop(void) { return m_stop; }
   virtual void      start(long _start) { m_start = _start; }
   virtual void      stop(long _stop) { m_stop = _stop; }

   // Model attributes
   virtual ulong     numParams(void) { return m_num_params; }
   virtual bool      upDateable(void) { return m_updateable; }
   virtual bool      closedForm(void) { return m_closedform; }
   virtual VolatilityUpdater* volUpdater(void) { return m_volupdater; }

   // --- Core Methods ---

   // Compute variance for a specific index (Default interface)
   virtual double    upDate(ulong _index, vector& parameters, vector& resids, vector& sigma2, vector& backcast, vector& varbounds)
     { return EMPTY_VALUE; }

   // Construct loose bounds to ensure variance remains positive/stable during optimization
   virtual matrix    varianceBounds(vector& resids, double power = 2.)
     {
      ulong nobs = resids.Size();
      ulong tau = MathMin(75, nobs);
      vector w = vector::Zeros(tau);
      for(ulong i = 0; i < tau; w[i] = pow(0.94, double(i)), ++i);
      w = w / w.Sum();
      double initial_value = w.Dot(pow(np::sliceVector(resids, 0, long(tau)), 2.));

      vector varbound = vector::Zeros(nobs);
      ewma_recursion(0.94, resids, varbound, resids.Size(), initial_value);

      matrix varbounds = matrix::Zeros(varbound.Size(), 2);
      varbounds.Col(varbound / 1.e6, 0); // Define lower bound
      varbounds.Col(varbound * 1.e6, 1); // Define upper bound

      double var = resids.Var();
      double min_upper_bound = 1. + (pow(resids, 2.)).Max();
      double lower_bound = var / 1.e8;
      double upper_bound = 1.e7 * (1. + (pow(resids, 2.)).Max());

      vector col0 = varbounds.Col(0);
      vector col1 = varbounds.Col(1);

      // Clip values to prevent extreme numerical overflow/underflow
      if(!col0.Clip(lower_bound, DBL_MAX) || !col1.Clip(min_upper_bound, upper_bound) || !varbounds.Col(col0, 0) || !varbounds.Col(col1, 1))
        {
         Print(__FUNCTION__, " error ", __LINE__, " ", GetLastError());
         return matrix::Zeros(0, 0);
        }

      if(power != 2.)
         varbounds = pow(varbounds, power / 2.);
      return varbounds;
     }

   // --- Initialization and Estimation Helpers ---
   virtual vector    startingValues(vector& resids) { return vector::Zeros(0); }

   // Backcast: Estimate initial variance before the start of the observed data
   virtual vector    backCast(vector& resids)
     {
      ulong tau = MathMin(75, resids.Size());
      vector w = vector::Zeros(tau);
      for(ulong i = 0; i < tau; w[i] = pow(0.94, double(i)), ++i);
      w = w / w.Sum();
      vector out(1);
      out[0] = (pow(np::sliceVector(resids, 0, long(tau)), 2.0) * w).Sum();
      return out;
     }

   virtual vector    backCastTransform(vector& backcast)
     {
      if(backcast.Min() < 0)
        {
         Print(__FUNCTION__, " User backcast value must be strictly positive ");
         return vector::Zeros(0);
        }
      return backcast;
     }

   virtual matrix    bounds(vector& resids) { return matrix::Zeros(0, 0); }
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) { return vector::Zeros(0); }
   virtual Constraints constraints(void) { return Constraints(); }

   // --- Polymorphic Forecast Engine ---
   // Routes the forecast request to analytic, simulation, or bootstrap logic
   VarianceForecast  forecast(vector &parameters, vector& resids, vector& _backcast, matrix& var_bounds, BootstrapRng &rng, int seed = 0, ulong _start = 0, ulong horizon = 1, ENUM_FORECAST_METHOD ForecastingMethod = FORECAST_ANALYTIC, ulong simulations = 1000)
     {
      VarianceForecast out;
      if(!horizon)
        {
         Print(__FUNCTION__, " horizon must be >= 1");
         return out;
        }
      if(!_check_forecasting_method(ForecastingMethod, horizon))
        { Print(__FUNCTION__, " Method not supported for this model"); return out; }

      if(!_start)
         _start = resids.Size() - 1;

      switch(ForecastingMethod)
        {
         case FORECAST_ANALYTIC:
            out = _analyticforecast(parameters, resids, _backcast, var_bounds, _start, horizon);
            break;
         case FORECAST_SIMULATION:
            if(rng.is_initialized())
               out = _simulationforecast(parameters, resids, _backcast, var_bounds, _start, horizon, simulations, rng);
            else
               Print(__FUNCTION__, " ERROR: BOOTSTRAP OBJECT IS NULL ");
            break;
         case FORECAST_BOOTSTRAP:
            if(_start < 10)
               Print(__FUNCTION__, " Bootstrap requires > 10 obs");
            else
               if(double(horizon / _start) > 0.2)
                  Print(__FUNCTION__, " Ratio of horizon-to-start < 20% required. ");
               else
                  out = _bootstrapforecast(parameters, resids, _backcast, var_bounds, _start, horizon, simulations, seed);
            break;
        }
      return out;
     }

   virtual matrix    simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) { return matrix::Zeros(0, 0); }
   virtual string    parameterNames(void) { return NULL; }
   virtual vector    garch_spec(void) { return vector::Zeros(4); }
  };
//+------------------------------------------------------------------+
//|  Constant volatility process                                     |
//+------------------------------------------------------------------+
class CConstantVariance: public CVolatilityProcess
  {
protected:
   //--- Verify the requested forecasting method as valid for the specification
   virtual bool      _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override
     {
      return true;
     }
   //--- Analytic multi-step volatility forecasts from the model
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override
     {
      VarianceForecast out;
      long t = (long)resids.Size();
      matrix forecasts = matrix::Zeros(t-_start,horizon);
      forecasts.Fill(parameters[0]);
      out.forecasts = forecasts;
      return out;
     }
   //--- Simulation-based volatility forecasts from the model
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      VarianceForecast out;
      long t = (long)resids.Size();
      matrix forecasts = matrix::Zeros(t-_start,horizon);
      ArrayResize(out.forecastpaths,int(t-_start));
      ArrayResize(out.shocks,int(t-_start));
      forecasts.Fill(parameters[0]);
      for(long i = 0; i<long(t-_start); ++i)
        {
         out.shocks[i] = sqrt(parameters[0])*rng.rng(simulations,horizon);
         out.forecastpaths[i].Fill(parameters[0]);
        }
      return out;
     }

public:
                     CConstantVariance(void)
     {
      m_initialized = CVolatilityProcess::_initialize(VOL_CONST,true,true,1,"Constant Variance");
     }
                     CConstantVariance(int seed=0,ulong nbootstraps=100)
     {
      m_initialized = CVolatilityProcess::_initialize(VOL_CONST,true,true,1,"Constant Variance",seed,0,0,0,0.0,0,-1,nbootstraps);
     }
                    ~CConstantVariance(void)
     {
     }
   //--- Compute the variance for the ARCH model
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) override
     {
      _sigma2.Fill(parameters[0]);
      return _sigma2;
     }
   //--- Returns starting values for the ARCH model
   virtual vector    startingValues(vector& resids) override
     {
      vector out(1);
      out[0] = resids.Var();
      return out;
     }
   //--- Construct parameter constraints arrays for parameter estimation
   virtual Constraints constraints(void) override
     {
      Constraints out;
      out._one = matrix::Ones(1,1);
      out._two = vector::Zeros(1);
      return out;
     }
   //--- Transformation to apply to user-provided backcast values
   virtual vector    backCastTransform(vector& backcast) override
     {
      backcast=backCastTransform(backcast);
      return backcast;
     }
   //--- Construct values for backcasting to start the recursion
   virtual vector    backCast(vector& resids) override
     {
      vector out(1);
      out[0] = resids.Var();
      return out;
     }
   //--- Returns bounds for parameters
   virtual matrix    bounds(vector& resids) override
     {
      double v = resids.Var();
      matrix out = matrix::Zeros(1,2);
      out[0,0] = v/100000.;
      out[0,1] = 10.0 * (v + pow(resids.Mean(), 2.0));
      return out.Transpose();
     }
   //--- Returns bounds for parameters
   virtual matrix       simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      vector errv = rng.rng(_nobs+burn);
      if(!errv.Size())
        {
         matrix container = rng.rng(_nobs+burn,1);
         errv = container.Col(0);
        }
      vector sigma2 = vector::Ones(_nobs+burn) * parameters[0];
      vector data =  sqrt(sigma2) * errv;
      sigma2 = np::sliceVector(sigma2,long(burn));
      data = np::sliceVector(data,long(burn));
      matrix out(data.Size(),2);
      out.Col(data,0);
      out.Col(sigma2,1);
      return out;
     }
   //--- Names of model parameters
   virtual string        parameterNames(void) override
     {
      return "sigma2";
     }
  };
//+------------------------------------------------------------------+
//|  GARCH and related model estimation                               |
//|  Base class for generalized autoregressive conditional           |
//|  heteroskedasticity processes                                     |
//+------------------------------------------------------------------+
class CGarchProcess: public CVolatilityProcess
  {
protected:
   ulong               m_p;    // Number of GARCH lags (variance lags)
   ulong               m_o;    // Number of asymmetry terms (leverage effect)
   ulong               m_q;    // Number of ARCH lags (residual squared lags)
   double              m_power; // Exponent applied to residuals (usually 2 for GARCH)

   //--- Packs three input vectors into all possible combinations (Cartesian product)
   //--- Used for generating grid search candidates or parameter combination matrices
   matrix              _product(vector& one, vector& two, vector& three)
     {
      ulong size = one.Size();
      ulong rows = 3;
      // Creates a matrix of size (size^3) x 3
      matrix out(ulong(pow(size, rows)), rows);
      ulong shift = 0;

      // Triple nested loop to populate every combination
      for(ulong i = 0; i < size; ++i)
        {
         for(ulong j = 0; j < size; ++j)
           {
            for(ulong k = 0; k < size; ++k)
              {
               out[shift, 0] = one[i];
               out[shift, 1] = two[j];
               out[shift, 2] = three[k];
               ++shift;
              }
           }
        }
      return out;
     }

   //--- Verify if the requested forecasting method is supported by this GARCH variant
   virtual bool        _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override
     {
      // Single-step forecasts are always supported
      if(horizon == 1)
         return true;

      // Analytic forecasts require a power of 2 (standard GARCH) for multi-step recursion.
      // Models like EGARCH (power=1) or GJR-GARCH may require simulation for multi-step.
      if(method == FORECAST_ANALYTIC && m_power != 2.)
        {
         Print(__FUNCTION__, " Analytic forecasts not available for horizon > 1 when power != 2");
         return false;
        }
      return true;
     }
   //--- Analytic multi-step volatility forecasts from the model
   //--- Computes recursive variance projections using the GARCH/GJR-GARCH structure
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override
     {
      VarianceForecast out;
      ulong t = resids.Size();
      vector sigma2;
      matrix forecasts;

      // Compute the immediate next-step forecast (h=1)
      _onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sigma2, forecasts);
      if(horizon == 1)
        {
         out.forecasts = forecasts;
         return out;
        }

      // Extract GARCH model components from parameter vector
      double omega = parameters[0];
      vector alpha = np::sliceVector(parameters, 1, long(m_p + 1));        // ARCH terms
      vector gamma = np::sliceVector(parameters, long(m_p + 1), long(m_p + m_o + 1)); // Leverage/Asymmetry terms
      vector beta = np::sliceVector(parameters, long(m_p + m_o + 1));      // GARCH (persistence) terms

      // Determine maximum lag depth required for recursion
      ulong m = MathMax(MathMax(m_p, m_o), m_q);
      vector _resids, _asym_resids, _sigma2, temp;
      _resids = _asym_resids = _sigma2 = vector::Zeros(m + horizon);

      // Iterate through historical data to set up the recursive projection start point
      for(ulong i = ulong(_start); i < t; ++i)
        {
         // If historical data exists for the full lag depth 'm'
         if(i - m + 1 >= 0)
           {
            temp = np::sliceVector(resids, long(i - m + 1), long(i + 1));
            np::vectorCopy(_resids, temp, 0, long(m));
            temp = np::sliceVector(_resids, 0, long(m));
            temp *= np::whereVectorIsLt(temp, 0.0); // Identify negative shocks for asymmetry
            np::vectorCopy(_asym_resids, temp, 0, long(m));
            temp = np::sliceVector(sigma2, long(i - m + 1), long(i + 1));
            np::vectorCopy(_sigma2, temp, 0, long(m));
           }
         // If we are at the start of the series (padding with backcast values)
         else
           {
            np::fillVector(_resids, sqrt(backcast[0]), 0, long(m - i - 1));
            temp = np::sliceVector(resids, 0, long(i + 1));
            np::vectorCopy(_resids, temp, long(m - i - 1), long(m));
            _asym_resids = _resids * (np::whereVectorIsLt(_resids, 0.0));
            np::fillVector(_asym_resids, sqrt(0.5 * backcast[0]), 0, long(m - i - 1));
            np::fillVector(_sigma2, backcast[0], 0, long(m));
            temp = np::sliceVector(sigma2, 0, long(i + 1));
            np::vectorCopy(_sigma2, temp, long(m - i - 10), long(m));
           }

         // Perform recursive forecasting for the requested horizon
         for(ulong h = 0; h < horizon; ++h)
           {
            ulong floc = i - ulong(_start);
            forecasts[floc, h] = omega;
            ulong sloc = h + m - 1;

            // Add ARCH components: alpha * eps^2_{t-j}
            for(ulong j = 0; j < m_p; ++j)
               forecasts[floc, h] += alpha[j] * pow(_resids[sloc - j], 2.);
            // Add GJR/Asymmetry components: gamma * I_{t-j} * eps^2_{t-j}
            for(ulong j = 0; j < m_o; ++j)
               forecasts[floc, h] += (gamma[j] * pow(_asym_resids[sloc - j], 2.));
            // Add GARCH persistence components: beta * sigma^2_{t-j}
            for(ulong j = 0; j < m_q; ++j)
               forecasts[floc, h] += beta[j] * _sigma2[sloc - j];

            // Update buffers with the new forecast to use for next-step recursion
            _resids[h + m] = sqrt(forecasts[floc, h]);
            _asym_resids[h + m] = sqrt(0.5 * forecasts[floc, h]);
            _sigma2[h + m] = forecasts[floc, h];
           }
        }

      out.forecasts = forecasts;
      return out;
     }
   //--- Simulates paths for forecasting operation using Monte Carlo methods
   //--- Projects future volatility by iteratively calculating shocks and updating variance
   virtual bool      _simulatepaths(ulong m, vector& parameters, ulong horizon, matrix& std_shocks, matrix& scaled_forecast_paths, matrix& scaled_shock, matrix& asym_scaled_shock, vector& mean_fpaths, matrix& fpaths, matrix& shocks)
     {
      // Extract model parameters: Omega (constant), Alpha (ARCH), Gamma (Asymmetry/Leverage), Beta (GARCH)
      double omega = parameters[0];
      vector alpha = ((m_p + 1) > 1) ? np::sliceVector(parameters, 1, long(m_p + 1)) : EMPTY_VECTOR;
      vector gamma = ((m_p + 1) < (m_p + m_o + 1)) ? np::sliceVector(parameters, long(m_p + 1), long(m_p + m_o + 1)) : EMPTY_VECTOR;
      vector beta = (parameters.Size() > (m_p + m_o + 1)) ? np::sliceVector(parameters, long(m_p + m_o + 1)) : EMPTY_VECTOR;

      matrix shock = scaled_forecast_paths;
      vector temp;

      // Iterate through the forecast horizon
      for(ulong h = 0; h < horizon; ++h)
        {
         ulong loc = h + m - 1;
         temp = vector::Zeros(scaled_forecast_paths.Rows());
         temp.Fill(omega);
         scaled_forecast_paths.Col(temp, h + m); // Set baseline to constant omega

         // Add ARCH contributions: alpha * epsilon^2
         for(ulong j = 0; j < m_p; ++j)
           {
            temp = scaled_forecast_paths.Col(h + m);
            temp += alpha[j] * scaled_shock.Col(loc - j);
            scaled_forecast_paths.Col(temp, h + m);
           }
         // Add Asymmetry/GJR-GARCH contributions: gamma * (indicator < 0) * epsilon^2
         for(ulong j = 0; j < m_o; ++j)
           {
            temp = scaled_forecast_paths.Col(h + m);
            temp += gamma[j] * asym_scaled_shock.Col(loc - j);
            scaled_forecast_paths.Col(temp, h + m);
           }
         // Add GARCH persistence contributions: beta * sigma^2
         for(ulong j = 0; j < m_q; ++j)
           {
            temp = scaled_forecast_paths.Col(h + m);
            temp += beta[j] * scaled_forecast_paths.Col(loc - j);
            scaled_forecast_paths.Col(temp, h + m);
           }

         // Generate the innovation (shock) for the current step based on projected volatility
         temp = std_shocks.Col(h) * pow(scaled_forecast_paths.Col(h + m), 1. / m_power);
         shock.Col(temp, h + m);

         // Calculate squared shock and flag asymmetry (negative shocks) for next recursion step
         vector lt_zero = np::whereVectorIsLt(shock.Col(h + m), 0.);
         scaled_shock.Col(pow(fabs(shock.Col(h + m)), m_power), h + m);
         asym_scaled_shock.Col(scaled_shock.Col(h + m) * lt_zero, h + m);
        }

      // Prepare final outputs: trim buffers, convert to variance scale, and calculate path means
      fpaths = np::sliceMatrixCols(scaled_forecast_paths, long(m));
      fpaths = pow(fpaths, 2. / m_power); // Standardize back to variance scale
      mean_fpaths = fpaths.Mean(0);       // Average of all simulation paths
      shocks = np::sliceMatrixCols(shock, long(m));

      return true;
     }
   //--- Simulation-based volatility forecasts from the model
   //--- Performs Monte Carlo simulations to project future variance paths for GARCH/GJR-GARCH models
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      VarianceForecast out;
      vector sigma2;
      matrix forecasts;

      // Compute base variance (one-step) to initialize simulation context
      _onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sigma2, forecasts);

      long t = (long)resids.Size();
      ArrayResize(out.forecastpaths, int(t - _start));
      ArrayResize(out.shocks, int(t - _start));

      ulong m = MathMax(MathMax(m_p, m_o), m_q); // Determine maximum lag depth

      // Pre-allocate matrices for simulation state (paths, squared residuals, and asymmetry/leverage components)
      matrix scaled_forecast_paths = matrix::Zeros(simulations, m + horizon);
      matrix scaled_shock = matrix::Zeros(simulations, m + horizon);
      matrix asym_scaled_shock = matrix::Zeros(simulations, m + horizon);

      matrix temp_m;
      vector temp_v, mask;
      ulong count = 0;

      // Iterate over each historical observation to initialize simulation paths
      for(ulong i = ulong(_start); i < ulong(t); ++i)
        {
         matrix std_shocks = rng.rng(simulations, horizon);

         // Initialization for early periods (using backcast values if historical data is insufficient)
         if(i - m < 0)
           {
            np::fillMatrix(scaled_forecast_paths, pow(backcast[0], m_power / 2.), BEGIN, END, STEP, BEGIN, long(m));
            np::fillMatrix(scaled_shock, pow(backcast[0], m_power / 2.), BEGIN, END, STEP, BEGIN, long(m));
            np::fillMatrix(asym_scaled_shock, 0.5 * pow(backcast[0], m_power / 2.), BEGIN, END, STEP, BEGIN, long(m));

            count = i + 1;
            // Populate initial buffers with historical sigma and residuals
            temp_v = np::sliceVector(sigma2, 0, long(count));
            temp_v = pow(temp_v, m_power / 2.);
            temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_forecast_paths.Rows());
            np::matrixCopy(scaled_forecast_paths, temp_m, BEGIN, END, STEP, long(m - count), long(m));

            temp_v = np::sliceVector(resids, 0, long(count));
            temp_v = pow(fabs(temp_v), m_power);
            temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_shock.Rows());
            np::matrixCopy(scaled_shock, temp_m, BEGIN, END, STEP, long(m - count), long(m));

            temp_v = np::sliceVector(resids, 0, long(count));
            temp_v = pow(fabs(temp_v), m_power) * np::whereVectorIsLt(temp_v, 0.0);
            temp_m = np::repeat_vector_as_rows_cols(temp_v, asym_scaled_shock.Rows());
            np::matrixCopy(asym_scaled_shock, temp_m, BEGIN, END, STEP, long(m - count), long(m));
           }
         // Initialization for standard periods using historical data slices
         else
           {
            temp_v = np::sliceVector(sigma2, long(i - m + 1), long(i + 1));
            temp_v = pow(temp_v, m_power / 2.);
            temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_forecast_paths.Rows());
            np::matrixCopy(scaled_forecast_paths, temp_m, BEGIN, END, STEP, BEGIN, long(m));

            temp_v = np::sliceVector(resids, long(i - m + 1), long(i + 1));
            mask = np::whereVectorIsLt(temp_v, 0.0);
            temp_v = pow(fabs(temp_v), m_power);
            temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_shock.Rows());
            np::matrixCopy(scaled_shock, temp_m, BEGIN, END, STEP, BEGIN, long(m));

            temp_m = np::sliceMatrixCols(scaled_shock, 0, long(m));
            temp_v = np::sliceVector(resids, long(i - m + 1), long(i + 1));
            mask = np::whereVectorIsLt(temp_v, 0.0);
            temp_m = np::multiply(temp_m, mask);
            np::matrixCopy(asym_scaled_shock, temp_m, BEGIN, END, STEP, BEGIN, long(m));
           }

         // Generate paths via internal simulation method
         vector f;
         matrix p, s;
         _simulatepaths(m, parameters, horizon, std_shocks, scaled_forecast_paths, scaled_shock, asym_scaled_shock, f, p, s);

         // Store results for the current path
         ulong loc = i - long(_start);
         out.forecasts.Row(f, loc);
         out.forecastpaths[loc] = p;
         out.shocks[loc] = s;
        }
      return out;
     }
   //--- Initialize class members and configure GARCH model structure
   virtual bool      _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 1, ulong o = 0, ulong q = 1, double power = 2., long start = 0, long stop = -1, ulong bootstrap_obs = 100) override
     {
      // Assign basic model configuration
      m_model = volmodel;
      m_updateable = updateable;
      m_num_params = nparams;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;

      // Synchronize random number generator seed if not already initialized
      if(m_normal.randomState() != uint(m_seed))
         m_normal.initialize(vector::Zeros(0), seed);

      // Configure model lag orders and power transformation parameter
      m_p = p; // ARCH lags
      m_o = o; // Asymmetry/Leverage lags
      m_q = q; // GARCH (persistence) lags
      m_power = power;

      // Validate power parameter (must be positive, e.g., 2.0 for standard GARCH)
      if(m_power <= 0.0)
        {
         Print(__FUNCTION__, " invalid value for power variable ");
         return false;
        }

      // Set default lag values if not explicitly provided, based on the specific model type
      switch(m_model)
        {
         case VOL_ARCH:
         case VOL_AVARCH:
            if(!m_p)
               m_p = 1;
            break;
         case VOL_GARCH:
         case VOL_AVGARCH:
            if(!m_p)
               m_p = 1;
            if(!m_q)
               m_q = 1;
            break;
         case VOL_GJR_GARCH:
         case VOL_TARCH:
            if(!m_p)
               m_p = 1;
            if(!m_q)
               m_q = 1;
            if(!m_o)
               m_o = 1;
            break;
        }

      // Recalculate total parameter count: 1 (omega) + p + o + q
      m_num_params = 1 + m_p + m_o + m_q;

      return (true);
     }
public:
   // Constructor: Default initialization for a standard GARCH(1,1) model
                     CGarchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_GARCH, true, false, 0, "GARCH", 0, 1, 0, 1, 2.0);
     }

   // Constructor: Parameterized initialization for custom lag structures and seeding
                     CGarchProcess(ulong p, ulong q, int seed = 0, ulong nbootstraps = 100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_GARCH, true, false, 0, "GARCH", seed, p, 0, q, 2.0, 0, -1, nbootstraps);
     }

   // Copy Constructor: Creates a new instance by copying properties from an existing GARCH process
                     CGarchProcess(CGarchProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closedForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
      m_q = other.get_q();
      m_p = other.get_p();
      m_o = other.get_o();
      m_power = other.get_power();
     }

   // Destructor
                    ~CGarchProcess(void)
     {
     }

   // Assignment Operator: Copies properties from another GARCH process to this instance
   void              operator=(CGarchProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closedForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
      m_q = other.get_q();
      m_p = other.get_p();
      m_o = other.get_o();
      m_power = other.get_power();
     }

   //--- Getters for model specifications (lags and power)
   virtual ulong     get_p(void) { return m_p; }
   virtual ulong     get_o(void) { return m_o; }
   virtual ulong     get_q(void) { return m_q; }
   virtual double    get_power(void) { return m_power; }

   //--- Construct boundaries for conditional variance estimation
   //--- Overrides base class to use the specific m_power exponent defined for this GARCH instance
   virtual matrix    varianceBounds(vector& resids, double power = 2.) override
     {
      return CVolatilityProcess::varianceBounds(resids, m_power);
     }
   //--- Compute the variance for the GARCH model using recursive updates
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) override
     {
      // Transform residuals based on the model power (e.g., standard GARCH uses power=2)
      vector fresids = pow(MathAbs(resids), m_power);

      // Capture sign of residuals to account for asymmetry (leverage effects) in GJR-GARCH models
      vector sresids(resids.Size());
      for(ulong i = 0; i < sresids.Size(); sresids[i] = np::sign(resids[i]), ++i);

      // Run recursive GARCH filtering to estimate conditional variance series
      _sigma2 = garch_recursion(parameters, fresids, sresids, _sigma2, m_p, m_o, m_q, resids.Size(), _backcast[0], varbounds);

      // Scale back to original variance units (inverse power transformation)
      double invp = 2. / m_power;
      return pow(_sigma2, invp);
     }

   //--- Returns starting values for the GARCH model estimation
   //--- Performs a grid search over potential starting parameter combinations to improve solver convergence
   virtual vector    startingValues(vector& resids) override
     {
      // Define candidate grid values for alphas, gammas, and persistence (betas)
      vector alphas = {0.01, 0.05, 0.1, 0.2};
      vector gammas = alphas;
      vector abg = {0.5, 0.7, 0.9, 0.98};

      // Generate Cartesian product of search candidates
      matrix abgs = _product(alphas, gammas, abg);

      // Scale target variance based on data properties to ensure reasonable starting scale
      double target = (pow(fabs(resids), m_power)).Mean();
      double scale = (pow(resids, 2.)).Mean() / pow(target, 2. / m_power);
      target *= pow(scale, m_power / 2.);

      vector svs[];
      ArrayResize(svs, int(abgs.Rows()));
      matrix vb = varianceBounds(resids);
      vector llfs = vector::Zeros(abgs.Rows());
      vector bc = backCast(resids);

      // Evaluate Log-Likelihood for each combination in the grid
      for(uint i = 0; i < svs.Size(); ++i)
        {
         vector row = abgs.Row(i);
         // Construct starting vector: constant omega + lags
         vector sv = (1. - row[2]) * target * vector::Ones(m_p + m_o + m_q + 1);

         // Distribute coefficients across ARCH (p), Leverage (o), and GARCH (q) components
         if(m_p > 0)
           {
            np::fillVector(sv, row[0] / double(m_p), long(1), long(1 + m_p));
            row[2] -= row[0];
           }
         if(m_o > 0)
           {
            np::fillVector(sv, row[1] / double(m_o), long(1 + m_p), long(1 + m_p + m_o));
            row[2] -= row[1] / 2.; // Adjust persistence for asymmetric impact
           }
         if(m_q > 0)
            np::fillVector(sv, row[2] / double(m_q), long(1 + m_p + m_o), long(1 + m_p + m_o + m_q));

         svs[i] = sv;
         llfs[i] = _gaussianloglikelihood(sv, resids, bc, vb);
        }

      // Select the parameter set that yielded the highest Log-Likelihood
      ulong loc = llfs.ArgMax();
      return svs[loc];
     }
   //--- Construct parameter constraints arrays for parameter estimation
   //--- Enforces positivity constraints and stationarity: sum(alpha + gamma/2 + beta) < 1
   virtual Constraints constraints(void) override
     {
      Constraints out;
      ulong k_arch = m_p + m_o + m_q;
      out._one = matrix::Zeros(k_arch + 2, k_arch + 1);

      // Enforce positivity: omega > 0 and all coefficients > 0
      for(ulong i = 0; i < (k_arch + 1); ++i)
         out._one[i, i] = 1.;

      // Additional constraints for asymmetry terms (if applicable)
      for(ulong i = 0; i < (m_o); ++i)
         if(i < m_p)
            out._one[i + m_p + 1, i + 1] = 1.;

      // Enforce stationarity: sum(alpha_i) + 0.5 * sum(gamma_i) + sum(beta_i) < 1
      // Implemented as linear inequality: - (sum(alpha) + 0.5*sum(gamma) + sum(beta)) >= -1
      for(ulong i = 1; i < (k_arch + 1); ++i)
         out._one[k_arch + 1, i] = -1.;

      // Apply 0.5 multiplier for leverage (asymmetry) terms per GJR-GARCH definition
      for(ulong i = m_p + 1; i < (m_p + m_o + 1); ++i)
         out._one[k_arch + 1, i] = -0.5;

      out._two = vector::Zeros(k_arch + 2);
      out._two[k_arch + 1] = -1.; // Sets the stationarity limit
      return out;
     }

   //--- Transformation to apply to user-provided backcast values
   //--- Converts variance-scale backcast to the model's power-scale (e.g., standard deviation for power=2)
   virtual vector    backCastTransform(vector& backcast) override
     {
      backcast = CVolatilityProcess::backCastTransform(backcast);
      if(backcast.Size())
         return pow(sqrt(backcast), m_power);
      else
         return vector::Zeros(0);
     }

   //--- Construct values for backcasting to start the recursion
   //--- Uses an Exponentially Weighted Moving Average (EWMA) of past residuals for initialization
   virtual vector    backCast(vector& resids) override
     {
      ulong tau = MathMin(75, resids.Size());
      // Create decay weights for the EWMA
      vector out = np::arange(tau);
      for(ulong i = 0; i < out.Size(); out[i] = pow(0.94, out[i]), ++i);
      out = out / out.Sum();

      // Calculate weighted average of squared (or power-transformed) residuals
      vector temp = np::sliceVector(resids, 0, long(tau));
      out = pow(MathAbs(temp), m_power) * out;

      vector bc(1);
      bc[0] = out.Sum();
      return bc;
     }
   //--- Returns parameter bounds for optimization
   //--- Defines the allowable range [min, max] for each GARCH parameter
   virtual matrix    bounds(vector& resids) override
     {
      // Calculate the mean of transformed residuals as a scale factor for variance
      double v = (pow(MathAbs(resids), m_power)).Mean();

      // Create a bounds matrix: rows = parameters (1 omega + p ARCH + o leverage + q GARCH), 2 columns (min, max)
      matrix out = matrix::Zeros(1 + m_p + m_o + m_q, 2);

      // Bounds for constant term (omega): near zero to 10x the mean variance
      out[0, 0] = 1e-8 * v;
      out[0, 1] = 10.0 * v;

      ulong from = 1;

      // Bounds for ARCH coefficients (alpha): range [0, 1]
      for(ulong i = 0; i < m_p; ++i)
        {
         out[from, 0] = 0.0;
         out[from++, 1] = 1.;
        }

      // Bounds for Leverage/Asymmetry coefficients (gamma):
      // For GJR-GARCH, if leverage term overlaps with ARCH lag, allow negative values (range [-1, 2])
      // otherwise, standard positive constraint (range [0, 2])
      for(ulong i = 0; i < m_o; ++i)
        {
         if(i < m_p)
           {
            out[from, 0] = -1.0;
            out[from++, 1] = 2.;
           }
         else
           {
            out[from, 0] = 0.0;
            out[from++, 1] = 2.;
           }
        }

      // Bounds for GARCH persistence coefficients (beta): range [0, 1]
      for(ulong i = 0; i < m_q; ++i)
        {
         out[from, 0] = 0.0;
         out[from++, 1] = 1.;
        }

      // Return transposed matrix to match the expected format [2 x N] for the optimizer
      return out.Transpose();
     }
   //--- Names of model parameters
   virtual string       parameterNames(void) override
     {
      return common_names(m_p,m_o,m_q);
     }
   //--- Simulate data from the GARCH/GJR-GARCH model
   virtual matrix    simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      // Generate random innovations (shocks) for the full period including burn-in
      vector errors = rng.rng(_nobs + burn);

      // Determine the initial long-run variance if not provided by the user
      if(initial_value == NULL)
        {
         // Calculate persistence = sum of ARCH + (0.5 * Leverage) + GARCH coefficients
         vector scale = vector::Ones(parameters.Size());
         for(ulong i = m_p + 1; i < (m_p + m_o + 1); scale[i] = 0.5, ++i);
         vector temp = np::sliceVector(parameters, 1) * np::sliceVector(scale, 1);
         double persistence = temp.Sum();

         // Unconditional Variance Formula: Omega / (1 - persistence)
         if((1. - persistence) > 0)
            initial_value = parameters[0] / (1. - persistence);
         else
           {
            Print(__FUNCTION__, " InitialValueWarning: Model parameters may imply non-stationarity");
            initial_value = parameters[0];
           }
        }

      vector sigma2, data, fsigma, fdata;
      sigma2 = data = fsigma = fdata = vector::Zeros(_nobs + burn);
      ulong max_lag = MathMax(MathMax(m_p, m_o), m_q);

      // Initialize pre-sample values to avoid cold-start bias
      np::fillVector(fsigma, initial_value, 0, long(max_lag));
      double dv = pow(initial_value, 2. / m_power);
      np::fillVector(sigma2, dv, 0, long(max_lag));

      for(ulong i = 0; i < max_lag; ++i)
        {
         data[i] = sqrt(sigma2[i]) * errors[i];
         fdata[i] = pow(fabs(data[i]), m_power);
        }

      // Recursive simulation of the volatility process
      for(ulong t = max_lag; t < (_nobs + burn); ++t)
        {
         ulong loc = 0;
         fsigma[t] = parameters[loc]; // Intercept (Omega)
         loc += 1;

         // ARCH terms: Sum of lagged squared residuals
         for(ulong j = 0; j < m_p; ++j)
           {
            fsigma[t] += parameters[loc] * fdata[t - 1 - j];
            loc += 1;
           }
         // Leverage (Asymmetry) terms: Gamma * I(res < 0) * res^2
         for(ulong j = 0; j < m_o; ++j)
           {
            fsigma[t] += parameters[loc] * fdata[t - 1 - j] * ((data[t - 1 - j] < 0.) ? 1. : 0.0);
            loc += 1;
           }
         // GARCH (Persistence) terms: Beta * sigma^2
         for(ulong j = 0; j < m_q; ++j)
           {
            fsigma[t] += parameters[loc] * fdata[t - 1 - j];
            loc += 1;
           }

         // Transform internal power scale back to variance
         sigma2[t] = pow(fsigma[t], 2. / m_power);
         // Generate observation based on current volatility
         data[t] = errors[t] * sqrt(sigma2[t]);
         fdata[t] = pow(MathAbs(data[t]), m_power);
        }

      // Discard burn-in period to ensure the series has reached its stationary distribution
      data = np::sliceVector(data, long(burn));
      sigma2 = np::sliceVector(sigma2, long(burn));

      // Construct output matrix [data, variance]
      matrix out(data.Size(), 2);
      out.Col(data, 0);
      out.Col(sigma2, 1);

      return out;
     }
   //--- Get Garch specfication
   virtual vector    garch_spec(void) override
     {
      vector specs(4);

      specs[0] = double(m_p);
      specs[1] = double(m_o);
      specs[2] = double(m_q);
      specs[3] = double(m_power);

      return specs;
     }
  };
//+------------------------------------------------------------------+
//| ARCH process                                                     |
//| Implementation of the Autoregressive Conditional Heteroskedasticity |
//| model, inheriting from the generalized GARCH framework.          |
//+------------------------------------------------------------------+
class CArchProcess: public CGarchProcess
  {
public:
   // Constructor: Default initialization for a standard ARCH(1) model
                     CArchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_ARCH, true, false, 0, "ARCH");
     }

   // Constructor: Custom initialization for ARCH(p) model with seed and bootstrap settings
                     CArchProcess(ulong p, int seed = 0, ulong nbootstraps = 100)
     {
      // ARCH models have no asymmetry (o=0) or GARCH (q=0) terms
      m_initialized = CGarchProcess::_initialize(VOL_ARCH, true, false, 0, "ARCH", seed, p, 0, 0, 2.0, 0, -1, nbootstraps);
     }

   // Destructor
                    ~CArchProcess(void)
     {
     }

   //--- Returns starting values for the ARCH model optimization
   //--- Performs a targeted grid search to find the initial alpha parameters that maximize the log-likelihood
   virtual vector    startingValues(vector& resids) override
     {
      // Create a grid of candidate ARCH persistence coefficients (alphas) to test
      vector alphas = np::arange(ulong(17), double(0.1), double(0.05));
      vector svs[17];
      vector bc = backCast(resids);
      vector llfs = alphas;
      matrix vb = varianceBounds(resids);
      vector sv, temp;

      // Iterate through candidates to find the best initial guess
      for(uint i = 0; i < svs.Size(); ++i)
        {
         // sv[0] is omega (constant), sv[1...p] are the ARCH coefficients
         // Initial omega is set based on the sample variance and assumed persistence
         sv = (1.0 - alphas[i]) * resids.Var() * vector::Ones(m_p + 1);

         // Evenly distribute the persistence coefficient (alpha) across all ARCH lags
         np::fillVector(sv, alphas[i] / double(m_p), 1);

         svs[i] = sv;
         // Evaluate log-likelihood for this parameter combination
         llfs[i] = _gaussianloglikelihood(sv, resids, bc, vb);
        }

      // Select the index that produced the highest log-likelihood
      ulong loc = llfs.ArgMax();
      return svs[loc];
     }
  };
//+------------------------------------------------------------------+
//| GJR-GARCH process                                                |
//| Models asymmetric volatility (leverage effect) using GARCH-style |
//| recursion with quadratic variance (power=2.0).                   |
//+------------------------------------------------------------------+
class CGjrGarchProcess: public CGarchProcess
  {
public:
   // Default GJR-GARCH(1,1,1) constructor
                     CGjrGarchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_GJR_GARCH, true, false, 0, "GJR-GARCH");
     }
   // Custom GJR-GARCH(p,o,q) constructor
                     CGjrGarchProcess(ulong p, ulong o, ulong q, int seed = 0, ulong nboots = 100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_GJR_GARCH, true, false, 0, "GJR-GARCH", seed, p, o, q, 2.0, 0, -1, nboots);
     }
  };

//+------------------------------------------------------------------+
//| AVARCH process                                                   |
//| Absolute Value ARCH process (power=1.0), modeling conditional   |
//| standard deviation instead of variance.                          |
//+------------------------------------------------------------------+
class CAvarchProcess: public CGarchProcess
  {
public:
   // Default AVARCH(1) constructor
                     CAvarchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_AVARCH, true, false, 0, "AVARCH", 0, 1, 0, 0, 1.0);
     }
   // Custom AVARCH(p) constructor
                     CAvarchProcess(ulong p, int seed = 0, ulong nboots = 100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_AVARCH, true, false, 0, "AVARCH", seed, p, 0, 0, 1.0, 0, -1, nboots);
     }
  };

//+------------------------------------------------------------------+
//| AVGARCH process                                                  |
//| Absolute Value GARCH process (power=1.0), extending AVARCH with   |
//| autoregressive variance components.                              |
//+------------------------------------------------------------------+
class CAvgarchProcess: public CGarchProcess
  {
public:
   // Default AVGARCH(1,1) constructor
                     CAvgarchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_AVGARCH, true, false, 0, "AVGARCH", 0, 1, 0, 1, 1.0);
     }
   // Custom AVGARCH(p,q) constructor
                     CAvgarchProcess(ulong p, ulong q, int seed = 0, ulong nboots = 100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_AVGARCH, true, false, 0, "AVGARCH", seed, p, 0, q, 1.0, 0, -1, nboots);
     }
  };

//+------------------------------------------------------------------+
//| TARCH process                                                    |
//| Threshold ARCH process (power=1.0), allowing for asymmetric      |
//| impacts on conditional standard deviation.                       |
//+------------------------------------------------------------------+
class CTarchProcess: public CGarchProcess
  {
public:
   // Default TARCH(1,1,1) constructor
                     CTarchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_TARCH, true, false, 0, "TARCH", 0, 1, 1, 1, 1.0);
     }
   // Custom TARCH(p,o,q) constructor
                     CTarchProcess(ulong p, ulong o, ulong q, int seed = 0, ulong nboots = 100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_TARCH, true, false, 0, "TARCH", seed, p, o, q, 1.0, 0, -1, nboots);
     }
  };
//+------------------------------------------------------------------+
//| HARCH process                                                    |
//| Heterogeneous Autoregressive Conditional Heteroskedasticity      |
//| Models volatility as a sum of lags over different frequencies    |
//+------------------------------------------------------------------+
class CHarchProcess: public CVolatilityProcess
  {
protected:
   //--- Custom structure to hold Common Forecast Components (CFC)
   //--- This acts as a data container for efficient recursive volatility projection
   struct CFC
     {
      double               Const;   // Omega (long-term variance component)
      vector               Arch;    // Flattened ARCH coefficients vector
      matrix               Resids2; // Matrix of squared residuals (historical + forecast buffer)

      // Default constructor initializing empty structures
                     CFC(void)
        {
         Const = 0;
         Arch = vector::Zeros(0);
         Resids2 = matrix::Zeros(0, 0);
        }

      // Parameterized constructor for pre-populated components
                     CFC(double c, vector& a, matrix& r)
        {
         Const = c;
         Arch = a;
         Resids2 = r;
        }

      // Copy constructor to handle structure duplication
                     CFC(CFC& other)
        {
         Const = other.Const;
         Arch = other.Arch;
         Resids2 = other.Resids2;
        }

      // Assignment operator to facilitate component passing
      void           operator=(CFC& other)
        {
         Const = other.Const;
         Arch = other.Arch;
         Resids2 = other.Resids2;
        }
     };

   ulong               m_num_lags; // Total number of frequency-based lag groups defined
   vector              m_lags;     // Vector containing the specific lag lengths for each group

   //--- Model transformation: Converts HARCH parameters (coefficients per frequency group)
   //--- into a standard ARCH(p) representation (one coefficient per individual lag)
   vector              _harch_to_arch(vector& parameters)
     {
      // ARCH representation requires one parameter per lag up to the maximum lag defined
      vector arch_params = vector::Zeros(1 + int(m_lags.Max()));
      arch_params[0] = parameters[0]; // Copy the constant (omega)

      double lag = 0;
      double param = 0;

      // Distribute the HARCH parameter weight evenly across all lags within that group
      for(ulong i = 1; i < arch_params.Size(); ++i)
        {
         param = parameters[i];
         lag = m_lags[i - 1];
         // Uniform distribution of the impact: alpha_i = weight_group / length_group
         for(ulong j = 1; j < ulong(lag + 1.); arch_params[j] += param / lag, ++j);
        }
      return arch_params;
     }

   //--- Verify the requested forecasting method as valid for the specification
   //--- HARCH models support both analytic and simulation methods natively
   virtual bool        _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override
     {
      return true;
     }

   //--- Intermediate forecast components: Prepares the history matrix for forecasting
   //--- Aligns squared residuals and backcast values to compute future conditional variance
   CFC                 _common_forecast_components(vector& parameters, vector& resids, vector& backcast, ulong horizon)
     {
      CFC out;
      // Transform HARCH parameters to an ARCH format for vectorization
      vector ap = _harch_to_arch(parameters);
      long t = long(resids.Size());
      long m = long(m_lags.Max());

      // Initialize the matrix with dimensions: Time x (MaxLag + Horizon)
      out.Resids2 = matrix::Zeros(t, m + horizon);

      // Seed the forecast horizon part of the matrix with initial backcast values
      for(long i = m; i < (long)out.Resids2.Rows(); ++i)
         for(long j = m; j < long(out.Resids2.Cols()); ++j)
            out.Resids2[i, j] = backcast[0];

      // Populate historical squared residuals into the rolling buffer matrix
      vector sqresids = pow(resids, 2.);
      for(long i = 0; i < m; ++i)
         for(long j = m - i - 1, k = 0; j < long(out.Resids2.Rows()) && k < long(t - (m - i - 1)); ++k, ++j)
            out.Resids2[j, i] = sqresids[k];

      // Store the constant and the ARCH terms in the CFC structure for use in _forecast
      out.Const = ap[0];
      out.Arch = np::sliceVector(ap, 1);
      return out;
     }
   //--- Analytic multi-step volatility forecasts from the HARCH model
   //--- Computes expected future variance using the linear projection of squared residuals
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override
     {
      // Retrieve the pre-calculated model constants and historical residual matrix
      CFC cfc = _common_forecast_components(parameters, resids, backcast, horizon);
      long m = long(m_lags.Max());

      // Trim history to the requested start point
      cfc.Resids2 = np::sliceMatrixRows(cfc.Resids2, _start);

      // Reverse the ARCH coefficient vector to align with the rolling residual buffer (t-1, t-2, ...)
      vector arch_rev = cfc.Arch;
      if(!np::reverseVector(arch_rev))
         return VarianceForecast();

      // Iterative projection of future variance components
      for(long i = 0; i < long(horizon); ++i)
        {
         // Extract the relevant lag window (t-m to t-1)
         matrix cc = np::sliceMatrixCols(cfc.Resids2, i, m + i);
         // Forecast: Omega + Sum(Alpha_j * Resids^2_{t-j})
         if(!cfc.Resids2.Col(cfc.Const + cc.MatMul(arch_rev), m + i))
           {
            Print(__FUNCTION__, ": Failed Columns assignment cfc.Resids2 ", GetLastError());
            return VarianceForecast();
           }
        }

      // Extract the projected variance paths (excluding initial lags)
      matrix out = np::sliceMatrixCols(cfc.Resids2, m);
      matrix empty[];
      return VarianceForecast(out, empty, empty);
     }

   //--- Simulation-based volatility forecasts from the HARCH model
   //--- Projects future variance paths using Monte Carlo integration
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) override
     {
      CFC cfc = _common_forecast_components(parameters, resids, backcast, horizon);
      long t = long(resids.Size());
      long m = long(m_lags.Max());

      matrix paths[], shocks[];
      matrix temp_resids2 = matrix::Zeros(simulations, m + horizon);

      ArrayResize(paths, int(t - _start));
      ArrayResize(shocks, int(t - _start));
      for(uint i = 0; i < paths.Size(); paths[i] = shocks[i] = matrix::Zeros(simulations, horizon), ++i);

      // Reverse ARCH coefficients for matrix multiplication alignment
      vector arch_rev = np::sliceVector(cfc.Arch, BEGIN_REVERSE, END_REVERSE, -1);

      // Run simulations for each historical point
      for(long i = _start; i < t; ++i)
        {
         matrix stdshocks = rng.rng(simulations, horizon);
         matrix rd2 = np::sliceMatrixRows(cfc.Resids2, i, i + 1);
         temp_resids2 = rd2;

         long path_loc = i - _start;

         // Perform recursive step-ahead simulation
         for(long j = 0; j < long(horizon); ++j)
           {
            rd2 = np::sliceMatrixCols(temp_resids2, j, m + j);

            // Update variance: Forecast = Omega + Sum(Alpha * PastSqResids)
            if(!paths[path_loc].Col(cfc.Const + rd2.MatMul(arch_rev), j))
              {
               Print(__FUNCTION__, ": Column assignment error paths", GetLastError());
               return VarianceForecast();
              }
            // Generate stochastic innovation based on the forecasted variance
            if(!shocks[path_loc].Col(stdshocks.Col(j) * sqrt(paths[path_loc].Col(j)), j))
              {
               Print(__FUNCTION__, ": Column assignment error shocks", GetLastError());
               return VarianceForecast();
              }
            // Feed back the squared innovation into the residual buffer for next-step recursion
            if(!temp_resids2.Col(pow(shocks[path_loc].Col(j), 2.), m + j))
              {
               Print(__FUNCTION__, ": Column assignment error temp_resids2", GetLastError());
               return VarianceForecast();
              }
           }
        }

      // Aggregate simulation results to find the mean volatility path
      matrix out(paths.Size(), paths[0].Cols());
      for(uint i = 0; i < paths.Size(); ++i)
         if(!out.Row(paths[i].Mean(1), i))
           {
            Print(__FUNCTION__, ": Row assignment error out", GetLastError());
            return VarianceForecast();
           }

      return VarianceForecast(out, paths, shocks);
     }
   //--- Initialize class members and configure HARCH model structure
   //--- HARCH models volatility as a weighted sum of returns at different frequency horizons
   virtual bool      _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100) override
     {
      m_updateable = updateable;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;
      m_model = volmodel;

      // Ensure the lag vector is populated; defaults to a single lag (standard ARCH(1)) if empty
      if(!m_lags.Size())
        {
         m_lags = vector::Ones(1);
         for(ulong i = 0; i < m_lags.Size(); m_lags[i] = 1. + double(i), ++i);
        }

      // Validate that all defined lags represent positive integers
      if(int(m_lags.Max()) < 1)
        {
         Print(__FUNCTION__, ": Invalid input. All lags should resolve to > 0 integers ");
         return false;
        }

      // Set internal counters for lag depth and total parameters (Constant + 1 per lag group)
      m_num_lags = m_lags.Size();
      m_num_params = m_num_lags + 1;

      // Initialize the random number generator for simulations
      return m_normal.initialize(vector::Zeros(0), seed);
     }

public:
   // Constructor: Default initialization for a basic HARCH process
                     CHarchProcess(void) : m_lags(vector::Ones(1))
     {
      m_initialized = _initialize(VOL_HARCH, true, false, 0, "HARCH");
     }

   // Constructor: Initializing with a custom vector of lag frequency groups
                     CHarchProcess(vector& lags)
     {
      // If a single value is provided, treat it as the depth for a uniform lag range
      if(lags.Size() == 1)
        {
         m_lags = vector::Ones(ulong(lags[0]));
         for(ulong i = 0; i < m_lags.Size(); m_lags[i] = lags[0] + double(i), ++i);
        }
      else
         m_lags = lags; // Use the explicitly provided set of lags

      m_initialized = _initialize(VOL_HARCH, true, false, 0, "HARCH");
     }
   //--- Returns bounds for parameters
   //--- Sets search intervals [min, max] for the HARCH parameters
   virtual matrix    bounds(vector& resids) override
     {
      // Bounds matrix: rows = parameters (1 omega + number of lag groups), 2 columns (min, max)
      matrix out = matrix::Zeros(m_lags.Size() + 1, 2);

      // Estimate sample variance to provide a logical upper bound for the constant (omega)
      vector squared = pow(resids, 2.);
      out[0, 1] = 10. * squared.Mean();

      // Set upper bounds for all ARCH coefficients (alpha_i) to 1.0 (enforcing local stability)
      for(ulong i = 1; i < out.Rows(); ++i)
         out[i, 1] = 1.0;

      // Transpose to match [2 x N] structure required by the optimization algorithm
      return out.Transpose();
     }

   //--- Construct parameter constraints arrays for parameter estimation
   //--- Enforces positivity (omega, alpha > 0) and stationarity: sum(alpha_i) < 1
   virtual Constraints constraints(void) override
     {
      // a[m_num_lags+1, m_num_lags+1] matrix for linear inequalities: A * params <= b
      matrix a = matrix::Zeros(m_num_lags + 2, m_num_lags + 1);

      // Enforce positivity: omega >= 0 and all alpha_i >= 0
      // Implemented as -omega <= 0 and -alpha_i <= 0, or identity matrix for lower bounds
      for(ulong i = 0; i < m_num_lags + 1; a[i, i] = 1., ++i);

      // Enforce stationarity: sum(alpha_i) < 1
      // Linear form: - (alpha_1 + alpha_2 + ... + alpha_k) >= -1  =>  (alpha_1 + ... + alpha_k) <= 1
      for(ulong j = 1; j < a.Cols(); ++j)
         a[m_num_lags + 1, j] = -1.;

      vector b = vector::Zeros(m_num_lags + 2);
      b[m_num_lags + 1] = -1.; // Constraint threshold

      Constraints out;
      out._one = a;
      out._two = b;
      return out;
     }

   //--- Compute the conditional variance for the HARCH model
   //--- Filters the residuals through the heterogeneous lag structure to produce a variance series
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& sigma2, vector& backcast, matrix& varbounds) override
     {
      vector lags = m_lags;
      int nobs = int(resids.Size());

      // Perform recursive filtering based on frequency-specific lag groups
      harch_recursion(parameters, resids, sigma2, lags, nobs, backcast[0], varbounds);

      return sigma2;
     }
   //--- Simulate synthetic data from the HARCH model
   //--- Generates a series of returns and their corresponding conditional variances
   virtual matrix    simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      // Generate random innovations (shocks) for the total period (including burn-in)
      vector errs = rng.rng(_nobs + burn);
      if(!errs.Size())
        {
         matrix container = rng.rng(_nobs + burn, 1);
         errs = container.Col(0);
        }

      // Calculate unconditional long-run variance to seed the process if not provided
      if(!initial_value)
        {
         // Variance = Omega / (1 - sum of all ARCH parameters)
         vector vp = np::sliceVector(parameters, 1);
         double psum = vp.Sum();
         if(1.0 - psum > 0)
            initial_value = parameters[0] / (1.0 - psum);
         else
           {
            Print(__FUNCTION__, ": Initial value WARNING: Parameters imply non-stationarity");
            initial_value = parameters[0];
           }
        }

      vector sigma2 = vector::Zeros(_nobs + burn);
      vector data = sigma2;
      long maxlag = long(m_lags.Max());

      // Initialize pre-sample buffer to prevent cold-start effects
      for(long i = 0; i < maxlag; ++i)
        {
         sigma2[i] = initial_value;
         data[i] = sqrt(initial_value);
        }

      // Recursive simulation of the HARCH(p) process
      double param = 0;
      for(long t = maxlag; t < long(_nobs + burn); ++t)
        {
         sigma2[t] = parameters[0]; // Intercept term (Omega)

         // Aggregate volatility contributions from different frequency lag groups
         for(ulong i = 0; i < m_lags.Size(); ++i)
           {
            // Distribute group coefficient across the specified lag window length
            param = parameters[1 + i] / double(m_lags[i]);
            for(ulong j = 0; j < ulong(m_lags[i]); ++j)
               sigma2[t] += param * pow(data[t - 1 - j], 2.);
           }

         // Update observation based on conditional variance
         data[t] = errs[t] * sqrt(sigma2[t]);
        }

      // Remove burn-in period to ensure stationarity
      data = np::sliceVector(data, long(burn));
      sigma2 = np::sliceVector(sigma2, long(burn));

      // Return resulting data as a [returns, variance] matrix
      matrix out(data.Size(), 2);
      out.Col(data, 0);
      out.Col(sigma2, 1);

      return out;
     }
   //--- Returns starting values for the HARCH model estimation
   //--- Provides an initial guess based on assumed persistence to help the optimizer converge
   virtual vector    startingValues(vector& resids) override
     {
      // Assume 90% persistence as a standard starting heuristic for long-memory models
      double alpha = 0.9;

      // Initialize omega (constant) based on sample variance, and split alpha persistence equally across all lag groups
      vector sv = (1. - alpha) * resids.Var() * vector::Ones(m_num_lags + 1);

      // Distribute the alpha weight evenly across the specified number of lag frequency groups
      for(ulong i = 1; i < sv.Size(); sv[i] = alpha / double(m_num_lags), ++i);

      return sv;
     }

   //--- Names of model parameters for logging and identification
   //--- Returns a comma-separated string listing omega followed by each frequency group label
   virtual string        parameterNames(void) override
     {
      string names = "omega";
      for(ulong i = 0; i < m_num_lags; ++i)
         // Format: alpha[lag_length]
         StringConcatenate(names, ",", StringFormat("%s%d%s", "alpha[", int(m_lags[i]), "]"));

      return names;
     }
  };
//+------------------------------------------------------------------+
//| FIGARCH process                                                  |
//| Fractionally Integrated GARCH: Models long-memory volatility     |
//| behavior using fractional differencing of the variance process.  |
//+------------------------------------------------------------------+
class CFiGarchProcess: public CVolatilityProcess
  {
protected:
   ulong               m_p;            // Order of the ARCH component
   ulong               m_q;            // Order of the GARCH component
   double              m_power;        // Power transformation (e.g., 2.0 for variance)
   ulong               m_truncation;   // Truncation lag for the infinite MA representation

   //--- Dynamically generate the model name string based on power and structure
   string              _genetrateName(void)
     {
      // Standard variance-based models
      if(m_power == 2.)
        {
         if(m_q == 0)
            return "FIARCH";
         else
            return "FIGARCH";
        }
      // Absolute value-based models (e.g., power=1)
      else
         if(m_power == 1.)
           {
            if(m_q == 0)
               return "FIAVARCH";
            else
               return "FIAVGARCH";
           }
         // Generalized Power models
         else
           {
            if(m_q == 0)
               return StringFormat("Power FIARCH (power: {%.1f})", m_power);
            else
               return StringFormat("Power FIGARCH (power: {%.1f})", m_power);
           }
     }

   //--- Verify if the requested forecasting method is mathematically valid
   //--- Analytic methods require specific power=2.0 structures for multi-step persistence
   virtual bool        _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override
     {
      // Single-step forecasts are always valid
      if(horizon == 1)
         return true;

      // Analytic multi-step forecasts are only supported for standard variance (power=2)
      if(method == FORECAST_ANALYTIC && m_power != 2.)
        {
         Print(__FUNCTION__, ": Analytic forecasts not available for horizon > 1 when power != 2");
         return false;
        }
      return true;
     }
   //--- Analytic multi-step volatility forecasts from the model
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon)
     {
      vector sg;
      matrix fc;
      //--- Establish the initial one-step look-ahead matrix boundaries to seed the multi-step recursion
      if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sg, fc))
         return VarianceForecast();

      //--- Short-circuit execution if only dealing with a single look-ahead window
      if(horizon == 1)
         return VarianceForecast(fc, EMPTY_MATRIX_ARRAY, EMPTY_MATRIX_ARRAY);

      //--- Extract model-specific parameters (skipping omega at parameters[0])
      vector params = np::sliceVector(parameters, 1);

      //--- Pre-calculate the expansion weights up to the allocated truncation lag
      vector lam = figarch_weights(params, m_p, m_q, m_truncation);
      vector lam_rev = lam;

      //--- Reverse the weights vector to align chronological steps for a vector dot product (convolution)
      if(!np::reverseVector(lam))
         return VarianceForecast();

      long t = long(resids.Size());
      double omega = parameters[0];
      double beta = (m_q) ? parameters[parameters.Size() - 1] : 0.0;
      double omega_tilde = omega / (1. - beta); //--- Long-run intercept scaling

      //--- Allocate a unified workspace vector to hold past realized shocks and future forecast values sequentially
      vector temp_forecasts = vector::Zeros(m_truncation + horizon);
      vector resids2 = pow(resids, 2.0); //--- Pre-square raw structural innovations

      //--- Loop across the time-series sliding calculation window starting at '_start'
      for(long i = _start; i < t; ++i)
        {
         long fcast_loc = i - _start;
         long available = long(i + 1 - fmax(0, i - m_truncation + 1));

         //--- Load historical squared innovations into the sliding workspace vector
         for(long j = long(m_truncation - available), k = (long)fmax(0, i - m_truncation + 1); j < long(m_truncation) && k < (i + 1); ++k, ++j)
            temp_forecasts[j] = resids2[k];

         //--- Hand off backcast estimates if historical data depth does not span the truncation requirements
         if(available < long(m_truncation))
            if(!np::fillVector(temp_forecasts, backcast[0], 0, long(m_truncation) - available))
               return VarianceForecast();

         //--- --- Multi-Step Analytic Loop ---
         //--- Step sequentially into the future horizon, resolving expectations via a convolution of reversed weights
         for(ulong h = 0; h < horizon; ++h)
           {
            vector lagged_forecasts = np::sliceVector(temp_forecasts, long(h), long(m_truncation + h));
            temp_forecasts[m_truncation + h] = omega_tilde + lam_rev.Dot(lagged_forecasts);
           }

         //--- Transfer the computed future horizons into the destination forecast result matrix row
         if(!fc.Row(np::sliceVector(temp_forecasts, long(m_truncation)), fcast_loc))
            return VarianceForecast();
        }
      return VarianceForecast(fc, EMPTY_MATRIX_ARRAY, EMPTY_MATRIX_ARRAY);
     }
   //--- Simulation-based volatility forecasts from the model
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      vector sig;
      matrix fc;
      //--- Baseline one-step layout verification
      if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sig, fc))
         return VarianceForecast();

      ulong t = resids.Size();
      matrix paths[], shocks[];
      ArrayResize(paths, int(t - _start));
      ArrayResize(shocks, int(t - _start));

      vector params = np::sliceVector(parameters, 1);
      vector lam = figarch_weights(params, m_p, m_q, m_truncation);
      vector lam_rev = lam;
      if(!np::reverseVector(lam_rev))
         return VarianceForecast();

      double omega = parameters[0];
      double beta = (m_q) ? parameters[parameters.Size() - 1] : 0.0;
      double omega_tilde = omega / (1. - beta);

      //--- Allocate parallel simulation matrices [simulations rows x (truncation + horizon) columns]
      matrix fpath = matrix::Zeros(simulations, m_truncation + horizon);
      vector fresids = pow(fabs(resids), 2.);

      //--- Loop chronologically across the forecasting evaluation footprint
      for(long i = _start; i < long(t); ++i)
        {
         //--- Extract structured random standard normal variables from the random number generator
         matrix std_shocks = rng.rng(simulations, horizon);
         long available = i + 1 - fmax(0, i - long(m_truncation) + 1);
         vector frd = np::sliceVector(fresids, fmax(0, i + 1 - m_truncation), i + 1);

         //--- Inject identical historical realized shocks across all path simulations
         for(long ii = long(m_truncation - available); ii < long(m_truncation); ++ii)
            fpath.Col(frd, ii);

         //--- Pad empty history with backcast values where data boundary overflows
         if(available < long(m_truncation))
            for(long ii = 0; ii < long(m_truncation - available); ++ii)
               fpath.Col(backcast, ii);

         //--- --- Horizon Path Integration Simulation ---
         for(ulong h = 0; h < horizon; ++h)
           {
            //--- Matrix column slice retrieves conditional tracking variables across all simulation paths simultaneously
            matrix lagged_forecasts = np::sliceMatrixCols(fpath, long(h), long(m_truncation + h));
            vector temp = omega_tilde + lagged_forecasts.MatMul(lam_rev);

            //--- Invert the structural transformation power variable back into basic variance space
            vector sigma2 = pow(temp, 2. / m_power);
            long path_loc = i - _start;

            //--- Generate dynamic synthesized shocks and variance values for this horizon index
            shocks[path_loc].Col(std_shocks.Col(h) * sqrt(sigma2), h);
            paths[path_loc].Col(sigma2, h);

            //--- Expected forecast point equals the mathematical mean across all simulated iterations
            fc[path_loc, h] = sigma2.Mean();

            //--- Feed back the generated structural shock into the path matrix for the next recursive horizon step
            fpath.Col(pow(fabs(shocks[path_loc].Col(h)), m_power), m_truncation + h);
           }
        }
      return VarianceForecast(fc, paths, shocks);
     }
   //--- Initialize class members
   virtual bool      _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100) override
     {
      //--- Synchronize configuration parameters to localized private class properties
      m_updateable = updateable;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;
      m_model = volmodel;
      m_p = p;
      m_q = q;
      m_truncation = o; //--- Maps the truncation lag dimension
      m_power = power;  //--- Declares structural model exponent (e.g. 2.0 for FIGARCH)
      m_num_params = 2 + m_p + m_q; //--- Core total parameter footprint (omega + phi + beta + d)

      //--- Initialize the distribution object wrapper
      return m_normal.initialize(vector::Zeros(0), seed);
     }
public:
   // Constructor: Default initialization for a standard FIGARCH(1,d,1) model
   // Uses default truncation of 1000 lags to approximate fractional integration
                     CFiGarchProcess(void)
     {
      m_initialized = _initialize(VOL_FIGARCH, true, false, 4, "FIGARCH", 0, 1, 1000, 1, 2.0);
     }

   // Constructor: Custom initialization for flexible FIGARCH(p,d,q) configurations
   // Allows user to define model orders, power transformation, and truncation depth
                     CFiGarchProcess(ulong p, ulong q, double power, ulong truncation, int seed, long start, long stop, ulong boot)
     {
      // Passes parameters to the internal _initialize method
      // Note: 'truncation' maps to the 'o' (asymmetry/leverage-style) parameter slot in the base process
      m_initialized = _initialize(VOL_FIGARCH, true, false, 0, "FIGARCH", seed, p, truncation, q, power, start, stop, boot);
     }
   //--- Returns starting values for the ARCH model using a grid-search heuristic initialization
   virtual vector    startingValues(vector& resids) override
     {
      //--- Establish parameter boundary array configurations for the grid search optimization initialization
      vector ds = {0.2, 0.5, 0.7};
      vector phi_ratio, beta_ratio;

      if(m_p)
        {
         phi_ratio.Resize(3);
         phi_ratio[0] = 0.2;
         phi_ratio[1] = 0.5;
         phi_ratio[2] = 0.8;
        }
      else
        {
         phi_ratio.Resize(1);
         phi_ratio[0] = 0;
        }

      if(m_q)
        {
         beta_ratio.Resize(3);
         beta_ratio[0] = 0.1;
         beta_ratio[1] = 0.5;
         beta_ratio[2] = 0.9;
        }
      else
        {
         beta_ratio.Resize(1);
         beta_ratio[0] = 0;
        }

      //--- Establish baseline variance expectations derived directly from residuals footprint
      double target = pow(fabs(resids), m_power).Mean();
      double scale  = pow(resids, 2.).Mean() / (pow(target, 2.0 / m_power));
      target *= pow(scale, m_power / 2.);

      //--- Map space to permute all parameter combinations
      matrix all_starting;
      all_starting.Resize(ds.Size() * phi_ratio.Size() * beta_ratio.Size(), 4);
      double beta, omega, d, phi, br, pr;
      vector temp = vector::Zeros(3);
      vector lam;
      ulong row = 0;

      //--- --- Grid Execution Matrix Search Loop ---
      for(ulong i = 0; i < ds.Size(); ++i)
        {
         d = ds[i];
         for(ulong j = 0; j < phi_ratio.Size(); ++j)
           {
            pr = phi_ratio[j];
            phi = (1. - d) / 2. * pr; //--- Force stationary constraints
            for(ulong k = 0; k < beta_ratio.Size(); ++k)
              {
               br = beta_ratio[k];
               beta = (d + phi) * br;

               temp[0] = phi;
               temp[1] = d;
               temp[2] = beta;
               lam = figarch_weights(temp, 1, 1, int(m_truncation));

               //--- Calculate corresponding consistent omega constant for the localized parameter combination
               omega = (1. - beta) * target * (1. - lam.Sum());

               all_starting[row, 0] = omega;
               all_starting[row, 1] = phi;
               all_starting[row, 2] = d;
               all_starting[row++, 3] = beta;
              }
           }
        }

      //--- Extract unique configurations to minimize processing waste
      matrix sv = np::unique(all_starting, 1);

      //--- Strip unused parameter tracking columns if lag settings are set to zero (p=0 or q=0)
      if(!m_q)
         sv = np::sliceMatrixCols(sv, 0, long(sv.Cols() - 1));
      if(!m_p)
        {
         long cols[];
         ArrayResize(cols, int(sv.Cols() - 1));
         cols[0] = 0;
         for(long i = 1; i < long(cols.Size()); cols[i] = i + 1, ++i);
         sv = np::selectMatrixCols(sv, cols);
        }

      matrix vb = varianceBounds(resids);
      vector bc = backCast(resids);
      vector llfs = vector::Zeros(sv.Rows());

      //--- Score every starting candidate against the target log-likelihood function
      for(ulong i = 0; i < sv.Rows(); ++i)
        {
         vector row = sv.Row(i);
         llfs[i] = _gaussianloglikelihood(row, resids, bc, vb);
        }

      //--- Return the parameter configuration vector that optimized/maximized the log-likelihood (ArgMax)
      return sv.Row(llfs.ArgMax());
     }
   //--- Construct values for backcasting to start the recursion
   virtual vector    backCast(vector& resids) override
     {
      //--- Limit calculation depth up to a safe window threshold (e.g. 75 observations)
      ulong tau = MathMin(75, resids.Size());
      vector w = vector::Zeros(tau);

      //--- Generate an exponentially decaying weights structure (RiskMetrics style, lambda = 0.94)
      for(ulong i = 0; i < w.Size(); w[i] = pow(0.94, double(i)), ++i);
      w = w / w.Sum(); //--- Normalize weight values to sum to 1.0

      //--- Extract the start of the series, apply transformation power, and scale by decay weights
      vector sres = np::sliceVector(resids, 0, long(tau));
      sres = pow(fabs(sres), m_power) * w;

      vector bc(1);
      bc[0] = sres.Sum(); //--- Aggregate into scalar backcast representation
      return bc;
     }

   //--- Transformation to apply to user-provided backcast values
   virtual vector    backCastTransform(vector& backcast) override
     {
      //--- Defer to core base class method transformations first
      backcast = CVolatilityProcess::backCastTransform(backcast);
      //--- Apply specific fractional variance power adjustments to the processed vector values
      vector _backcast = pow(sqrt(backcast), m_power);
      return _backcast;
     }

   //--- Construct loose bounds for conditional variances
   virtual matrix    varianceBounds(vector& resids, double power = 2.) override
     {
      //--- Pass computation parameters over to internal system tracker
      return CVolatilityProcess::varianceBounds(resids, m_power);
     }

   //--- Returns boundary matrix structures for parameters optimization constraints
   virtual matrix    bounds(vector& resids) override
     {
      double eps_half = sqrt(2.220446049250313e-16); //--- System machine precision constant limit
      double v = (pow(fabs(resids), m_power)).Mean();
      matrix bnds = matrix::Zeros(m_p + m_q + 2, 2);

      bnds[0, 1] = 10. * v; //--- Enforce an upper bound ceiling on baseline intercept omega constant
      for(ulong i = 1; i < m_p + 1; bnds[i, 1] = 0.5, ++i); //--- Limit short term memory parameters
      for(ulong i = m_p + 1; i < bnds.Rows(); bnds[i, 1] = 1. - eps_half, ++i); //--- Bound persistence variables below 1.0

      return bnds.Transpose();
     }

   //--- Compute the variance for the ARCH model
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds)
     {
      //--- Transform base residuals into process space based on model's power coefficient
      vector fresids = pow(fabs(resids), m_power);
      int nobs = (int)resids.Size();

      //--- Execute the fractional integration recursion engine to write conditional variance calculations directly into _sigma2
      figarch_recursion(parameters, fresids, _sigma2, int(m_p), int(m_q), nobs, m_truncation, _backcast[0], varbounds);

      //--- Apply inverse power mapping to re-align output to standard variance units
      double inv_power = 2. / m_power;
      _sigma2 = pow(_sigma2, inv_power);
      return _sigma2;
     }

   //--- Construct parameter constraints arrays for parameter estimation (Boundary bounding vectors)
   virtual Constraints constraints(void) override
     {
      //--- Initialize hard-coded inequality constraint matrices (A) for OLS/ML parameter bounding
      matrix a =
        {
           {1, 0, 0, 0},
           {0, 1, 0, 0},
           {0, -2, -1, 0},
           {0, 0, 1, 0},
           {0, 0, -1, 0},
           {0, 0, 0, 1},
           {0, 1, 1, -1},
        };
      matrix aa =
        {
           {1, 0, 0},
           {0, 1, 0},
           {0, -1, 0},
           {0, 0, 1},
           {0, 1, -1},
        };
      matrix aaa =
        {
           {1, 0, 0},
           {0, 1, 0},
           {0, -2, -1},
           {0, 0, 1},
           {0, 0, -1}
        };

      //--- Associated boundary mapping verification vectors (B) where: A * parameters >= B
      vector b = {0, 0, -1, 0, -1, 0, 0};
      vector bb = {0, 0, -1, 0, 0};
      vector bbb = {0, 0, -1, 0, -1};

      Constraints constr;
      constr._one = a;
      constr._two = b;

      //--- Adjust structural linear constraint layout mapping depending on specified lag dimension flags (p=0 or q=0)
      if(!m_q)
        {
         constr._one = aaa;
         constr._two = bbb;
        }
      if(!m_p)
        {
         constr._one = aa;
         constr._two = bb;
        }
      return constr;
     }

   //--- Simulate data paths from the model structure
   virtual matrix      simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      vector params = np::sliceVector(parameters, 1);

      //--- Extract long memory hyperbolic expansion decay parameter array
      vector lam = figarch_weights(params, m_p, m_q, m_truncation);
      vector lam_rev = lam;
      if(!np::reverseVector(lam_rev))
         return matrix::Zeros(0, 0);

      //--- Pre-allocate standardized random error distribution array including the requested burn-in footprint
      vector errs = rng.rng(ulong(m_truncation + _nobs + burn));

      //--- Calculate stationary unconditional expectation values if no initial value parameter was provided
      if(initial_value == NULL)
        {
         double persistence = lam.Sum();
         double beta = (m_q) ? parameters[parameters.Size() - 1] : 0.0;
         initial_value = parameters[0];
         if(beta < 1.)
            initial_value /= 1. - beta;
         if(persistence < 1.)
            initial_value /= 1. - persistence;
         if(persistence >= 1. || beta >= 1.)
            Print(__FUNCTION__": WARNING: Initial value warning (non-stationary bounds detected)");
        }

      if(initial_value == NULL)
        {
         Print(__FUNCTION__, ": Initial value is NULL");
         return matrix::Zeros(0, 0);
        }

      //--- Allocate unified arrays tracking simulated states over time
      vector sigma2 = vector::Zeros(ulong(m_truncation + _nobs + burn));
      vector data = sigma2;
      vector fsigma = data;
      vector fdata = fsigma;

      //--- Initialize index blocks spanning pre-sample truncation window constraints
      if(!np::vectorFill(fsigma, initial_value, 0, long(m_truncation)) ||
         !np::vectorFill(sigma2, pow(initial_value, 2. / m_power), 0, long(m_truncation)))
        {
         Print(__FUNCTION__, ": vector fill error");
         return matrix::Zeros(0, 0);
        }

      data = sqrt(sigma2) * errs;
      fdata = pow(fabs(data), m_power);

      double omega = parameters[0];
      double beta = (m_q) ? parameters[parameters.Size() - 1] : 0;
      double omega_tilde = 0;
      if(beta < 1.)
         omega_tilde = omega / (1. - beta);
      else
        {
         Print(__FUNCTION__, ": Beta >= 1.0, using omega as intercept since long-run variance is ill-defined.");
         omega_tilde = omega;
        }

      //--- --- Simulation Time March Loop ---
      //--- Recursively project long-memory volatility evolution paths forward chronologically
      for(long t = long(m_truncation); t < long(m_truncation + _nobs + burn); ++t)
        {
         //--- Convolve the reversed weights against the sliding window vector slice of previous absolute values
         fsigma[t] = omega_tilde + lam_rev.Dot(np::sliceVector(fdata, long(t - m_truncation), t));
         sigma2[t] = pow(fsigma[t], 2.0 / m_power); //--- Convert out of process shape parameters
         data[t] = errs[t] * sqrt(sigma2[t]);        //--- Generate synthetic asset return realization
         fdata[t] = pow(fabs(data[t]), m_power);     //--- Cache internal tracking absolute shape value
        }

      //--- Slice out and discard the initial initialization burn-in sequence block
      data = np::sliceVector(data, long(m_truncation + burn));
      sigma2 = np::sliceVector(sigma2, long(m_truncation + burn));

      //--- Package the output data streams into a 2-column matrix structure [Returns, Variance]
      matrix out = matrix::Zeros(data.Size(), 2);
      if(!out.Col(data, 0) || !out.Col(sigma2, 1))
        {
         Print(__FUNCTION__, ": Column insertion error out ", GetLastError());
         return matrix::Zeros(0, 0);
        }
      return out;
     }

   //--- Names of model parameters
   virtual string       parameterNames(void) override
     {
      string names = "omega";
      if(m_p)
         StringConcatenate(names,",","phi");
      if(m_q)
         StringConcatenate(names,",","beta");
      return names;
     }
  };
//+------------------------------------------------------------------+
//| EGARCH process                                                   |
//| Models volatility as a sum of lags over different frequencies    |
//+------------------------------------------------------------------+
class CEgarchProcess: public CVolatilityProcess
  {
protected:
   ulong               m_p;            // Order of the ARCH component
   ulong               m_o;            // Order of Assymetry component
   ulong               m_q;            // Order of the GARCH component
   double              m_power;        // Power transformation 
   vector              m_vectors[];    // Local container

   matrix            _product(vector& one, vector& two, vector& three)
     {
      ulong rows = 3;
      matrix out(one.Size()*two.Size()*three.Size(),3);
      ulong shift = 0;
      for(ulong i = 0; i<one.Size(); ++i)
        {
         for(ulong j = 0; j<two.Size(); ++j)
           {
            for(ulong k = 0; k<three.Size(); ++k)
              {
               out[shift,0] = one[i];
               out[shift,1] = two[j];
               out[shift,2] = three[k];
               ++shift;
              }
           }
        }
      return out;
     }

   //--- Verify if the requested forecasting method is mathematically valid
   //--- Analytic methods require specific power=2.0 structures for multi-step persistence
   virtual bool        _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override
     {
      // Analytic forecasts are only supported for one step ahead forecasts
      if(method == FORECAST_ANALYTIC && horizon > 1)
        {
         Print(__FUNCTION__, ": Analytic forecasts not available for horizon > 1 ");
         return false;
        }
      return true;
     }
   //--- Analytic multi-step volatility forecasts from the model
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon)
     {
      vector sg;
      matrix fc;
      //--- Establish the initial one-step look-ahead matrix boundaries to seed the multi-step recursion
      if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sg, fc))
         return VarianceForecast();

      return VarianceForecast(fc, EMPTY_MATRIX_ARRAY, EMPTY_MATRIX_ARRAY);
     }
   //--- Simulation-based volatility forecasts from the model
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      vector sig;
      matrix fc;
      //--- Baseline one-step layout verification
      if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sig, fc))
         return VarianceForecast();

      ulong t = resids.Size();
      ulong m = fmax(fmax(m_p,m_o),m_q);

      vector Insigma2 = log(sig);
      vector e = resids/sqrt(sig);

      matrix Insigma2_mat = matrix::Zeros(t,m);
      if(backcast.Size()==m)
        {
         for(ulong i = 0; i<Insigma2_mat.Rows(); ++i)
            Insigma2_mat.Row(backcast,i);
        }
      else
         if(backcast.Size() == 1)
            Insigma2_mat.Fill(backcast[0]);

      matrix e_mat = matrix::Zeros(t,m);
      matrix abs_e_mat = matrix::Zeros(t,m);
      abs_e_mat.Fill(sqrt(2./M_PI));

      for(ulong i = 0; i<m; ++i)
        {
         for(ulong j = (m - i - 1), k = 0; j<t && k<(t - (m - 1)); ++j, ++k)
           {
            Insigma2_mat[j,i] = Insigma2[k];
            e_mat[j,i] = e[k];
            abs_e_mat[j,i] = fabs(e[k]);
           }

        }
      matrix paths[], shocks[];

      ArrayResize(paths, int(t - _start));
      ArrayResize(shocks, int(t - _start));

      for(uint i = 0; i<paths.Size(); ++i)
        {
         paths[i] = matrix::Zeros(simulations,horizon);
         shocks[i] = paths[i];
        }

      double sqrt2pi = sqrt(2./M_PI);

      matrix _lnsigma2 = matrix::Zeros(simulations,m+horizon);
      matrix _e = _lnsigma2;
      matrix _abs_e = _e;
      matrix std_shocks;
      ulong loc;
      for(ulong i = ulong(_start); i<t; ++i)
        {
         std_shocks = rng.rng(simulations,horizon);
         for(ulong k  = 0; k < simulations; ++k)
           {
            for(ulong j = 0; j < m; ++j)
              {
               _lnsigma2[k,j] = Insigma2_mat[i,j];
               _e[k,j] = e_mat[i,j];
               _abs_e[k,j] = abs_e_mat[i,j];
              }
            for(ulong l = m; l < m+horizon; ++l)
              {
               for(ulong n = 0; n<horizon; ++n)
                 {
                  _e[k,l] = std_shocks[k,n];
                  _abs_e[k,l] = fabs(std_shocks[k,n]);
                 }
              }
           }
         for(ulong j = 0; j < horizon; ++j)
           {
            loc  = 0;
            for(ulong jj = 0; jj< _lnsigma2.Rows(); ++jj)
               _lnsigma2[jj,m + j] = parameters[loc];
            loc += 1;
            for(ulong k = 0; k < m_p; ++k)
              {
               for(ulong jj = 0; jj< _lnsigma2.Rows(); ++jj)
                  _lnsigma2[jj,m + j] += parameters[loc] * (_abs_e[jj, m + j - 1 - k] - sqrt2pi);
               loc += 1;
              }
            for(ulong k = 0; k < m_o; ++k)
              {
               for(ulong jj = 0; jj< _lnsigma2.Rows(); ++jj)
                  _lnsigma2[jj,m + j] += parameters[loc] * _e[jj, m + j - 1 - k];
               loc += 1;
              }

            for(ulong k = 0; k < m_q; ++k)
              {
               for(ulong jj = 0; jj< _lnsigma2.Rows(); ++jj)
                  _lnsigma2[jj,m + j] += parameters[loc] * _lnsigma2[jj, m + j - 1 - k];
               loc += 1;
              }
           }
         loc = i - ulong(_start);
         paths[loc] = exp(np::sliceMatrixCols(_lnsigma2,long(m)));
         shocks[loc] = sqrt(paths[loc]) * std_shocks;
        }

      fc = matrix::Zeros(paths.Size(),paths[0].Cols());
      for(ulong i = 0; i < fc.Rows(); ++i)
         if(!fc.Row(paths[i].Mean(0),i))
           {
            Print(__FUNCTION__,": Row insertion error ", GetLastError());
            return VarianceForecast();
           }

      return VarianceForecast(fc,paths,shocks);
     }
   //--- Initialize class members
   virtual bool      _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100) override
     {
      //--- Synchronize configuration parameters to localized private class properties
      m_updateable = updateable;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;
      m_model = volmodel;
      m_p = p;
      m_o = o;
      m_q = q;
      m_power = power;
      m_num_params = 1 + m_o + m_p + m_q;
      m_name = (m_q>0)?"EGARCH":"EARCH";

      //--- Initialize the distribution object wrapper
      return m_normal.initialize(vector::Zeros(0), seed);
     }
public:
                     CEgarchProcess(void)
     {
      m_initialized = _initialize(VOL_EGARCH,true,false,0,NULL,0,1,0,1,2.0,0,-1,100);
     }
                     CEgarchProcess(ulong p, ulong o, ulong q,double power, int seed, long start, long stop, ulong boot)
     {
      m_initialized = _initialize(VOL_EGARCH,true,false,0,NULL,seed,p,o,q,2.0,start,stop,boot);
     }

   //--- Returns starting values for the ARCH model using a grid-search heuristic initialization
   virtual vector    startingValues(vector& resids) override
     {

      vector alphas = {0.01, 0.05, 0.1, 0.2};
      vector gammas = {-0.1, 0.0, 0.1};
      vector betas = {0.5, 0.7, 0.9, 0.98};
      matrix agbs = _product(alphas,gammas,betas);

      double target = log(pow(resids,2).Mean());
      
      vector svs[];
      ArrayResize(svs, int(agbs.Rows()));
      matrix vb = varianceBounds(resids);
      vector llfs = vector::Zeros(agbs.Rows());
      vector bc = backCast(resids);

      // Evaluate Log-Likelihood for each combination in the grid
      for(uint i = 0; i < svs.Size(); ++i)
        {
         vector row = agbs.Row(i);
         // Construct starting vector: constant omega + lags
         vector sv = (1. - row[2]) * target * vector::Ones(m_p + m_o + m_q + 1);

         // Distribute coefficients across ARCH (p), Leverage (o), and GARCH (q) components
         if(m_p > 0)
            np::fillVector(sv, row[0] / double(m_p), long(1), long(1 + m_p));
         if(m_o > 0)
            np::fillVector(sv, row[1] / double(m_o), long(1 + m_p), long(1 + m_p + m_o));
         if(m_q > 0)
            np::fillVector(sv, row[2] / double(m_q), long(1 + m_p + m_o), long(1 + m_p + m_o + m_q));

         svs[i] = sv;
         llfs[i] = _gaussianloglikelihood(sv, resids, bc, vb);
        }

      // Select the parameter set that yielded the highest Log-Likelihood
      ulong loc = llfs.ArgMax();
      return svs[loc];
     }
   //--- Construct values for backcasting to start the recursion
   virtual vector    backCast(vector& resids) override
     {
      return log(CVolatilityProcess::backCast(resids));
     }

   //--- Transformation to apply to user-provided backcast values
   virtual vector    backCastTransform(vector& backcast) override
     {
      //--- Defer to core base class method transformations first
      backcast = CVolatilityProcess::backCastTransform(backcast);
      return log(backcast);
     }

   //--- Construct loose bounds for conditional variances
   virtual matrix    varianceBounds(vector& resids, double power = 2.) override
     {
      //--- Pass computation parameters over to internal system tracker
      return CVolatilityProcess::varianceBounds(resids, m_power);
     }

   //--- Returns boundary matrix structures for parameters optimization constraints
   virtual matrix    bounds(vector& resids) override
     {
      double log_const = log(10000.0);
      double v = (pow(fabs(resids), m_power)).Mean();
      matrix bnds = matrix::Zeros(1 + m_p + m_o + m_q, 2);

      bnds[0, 0] = log(v) - log_const; 
      bnds[0, 1] = log(v) + log_const;
      for(ulong i = 1; i < m_p + m_o + 1; bnds[i, 0] = -1*double("inf"), bnds[i,1] = double("inf"), ++i); //--- Limit short term memory parameters
      for(ulong i = m_p + m_o + 1; i < bnds.Rows(); bnds[i, 1] = double(m_q), ++i); 

      return bnds.Transpose();
     }

   //--- Compute the variance for the ARCH model
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds)
     {
      //--- Transform base residuals into process space based on model's power coefficient
      //vector lnsigma2, std_resids, abs_std_resids;
      int nobs = (int)resids.Size();
      
      if(!m_vectors.Size() || m_vectors[0].Size() != resids.Size())
       {
        ArrayResize(m_vectors,3);
        m_vectors[0] = m_vectors[1] = m_vectors[2] = vector::Zeros(nobs);
       }
       
      egarch_recursion(parameters,resids,_sigma2,m_p,m_o,m_q,ulong(nobs),_backcast[0],varbounds,m_vectors[0],m_vectors[1],m_vectors[2]);
      
      return _sigma2;
     }

   //--- Construct parameter constraints arrays for parameter estimation (Boundary bounding vectors)
   virtual Constraints constraints(void) override
     {
      matrix a = matrix::Zeros(1,m_p+m_o+m_q+1);
      vector b = vector::Zeros(1);
      b[0] = -1.0;
      for(ulong i = m_p+m_o+1; i<a.Cols(); a[0,i] = -1.0, ++i);
      
      Constraints constr;
      constr._one = a;
      constr._two = b;
      
      return constr;
     }

   //--- Simulate data paths from the model structure
   virtual matrix      simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      vector params;// = np::sliceVector(parameters, 1);
      vector errs = rng.rng(_nobs+burn);
      double beta_sum = 0.0;
      if(initial_value == NULL)
       {
        if(m_q)
         {
          vector prm = np::sliceVector(parameters,long(m_p+m_o+1));
          beta_sum = prm.Sum();
         }
        else 
         beta_sum = 0.0;
        if(beta_sum<1.)
         initial_value = parameters[0]/(1.0/beta_sum);
        else
         initial_value = parameters[0];
       }

      
      //--- Allocate unified arrays tracking simulated states over time
      vector sigma2 = vector::Zeros(ulong(_nobs + burn));
      vector data = sigma2;
      vector lnsigma2 = data;
      vector abserrors = fabs(errs);
      double norm_const = sqrt(2./M_PI);
      ulong max_lag = fmax(fmax(m_p,m_o),m_q);

      //--- Initialize index blocks spanning pre-sample truncation window constraints
      if(!np::vectorFill(lnsigma2, initial_value, 0, long(max_lag)) ||
         !np::vectorFill(sigma2, exp(initial_value), 0, long(max_lag)))
        {
         Print(__FUNCTION__, ": vector fill error");
         return matrix::Zeros(0, 0);
        }

      data = sqrt(sigma2) * errs;

      double omega = parameters[0];
      double beta = (m_q) ? parameters[parameters.Size() - 1] : 0;
      double omega_tilde = 0;
      if(beta < 1.)
         omega_tilde = omega / (1. - beta);
      else
        {
         Print(__FUNCTION__, ": Beta >= 1.0, using omega as intercept since long-run variance is ill-defined.");
         omega_tilde = omega;
        }

      //--- --- Simulation Time March Loop ---
      ulong loc = 0;
      for(long t = long(max_lag); t < long(_nobs + burn); ++t)
        {
         loc = 0;
         lnsigma2[t] = parameters[loc];
         loc+=1;
         for(ulong j = 0; j<m_p; ++j)
           {
            lnsigma2[t] += parameters[loc] * (abserrors[t - 1 - j] - norm_const);
            loc+=1;
           }
          for(ulong j = 0; j<m_o; ++j)
           {
            lnsigma2[t] += parameters[loc] * (errs[t - 1 - j]);
            loc+=1;
           }
          for(ulong j = 0; j<m_q; ++j)
           {
            lnsigma2[t] += parameters[loc] * (lnsigma2[t - 1 - j]);
            loc+=1;
           }
         }

      //--- Slice out and discard the initial initialization burn-in sequence block
      data = errs*sqrt(sigma2);
      data = np::sliceVector(data, long(burn));
      sigma2 = exp(lnsigma2);
      sigma2 = np::sliceVector(sigma2, long(burn));

      //--- Package the output data streams into a 2-column matrix structure [Returns, Variance]
      matrix out = matrix::Zeros(data.Size(), 2);
      if(!out.Col(data, 0) || !out.Col(sigma2, 1))
        {
         Print(__FUNCTION__, ": Column insertion error out ", GetLastError());
         return matrix::Zeros(0, 0);
        }
      return out;
     }

   //--- Names of model parameters
   virtual string       parameterNames(void) override
     {
      return common_names(m_p,m_o,m_q);
     }
  };
//+------------------------------------------------------------------+