//+------------------------------------------------------------------+
//|                                                         mean.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"volatility.mqh"
#ifdef __SLSQP__
#include"..\..\slsqp.mqh"
#else
#include<Math\Alglib\delegatefunctions.mqh>
#include<Math\Alglib\optimization.mqh>
#endif
//+------------------------------------------------------------------+
//| Abstract base class for mean models                              |
//+------------------------------------------------------------------+
class CArchModel: public CObject
  {
protected:
   //--- Class Members
   ArchParameters    m_model_spec;       // Structure containing execution specifications
   vector            m_params;           // Optimized joint parameter values array
   vector            m_fit_indices;      // Index flags mapping timestamps during optimization
   vector            m_backcast;         // Initial pre-sample variance conditions vector
   vector            m_fit_y;            // In-sample target timeline observation variables
   matrix            m_fit_x;            // Independent exogenous variables design matrix
   matrix            m_fit_regressors;   // Cached design matrix mapping lags/regressors
   double            m_scale;            // Internal adjustment scaling multiplier
   ulong             m_nobs;             // Number of total parsed dataset elements
   ulong             m_holdback;         // Truncated initial sequence element offset
   ulong             m_maxlags;          // Maximum depth threshold across mean/volatility processes
   ulong             m_num_params;       // Total dedicated mean-model estimation parameter counts
   bool              m_converged;        // Boolean indicating successful optimizer exit metrics
   bool              m_rescale;          // Operational flag forcing internal auto-scaling routines
   bool              m_constant;         // Structural toggle introducing intersect intercepts
   bool              m_rotated;          // Indicator mapping variance targeting or HAR adjustments
   string            m_name;             // Custom textual identifier for model variants
   matrix            m_regressors;       // Dynamic matrix tracking user external series structures
   matrix            m_lags;             // Historical lag tracking layout metrics
   matrix            m_var_bounds;       // Linear constraint boundaries array for variance limits
   bool              m_delete_vp;        // Lifecycle memory management tracking for Volatility object
   bool              m_delete_dist;      // Lifecycle memory management tracking for Distribution object

   CDistribution      *m_distribution;   // Object pointer managing target error innovations
   CVolatilityProcess *m_vp;             // Object pointer managing conditional variance steps

   bool                m_initialized;    // System lock indicating model verification validity

   //--- Calculate epsilon vector for numerical optimization step intervals
   vector            _get_epsilon(vector &x, double s, ulong n, double epsilon = EMPTY_VALUE)
     {
      vector out;
      if(epsilon == EMPTY_VALUE)
         // Machine epsilon scaling method based on parameter magnitude bounds
         out = pow(DBL_EPSILON,1./s)*np::maximum(MathAbs(x),0.1);
      else
        {
         out = vector::Zeros(n);
         out.Fill(epsilon);
        }
      return out;
     }

   //--- Approximates the partial derivatives (Jacobian/Gradient vector) of the objective function
   matrix            approx_fprime(vector &x,vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false,bool centered = false, double epsilon = EMPTY_VALUE)
     {
      ulong n = x.Size();
      vector f0 = _loglikelihood(x,sigma2,backcast,varbounds,individual);
      matrix grad = matrix::Zeros(n,f0.Size());
      vector ei = vector::Zeros(n);
      vector eps,row;

      if(!centered) // Standard Forward Finite-Difference estimation
        {
         eps = _get_epsilon(x,2.,n,epsilon);
         for(ulong i = 0; i<n; ++i)
           {
            ei[i] = eps[i];
            row = (_loglikelihood(x+ei,sigma2,backcast,varbounds,individual) - f0)/eps[i];
            ei[i] = 0.0;
            grad.Row(row,i);
           }
        }
      else // Higher accuracy O(h^2) Central Finite-Difference estimation
        {
         eps = _get_epsilon(x,3,n,epsilon)/2.;
         for(ulong i = 0; i<n; ++i)
           {
            ei[i] = eps[i];
            row = (_loglikelihood(x+ei,sigma2,backcast,varbounds,individual) - _loglikelihood(x-ei,sigma2,backcast,varbounds,individual))/(2. * eps[i]);
            ei[i] = 0.0;
            grad.Row(row,i);
           }
        }
      return grad.Transpose();
     }

   //--- Approximates the structural Hessian matrix via localized finite-difference increments
   matrix            approx_hess(vector &x,vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false,double epsilon = EMPTY_VALUE)
     {
      ulong n = x.Size();
      vector h = _get_epsilon(x,4,n,epsilon);
      matrix ee;
      ee.Diag(h); // Diagonal perturbation step sizing arrangement
      matrix hess = h.Outer(h);

      // Compute symmetric multi-directional cross-partial evaluations
      for(ulong i = 0; i<n; ++i)
        {
         for(ulong j = i; j<n; ++j)
           {
            vector temp = (_loglikelihood(x+ee.Row(i)+ee.Row(j),sigma2,backcast,varbounds,individual) -
                           _loglikelihood(x+ee.Row(i)-ee.Row(j),sigma2,backcast,varbounds,individual) -
                           (_loglikelihood(x-ee.Row(i)+ee.Row(j),sigma2,backcast,varbounds,individual) -
                            _loglikelihood(x-ee.Row(i)-ee.Row(j),sigma2,backcast,varbounds,individual)))/(4.*hess[i,j]);
            hess[i,j] = temp[0];
            hess[j,i] = hess[i,j]; // Enforce mathematical matrix symmetry
           }
        }
      return hess;
     }

   //--- Rescales the residuals vector to defend optimizer from underflow/overflow bounds
   void              _checkscale(vector& resids)
     {
      double ogscale = resids.Var();
      double scale = ogscale;
      double rescale = 1.;

      // Loop iteratively to determine decimal step shifting scale bounds
      while((scale>=1000.|| scale<0.1) && scale>0.0)
        {
         if(scale<1.0)
            rescale*=10.0;
         else
            rescale/=10.0;
         scale = ogscale*pow(rescale,2.0);
        }
      if(rescale == 1.0)
         return;
      if(!m_rescale)
        {
         Print(__FUNCTION__, " y is poorly scaled, which may affect convergence of the optimizer when estimating the model parameters.\n","Calculated scale is ",ogscale," . Parameter estimation works better when this value is between 1 and 1000. The recommended rescaling is ",rescale," * y  ");
         return;
        }
      m_scale = rescale;
     }

   //--- Virtual hook: triggered natively when localized array tracking scales change
   virtual void      _scalechanged(void)
     {
      return;
     }

   //--- Native placeholder computing the structural Coefficient of Determination ($R^2$)
   virtual double    _r2(vector& params)
     {
      return EMPTY_VALUE;
     }

   //--- Virtual hook: calculates structural analytical start guesses bypassing ARCH processes
   virtual vector    _fit_no_arch_normal_errors_params(void)
     {
      return vector::Zeros(0);
     }

   //--- Evaluates a baseline static Gaussian Log-Likelihood sequence
   double            _static_gaussian_loglikelihood(vector &resids_)
     {
      ulong nobs_ = resids_.Size();
      double sigma2 = resids_.Dot(resids_)/double(nobs_);
      double ll = -0.5 * double(nobs_) * log(2.0 * M_PI);
      ll -= 0.5*double(nobs_)*log(sigma2);
      ll -= 0.5*double(nobs_);
      return ll;
     }

   //--- Unpacks the unified parameter array into mean, volatility, and error component slices
   bool              _parse_parameters(vector &x, vector& out[])
     {
      ulong km = m_num_params;
      ulong kv = m_vp.numParams();
      ArrayResize(out,3);

      // Slice array segments explicitly based on model configuration widths
      out[0] = (km)?np::sliceVector(x,0,long(km)):vector::Zeros(0);
      out[1] = (kv)?np::sliceVector(x,long(km),long(km+kv)):vector::Zeros(0);
      out[2] = ((km+kv)<x.Size())?np::sliceVector(x,long(km+kv)):vector::Zeros(0);
      return true;
     }

   //--- Computes Negative Log-Likelihood aggregates across all linked subsystems
   vector            _loglikelihood(vector& parameters, vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false)
     {
      vector prs[];
      _parse_parameters(parameters,prs);

      // Extract sequential system steps
      vector resids = resids(prs[0],EMPTY_VECTOR,EMPTY_MATRIX);
      sigma2 = m_vp.computeVariance(prs[1],resids,sigma2,backcast,varbounds);
      vector llf = m_distribution.loglikelihood(prs[2],resids,sigma2,individual);

      return (-1.*llf); // Inverse returns optimization-ready minimized objective shape
     }

   //--- Performs localized dataset boundary shifting adjustments
   virtual void       _adjust_sample(long& first, long& last)
     {
      return;
     }

   //--- Safely flushes dynamic sub-allocations on execution termination sequences
   virtual void      deinitialize(void)
     {
      if(m_delete_vp && CheckPointer(m_vp) == POINTER_DYNAMIC)
         delete m_vp;
      m_vp = NULL;
      if(m_delete_dist && CheckPointer(m_distribution) == POINTER_DYNAMIC)
         delete m_distribution;
      m_distribution = NULL;
     }

public:
   //+------------------------------------------------------------------+
   //| Constructor & Getters                                            |
   //+------------------------------------------------------------------+
                     CArchModel(void):m_vp(NULL),
                     m_initialized(false),
                     m_distribution(NULL),
                     m_nobs(0),
                     m_holdback(0),
                     m_converged(false),
                     m_rescale(false),
                     m_constant(true),
                     m_delete_vp(true),
                     m_delete_dist(true),
                     m_scale(1.0),
                     m_name(NULL),
                     m_fit_y(vector::Zeros(0)),
                     m_params(vector::Zeros(0)),
                     m_fit_indices(vector::Zeros(0)),
                     m_fit_regressors(matrix::Zeros(0,0)),
                     m_fit_x(matrix::Zeros(0,0)),
                     m_regressors(matrix::Zeros(0,0))
     {
     }

   ulong              nobs(void) { return m_nobs; }
   ulong              hold_back(void) { return m_holdback; }
   bool               converged(void) { return m_converged;}
   bool               rescaled(void) { return m_rescale; }
   double             get_scale(void) { return m_scale; }
   vector             get_y(void) { return m_model_spec.observations; }
   vector             get_params(void) { return m_params; }
   ArchParameters     get_specification(void) { return m_model_spec; }

   //+------------------------------------------------------------------+
   //| File Input/Output Binary Serialization Handlers                  |
   //+------------------------------------------------------------------+
   virtual bool      Save(const int file_handle) override
     {
      if(file_handle!=INVALID_HANDLE)
        {
         // Write validation marker block
         if(FileWriteLong(file_handle,-1)==sizeof(long))
           {
            // Serialize Parameter structures sizes and contents
            if(FileWriteInteger(file_handle,int(m_params.Size()),INT_VALUE)!=INT_VALUE)
               return(false);
            for(ulong i = 0; i<m_params.Size(); ++i)
              {
               if(FileWriteDouble(file_handle,m_params[i])!=sizeof(double))
                 {
                  Print(__FUNCTION__, " file write error ", GetLastError());
                  return false;
                 }
              }

            // Serialize Mean Model configurations
            if(FileWriteInteger(file_handle,int(m_model_spec.mean_model_type),INT_VALUE)!=INT_VALUE ||
               FileWriteLong(file_handle,long(m_model_spec.holdout_size))!=sizeof(long) ||
               FileWriteInteger(file_handle,int(m_model_spec.is_rescale_enabled),INT_VALUE)!=INT_VALUE ||
               FileWriteInteger(file_handle,int(m_model_spec.include_constant),INT_VALUE)!=INT_VALUE ||
               FileWriteInteger(file_handle,int(m_model_spec.use_har_rotation),INT_VALUE)!=INT_VALUE)
               return(false);

            if(FileWriteDouble(file_handle,m_model_spec.scaling_factor)!=sizeof(double))
               return(false);
            if(FileWriteLong(file_handle,long(m_model_spec.mean_lags.Size()))!=sizeof(long))
               return(false);

            for(ulong i = 0; i<m_model_spec.mean_lags.Size(); ++i)
              {
               if(FileWriteDouble(file_handle,m_model_spec.mean_lags[i])!=sizeof(double))
                 {
                  Print(__FUNCTION__, " file write error ", GetLastError());
                  return false;
                 }
              }

            // Serialize Volatility specification tracking indices
            if(FileWriteInteger(file_handle,int(m_model_spec.vol_model_type),INT_VALUE)!=INT_VALUE)
               return(false);
            if(FileWriteInteger(file_handle,int(m_model_spec.vol_rng_seed),INT_VALUE)!=INT_VALUE ||
               FileWriteLong(file_handle,m_model_spec.sample_start_idx)!=sizeof(long) ||
               FileWriteLong(file_handle,m_model_spec.sample_end_idx)!=sizeof(long) ||
               FileWriteLong(file_handle,long(m_model_spec.min_bootstrap_sims))!=INT_VALUE ||
               FileWriteLong(file_handle,long(m_model_spec.garch_p))!=sizeof(long) ||
               FileWriteLong(file_handle,long(m_model_spec.garch_o))!=sizeof(long) ||
               FileWriteLong(file_handle,long(m_model_spec.garch_q))!=sizeof(long))
               return(false);

            if(FileWriteDouble(file_handle,m_model_spec.vol_power)!=sizeof(double))
               return false;
            if(FileWriteInteger(file_handle,int(m_model_spec.dist_type),INT_VALUE)!=INT_VALUE)
               return(false);
            if(FileWriteInteger(file_handle,int(m_model_spec.dist_rng_seed),INT_VALUE)!=INT_VALUE)
               return(false);
            if(FileWriteLong(file_handle,long(m_model_spec.observations.Size()))!=sizeof(long))
               return(false);

            for(ulong i = 0; i<m_model_spec.observations.Size(); ++i)
              {
               if(FileWriteDouble(file_handle,m_model_spec.observations[i])!=sizeof(double))
                 {
                  Print(__FUNCTION__, " file write error ", GetLastError());
                  return false;
                 }
              }

            // Serialize multi-dimensional Exogenous Matrix arrays
            if(FileWriteLong(file_handle,long(m_model_spec.exog_data.Rows()))!=sizeof(long))
               return(false);
            if(FileWriteLong(file_handle,long(m_model_spec.exog_data.Cols()))!=sizeof(long))
               return(false);
            for(ulong i = 0; i<m_model_spec.exog_data.Rows(); ++i)
              {
               for(ulong j = 0; j<m_model_spec.exog_data.Cols(); ++j)
                 {
                  if(FileWriteDouble(file_handle,m_model_spec.exog_data[i,j])!=sizeof(double))
                    {
                     Print(__FUNCTION__, " file write error ", GetLastError());
                     return false;
                    }
                 }
              }
            return true;
           }
        }
      return false;
     }

   virtual bool      Load(const int file_handle) override
     {
      if(file_handle!=INVALID_HANDLE)
        {
         long ff = FileReadLong(file_handle);
         if(ff == -1)
           {
            ResetLastError();
            ulong size = ulong(FileReadLong(file_handle));
            if(size)
              {
               m_params.Resize(size);
               m_converged = true;
               for(ulong i = 0; i<m_params.Size(); ++i)
                  m_params[i]=FileReadDouble(file_handle);
               if(GetLastError())
                 {
                  Print(__FUNCTION__, " possible read error ", GetLastError());
                  return false;
                 }
              }
            else
              {
               m_converged = false;
               m_params = vector::Zeros(0);
              }

            // Unpack deserialized variables back into configuration structure
            m_model_spec.mean_model_type = ENUM_MEAN_MODEL(FileReadInteger(file_handle,INT_VALUE));
            m_model_spec.holdout_size = ulong(FileReadLong(file_handle));
            m_model_spec.is_rescale_enabled = bool(FileReadInteger(file_handle,INT_VALUE));
            m_model_spec.include_constant = bool(FileReadInteger(file_handle,INT_VALUE));
            m_model_spec.use_har_rotation = bool(FileReadInteger(file_handle,INT_VALUE));

            if(GetLastError())
              {
               Print(__FUNCTION__, " possible read error ", GetLastError());
               return false;
              }

            m_model_spec.scaling_factor=FileReadDouble(file_handle);
            size = ulong(FileReadLong(file_handle));
            if(size)
              {
               m_model_spec.mean_lags.Resize(size);
               for(ulong i = 0; i<m_model_spec.mean_lags.Size(); ++i)
                  m_model_spec.mean_lags[i]=FileReadDouble(file_handle);
              }
            else
              {
               m_model_spec.mean_lags = vector::Zeros(0);
              }

            if(GetLastError())
              {
               Print(__FUNCTION__, " possible read error ", GetLastError());
               return false;
              }

            m_model_spec.vol_model_type = ENUM_VOLATILITY_MODEL(FileReadInteger(file_handle,INT_VALUE));
            m_model_spec.vol_rng_seed = (FileReadInteger(file_handle,INT_VALUE));
            m_model_spec.sample_start_idx = long(FileReadLong(file_handle));
            m_model_spec.sample_end_idx = long(FileReadLong(file_handle));
            m_model_spec.min_bootstrap_sims = ulong(FileReadLong(file_handle));
            m_model_spec.garch_p = ulong(FileReadLong(file_handle));
            m_model_spec.garch_o = ulong(FileReadLong(file_handle));
            m_model_spec.garch_q = ulong(FileReadLong(file_handle));

            m_model_spec.vol_power = FileReadDouble(file_handle);
            m_model_spec.dist_type = ENUM_DISTRIBUTION_MODEL(FileReadInteger(file_handle,INT_VALUE));
            m_model_spec.dist_rng_seed = FileReadInteger(file_handle,INT_VALUE);

            if(GetLastError())
              {
               Print(__FUNCTION__, " possible read error ", GetLastError());
               return false;
              }

            size = ulong(FileReadLong(file_handle));
            if(size)
              {
               m_model_spec.observations.Resize(size);
               for(ulong i = 0; i<m_model_spec.observations.Size(); ++i)
                  m_model_spec.observations[i] = FileReadDouble(file_handle);
              }
            else
               m_model_spec.observations = vector::Zeros(0);

            long rws = FileReadLong(file_handle);
            long cls = FileReadLong(file_handle);
            m_model_spec.exog_data = matrix::Zeros(rws,cls);
            if(m_model_spec.exog_data.Rows())
              {
               for(ulong i = 0; i<m_model_spec.exog_data.Rows(); ++i)
                 {
                  for(ulong j = 0; j<m_model_spec.exog_data.Cols(); ++j)
                     m_model_spec.exog_data[i,j] = FileReadDouble(file_handle);
                 }
              }

            if(GetLastError())
              {
               Print(__FUNCTION__, " possible read error ", GetLastError());
               return false;
              }
            return true;
           }
        }
      return false;
     }

   virtual bool       is_initialized(void) { return m_initialized; }
                    ~CArchModel(void)     { deinitialize(); }

   //+------------------------------------------------------------------+
   //| Core Volatility and Distribution Access Pointers                 |
   //+------------------------------------------------------------------+
   CVolatilityProcess *volatility(void)          { return m_vp; }
   CDistribution      *distribution(void)        { return m_distribution; }
   ENUM_MEAN_MODEL     mean_model_type(void)     { return m_model_spec.mean_model_type; }
   ENUM_VOLATILITY_MODEL volatility_process_type(void) { return m_vp.volatilityprocess(); }
   ENUM_DISTRIBUTION_MODEL distribution_type(void)    { return m_distribution.distribution_type(); }

   //--- Wrapper calling private log-likelihood objective processes
   vector            objective(vector& parameters, vector& sigma2, vector &backcast, matrix& varbounds, bool individual=false)
     {
      return _loglikelihood(parameters,sigma2,backcast,varbounds,individual);
     }

   //--- Wrapper tracking optimization Jacobian matrices
   matrix            jacobian(vector& parameters, vector& sigma2, vector &backcast, matrix& varbounds, bool individual=false,bool centered=false, double eps=EMPTY_VALUE)
     {
      return approx_fprime(parameters,sigma2,backcast,varbounds,individual,centered,eps);
     }

   virtual vector    starting_values(void)
     {
      return _fit_no_arch_normal_errors_params();
     }

   virtual vector    resids(vector& params, vector& y, matrix& regressors)
     {
      return vector::Zeros(0);
     }

   //--- Calculates the Parameter Covariance Matrix using robust Sandwich Estimation ($H^{-1} J H^{-1}$)
   matrix            compute_param_cov(vector& params, vector& backcast, bool robust = true)
     {
      vector sv_ = starting_values();
      vector fy = vector::Zeros(0);
      matrix fr = matrix::Zeros(0,0);
      vector resids_ = resids(sv_,fy,fr);
      matrix vb = m_vp.varianceBounds(resids_);
      ulong nobs = resids_.Size();

      if(!backcast.Size() && !m_backcast.Size())
        {
         backcast = m_vp.backCast(resids_);
         m_backcast = backcast;
        }
      else
         if(!backcast.Size())
            backcast = m_backcast;

      vector sgma2=vector::Zeros(nobs);
      matrix hess = approx_hess(params,sgma2,backcast,vb,false);
      hess/=double(nobs);
      matrix inv_hess = hess.Inv(); // Invert the Hessian ($H^{-1}$)

      if(robust) // Quasi-Maximum Likelihood (QMLE) Robust Sandwich Variance Covariance calculation
        {
         matrix scores = approx_fprime(params,sgma2,backcast,vb,true);
         matrix score_cov = scores.Transpose().Cov(); // Inner Outer-Product Jacobian ($J$)
         return (inv_hess.MatMul(score_cov).MatMul(inv_hess))/double(nobs);
        }
      else
         return inv_hess/double(nobs);
     }

   virtual Constraints          constraints(void)
     {
      Constraints out;
      return out;
     }

   virtual matrix                bounds(void)
     {
      vector temp = {AL_NEGINF,AL_POSINF};
      return np::repeat_vector_as_rows_cols(temp,m_num_params,false);
     }

   virtual string    name(void)
     {
      return m_name;
     }

   virtual ArchForecast forecast(vector& params,matrix& x[],ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
     {
      ArchForecast out;
      return out;
     }
  };
//+------------------------------------------------------------------+
//| Arch model fixed result struct                                   |
//| Holds the optimization outputs and diagnostic metrics for a      |
//| finalized ARCH/GARCH execution state.                             |
//+------------------------------------------------------------------+
struct ArchModelFixedResult
  {
public:
   double            loglikelihood;          // Maximized Joint Log-Likelihood objective value
   vector            params;                 // Final estimated parameter values vector
   vector            conditional_volatility; // Modeled time-varying conditional standard deviation (\sigma_t)
   ulong             nobs;                   // Total sample size count used in fitting
   vector            resid;                  // Extracted mean model residual errors (\epsilon_t)

   //--- Default Constructor (Initializes to empty state)
                     ArchModelFixedResult(void)
     {
      loglikelihood = EMPTY_VALUE;
      nobs = 0;
      params = resid = conditional_volatility = vector::Zeros(0);
     }

   //--- Copy Constructor
                     ArchModelFixedResult(ArchModelFixedResult &other)
     {
      loglikelihood = other.loglikelihood;
      nobs = other.nobs;
      params = other.params;
      resid = other.resid;
      conditional_volatility = other.conditional_volatility;
     }

   //--- Parameterized Constructor
                     ArchModelFixedResult(vector &_params, vector &_resid, vector &_volatility, vector &_depvar, double _loglikelihood)
     {
      loglikelihood = _loglikelihood;
      params = _params;
      resid = _resid;
      conditional_volatility = _volatility;
      nobs = _resid.Size(); // Dynamically infers observation size from residual length
     }

   //--- Destructor
                    ~ArchModelFixedResult(void)
     {
     }

   //--- Overloaded Assignment Operator
   void              operator=(ArchModelFixedResult &other)
     {
      loglikelihood = other.loglikelihood;
      nobs = other.nobs;
      params = other.params;
      resid = other.resid;
      conditional_volatility = other.conditional_volatility;
     }

   //--- Akaike Information Criterion (AIC): Penalizes model complexity linearly
   double            aic(void)
     {
      return -2 * loglikelihood + 2 * double(num_params());
     }

   //--- Returns the total parameter count estimated across the complete system
   ulong             num_params(void)
     {
      return params.Size();
     }

   //--- Schwarz/Bayesian Information Criterion (BIC): Penalizes complexity logarithmic to sample size
   double            bic(void)
     {
      return -2.*loglikelihood+log(nobs)*double(num_params());
     }

   //--- Returns standardized residuals (z_t = \epsilon_t / \sigma_t) for distribution testing
   vector            std_resid(void)
     {
      return resid/conditional_volatility;
     }
  };
//+------------------------------------------------------------------+
//| Arch model result struct                                         |
//| Inherits base properties from ArchModelFixedResult to add        |
//| statistical inference methods and parameter covariance parsing.   |
//+------------------------------------------------------------------+
struct ArchModelResult: public ArchModelFixedResult
  {
public:
   long              fit_indices[2];         // [0] = Start index, [1] = Stop index mapping the target data slice
   matrix            param_cov;              // Calculated parameter Variance-Covariance matrix
   double            r2;                     // Unadjusted Coefficient of Determination (R-squared)
   ENUM_COVAR_TYPE   cov_type;               // Strategy flag applied during variance calculations (Standard vs Robust)

   //--- Default Constructor (Sets defaults for invalid metrics)
                     ArchModelResult(void)
     {
      fit_indices[0] = fit_indices[1] = 0;
      param_cov = matrix::Zeros(0,0);
      r2 = EMPTY_VALUE;
      cov_type = WRONG_VALUE;
     }

   //--- Parameterized Constructor
                     ArchModelResult(vector& _params, matrix &_param_cov, double _r2, vector& _resid, vector & _volatility, ENUM_COVAR_TYPE _cov_type, double _loglikelihood, long fit_start, long fit_stop)
     {
      loglikelihood = _loglikelihood;
      params = _params;
      param_cov = _param_cov;
      resid = _resid;
      conditional_volatility = _volatility;
      nobs = _resid.Size(); // Tracks observation sample depth internally

      fit_indices[0] = fit_start;
      fit_indices[1] = fit_stop;
      r2 = _r2;
      cov_type = _cov_type;
     }

   //--- Destructor
                    ~ArchModelResult(void)
     {
     }

   //--- Copy Constructor
                     ArchModelResult(ArchModelResult &other)
     {
      loglikelihood = other.loglikelihood;
      nobs = other.nobs;
      params = other.params;
      resid = other.resid;
      conditional_volatility = other.conditional_volatility;

      fit_indices[0] = other.fit_indices[0];
      fit_indices[1] = other.fit_indices[1];
      r2 = other.r2;
      cov_type = other.cov_type;
      param_cov = other.param_cov;
     }

   //--- Overloaded Assignment Operator
   void              operator=(ArchModelResult &other)
     {
      loglikelihood = other.loglikelihood;
      nobs = other.nobs;
      params = other.params;
      resid = other.resid;
      conditional_volatility = other.conditional_volatility;

      fit_indices[0] = other.fit_indices[0];
      fit_indices[1] = other.fit_indices[1];
      r2 = other.r2;
      cov_type = other.cov_type;
      param_cov = other.param_cov;
     }

   //--- Dynamically evaluates or updates parameter covariance states using standard or robust methods
   matrix            parameter_covariance(CArchModel &archmodel)
     {
      vector fitted_params = archmodel.get_params();
      if(cov_type == COVAR_ROBUST)
         param_cov = archmodel.compute_param_cov(fitted_params,EMPTY_VECTOR); // Triggers QMLE Sandwich Estimation
      else
         param_cov = archmodel.compute_param_cov(fitted_params,EMPTY_VECTOR,false); // Standard Hessian Inverse
      return param_cov;
     }

   //--- Computes asymptotic symmetric confidence intervals based on Inverse Normal critical thresholds
   matrix            conf_int(double alpha = 0.05)
     {
      // Calculate two-tailed critical value Z_{1 - \alpha/2}
      double cv = CNormalDistr::InvNormalCDF(1.0-alpha/2.0);
      vector se = std_err();
      matrix out = matrix::Zeros(params.Size(),2);

      out.Col(params - cv * se, 0); // Lower Bound: Estimate - (Z * SE)
      out.Col(params + cv * se, 1); // Upper Bound: Estimate + (Z * SE)

      return out;
     }

   //--- Adjusted R-squared: factors parameter-to-sample loss penalty adjustments
   double            rsquared_adj(void)
     {
      return 1.0 - ((1.0-r2)*double(nobs-1)/double(nobs-num_params()));
     }

   //--- Isolates individual parameter standard errors by extracting the square root of the covariance diagonal
   vector            std_err(void)
     {
      return sqrt(param_cov.Diag());
     }

   //--- Computes two-tailed survival probability p-values under the Standard Normal curve
   vector            pvalues(void)
     {
      vector pvals = tvalues();
      for(ulong i = 0; i<pvals.Size(); ++i)
         // Mirror distributions for lower tail and double for a two-tailed probability assessment
         pvals[i] = CNormalDistr::NormalCDF(-1.0*MathAbs(pvals[i]))*2.0;

      return pvals;
     }

   //--- Calculates asymptotic t-statistics (Wald ratios) testing the null hypothesis H_0: \beta_i = 0
   vector            tvalues(void)
     {
      return params/std_err();
     }
  };
//+------------------------------------------------------------------+
//| The Lagrange multiplier (LM) test for ARCH effects               |
//| Computes Engle's LM test statistic ($T \times R^2$) to test for  |
//| conditional heteroskedasticity in series observations/residuals. |
//+------------------------------------------------------------------+
WaldTestStatistic archlmtest(vector& residuals, vector& conditional_volatility, ulong lags, bool standardized = false)
  {
   WaldTestStatistic out;
   vector resids = residuals;

//--- Standardize residuals if requested ($z_t = \epsilon_t / \sigma_t$)
   if(standardized)
      resids = resids/conditional_volatility;

//--- Data Integrity: Scrub and filter out missing values (NaNs)
   if(resids.HasNan())
     {
      vector nresids = vector::Zeros(resids.Size() - resids.HasNan());
      for(ulong i = 0, k = 0; i<resids.Size(); ++i)
         if(MathClassify(resids[i]) == FP_NAN)
            continue;
         else
            nresids[k++] = resids[i];
      resids = nresids;
     }

   int nobs = (int)resids.Size();
   vector resid2 = MathPow(resids,2.0); // Convert series into squared innovations ($\epsilon_t^2$)

//--- Dynamic Lag Selection: Default to Schwert rule-of-thumb if lags = 0
   if(!lags)
      lags = ulong(ceil(12.0*pow(nobs/100.0,1.0/4.0)));

//--- Ensure lags fit within mathematical bounds relative to sample size
   lags = MathMax(MathMin(resids.Size()/2 - 1, lags), 1);

//--- Validation: Ensure adequate sample size exists for secondary regression
   if(resid2.Size()<3)
     {
      Print(__FUNCTION__, " Test requires at least 3 non-nan observations ");
      return out;
     }

//--- Reshape squared series into a 1-column matrix for matrix lag operations
   matrix matres = matrix::Zeros(resid2.Size(),1);
   matres.Col(resid2,0);
   matrix lag[];

//--- Construct the design lag matrices via lagmat (TRIM_BOTH drops pre-sample rows)
   if(!lagmat(matres,lag,lags,TRIM_BOTH,ORIGINAL_SEP))
     {
      Print(__FUNCTION__, " lagmat failed ");
      return out;
     }

//--- Append an intercept/constant column to the independent regressor matrix (lag[0])
   matrix lagout;
   if(!addtrend(lag[0],lagout,TREND_CONST_ONLY,true,HAS_CONST_SKIP))
     {
      Print(__FUNCTION__, " addtrend failed ");
      return out;
     }
   lag[0] = lagout;

//--- Fit the auxiliary OLS regression: $\epsilon_t^2 = \alpha_0 + \alpha_1 \epsilon_{t-1}^2 + ... + \alpha_p \epsilon_{t-p}^2 + v_t$
   OLS ols;
   if(!ols.Fit(lag[1].Col(0),lag[0]))
     {
      Print(__FUNCTION__, " OLS fitting failed ");
      return out;
     }

//--- Calculate Engle's LM test statistic: Sample Size ($T$) $\times$ R-squared ($R^2$)
   out.stat = nobs*ols.Rsqe();

//--- Document testing hypothesis strings
   if(standardized)
     {
      out.null = "Standardized residuals are homoskedastic";
      out.alternative = "Standardized residuals are conditionally heteroskedastic";
     }
   else
     {
      out.null = "Residuals are homoskedastic";
      out.alternative = "Residuals are conditionally heteroskedastic";
     }

//--- Save degrees of freedom (corresponds asymptotically to a $\chi^2$ distribution with $p$ lags)
   out.df = lags;
   return out;
  }
//+------------------------------------------------------------------+
//| Forecast helper function                                         |
//| Pads the front of a forecast matrix array with NaN entries to    |
//| align historical timeline intervals or pre-sample periods.       |
//+------------------------------------------------------------------+
bool _forecast_pad(int count, matrix &forecasts[], matrix &out[])
  {
//--- Expand target output array to hold both padding frames and raw forecasts
   ArrayResize(out, count + (int)forecasts.Size());

//--- Construct a structural matrix populated cleanly with NaN values
   matrix fill = matrix::Zeros(forecasts[0].Rows(), forecasts[0].Cols());
   fill.Fill(double("nan"));

//--- Phase 1: Inject empty padding matrices up to the specified boundary offset
   for(int i = 0; i < count; ++i)
      out[i] = fill;

//--- Phase 2: Copy raw forecast entries into the remaining indexed positions
   for(int i = count; i < int(out.Size()); ++i)
      out[i] = forecasts[i - count];

   return true;
  }
//+------------------------------------------------------------------+
//| Forecasting helper function (Matrix Overload)                    |
//| Prepend empty rows (NaN) to a matrix to align structural lengths  |
//| across historical timeline matrices.                             |
//+------------------------------------------------------------------+
matrix _forecast_pad(int count, matrix &forecasts)
  {
//--- Initialize a destination matrix expanded by the offset row count
   matrix fill = matrix::Zeros(count+forecasts.Rows(), forecasts.Cols());
   fill.Fill(double("nan")); // Default fill everything with NaN

//--- Shift and copy original matrix rows sequentially into the padded matrix
   for(ulong i = (ulong)count; i<fill.Rows(); ++i)
      fill.Row(forecasts.Row(ulong(i-count)),i);

   return fill;
  }

//+------------------------------------------------------------------+
//| Forecasting helper function (Vector Overload)                    |
//| Prepend empty values (NaN) to a vector to match data streams      |
//| and account for pre-sample estimation gaps.                       |
//+------------------------------------------------------------------+
vector _forecast_pad(int count, vector &forecasts)
  {
//--- Initialize a destination vector expanded by the offset size count
   vector fill = vector::Zeros(count+forecasts.Size());
   fill.Fill(double("nan")); // Default fill everything with NaN

//--- Shift and copy original vector elements into the padded vector positions
   for(ulong i = (ulong)count; i<fill.Size(); ++i)
      fill[i] = forecasts[ulong(i-count)];

   return fill;
  }
//+------------------------------------------------------------------+
//| Generate mean forecasts from an AR-X model                      |
//| Operates recursively over an out-of-sample forecast horizon by    |
//| rolling autoregressive lags forward and adding exogenous impacts. |
//+------------------------------------------------------------------+
matrix _ar_forecast(vector& y, ulong horizon, ulong start_index, double constant, vector & arp, matrix &x[], vector &exogp)
  {
   ulong t = y.Size();
   ulong p = arp.Size();
   ulong first,last;
   vector col;

//--- Initialize destination forecast canvas matrix
//--- Columns span the historic autoregressive lags [0...p-1] plus the forward horizon [p...p+horizon-1]
   matrix fcasts = matrix::Zeros(t - start_index,p+horizon);

//--- Seed the first 'p' columns of the matrix with actual historical lag observations
   for(ulong i = 0; i<p; ++i)
     {
      first = start_index - p + i + 1;
      last = t - p + i + 1;
      col = np::sliceVector(y,long(first),long(last));
      fcasts.Col(col,i);
     }

//--- Reverse AR parameters to naturally match chronological history via dot-product matrix multiplications
   vector arp_rev  = arp;
   np::reverseVector(arp_rev);
   matrix mlt;
   vector mat;

//--- Main loop tracking forward time-steps sequentially over the forecast horizon
   for(ulong i = p; i<horizon+p; ++i)
     {
      //--- Step 1: Extract the rolling sliding window containing the preceding 'p' lags (mix of historical and generated forecasts)
      mlt = ((i-p)!=i)?np::sliceMatrixCols(fcasts,long(i-p),long(i)):matrix::Zeros(0,0);

      //--- Step 2: Calculate the standard Autoregressive linear combination ($AR(p) = \phi_1 y_{t-1} + ... + \phi_p y_{t-p}$)
      mat = (arp_rev.Size()&&mlt.Cols())?mlt.MatMul(arp_rev):vector::Zeros(fcasts.Rows());

      //--- Step 3: Accumulate structural intercept constants into the conditional mean column path
      fcasts.Col(constant+mat,i);

      //--- Step 4: Map exogenous features ($\beta \cdot X_t$) if external matrix trackers are explicitly provided
      if(x.Size())
        {
         mlt = matrix::Zeros(x[0].Rows(),x.Size());
         for(uint k = 0; k<x.Size(); ++k)
            mlt.Col(x[k].Col(i-p),k); // Col(i-p) isolates the precise forward index for feature step $h$

         //--- Append the combined dot-product exogenous additions back into current forecast matrix coordinate arrays
         col = fcasts.Col(i);
         col += mlt.MatMul(exogp);
         fcasts.Col(col,i);
        }
     }

//--- Clean-up and trim away historical pre-sample columns, leaving only pure calculated out-of-sample forecasts
   fcasts = np::sliceMatrixCols(fcasts,long(p));
   return fcasts;
  }

//+------------------------------------------------------------------+
//| Convert AR coefficients to an Impulse Response Function (IRF)    |
//| Recursively traces the long-term impact of a single structural    |
//| shock unit ($\epsilon_0 = 1$) forward through time step horizons. |
//+------------------------------------------------------------------+
vector _ar_to_impulse(ulong steps, vector& params)
  {
   ulong p = params.Size();
   vector impulse = vector::Zeros(steps);

//--- Seed the shock at time zero ($t_0 = 1$)
   impulse[0] = 1.0;

//--- If there are no AR parameters, the impact is purely transitory (white noise)
   if(p == 0)
      return impulse;

   long k, st;
   vector rparams, imp;

//--- Recursively project the shock propagation across subsequent time intervals
   for(long i = 1; i<long(steps); ++i)
     {
      //--- Dynamically bound indices to avoid reading beyond available AR parameters or past history
      k = (long)fmin(long(p) - 1, i - 1);
      st = (long)fmax(i - long(p), 0);

      //--- Isolate historical impulse factors ($y_{t-1}, y_{t-2}, ...$)
      imp  = np::sliceVector(impulse, st, i);

      //--- Reverse the applicable slice of AR parameters to align correctly chronologically
      rparams = np::sliceVector(params, k, END_REVERSE, STEP_REVERSE);

      //--- Compute the dot product to dynamically determine the current step's decay or amplification
      if(imp.Size()==rparams.Size() && imp.Size())
         impulse[i] = imp.Dot(rparams);
     }

   return impulse;
  }
//+------------------------------------------------------------------+
//| Heterogeneous Autoregression (HAR)                               |
//+------------------------------------------------------------------+
class HARX: public CArchModel
  {
protected:
#ifdef __SLSQP__
   //+------------------------------------------------------------------+
   //| Objective function for SLSQP optimizer                           |
   //| Extends CFunctor to provide a standard interface for evaluating   |
   //| the log-likelihood objective value and its respective gradients. |
   //+------------------------------------------------------------------+
   class CObjectiveF: public CFunctor
     {
   private:
      vector            m_sigma, m_bc;       // Variance components and parameter constraint vectors
      matrix            m_vb;                // Basis matrix/regressor blocks for variance computations
      HARX*          m_obj;               // Pointer to the underlying HARX model instance
      bool              m_centered, m_individual; // Flags controlling structural type of gradient calculation
      double            m_epsilon;           // Step size (perturbation value) used for numerical differentiation
      vector            obj_result;          // Cached tracking vector containing objective evaluation scores
      matrix            obj_gradient;        // Cached tracking matrix containing the calculated Jacobian/Gradient row blocks
   public:
      //--- Constructor
                     CObjectiveF(void) {}

      //--- Destructor (Handles explicit resource resets and safe pointer detachment)
                    ~CObjectiveF(void)
        {
         m_obj = NULL;
         obj_result = vector::Zeros(1);
         obj_gradient = matrix::Zeros(1,2);
        }

      //--- Sets internal parameter constants of the objective function
      void              setFuncParams(vector& s, vector& bc, matrix &vb, HARX &obj, bool individual=false, bool center = false, double ep=EMPTY_VALUE)
        {
         m_sigma  = s;
         m_bc = bc;
         m_vb = vb;
         m_obj = GetPointer(obj); // Safe encapsulation using structural runtime pointers
         m_individual = individual;
         m_centered = center;
         m_epsilon = ep;
         m_clip = true;                               // Activates parameter clipping boundaries to prevent parameter explosion
         m_grad_options.method = GRAD_POINT_CALLABLE; // Tells the solver to utilize our explicit analytical/numerical override
        }

      //--- Objective function wrapper: Returns scalar evaluation ($f(x)$) for optimizer line search steps
      virtual double    orig_fun(vector& x) override
        {
         obj_result = m_obj.objective(x,m_sigma,m_bc,m_vb);
         return obj_result[0];
        }

      //--- Gradient of objective function wrapper: Extracts and returns the local gradient vector ($\nabla f(x)$)
      virtual vector grad_fun(vector& x) override
        {
         obj_gradient = m_obj.jacobian(x,m_sigma,m_bc,m_vb,m_individual,m_centered,m_epsilon);
         return obj_gradient.Row(0); // Isolates the active row mapping the target scalar parameter bounds
        }
     };

   //+------------------------------------------------------------------+
   //| Inequality constraints class                                     |
   //| Formalizes boundary conditions structured as: $b - A \cdot x \ge 0$|
   //+------------------------------------------------------------------+
   class CIneqConstraints: public CConstraints
     {
   private:
      matrix            a;                   // Linear constraint coefficient matrix ($A$)
      vector            b;                   // Boundary threshold vector constants ($b$)
      vector            obj_result;          // Cached tracking array mapping constraint slack values
   public:
      //--- Constructor
                     CIneqConstraints(void)
        {
        }

      //--- Destructor
                    ~CIneqConstraints(void)
        {
        }

      //--- Registers the linear parameters defining the system boundaries
      void              setConstraints(matrix& a_, vector& b_)
        {
         a = a_;
         b = b_;
         m_m = int(a.Rows());                 // Tracks total count of structural constraint equations
         m_n = int(a.Cols());                 // Tracks dimensional parameter inputs expected
         obj_result =  vector::Zeros(m_n);
         m_options.method = GRAD_POINT_2;     // Configures numerical finite-differencing strategy for boundaries
        }

      //--- Evaluates constraint spaces ($c(x) = b - A \cdot x$); positive values represent valid spaces
      virtual vector    orig_fun(vector& x) override
        {
         if(!a.Rows())
            obj_result = vector::Zeros(0);    // Handles unconstrained or open boundary executions safely
         else
            obj_result = b - a.MatMul(x);     // Computes the linear coordinate residuals against boundary thresholds
         return obj_result;
        }
     };
#else
   //+------------------------------------------------------------------+
   //| Objective function adapter for Alglib Optimizer                  |
   //| Inherits from CNDimensional_Func to map multi-dimensional double  |
   //| arrays/rows to the underlying HARX model's math vector objective. |
   //+------------------------------------------------------------------+
   class CObjectiveFunc: public CNDimensional_Func
     {
   private:
      vector            m_sigma, m_bc;       // Variance properties and optimization boundary limits
      matrix            m_vb;                // Linear basis data partitions/regressors
      HARX*          m_obj;               // Structural tracking link pointing back to the core HARX model
      bool              m_centered, m_individual; // Spatial properties configuration flags for differentiation
      double            m_epsilon;           // Perturbation constraint index used during gradient calculations
   public:
      //--- Constructor
                     CObjectiveFunc(void) {}

      //--- Destructor (Safely detaches the internal model pointer reference)
                    ~CObjectiveFunc(void)
        {
         m_obj = NULL;
        }

      //--- Configures parameter constants required by the target objective calculation system
      void              SetParams(vector& s, vector& bc, matrix &vb, HARX &obj, bool individual=false, bool center = false, double ep=EMPTY_VALUE)
        {
         m_sigma  = s;
         m_bc = bc;
         m_vb = vb;
         m_obj = GetPointer(obj);
         m_individual = individual;
         m_centered = center;
         m_epsilon = ep;
        }

      //--- Evaluates objective function via raw double array conversion interfaces
      virtual void      Func(double &x[], double &func, CObject &obj) override
        {
         vector temp;
         temp.Assign(x); // Convert native Alglib parameter array layout into a native math vector
         vector r = m_obj.objective(temp, m_sigma, m_bc, m_vb);
         func = r[0];    // Isolates the primary log-likelihood target scalar metric
        }

      //--- Overloaded alternative evaluating objective space via Alglib's CRowDouble format
      virtual void      Func(CRowDouble &x, double &func, CObject &obj) override
        {
         vector temp = x.ToVector(); // Safe translation layer from specialized Alglib types to math vector profiles
         vector r = m_obj.objective(temp, m_sigma, m_bc, m_vb);
         func = r[0];
        }
     };

   //+------------------------------------------------------------------+
   //| Jacobian fvec vector implementation for Alglib                   |
   //| Inherits from CNDimensional_Jac to handle simultaneous evaluation |
   //| of the residual vector ($f_{vec}$) and its first-order derivative matrix. |
   //+------------------------------------------------------------------+
   class CObjectiveJacFVec: public CNDimensional_Jac
     {
   private:
      vector            m_sigma, m_bc;       // Parameter limits and target variance matrices
      matrix            m_vb;                // Basis matrix array partitions
      HARX*          m_obj;               // Structural pointer pointing to the active execution engine
      bool              m_centered, m_individual; // Gradient computation alignment strategies
      double            m_epsilon;           // Scaling value used for numeric parameter mutations
   public:
      //--- Constructor
                     CObjectiveJacFVec(void) {}

      //--- Destructor
                    ~CObjectiveJacFVec(void)
        {
         m_obj = NULL;
        }

      //--- Binds optimization parameters to internal tracker fields
      void              SetParams(vector& s, vector& bc, matrix &vb, HARX &obj, bool individual=false, bool center = false, double ep=EMPTY_VALUE)
        {
         m_sigma  = s;
         m_bc = bc;
         m_vb = vb;
         m_obj = GetPointer(obj);
         m_individual = individual;
         m_centered = center;
         m_epsilon = ep;
        }

      //--- Evaluates both individual residual rows ($f_i$) and the Jacobian matrix simultaneously using native arrays
      virtual void      Jac(double &x[], double &fi[], CMatrixDouble &jac, CObject &obj) override
        {
         vector temp;
         temp.Assign(x);

         //--- Evaluate objective baseline at the current spatial parameters coordinate
         vector r = m_obj.objective(temp, m_sigma, m_bc, m_vb);
         fi[0] = r[0];

         //--- Calculate the system Jacobian matrix ($\mathbf{J}_{i,j} = \frac{\partial f_i}{\partial x_j}$)
         matrix mjac = m_obj.jacobian(temp, m_sigma, m_bc, m_vb, m_individual, m_centered, m_epsilon);
         jac = mjac; // Implicit casting assignment to Alglib's internal CMatrixDouble format
        }

      //--- Overloaded variant mapping spatial dimensions directly from Alglib CRowDouble components
      virtual void      Jac(CRowDouble &x, CRowDouble &fi, CMatrixDouble &jac, CObject &obj) override
        {
         vector xv = x.ToVector();

         //--- Evaluate the multidimensional residual mapping vector block
         vector temp = m_obj.objective(xv, m_sigma, m_bc, m_vb);
         fi = temp; // Directly overwrite the multi-row output reference structure

         //--- Compute the matrix partial derivatives layout block
         matrix mjac = m_obj.jacobian(xv, m_sigma, m_bc, m_vb, m_individual, m_centered, m_epsilon);
         jac = mjac;
        }
     };
#endif
   int               m_extra_simulation_params;
   //+------------------------------------------------------------------+
   //| Calculate the Coefficient of Determination (R-squared)           |
   //| Computes $R^2 = 1 - \frac{SSR}{TSS}$ to determine the percentage |
   //| of explained variance captured by the structural model.           |
   //+------------------------------------------------------------------+
   virtual double    _r2(vector &params) override
     {
      vector y = m_fit_y;
      bool constant = false;
      matrix x = m_fit_regressors;

      //--- Check for an explicit or implicit structural model intercept constant
      if(x.Cols())
         constant = (m_constant||implicit_constant(x));

      //--- Center the dependent vector around its mean if a constant exists
      if(constant)
        {
         //--- If the matrix is intercept-only, structural variance cannot be explained ($R^2 = 0$)
         if(x.Cols() ==1)
            return 0.0;
         y = y - y.Mean(); // Transforms the vector to compute centered Sum of Squares
        }

      //--- Calculate Total Sum of Squares ($TSS = \sum (y_i - \bar{y})^2$) via vector dot-product
      double tss = y.Dot(y);
      if(tss<=0.)
         return AL_NaN; // Guard against division by zero errors on flat data streams

      //--- Extract model tracking errors / residuals ($\epsilon$)
      vector e = resids(params,EMPTY_VECTOR,EMPTY_MATRIX);

      //--- Return final metric: $R^2 = 1 - \frac{Residual\ Sum\ of\ Squares\ (e \cdot e)}{Total\ Sum\ of\ Squares\ (TSS)}$
      return 1.-(e.Dot(e))/tss;
     }

   //+------------------------------------------------------------------+
   //| Dynamic sample framing array adjuster                            |
   //| Re-calibrates active training timeline windows, accounting for    |
   //| designated lag pre-sample allocations or back-testing holdbacks. |
   //+------------------------------------------------------------------+
   bool              _adjust_sample(long first_obs, long last_obs)
     {
      //--- Offset start boundary forward to shift past the warm-up/holdback lag horizon
      first_obs+=long(m_holdback);
      if(last_obs<= first_obs)
        {
         Print(__FUNCTION__, " invalid indices specified ");
         return false;
        }

      //--- Cache target data window indices inside the fit index tracking vector
      vector f = {double(first_obs),double(last_obs)};
      m_fit_indices = f;

      //--- Slice and align the dependent and independent series to the active estimation block
      m_fit_y = adjust_y();
      m_fit_x = adjust_x();

      //--- Extract active explanatory feature rows (regressors) from the master collection matrix
      m_fit_regressors = (m_regressors.Rows())?np::sliceMatrixRows(m_regressors,long(m_fit_indices[0]),long(m_fit_indices[1])):matrix::Zeros(0,0);

      //--- Sync structural pointer/state indicators for the local model processing space
      m_vp.start(long(first_obs));
      m_vp.stop(long(last_obs));

      return true;
     }
   //+------------------------------------------------------------------+
   //| Adjust Dependent Variable Vector (y)                             |
   //| Slices the master observation stream to isolate the active        |
   //| sample timeframe mapped by the current model estimation window.  |
   //+------------------------------------------------------------------+
   vector            adjust_y(void)
     {
      //--- Extract and slice data if the master observations vector contains entries
      if(m_model_spec.observations.Size())
         return np::sliceVector(m_model_spec.observations,long(m_fit_indices[0]),long(m_fit_indices[1]));
      return vector::Zeros(0); // Return an empty vector frame if no observations exist
     }

   //+------------------------------------------------------------------+
   //| Adjust Exogenous Regressor Matrix (X)                            |
   //| Slices the master exogenous data matrix by either rows or        |
   //| columns to align tracking structures with the active sample window|
   //+------------------------------------------------------------------+
   matrix            adjust_x(bool rows = true)
     {
      //--- Proceed only if the master matrix container contains data rows
      if(m_model_spec.exog_data.Rows())
        {
         //--- Slice horizontally (default) to extract active temporal periods/rows
         if(rows)
            return np::sliceMatrixRows(m_model_spec.exog_data,long(m_fit_indices[0]),long(m_fit_indices[1]));
         //--- Otherwise slice vertically to extract specific structural features/columns
         else
            return np::sliceMatrixCols(m_model_spec.exog_data,long(m_fit_indices[0]),long(m_fit_indices[1]));
        }
      else
         return matrix::Zeros(0,0); // Return empty 0x0 matrix if no features exist
     }

   //+------------------------------------------------------------------+
   //| Compute Baseline Parameters under Homoskedastic Normal Errors     |
   //| Generates starting parameters for the conditional mean equation  |
   //| via OLS using the Moore-Penrose Pseudo-Inverse ($\mathbf{X}^+$). |
   //+------------------------------------------------------------------+
   vector            _fit_no_arch_normal_errors_params(void)
     {
      vector y = m_fit_y;
      ulong nobs = y.Size();

      //--- Data Check: Enforce degrees-of-freedom requirements ($N > K$)
      if(nobs<m_num_params)
        {
         Print(__FUNCTION__, " insufficient data ");
         return vector::Zeros(0);
        }

      matrix x = m_fit_regressors;
      if(x.Cols())
        {
         ResetLastError();

         //--- Compute Moore-Penrose Pseudo-Inverse: $\mathbf{X}^+ = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T$
         matrix xpinv = np::pinv(x);

         //--- Trapping singular matrices or calculation underflows
         if(xpinv.HasNan())
           {
            Print(__FUNCTION__, " PInv() returned invalid number ", GetLastError());
            return vector::Zeros(0);
           }

         //--- Calculate OLS coefficients: $\hat{\beta} = \mathbf{X}^+ \cdot \mathbf{y}$
         return np::matmul(xpinv,y);
        }
      else
         return vector::Zeros(0); // Return empty if there are no regressors in the system
     }
   //+------------------------------------------------------------------+
   //| Estimates baseline model parameters with homoskedastic errors   |
   //| Fits OLS parameters, extracts residuals, and builds information  |
   //| matrices for either standard (Classic) or White's Robust errors. |
   //+------------------------------------------------------------------+
   ArchModelResult   _fit_no_arch_normal_errors(ENUM_COVAR_TYPE cov_type = COVAR_ROBUST)
     {
      ArchModelResult out;
      matrix x = m_fit_regressors;
      vector y = m_fit_y;
      vector fitted, rparams;
      fitted = rparams = vector::Zeros(0);
      matrix xpxi = matrix::Zeros(0,0);

      //--- 1. Parameter Estimation via OLS
      if(x.Cols())
        {
         matrix xpinv = np::pinv(x); // Moore-Penrose Pseudo-Inverse
         rparams = np::matmul(xpinv,y); // $\hat{\beta} = \mathbf{X}^+ \cdot \mathbf{y}$
         matrix xt = (x.Transpose().MatMul(x))/double(y.Size()); // Normalized $\frac{\mathbf{X}^T\mathbf{X}}{T}$
         xpxi = xt.Inv(); // $(\mathbf{X}^T\mathbf{X})^{-1}$ tracking block
         fitted = x.MatMul(rparams); // Predicted conditional mean outputs
        }

      //--- 2. Extract Residual Errors ($\epsilon$) and Variance Vector ($\sigma^2$)
      vector e = np::minus(y,fitted);
      double sigma2 = e.Dot(e)/double(y.Size()); // Homoskedastic variance estimate
      vector params = rparams;

      //--- Append the variance estimate ($\sigma^2$) as the final structural parameter row
      rparams.Resize(params.Size()+1);
      rparams[params.Size()] = sigma2;

      //--- 3. Construct the Analytical Log-Likelihood Information Hessian Matrix
      matrix hessian = matrix::Zeros(m_num_params+1,m_num_params+1);
      np::matrixCopy(hessian,-1.*xpxi,0,long(m_num_params),1,0,long(m_num_params));
      hessian[hessian.Rows()-1,hessian.Cols()-1] = -1.; // Map scale parameter derivative threshold

      matrix param_cov = matrix::Zeros(0,0);
      string nm;

      //--- 4. Structural Covariance Type Resolution
      switch(cov_type)
        {
         //--- Standard Asymptotic Covariance Matrix Evaluation
         case COVAR_CLASSIC:
            param_cov = sigma2 *(-1.*hessian);
            param_cov[m_num_params,m_num_params] = 2.*pow(sigma2,2.); // Asymptotic variance of error variance
            param_cov/=double(y.Size()); // Scale matrix by sample depth ($T$)
            nm = "classic_ols";
            break;

         //--- White's Heteroskedasticity-Consistent (Sandwich) Covariance Estimator
         case COVAR_ROBUST:
           {
            // Initialize score matrix to hold individual structural gradients/moment conditions
            matrix scores = matrix::Zeros(y.Size(),m_num_params+1);
            for(ulong i =0; i<m_num_params; ++i)
               scores.Col(e*x.Col(i),i); // Mean equation parameter scores: $\epsilon_t \cdot \mathbf{X}_{t,i}$

            // Variance parameter score column: $\epsilon_t^2 - \sigma^2$
            scores.Col(pow(e,2.)-sigma2,scores.Cols()-1);

            // Outer product of scores matrix ($\mathbf{B} = \frac{\mathbf{G}^T\mathbf{G}}{T}$)
            matrix scorescov = scores.Transpose().MatMul(scores)/double(y.Size());

            // Core Sandwich Assembly: $\mathbf{A}^{-1} \mathbf{B} \mathbf{A}^{-1}$
            param_cov = hessian.MatMul(scorescov);
            param_cov = param_cov.MatMul(hessian);
            param_cov/=double(y.Size()); // Divide by overall sample length $T$
            nm = "white";
           }
         break;
        }

      //--- 5. Diagnostics, Vector Alignment, and Output Structure Packaging
      double r2 = _r2(params); // Compute R-squared metric

      // Prepare full-sized array allocations aligned with the complete master timeline
      vector resids_ = vector::Zeros(m_model_spec.observations.Size());
      resids_.Fill(AL_NaN); // Pad unestimated rows with NaN
      np::vectorCopy(resids_,e,long(m_fit_indices[0]),long(m_fit_indices[1])); // Splice in active window errors

      vector vol = vector::Zeros(resids_.Size());
      vol.Fill(AL_NaN);
      double ss = sqrt(sigma2); // Calculate constant standard deviation ($\sigma$)

      // Populate the volatility vector for the active estimation window slice
      for(ulong i = ulong(m_fit_indices[0]); i<ulong(m_fit_indices[1]); vol[i] = ss, ++i);

      // Compute final Gaussian Log-Likelihood objective score
      double ll = _static_gaussian_loglikelihood(e);

      // Save state variables to internal class properties
      m_params = out.params = rparams;
      m_converged = (m_params.Size()>0);
      out.param_cov = param_cov;
      out.r2 = r2;
      out.resid = resids_;
      out.conditional_volatility = vol;
      out.loglikelihood = ll;
      out.nobs = e.Size();
      out.fit_indices[0] = long(m_fit_indices[0]);
      out.fit_indices[1] = long(m_fit_indices[1]);

      return out;
     }
   //--- Should be called whenever the model is initialized or changed
   bool              _init_model(void)
     {
      ResetLastError();

      //--- 1. Process and Reformat Model Lag Specifications
      if(m_model_spec.mean_lags.Size())
        {
         if(!_reformat_lags())
            return false;

         // Determine the largest autoregressive lag index to protect pre-sample sequences
         m_maxlags = ulong(m_model_spec.mean_lags.Max());

         // Enforce a minimum data holdback sample size to prevent index-out-of-bounds errors during lag construction
         if(!m_holdback)
            m_holdback = m_maxlags;
         else
           {
            if(m_holdback<m_maxlags)
              {
               Print(__FUNCTION__, "  hold_back is less then the minimum number given the lags selected");
               m_holdback = m_maxlags;
              }
           }
        }

      //--- 2. Validate structural constraints across model options
      if(!_check_specification())
         return false;

      ulong nobs_orig = m_model_spec.observations.Size();
      vector reg_constant = vector::Zeros(0);
      matrix reg_lags = matrix::Zeros(0,0);

      //--- 3. Construct the Intercept Column vector (ones vector)
      if(m_constant)
         reg_constant = vector::Ones(nobs_orig);

      //--- 4. Build Autoregressive Lag Arrays via lagmat
      if(m_lags.Rows() && nobs_orig)
        {
         double maxlag = m_lags.Max();
         matrix ymat(m_model_spec.observations.Size(),1);
         ymat.Col(m_model_spec.observations,0);
         matrix lagarray[];

         // Build the pre-sample lag matrix pool (TRIM_FORWARD protects matrix dimensions)
         if(!lagmat(ymat,lagarray,ulong(maxlag),TRIM_FORWARD,ORIGINAL_EX))
            return false;

         // Extract and aggregate specific lag column matrices (averaging or slicing subsets)
         reg_lags = matrix::Zeros(nobs_orig,m_lags.Cols());
         for(ulong i = 0; i<m_lags.Cols(); ++i)
           {
            matrix temp = np::sliceMatrixCols(lagarray[0],long(m_lags[0,i]),long(m_lags[1,i]));
            vector vtemp = temp.Mean(1); // Compute mean across specified lag bounds if grouping is applied
            reg_lags.Col(vtemp,i);
           }
        }

      //--- 5. Initialize the Consolidated Design Matrix Container (m_regressors)
      // Total column allocations: Constant (0 or 1) + AR Lag columns + Exogenous Feature columns
      m_regressors = matrix::Zeros(nobs_orig,ulong(nobs_orig>0)*(ulong(m_constant)+reg_lags.Cols()+m_model_spec.exog_data.Cols()));

      // Slot in the model intercept constant vector block at column index 0
      if(m_constant && nobs_orig)
         m_regressors.Col(reg_constant,0);

      // Copy the constructed autoregressive series variables sequentially into adjacent channels
      if(reg_lags.Cols())
         for(ulong i = 0; i<reg_lags.Cols(); ++i)
            m_regressors.Col(reg_lags.Col(i),i+ulong(m_constant));

      // Append external exogenous design elements into the remaining trailing columns
      if(m_model_spec.exog_data.Cols())
         for(ulong i = 0; i<m_model_spec.exog_data.Cols(); ++i)
            m_regressors.Col(m_model_spec.exog_data.Col(i),i+ulong(m_constant)+reg_lags.Cols());

      //--- 6. Structural error boundary check before unlocking
      if(GetLastError()!=0)
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }
      return true;
     }
   //+------------------------------------------------------------------+
   //| Scale Change Event Listener Override                             |
   //| Automatically re-initializes and updates the structural design   |
   //| matrix layouts if data scales or internal parameters are altered. |
   //+------------------------------------------------------------------+
   virtual void      _scalechanged(void) override
     {
      //--- Trigger an inline model refresh to rebuild vectors and track matrix sizing
      if(!_init_model())
         Print(__FUNCTION__, " error ");
     }

   //+------------------------------------------------------------------+
   //| Validate Model Structural Specifications                         |
   //| Performs checking boundaries on data shapes to confirm exogenous |
   //| matrix lengths perfectly align with the core observation series. |
   //+------------------------------------------------------------------+
   bool              _check_specification(void)
     {
      //--- Evaluate matrix dimensional integrity if exogenous data blocks exist
      if(m_model_spec.exog_data.Rows())
        {
         //--- Dimensional Gate: Enforce structural equality along the time axis (Rows == N)
         if(m_model_spec.exog_data.Rows()!=m_model_spec.observations.Size())
           {
            Print(__FUNCTION__, " regressors shape does not match input data ");
            return false; // Rejects initialization due to broken sample alignments
           }
        }
      return true; // Pass specifications check
     }
   //+------------------------------------------------------------------+
   //| Reformat input lags to simplify matrix operations                |
   //| Standardizes a 1-row lag vector or 2-row interval array into a   |
   //| normalized indexing matrix and validates entry uniqueness via rank.|
   //+------------------------------------------------------------------+
   bool              _reformat_lags(void)
     {
      //--- If no lags are specified, bypass processing cleanly
      if(!m_lags.Rows())
         return true;

      matrix lags = m_lags;

      //--- Case 1: Input is a simple 1-Row array of specific lag milestones
      if(lags.Rows()==1)
        {
         // Validate data inputs; lag values must be positive integers
         if(lags.Min()<=0)
           {
            Print(__FUNCTION__, " invalid lags input ");
            return false;
           }

         // Extract row, remove duplicate instances, and sort ascendingly
         vector row = lags.Row(0);
         row = np::unique(row);

         // Convert to a 2-row matrix to represent continuous intervals [Start, Stop]
         matrix temp = matrix::Zeros(2,row.Size());
         if(!temp.Row(row,0) || !temp.Row(row,1))
            return false;

         //--- Process Rotated Lags: Maps cumulative/overlapping periods (e.g., HAR-Daily/Weekly/Monthly)
         if(m_rotated)
           {
            vector v = np::sliceVector(row,0,long(row.Size()-1));
            for(ulong i = 1; i<temp.Cols(); temp[0,i] = v[i-1], ++i);
            temp[0,0] = 0.; // Sets base offset for the initial lag index
           }
         //--- Process Standard Lags: Disjoint individual time steps
         else
            temp.Row(vector::Zeros(temp.Cols()),0); // Base remains zero for direct offset parsing

         m_lags = temp;
        }
      //--- Case 2: Input is a 2-Row matrix specifying range intervals [Start, Stop]
      else
        {
         vector mins = lags.Min(1);
         // Validate that all bounds are positive and upper boundaries exceed lower bounds
         if(lags.Min()<=0 || mins[1]<mins[0])
           {
            Print(__FUNCTION__, " invalid lags input ");
            return false;
           }

         // Normalize structure directions across rows and sort columns
         if(!np::reverseMatrixRows(lags) || !np::sort(lags,true))
            return false;

         //--- Validate Multicollinearity / Redundant overlapping intervals
         // Build a binary indicator profile matrix matching the max lag width
         matrix testmat = matrix::Zeros(lags.Cols(),(ulong)lags.Max());
         vector subtract = {1.,0.};
         matrix sub = np::repeat_vector_as_rows_cols(subtract,lags.Cols(),false);
         lags = lags-sub; // Zero-index adjustment subtraction layer

         // Fill matrix coordinates representing the active ranges
         for(ulong i = 0; i<lags.Cols(); ++i)
            np::matrixFill(testmat,1.0,long(i),long(i+1),STEP,long(lags[0,i]),long(lags[1,i]));

         // Check rank condition: Number of non-zero singular vectors must equal total column dimensions
         ulong rank = np::matrix_rank(testmat);
         if(rank != lags.Cols())
           {
            Print(__FUNCTION__," lags contains redundant entries ");
            return false; // Rejects collinear/overlapping definitions
           }

         if(m_rotated)
            Print(__FUNCTION__, " Rotation is not available for this model specification ");

         m_lags = lags;
        }
      return true;
     }
   //+------------------------------------------------------------------+
   //| Reformat and validate exogenous forecast parameters             |
   //| Enforces structural integrity on a 3-D container (array of       |
   //| matrices) containing multi-step exogenous values for forecasting. |
   //+------------------------------------------------------------------+
   bool              _reformat_forecast_x(ulong start, ulong horizon, matrix& x[], matrix& out[])
     {
      //--- 1. Verification Gate: Confirm presence/absence of exogenous structures aligns with specification
      if(!x.Size())
        {
         if(!m_model_spec.exog_data.Rows())
            return true; // No exogenous data expected or provided; bypass cleanly
         else
           {
            Print(__FUNCTION__, " error x function variable should not be empty since model specified with exogenous variables ");
            return false;
           }
        }
      else
         if(!m_model_spec.exog_data.Rows())
           {
            Print(__FUNCTION__, " x is not None but the model does not contain any exogenous variables.");
            return false;
           }

      //--- 2. Dimensional Gates: Check feature space and step horizon matches
      ulong nx = m_model_spec.exog_data.Cols();
      if(x.Size()!=uint(nx))
        {
         Print(__FUNCTION__, " the number of elements in x should match number of columns of m_model_spec.exog_data property ");
         return false;
        }
      if(x[0].Cols()!=horizon)
        {
         Print(__FUNCTION__, " the number of values passed does not match the horizon of the forecasts ");
         return false;
        }

      //--- 3. Validation: Verify active row depth conforms to the forecasting window index space
      if(x[0].Rows()!=(m_model_spec.observations.Size()-start))
        {
         Print(__FUNCTION__," the shape of x does not have the proper dimensions ");
         return false;
        }

      //--- 4. Memory Preparation and Row-Slicing Alignment
      ArrayResize(out,x.Size());

      // If the input matrices contain extra historical rows, slice starting from the active 'start' forecast index
      if(x[0].Rows()>(m_model_spec.observations.Size()-start))
        {
         for(uint i = 0; i<x.Size(); ++i)
            out[i] = np::sliceMatrixRows(x[i],long(start));
        }
      // Otherwise, the arrays are already perfectly trimmed; perform a direct copy swap
      else
        {
         for(uint i = 0; i<x.Size(); ++i)
            out[i] = x[i];
        }

      return true;
     }
   //+------------------------------------------------------------------+
   //| Model transformation: HAR parameters to standard AR parameters   |
   //| Distributes cumulative HAR lag weights (e.g., weekly or monthly   |
   //| averages) uniformly across standard daily AR lag structures.     |
   //+------------------------------------------------------------------+
   vector            _har_to_ar(vector &params)
     {
      //--- If no lags exist, isolate and return the intercept constant parameter block
      if(!m_maxlags)
         return np::sliceVector(params,0,long(m_constant));

      //--- Strip away the intercept (if present) to isolate the raw HAR coefficients
      vector har = np::sliceVector(params,long(m_constant));
      vector ar = vector::Zeros(m_maxlags);
      vector temp;

      //--- 1. Linear Mapping: Disentangle intervals and distribute weights
      for(ulong i = 0; i<m_lags.Cols(); ++i)
        {
         // Distribute the aggregated HAR weight uniformly across its interval span ($W_j = \phi_{HAR} / \Delta t$)
         for(ulong j = ulong(m_lags[0,i]); j<ulong(m_lags[1,i]); ++j)
            ar[j] += har[i]/(m_lags[1,i] - m_lags[0,i]);
        }

      //--- 2. Structure Reassembly: Re-attach the model intercept constant
      if(m_constant)
        {
         vector a(ar.Size()+1);
         a[0] = params[0]; // Set the intercept constant at position index 0

         // Shift original standard AR lag values back into chronological vector channels
         for(ulong i = 1; i<a.Size(); a[i] = ar[i-1], ++i);
         ar = a;
        }

      return ar;
     }
   //+------------------------------------------------------------------+
   //| Runs a recursive simulation loop to calculate the model mean     |
   //| Generates a structural time-series sequence by combining path     |
   //| initializations, lag dependencies, exogenous variables, and errors.|
   //+------------------------------------------------------------------+
   vector            _simulate_mean(vector &parameters, matrix &x, vector &errs, vector& initial_values, vector& cond_var)
     {
      ulong max_lag = (!m_lags.Rows())?0:(ulong)m_lags.Max();
      ulong nobs_and_burn = errs.Size();
      vector y =  vector::Zeros(nobs_and_burn);
      vector iv;
      ulong iv_size = initial_values.Size();

      //--- 1. Condition and Seed Pre-Sample Initialization Arrays
      if(iv_size && iv_size!=max_lag)
        {
         Print(__FUNCTION__, " initial_value has the wrong shape ");
         return vector::Zeros(0);
        }
      else
         if(!iv_size)
            iv.Resize(max_lag); // If omitted, defaults to an array of zeros
         else
            iv = initial_values;

      // Warm up the beginning of the time-series array with specified seed histories
      for(ulong i = 0; i<max_lag; y[i] = iv[i],++i);
      ulong kx = x.Cols();

      ulong ind;
      vector temp, col;

      //--- 2. Main Chronological Simulation Iteration Loop
      for(ulong t = max_lag; t<nobs_and_burn; ++t)
        {
         ind = 0; // Tracks parameter index mapping across structural boundaries

         // Phase A: Apply Intercept Model Constants ($\mu$)
         if(m_constant)
           {
            y[t] = parameters[ind];
            ind+=1;
           }

         // Phase B: Compute Autoregressive / HAR Interval Lag Components
         for(ulong lag = 0; lag<m_lags.Cols(); ++lag)
           {
            col = m_lags.Col(lag); // Extract the interval bounds [Start, Stop]

            // Slice the active window tracking period chronologically relative to step 't'
            temp  = np::sliceVector(y,long(t - (ulong)col[1]), long(t - ulong(col[0])));

            // Apply structural coefficient against the slice mean (supports standard and grouped/HAR components)
            y[t] += parameters[ind] * temp.Mean();
            ind+=1;
           }

         // Phase C: Incorporate Exogenous Regressor Linear Metrics ($\sum \beta_i X_{t,i}$)
         for(ulong i = 0; i<kx; ++i)
           {
            y[t] += parameters[ind] * x[t,i];
            ind += 1;
           }

         // Phase D: Inject Random/Conditional Error Innovations ($\epsilon_t$)
         y[t] += errs[t];
        }

      return y;
     }
   //+------------------------------------------------------------------+
   //| Fit model without parameters                                     |
   //| Bypasses the optimizer when zero parameters require estimations.  |
   //| Directly extracts raw residuals and profiles the variance series.|
   //+------------------------------------------------------------------+
   ArchModelResult   _fit_parameterless_model(ENUM_COVAR_TYPE cv, vector& backcast)
     {
      ArchModelResult out;

      vector y = m_fit_y;
      vector params = vector::Zeros(0);      // No parameters to optimize
      matrix param_cov = matrix::Zeros(0,0);  // Covariance matrix is empty
      ulong first_obs = (ulong)m_fit_indices[0];
      ulong last_obs = (ulong)m_fit_indices[1];

      //--- 1. Isolate and Extract Model Residuals
      vector resids_final = vector::Zeros(y.Size());
      if(!np::vectorCopy(resids_final, y, first_obs, last_obs))
        {
         Print(__FUNCTION__,  " failed vector copy ");
         return out;
        }
      
      //--- 2. Compute Target Variance Boundaries and Structural Conditions
      matrix var_bounds = m_vp.varianceBounds(y);
      vector vol = vector::Zeros(y.Size());
      m_vp.computeVariance(params, y, vol, backcast, var_bounds);
      // Transform raw variance values into volatility standard deviations ($\sigma_t = \sqrt{\sigma_t^2}$)
      vol = sqrt(vol);

      //--- 3. Slice and Align Volatility Elements inside the Timeline Window
      vector vol_final = vector::Zeros(y.Size());
      vol_final.Fill(AL_NaN); // Initialize array with NaN to pad unestimated spaces
      if(!np::vectorCopy(vol_final, vol, first_obs, last_obs))
        {
         Print(__FUNCTION__,  " failed vector copy ");
         return out;
        }

      //--- 4. Calculate Summary Diagnostics and Core Model Log-Likelihood
      double r2 = _r2(params);
      vector ll = -1.0 * _loglikelihood(params, pow(vol, 2) * vector::Ones(ulong(m_fit_indices[1] - m_fit_indices[0])), backcast, var_bounds);

      //--- 5. Package and Deliver Output Result Structures
      out.params = params;
      out.param_cov = param_cov;
      out.r2 = r2;
      out.resid = resids_final;
      out.conditional_volatility = vol_final;
      out.cov_type = cv;
      out.loglikelihood = ll[0];
      out.fit_indices[0] = long(m_fit_indices[0]);
      out.fit_indices[1] = long(m_fit_indices[1]);

      return out;
     }
   //+------------------------------------------------------------------+
   //| Hidden class initialization routine                              |
   //| Configures the complete execution framework by validating the    |
   //| mean specification, instantiating the chosen volatility process, |
   //| and setting up the target innovation error distribution model.   |
   //+------------------------------------------------------------------+
   bool              _initialize(ArchParameters &vol_dist_params)
     {
      //--- 1. Verification Gate: Enforce structural class and specification alignment
      if(m_model_spec.mean_model_type!=WRONG_VALUE && vol_dist_params.mean_model_type!=WRONG_VALUE && vol_dist_params.mean_model_type!=m_model_spec.mean_model_type)
        {
         m_initialized = false;
         Print(__FUNCTION__" Incorrect initialization of mean model. \nMake sure input parameters correspond with the correct class ");
         return false;
        }

      // Sync the specifications type if the target parameters structure is uninitialized
      if(m_model_spec.mean_model_type!=WRONG_VALUE && vol_dist_params.mean_model_type==WRONG_VALUE)
         vol_dist_params.mean_model_type = m_model_spec.mean_model_type;

      // Bind the incoming parameter structure to the internal master specification state
      m_model_spec = vol_dist_params;

      //--- 2. Mean Model Configuration State Machine
      switch(m_model_spec.mean_model_type)
        {
         case MEAN_CONSTANT:
            m_model_spec.mean_lags = vector::Zeros(0);
            m_name = "Constant Mean";
            m_model_spec.include_constant = true;
            m_model_spec.use_har_rotation = false;
            break;
         case MEAN_AR:
            m_model_spec.use_har_rotation = false;
            m_model_spec.exog_data = matrix::Zeros(0,0);
            m_name = "AR";
            break;
         case MEAN_ARX:
            m_model_spec.use_har_rotation = false;
            m_name ="AR-X";
            break;
         case MEAN_HAR:
            m_model_spec.exog_data = matrix::Zeros(0,0);
            m_name = "HAR";
            break;
         case MEAN_HARX:
            m_name = "HAR-X";
            break;
         case MEAN_ZERO:
            m_model_spec.exog_data=matrix::Zeros(0,0);
            m_model_spec.mean_lags = vector::Zeros(0);
            m_model_spec.include_constant = false;
            m_model_spec.use_har_rotation = false;
            m_name = "Zero Mean";
            break;
         default:
            m_name = "HAR-X";
            break;
        }

      //--- 3. Set Timeline Estimation Boundaries and Format Structural Lags
      m_fit_indices = vector::Zeros(2);
      m_fit_indices[1] = (m_model_spec.observations.Size())?double(m_model_spec.observations.Size()):-1.;

      if(m_model_spec.mean_lags.Size())
        {
         switch(m_model_spec.mean_model_type)
           {
            // Standard Autoregressive: Format structural arrays for point-to-point lag entries
            case MEAN_AR:
            case MEAN_ARX:
               m_lags = matrix::Zeros(2,m_model_spec.mean_lags.Size());
               m_lags.Row(m_model_spec.mean_lags,0);
               m_lags.Row(m_model_spec.mean_lags,1);
               break;
            // Heterogeneous Autoregressive: Initialize arrays to track multi-day windows
            case MEAN_HAR:
            case MEAN_HARX:
               m_lags = matrix::Zeros(1,m_model_spec.mean_lags.Size());
               m_lags.Row(m_model_spec.mean_lags,0);
               break;
           }
        }
      else
         m_lags = matrix::Zeros(0,0);

      // Sync local runtime toggles with the verified specification snapshot
      m_constant = m_model_spec.include_constant;
      m_rescale = m_model_spec.is_rescale_enabled;
      m_rotated = m_model_spec.use_har_rotation;
      m_holdback = m_model_spec.holdout_size;

      // Initialize core matrix calculations and allocate design variables
      m_initialized = _init_model();
      if(!m_initialized)
         return false;

      m_num_params = num_params();

      //--- 4. Dynamic Polymorphic Volatility Process Object Allocation
      if(CheckPointer(m_vp)==POINTER_DYNAMIC)
         delete m_vp; // Prevent memory leaks by freeing old dynamic instances
      m_vp = NULL;

      switch(vol_dist_params.vol_model_type)
        {
         case VOL_CONST:
            m_vp = new CConstantVariance(m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_ARCH:
            m_vp = new CArchProcess(m_model_spec.garch_p,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_GARCH:
            m_vp = new CGarchProcess(m_model_spec.garch_p,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_AVARCH:
            m_vp = new CAvarchProcess(m_model_spec.garch_p,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_AVGARCH:
            m_vp = new CAvgarchProcess(m_model_spec.garch_p,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_TARCH:
            m_vp = new CTarchProcess(m_model_spec.garch_p,m_model_spec.garch_o,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_GJR_GARCH:
            m_vp = new CGjrGarchProcess(m_model_spec.garch_p,m_model_spec.garch_o,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_HARCH:
            m_vp = new CHarchProcess(m_model_spec.harch_lags);
            break;
         case VOL_FIGARCH:
            m_vp = new CFiGarchProcess(m_model_spec.garch_p,m_model_spec.garch_q,m_model_spec.vol_power,m_model_spec.figarch_truncation,m_model_spec.vol_rng_seed,m_model_spec.sample_start_idx,m_model_spec.sample_end_idx,m_model_spec.min_bootstrap_sims);
            break;
         case VOL_EGARCH:
            m_vp = new CEgarchProcess(m_model_spec.garch_p,m_model_spec.garch_o,m_model_spec.garch_q,2.0,m_model_spec.vol_rng_seed,m_model_spec.sample_start_idx,m_model_spec.sample_end_idx,m_model_spec.min_bootstrap_sims);
            break;
         default:
            m_vp = new CConstantVariance(m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
        }

      //--- 5. Dynamic Innovation Error Distribution Object Allocation
      if(CheckPointer(m_distribution)==POINTER_DYNAMIC)
         delete m_distribution;
      m_distribution = NULL;

      switch(vol_dist_params.dist_type)
        {
         case DIST_NORMAL:
            m_distribution = new CNormal();
            break;
         case DIST_STUDENT:
            m_distribution = new CStudentsT();
            break;
         case DIST_SKEW_STUDENT:
            m_distribution = new CSkewStudent();
            break;
         case DIST_GEN_ERROR:
            m_distribution = new CGeneralizedError();
            break;
         default:
            m_distribution = new CNormal();
            break;
        }

      //--- 6. Combined Subsystem Initialization and Verification Audit
      m_initialized = (
                         m_distribution!=NULL  &&
                         m_vp!=NULL            &&
                         m_vp.is_initialized() &&
                         m_distribution.initialize(m_model_spec.dist_init_params,
                               m_model_spec.dist_rng_seed)
                      );

      //--- 7. Back-register Derived Subsystem Hyperparameters to Master Specifications
      vector final_vol_specs = m_vp.garch_spec();

      m_model_spec.garch_p = (ulong)final_vol_specs[0];
      m_model_spec.garch_o = (ulong)final_vol_specs[1];
      m_model_spec.garch_q = (ulong)final_vol_specs[2];
      m_model_spec.vol_power = final_vol_specs[3];

      return m_initialized;
     }
public:
                     HARX(void):m_extra_simulation_params(0)
     {
      m_model_spec.mean_model_type = MEAN_HARX;
     }
                     HARX(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_HARX;
      initialize(model_specification);
     }

                    ~HARX(void)
     {
      deinitialize();
     }
   //--- Initialization method checks to ensure model correctly specified
   bool                      initialize(ArchParameters& vol_dist_params)
     {
      return _initialize(vol_dist_params);
     }
   //--- Get exogenous variables if specified
   matrix                    get_x(void) { return m_model_spec.exog_data; }
   //--- Get regressors
   matrix                    get_regressors(void) { return m_regressors; }
   //--- Get observations
   vector                    get_y(void) { return m_model_spec.observations; }
   //--- Calculate residuals
   virtual vector            resids(vector &params, vector &y, matrix &regressors) override
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return EMPTY_VECTOR;
        }

      vector vec;
      if(!y.Size())
        {
         vec = m_fit_regressors.MatMul(params);
         return np::minus(m_fit_y,vec);
        }
      else
        {
         vec = regressors.MatMul(params);
         return np::minus(y,vec);
        }
     }
   //--- Get the number of model parameters
   virtual ulong             num_params(void)
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return 0;
        }
      return m_regressors.Cols();
     }
   //--- Manually set a previously initialized volatility process
   bool              set_volatility_process(CVolatilityProcess* &vp)
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return false;
        }
      if(CheckPointer(vp) == POINTER_INVALID)
        {
         Print(__FUNCTION__, " invalid pointer ");
         return false;
        }

      if(!vp.is_initialized())
        {
         Print(__FUNCTION__, " volatility process has not been initialized ");
         return false;
        }

      if(CheckPointer(m_vp) == POINTER_DYNAMIC)
         delete m_vp;

      m_vp = NULL;
      m_model_spec.vol_model_type = vp.volatilityprocess();
      m_vp = vp;
      m_delete_vp = false;

      return true;
     }
   //--- Specify a previously initialized error distribution
   bool              set_distribution(CDistribution* &dist)
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return false;
        }
      if(CheckPointer(dist) == POINTER_INVALID)
        {
         Print(__FUNCTION__, " invalid pointer ");
         return false;
        }

      if(!dist.is_initialized())
        {
         Print(__FUNCTION__, " distribution has not been initialized ");
         return false;
        }

      if(CheckPointer(m_distribution) == POINTER_DYNAMIC)
         delete m_distribution;

      m_distribution = NULL;
      m_model_spec.dist_type = dist.distribution_type();
      m_distribution = dist;
      m_delete_dist = false;

      return true;
     }
   //+------------------------------------------------------------------+
   //| Simulates data from a linear regression, AR or HAR models       |
   //| Generates synthetic realization paths by combining a structural   |
   //| mean timeline with dynamic volatility pipelines and error shocks. |
   //+------------------------------------------------------------------+
   virtual matrix            simulate(vector &params, ulong nobs, ulong burn, vector &initial_vals, matrix &x, vector &initial_vals_vol)
     {
      //--- 1. Verification Gate: Check class instance state
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return EMPTY_MATRIX;
        }

      ulong kx = x.Cols();
      ulong mc = 0;

      //--- 2. Dimensional Gates: Verify exogenous matrix contains burn-in space allocations
      if(kx && x.Rows()<nobs+burn)
        {
         Print(__FUNCTION__, " x must have nobs + burn rows");
         return matrix::Zeros(0,0);
        }

      // Count parameters assigned exclusively to mean equation processes
      if(m_model_spec.mean_lags.Size())
         mc = ulong(m_constant)+m_lags.Cols()+kx+m_extra_simulation_params;

      // Extract parameter allocation bounds across variance and error distributions
      ulong vc = m_vp.numParams();
      ulong dc = m_distribution.numParams();
      ulong numparams = mc+vc+dc;

      if(params.Size()!=numparams)
        {
         Print(__FUNCTION__, " params has the wrong number of elements.");
         return matrix::Zeros(0,0);
        }

      //--- 3. Parameter Sub-vector Segmentation
      vector allparams[];
      if(!_parse_parameters(params,allparams))
         return matrix::Zeros(0,0);

      //--- 4. Structural Volatility and Innovation Error Modeling
      // Initialize the specialized pseudorandom seed state engine based on target distribution rules
      BootstrapRng rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());

      // Compute joint volatility pathways and return a multi-layered simulation snapshot matrix
      matrix sim_data = m_vp.simulate(allparams[1],nobs+burn,rng,burn,initial_vals_vol.Size()?initial_vals_vol[0]:0.0);

      // Extract unconditioned shock innovations and truncate the pre-sample burn-in allocation block
      vector errors = sim_data.Row(0);
      errors = np::sliceVector(errors,long(burn));

      // Extract conditional variance, convert to volatility ($\sigma_t = \sqrt{\sigma_t^2}$), and drop burn-in entries
      vector vol = sqrt(sim_data.Row(1));
      vol = np::sliceVector(vol,long(burn));

      //--- 5. Recursive Mean Synthesis Loop
      vector y = _simulate_mean(allparams[0],x,errors,initial_vals,sim_data.Row(1));
      y = np::sliceVector(y,long(burn)); // Trim transient burn-in rows from simulated dependent track

      //--- 6. Construct Matrix Output Payload [Columns: Realized Series, Volatility, Error Shocks]
      matrix out = matrix::Zeros(y.Size(),3);
      if(!out.Col(y,0) || !out.Col(vol,1) || !out.Col(errors,2))
        {
         Print(__FUNCTION__, " ", GetLastError());
         return matrix::Zeros(0,0);
        }

      return out;
     }
   //--- Make forecasts from a model
   ArchForecast      forecast(ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
     {
      return forecast(EMPTY_VECTOR,EMPTY_MATRIX_ARRAY,horizon,start,method,simulations,seed);
     }
   //--- Forecast overload
   ArchForecast      forecast(matrix& x[],ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
     {
      return forecast(EMPTY_VECTOR,x,horizon,start,method,simulations,seed);
     }
   //+------------------------------------------------------------------+
   //| Multi-step conditional mean and volatility forecasting           |
   //| Generates multi-period analytic or simulation-based paths,       |
   //| handling rolling lag spaces and exogenous profile inputs ($X$).   |
   //+------------------------------------------------------------------+
   virtual ArchForecast forecast(vector& params, matrix& x[], ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
     {
      ArchForecast out;

      //--- 1. Verification Gate: Check object instance readiness
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return out;
        }

      // Automatically handle sample configuration boundaries if the active set is blank
      if(!m_fit_y.Size())
        {
         long from = 0;
         long to = long(m_model_spec.observations.Size());
         if(!_adjust_sample(from,to))
            return out;
        }

      //--- 2. Sample Alignment and Edge Case Boundary Assertions
      long earliest = long(m_fit_indices[0]);
      long default_start = long(m_fit_indices[1]);
      default_start = MathMax(0,default_start-1);
      long start_index = (start>-1 && start<=default_start)?start:default_start;

      // Enforce the data boundary floor: Reject attempts to look behind the max lag memory window
      if(start_index<(earliest-1))
        {
         Print(__FUNCTION__," Due to backcasting and/or data availability start cannot be less "\
               "than the index of the largest value in the right-hand-side "\
               "variables used to fit the first observation.");
         return out;
        }

      //--- 3. Parameter Deconstruction and Signal Parsing
      vector allparams[];
      if(!CArchModel::_parse_parameters(params.Size()?params:m_params,allparams))
        {
         Print(__FUNCTION__, " failed to parse the parameters ");
         return out;
        }

      // Extract structural residuals and calculate target variance backcast roots
      vector res = resids(allparams[0],vector::Zeros(0),matrix::Zeros(0,0));
      vector bc = m_vp.backCast(res);
      vector ys = np::sliceVector(m_model_spec.observations,long(earliest));
      matrix xs = np::sliceMatrixRows(m_regressors,long(earliest));
      vector full_res = resids(allparams[0],ys,xs);
      matrix vb = m_vp.varianceBounds(full_res,2.);

      // Initialize random distribution engines
      BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
      long vstart = MathMax(0,start_index - earliest);

      //--- 4. Generate Core Conditional Variance Forecast
      VarianceForecast vfcast = m_vp.forecast(allparams[1],full_res,bc,vb,Rng,seed,vstart,horizon,method,simulations);
      if(!vfcast.forecasts.Cols())
         return out;

      matrix var_fcasts = vfcast.forecasts;
      if(start_index<earliest)
         var_fcasts = _forecast_pad(int(earliest-start_index),var_fcasts);

      // Convert HAR multi-day weights into clean point-to-point standard daily AR lag arrays
      vector arp =  _har_to_ar(allparams[0]);
      ulong nexog = m_model_spec.exog_data.Cols();
      vector exog_p = (nexog)?np::sliceVector(allparams[0],-1*long(nexog)):vector::Zeros(0);
      double constant = m_constant?arp[0]:0.0;
      vector dynp = (arp.Size()>ulong(m_constant))?np::sliceVector(arp,long(m_constant)):vector::Zeros(0);

      // Align upcoming exogenous feature tracks ($X_{T+h}$)
      matrix expected_x[];
      if(!_reformat_forecast_x(ulong(start_index),horizon,x,expected_x))
        {
         Print(__FUNCTION__, " reformating error ");
         return out;
        }

      //--- 5. Generate Conditional Mean Forecast Pathways
      matrix mean_fcast = _ar_forecast(m_model_spec.observations,horizon,start_index,constant,dynp,expected_x,exog_p);
      if(!mean_fcast.Cols())
         return out;

      //--- 6. Construct Long-Run Cumulative Volatility Profiles
      // Maps localized shocks to future windows using an MA($\infty$) impulse map ($\psi_h$)
      vector impulse = _ar_to_impulse(horizon,dynp);
      matrix longrun_var_fcasts = var_fcasts;
      matrix tempm;
      vector tempv, lrf;

      // Accumulate squared coefficients along the impulse timeline: $\sigma^2_{T+h} = \sum \psi_j^2 \cdot \sigma^2_{T+h-j}$
      for(ulong i = 0; i<horizon; ++i)
        {
         tempm = np::sliceMatrixCols(var_fcasts,0,long(i+1));
         tempv = np::sliceVector(impulse,long(i),END_REVERSE,STEP_REVERSE);
         tempv = pow(tempv,2.);
         lrf = tempm.MatMul(tempv);
         longrun_var_fcasts.Col(lrf,i);
        }

      //--- 7. Path-Simulation Processing Layer [Bootstrap / Monte Carlo Simulation Paths]
      matrix variance_paths[], mean_paths[];
      matrix shocks[], long_run_variance_paths[];

      if(method == FORECAST_BOOTSTRAP || method == FORECAST_SIMULATION)
        {
         // Copy simulated variance paths and raw standardized shock structures
         ArrayResize(variance_paths,vfcast.forecastpaths.Size());
         for(uint i = 0; i<variance_paths.Size(); ++i)
            variance_paths[i] = vfcast.forecastpaths[i];

         ArrayResize(shocks,vfcast.shocks.Size());
         for(uint i = 0; i<shocks.Size(); ++i)
            shocks[i] = vfcast.shocks[i];

         if(start_index<earliest)
           {
            _forecast_pad(int(earliest-start_index),variance_paths,variance_paths);
            _forecast_pad(int(earliest-start_index),shocks,shocks);
           }

         ArrayResize(long_run_variance_paths,variance_paths.Size());
         for(uint i = 0; i<variance_paths.Size(); ++i)
            long_run_variance_paths[i] = vfcast.forecastpaths[i];

         matrix imps;
         vector _impulses;

         // Evaluate cumulative long-run variance profiles for simulated trials
         for(ulong i = 0; i<horizon; ++i)
           {
            _impulses = np::sliceVector(impulse,long(i),END_REVERSE,STEP_REVERSE);
            imps = matrix::Zeros(_impulses.Size(),1);
            _impulses = pow(_impulses,2.);
            imps.Col(_impulses,0);
            for(uint k = 0; k<variance_paths.Size(); ++k)
              {
               tempm = np::sliceMatrixCols(variance_paths[k],BEGIN,long(i+1));
               tempm = tempm.MatMul(imps);
               long_run_variance_paths[k].Col(tempm.Col(0),i);
              }
           }

         //--- 8. Recursive Mean Synthesis across Simulation Iterations
         ulong t = m_model_spec.observations.Size();
         ArrayResize(mean_paths,shocks.Size());
         for(uint i =0; i<mean_paths.Size(); ++i)
            mean_paths[i] = matrix::Zeros(shocks[0].Rows(),m_maxlags+horizon);

         vector dynp_rev = dynp;
         np::reverseVector(dynp_rev); // Match lag coefficients to reversed chronological rows

         for(ulong i = start_index; i<t; ++i)
           {
            ulong ploc = i - start_index;
            tempv = np::sliceVector(m_model_spec.observations,long(i-m_maxlags+1),long(i+1));
            tempm = np::repeat_vector_as_rows_cols(tempv,mean_paths[ploc].Rows());
            np::matrixCopy(mean_paths[ploc],tempm,BEGIN,END,STEP,BEGIN,long(m_maxlags));

            // Loop step-by-step through the horizon window, updating mean pathways dynamically
            for(ulong j = 0; j<horizon; ++j)
              {
               tempm = np::sliceMatrixCols(mean_paths[ploc],long(j),long(m_maxlags+j));
               tempv = constant + tempm.MatMul(dynp_rev) + shocks[ploc].Col(j);
               mean_paths[ploc].Col(tempv,m_maxlags+j);

               // Inject multi-step exogenous model parameters if assigned
               if(expected_x.Size())
                 {
                  tempv = vector::Zeros(mean_paths[ploc].Rows());
                  for(uint z = 0; z<expected_x.Size(); ++z)
                     tempv[z] = expected_x[z][ploc,j];
                  double dotproduct = tempv.Dot(exog_p);
                  for(ulong k = 0; k<mean_paths[ploc].Rows(); ++k)
                     mean_paths[ploc][k,m_maxlags+j] += dotproduct;
                 }
              }
           }

         // Strip off seed historical lag nodes to leave only pure out-of-sample forecast horizons
         for(uint i = 0; i<mean_paths.Size(); ++i)
            mean_paths[i] = np::sliceMatrixCols(mean_paths[i],long(m_maxlags));
        }

      //--- 9. Package final data arrays into the output container
      return ArchForecast(start_index,mean_fcast,longrun_var_fcasts,var_fcasts,mean_paths,shocks,long_run_variance_paths,variance_paths);
     }
   //+------------------------------------------------------------------+
   //| Construct a fixed-parameter model execution result               |
   //| Bypasses estimation optimization to run evaluation diagnostics  |
   //| directly using a pre-determined user-defined parameter array.    |
   //+------------------------------------------------------------------+
   ArchModelFixedResult fix(vector& params, long first_obs = 0, long last_obs = -1)
     {
      ArchModelFixedResult out;

      //--- 1. Verification Gate: Check class instance state readiness
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return out;
        }

      // Synchronize data sequence indexes with given sample limits
      if(!_adjust_sample(first_obs,last_obs))
         return out;

      //--- 2. Pre-Sample Backcast Conditioning and Target Variance Bounds
      vector sv = starting_values();
      vector res = resids(sv,EMPTY_VECTOR,EMPTY_MATRIX);
      vector sigma2 = vector::Zeros(res.Size());
      vector back_cast = m_vp.backCast(res);
      m_backcast = back_cast; // Cache the historical conditioning backcast value
      matrix vb = m_vp.varianceBounds(res);

      // Compute the overall objective negative log-likelihood score vector
      vector logl = _loglikelihood(params,sigma2,back_cast,vb)*(-1.0);

      //--- 3. Parameter Deconstruction and Signal Slicing
      vector parameters[];
      if(!_parse_parameters(params,parameters))
         return out;

      //--- 4. Direct Volatility Evolution Trajectory Generation
      vector vol = vector::Zeros(res.Size());
      m_vp.computeVariance(parameters[1],res,vol,back_cast,vb);

      // Convert raw variance states into volatility standard deviations ($\sigma_t = \sqrt{\sigma_t^2}$)
      vol = sqrt(vol);

      //--- 5. Synchronize Timeline Data Windows and Padding Assignments
      first_obs = long(m_fit_indices[0]);
      last_obs = long(m_fit_indices[1]);

      // Copy core model residuals into a NaN-padded structural timeline vector
      vector res_final = vector::Zeros(res.Size());
      res_final.Fill(AL_NaN);
      np::vectorCopy(res_final,res,long(first_obs),long(last_obs));

      // Copy estimated conditional volatilities into a matching NaN-padded container
      vector vol_final = vector::Zeros(vol.Size());
      vol_final.Fill(AL_NaN);
      np::vectorCopy(vol_final,vol,long(first_obs),long(last_obs));

      //--- 6. Package final variables into the output object payload
      out.params = params;
      out.resid = res_final;
      out.conditional_volatility = vol_final;
      out.loglikelihood = logl[0];

      return out;
     }
#ifdef __SLSQP__
   //--- Fit a model to data using SLSQP optimizer
   ArchModelResult   fit(double tol = 1e-6,uint maxits = 0,ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1)
     {
      vector starting_values,back_cast;
      return fit(starting_values,back_cast,tol,maxits,cov_type,first,last);
     }
   //+------------------------------------------------------------------+
   //| Fit model using numerical optimization                           |
   //| Orchestrates the optimization lifecycle by concatenating param   |
   //| segments, stacking dynamic boundaries, and passing everything    |
   //| into the SLSQP solver engine to extract optimal coefficients.    |
   //+------------------------------------------------------------------+
   ArchModelResult   fit(vector& startingvalues, vector& backcast, double tol = 1e-6, uint maxits = 0, ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1)
     {
      ArchModelResult out;

      //--- 1. Verification Gates: Validate initialization and data state
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return out;
        }
      if(!m_model_spec.observations.Size())
        {
         Print(__FUNCTION__," Cannot estimate model without data.");
         return out;
        }

      //--- 2. Parse Parameter Counts and Set Sample Boundaries
      vector offsets(3);
      offsets[0] = (double)m_num_params;            // Mean equation parameters
      offsets[1] = (double)m_vp.numParams();         // Volatility process parameters
      offsets[2] = (double)m_distribution.numParams(); // Error distribution shape parameters
      
      double total_params = offsets.Sum();

      if(last<0 || last<first)
         last = long(m_model_spec.observations.Size());

      if(!_adjust_sample(first,last))
         return out;

      // Check for closed-form shortcuts (e.g., constant volatility, normal distribution)
      bool hascloseform = (offsets[2] == 0.0 && m_vp.closedForm() && m_vp.volatilityprocess()==VOL_CONST);
      
      //--- 3. Data Scaling Check and Adjustment
      vector sv = starting_values();
      vector res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
      _checkscale(res);

      if(m_scale!=1.)
        {
         m_model_spec.observations = m_model_spec.observations*m_scale;
         _scalechanged();
         _adjust_sample(first,last);
         sv = starting_values();
         res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
        }
      
      //--- 4. Configure Pre-sample Backcast Values
      vector bcast;
      if(!backcast.Size())
         bcast = m_vp.backCast(res);
      else
         bcast = m_vp.backCastTransform(backcast);
      
      // Direct routing shortcuts for exceptional model layouts
      if(hascloseform)
         return _fit_no_arch_normal_errors(cov_type);
      if(total_params == 0.0)
         return _fit_parameterless_model(cov_type,bcast);
      // Initialize conditional variance pathways and bounds
      vector sigma2 = vector::Zeros(res.Size());
      m_backcast = bcast;
      vector sv_v = m_vp.startingValues(res);
      m_var_bounds = m_vp.varianceBounds(res);
      matrix vb = m_var_bounds;
      m_vp.computeVariance(sv_v,res,sigma2,bcast,vb);
      vector std_res = res/sqrt(sigma2);

      //--- 5. Collect and Vertically Stack Inequality Constraint Matrices
      Constraints cst[3];
      cst[0] = constraints();
      cst[1] = m_vp.constraints();
      cst[2] = m_distribution.constraints();

      vector num_cons(3);
      for(uint i = 0; i<cst.Size(); ++i)
         num_cons[i] = (double)cst[i]._one.Rows();

      matrix a = matrix::Zeros(ulong(num_cons.Sum()),ulong(total_params));
      vector b = vector::Zeros(ulong(num_cons.Sum()));
      long ren, cen, rst, ct;
      vector nc, of;

      for(uint i = 0; i<3; ++i)
        {
         nc = np::sliceVector(num_cons,0,long(i+1));
         of = np::sliceVector(offsets,0,long(i+1));
         ren = long(nc.Sum());
         cen = long(of.Sum());
         rst = ren - long(num_cons[i]);
         ct = cen - long(offsets[i]);
         if(ren-rst>0)
           {
            np::matrixCopy(a,cst[i]._one,rst,ren,STEP,ct,cen);
            np::vectorCopy(b,cst[i]._two,rst,ren);
           }
        }

      //--- 6. Gather and Horizontal Stitch Parameter Min/Max Bounds
      matrix bound[3];
      bound[0] = bounds();
      bound[1] = m_vp.bounds(res);
      bound[2] = m_distribution.bounds(res);
      matrix all_bounds(2,bound[0].Cols()+bound[1].Cols()+bound[2].Cols());

      long cols = 0;
      for(uint i = 0; i<bound.Size(); ++i)
        {
         if(bound[i].Cols())
           {
            cols+=(i)?long(bound[i-1].Cols()):0;
            np::matrixCopyCols(all_bounds,bound[i],cols,long(cols+bound[i].Cols()));
           }
        }

      //--- 7. Validate or Automatically Seed Initial Starting Parameters
      if(startingvalues.Size())
        {
         bool valid = (startingvalues.Size()==ulong(total_params));
         if(a.Rows()&&valid)
           {
            vector temp = b - a.MatMul(startingvalues);
            valid = (valid && temp.Max()<0.0);
           }
         for(uint i = 0; i<startingvalues.Size(); ++i)
            valid = (valid && all_bounds[0,i]<=startingvalues[i] && startingvalues[i]<=all_bounds[1,i]);
         if(!valid)
            startingvalues = EMPTY_VECTOR;
        }

      if(!startingvalues.Size())
        {
         sv = starting_values();
         ulong ss = sv.Size();
         vector temp = m_distribution.startingValues(std_res);
         sv.Resize(ss+sv_v.Size()+temp.Size());
         if(sv_v.Size())
            np::vectorCopy(sv,sv_v,long(ss),long(ss+sv_v.Size()));
         if(temp.Size())
            np::vectorCopy(sv,temp,long(ss+sv_v.Size()));
        }
      else
         sv = startingvalues;

      //--- 8. Prepare and Initialize Objective Functions and Constraints
      CObject             obj;
      CObjectiveF         obfunc;
      CIneqConstraints    obfunc_constraints;

      obfunc.setFuncParams(sigma2,bcast,vb,this);
      obfunc_constraints.setConstraints(a,b);

      all_bounds = all_bounds.Transpose();

      obfunc.setBounds(all_bounds);
      obfunc_constraints.setBounds(all_bounds);

      if(!obfunc.initialize(sv) || !obfunc_constraints.initialize(sv))
        {
         m_converged = false;
         return out;
        }

      //--- 9. Execute Numerical Non-Linear Optimization Loop (SLSQP Solver)
      CSlsqp slsqp_minim;
      slsqp_minim.SetInequalityConstraints(obfunc_constraints);

      if(maxits)
         slsqp_minim.SetMaxEval(int(maxits));
      if(fabs(tol))
        {
         slsqp_minim.SetAcc(tol);
         slsqp_minim.SetXtolRel(tol);
        }

      OptimizeResult minim_result = slsqp_minim.Minimize(obfunc);

      //--- 10. Process Optimization Solutions and Update Internal States
      if(minim_result.return_code<0)
        {
         Print(__FUNCTION__, " termination reason ", minim_result.return_code);
         m_converged = false;
         return out;
        }
      else
        {
         m_params = minim_result.solution;
         m_converged = true;
         //Print(__FUNCTION__, " termination reason ", minim_result.return_code);
        }

      //--- 11. Extract Convergence Data Metrics and Pack Final Results Struct
      vector allparams[];
      vector llh = _loglikelihood(m_params,sigma2,bcast,vb);

      _parse_parameters(m_params,allparams);
      res = resids(allparams[0],EMPTY_VECTOR,EMPTY_MATRIX);

      vector vol = vector::Zeros(res.Size());
      m_vp.computeVariance(allparams[1],res,vol,bcast,vb);
      vol = sqrt(vol);

      double r2 = _r2(allparams[0]);

      long f_obs = (long)m_fit_indices[0];
      long l_obs = (long)m_fit_indices[1];

      vector res_final, vol_final;
      res_final = vol_final = m_model_spec.observations;
      res_final.Fill(AL_NaN);
      vol_final = res_final;

      np::vectorCopy(res_final,res,f_obs,l_obs);
      np::vectorCopy(vol_final,vol,f_obs,l_obs);

      matrix pcov = compute_param_cov(m_params,bcast,(cov_type==COVAR_ROBUST));

      return ArchModelResult(m_params,pcov,r2,res_final,vol_final,cov_type,-1.0*llh[0],f_obs,l_obs);
     }
#else
   //--- Fit model to data using Alglib optimizer
   ArchModelResult   fit(double scaling = 1.0, uint maxits = 0,ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1, double tol = 1e-9, bool guardsmoothness=false, double gradient_test_step = 0.0)
     {
      return fit(EMPTY_VECTOR,EMPTY_VECTOR,scaling,maxits,cov_type,first,last,tol,guardsmoothness,gradient_test_step);
     }
   //+------------------------------------------------------------------+
   //| Fit model using ALGLIB Nonlinear Constrained Optimization (AUL)  |
   //| Orchestrates estimation via Augmented Lagrangian Multipliers,    |
   //| integrating smoothness safeguards, gradient test steps, and      |
   //| custom interactive Jacobian evaluation loops.                    |
   //+------------------------------------------------------------------+
   ArchModelResult   fit(vector& startingvalues, vector& backcast, double scaling = 1.0, uint maxits = 0, ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1, double tol = 1e-9, bool guardsmoothness=false, double gradient_test_step = 0.0)
     {
      ArchModelResult out;

      //--- 1. Verification Gates: Validate initialization and data state
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return out;
        }
      if(!m_model_spec.observations.Size())
        {
         Print(__FUNCTION__," Cannot estimate model without data.");
         return out;
        }

      //--- 2. Parse Parameter Counts and Set Sample Boundaries
      vector offsets(3);
      offsets[0] = (double)m_num_params;            // Mean model parameter count
      offsets[1] = (double)m_vp.numParams();         // Volatility process parameter count
      offsets[2] = (double)m_distribution.numParams(); // Error distribution shape parameter count

      double total_params = offsets.Sum();

      if(last<0 || last<first)
         last = long(m_model_spec.observations.Size());

      if(!_adjust_sample(first,last))
         return out;

      // Check for closed-form shortcuts (e.g., constant volatility, normal distribution)
      bool hascloseform = (offsets[2] == 0.0 && m_vp.closedForm() && m_vp.volatilityprocess()==VOL_CONST);

      //--- 3. Data Scaling Check and Adjustment
      vector sv = starting_values();
      vector res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
      _checkscale(res);

      if(m_scale!=1.)
        {
         m_model_spec.observations = m_model_spec.observations*m_scale;
         _scalechanged();
         _adjust_sample(first,last);
         sv = starting_values();
         res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
        }

      //--- 4. Configure Pre-sample Backcast Values
      vector bcast;
      if(!backcast.Size())
         bcast = m_vp.backCast(res);
      else
         bcast = m_vp.backCastTransform(backcast);

      // Direct routing shortcuts for exceptional model layouts
      if(hascloseform)
         return _fit_no_arch_normal_errors(cov_type);
      if(total_params == 0.0)
         return _fit_parameterless_model(cov_type,bcast);

      // Initialize conditional variance pathways and bounds
      vector sigma2 = vector::Zeros(res.Size());
      m_backcast = bcast;
      vector sv_v = m_vp.startingValues(res);
      m_var_bounds = m_vp.varianceBounds(res);
      matrix vb = m_var_bounds;
      m_vp.computeVariance(sv_v,res,sigma2,bcast,vb);
      vector std_res = res/sqrt(sigma2);

      //--- 5. Collect and Format Inequality Constraint Systems ($A \cdot x \ge b$)
      Constraints cst[3];
      cst[0] = constraints();
      cst[1] = m_vp.constraints();
      cst[2] = m_distribution.constraints();

      vector num_cons(3);
      for(uint i = 0; i<cst.Size(); ++i)
         num_cons[i] = (double)cst[i]._one.Rows();

      // Allocate space with an extra column to hold linear condition offsets natively
      matrix a = matrix::Zeros(ulong(num_cons.Sum()),ulong(total_params)+1);
      vector b = vector::Zeros(ulong(num_cons.Sum()));
      long ren, cen, rst, ct;
      vector nc, of;

      for(uint i = 0; i<3; ++i)
        {
         nc = np::sliceVector(num_cons,0,long(i+1));
         of = np::sliceVector(offsets,0,long(i+1));
         ren = long(nc.Sum());
         cen = long(of.Sum());
         rst = ren - long(num_cons[i]);
         ct = cen - long(offsets[i]);
         if(ren-rst>0)
           {
            np::matrixCopy(a,cst[i]._one,rst,ren,STEP,ct,cen);
            np::vectorCopy(b,cst[i]._two,rst,ren);
           }
        }

      // Convert constraints into standard ALGLIB linear format: [A | -b] . x >= 0
      a.Col(b,ulong(total_params));
      b = vector::Ones(a.Rows());

      //--- 6. Gather Boundary Matrices [Min / Max Box Constraints]
      matrix bound[3];
      bound[0] = bounds();
      bound[1] = m_vp.bounds(res);
      bound[2] = m_distribution.bounds(res);
      matrix all_bounds(2,bound[0].Cols()+bound[1].Cols()+bound[2].Cols());

      long cols = 0;
      for(uint i = 0; i<bound.Size(); ++i)
        {
         if(bound[i].Cols())
           {
            cols+=(i)?long(bound[i-1].Cols()):0;
            np::matrixCopyCols(all_bounds,bound[i],cols,long(cols+bound[i].Cols()));
           }
        }

      //--- 7. Validate Provided Starting Vector Against Linear and Box Constraints
      if(startingvalues.Size())
        {
         bool valid = (startingvalues.Size()==ulong(total_params));
         matrix aa = np::sliceMatrixCols(a,0,long(a.Cols()-1));
         vector bb = a.Col(a.Cols()-1);
         if(aa.Rows()&&valid)
           {
            vector temp = aa.MatMul(startingvalues) - bb;
            valid = (valid && temp.Min()>0.0);
           }
         for(uint i = 0; i<startingvalues.Size(); ++i)
            valid = (valid && all_bounds[0,i]<=startingvalues[i] && startingvalues[i]<=all_bounds[1,i]);
         if(!valid)
            startingvalues = EMPTY_VECTOR;
        }

      // Compute structural initial guesses if the vector remains empty
      if(!startingvalues.Size())
        {
         sv = starting_values();
         ulong ss = sv.Size();
         vector temp = m_distribution.startingValues(std_res);
         sv.Resize(ss+sv_v.Size()+temp.Size());
         if(sv_v.Size())
            np::vectorCopy(sv,sv_v,long(ss),long(ss+sv_v.Size()));
         if(temp.Size())
            np::vectorCopy(sv,temp,long(ss+sv_v.Size()));
        }
      else
         sv = startingvalues;

      //--- 8. Initialize ALGLIB Nonlinear Constrained State Machine (MinNLC)
      CMinNLCState optim_state;
      CRowDouble _initial_values = sv;

      CRowDouble _ones = vector::Ones(sv.Size());
      if(scaling>0.0)
         _ones*=scaling;

      CRowDouble bl = all_bounds.Row(0);
      CRowDouble bu = all_bounds.Row(1);
      CMatrixDouble _constraints = a;

      CMinNLC::MinNLCCreate(_initial_values.Size(),_initial_values,optim_state);
      CMinNLC::MinNLCSetCond(optim_state,tol,int(maxits));
      CMinNLC::MinNLCSetScale(optim_state,_ones);
      CMinNLC::MinNLCSetSTPMax(optim_state,0.0);
      CMinNLC::MinNLCSetAlgoAUL(optim_state,10,0); // Select Augmented Lagrangian Multipliers
      CMinNLC::MinNLCSetBC(optim_state,bl,bu);    // Bind simple box boundaries

      CRowInt intones;
      intones.Resize(_constraints.Rows());
      intones.Fill(1); // Identify all rows as inequality expressions (>= 0)
      CMinNLC::MinNLCSetLC(optim_state,_constraints,intones,intones.Size());

      // Configure algorithmic safety features
      if(guardsmoothness)
         CMinNLC::MinNLCOptGuardSmoothness(optim_state,guardsmoothness);
      if(gradient_test_step>0.0)
         CMinNLC::MinNLCOptGuardGradient(optim_state,gradient_test_step);

      CMinNLCReport rep;
      CRowDouble solution;

      // Set up objective function structural parameters
      CObjectiveJacFVec fvec;
      fvec.SetParams(sigma2,bcast,vb,this);

      CNDimensional_Rep f_rep;
      CObject object;

      if(!CAp::Assert(GetPointer(fvec)!=NULL,"ALGLIB: error in 'minnlcoptimize()' (fvec is null)"))
        {
         m_converged = false;
         return out;
        }

      //--- 9. Interactive Optimization Routine (Jacobian / Step-Evaluation Cycle)
      while(CMinNLC::MinNLCIteration(optim_state))
        {
         // Handle request to re-evaluate function vector and full Jacobian matrix
         if(optim_state.m_needfij)
           {
            fvec.Jac(optim_state.m_x,optim_state.m_fi,optim_state.m_j,object);
            continue;
           }
         // Handle request to report update steps
         if(optim_state.m_xupdated)
           {
            if(GetPointer(f_rep)!=NULL)
               f_rep.Rep(optim_state.m_x,optim_state.m_f,object);
            continue;
           }
         CAp::Assert(false,"ALGLIB: error in 'minnlcoptimize' (some derivatives were not provided?)");
         break;
        }

      //--- 10. Extract Optimization Results and Assess Convergence Status
      CMinNLC::MinNLCResults(optim_state,solution,rep);
      if(rep.m_terminationtype<0)
        {
         Print(__FUNCTION__," failed convergence: ", rep.m_terminationtype);
         m_converged = false;
         return out;
        }
      else
         m_converged = true;

      //--- 11. Extract Convergence Data Metrics and Pack Final Results Struct
      vector allparams[];
      m_params = solution.ToVector();

      // Calculate final target log-likelihood
      vector llh = _loglikelihood(m_params,sigma2,bcast,vb);

      _parse_parameters(m_params,allparams);
      res = resids(allparams[0],EMPTY_VECTOR,EMPTY_MATRIX);

      vector vol = vector::Zeros(res.Size());
      m_vp.computeVariance(allparams[1],res,vol,bcast,vb);
      vol = sqrt(vol);

      double r2 = _r2(allparams[0]);

      long f_obs = (long)m_fit_indices[0];
      long l_obs = (long)m_fit_indices[1];

      vector res_final, vol_final;
      res_final = vol_final = m_model_spec.observations;
      res_final.Fill(AL_NaN);
      vol_final = res_final;

      np::vectorCopy(res_final,res,f_obs,l_obs);
      np::vectorCopy(vol_final,vol,f_obs,l_obs);

      matrix pcov = compute_param_cov(m_params,bcast,(cov_type==COVAR_ROBUST));

      return ArchModelResult(m_params,pcov,r2,res_final,vol_final,cov_type,-1.0*llh[0],f_obs,l_obs);
     }
#endif
   //+------------------------------------------------------------------+
   //| Accessible loglikelihood evaluation routine                      |
   //| Replicates the exact environment setup used during estimation    |
   //| to compute a clean out-of-sample or spot likelihood metric.      |
   //+------------------------------------------------------------------+
   bool              eval_loglikelihood(vector& params, double &out_loglikelihood)
     {
      long first = 0;
      long last = -1;

      vector startingvalues = EMPTY_VECTOR;
      vector backcast = EMPTY_VECTOR;

      //--- 1. Verification Gates: Validate initialization and data state
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return false;
        }
      if(!m_model_spec.observations.Size())
        {
         Print(__FUNCTION__," Cannot estimate model without data.");
         return false;
        }

      //--- 2. Parse Parameter Counts and Set Sample Boundaries
      vector offsets(3);
      offsets[0] = (double)m_num_params;            // Mean model parameter count
      offsets[1] = (double)m_vp.numParams();         // Volatility process parameter count
      offsets[2] = (double)m_distribution.numParams(); // Error distribution shape parameter count

      double total_params = offsets.Sum();

      if(last<0 || last<first)
         last = long(m_model_spec.observations.Size());

      if(!_adjust_sample(first,last))
         return false;

      // Check for closed-form shortcuts (e.g., constant volatility, normal distribution)
      bool hascloseform = (offsets[2] == 0.0 && m_vp.closedForm() && m_vp.volatilityprocess()==VOL_CONST);

      //--- 3. Data Scaling Check and Adjustment
      vector sv = starting_values();
      vector res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
      _checkscale(res);

      if(m_scale!=1.)
        {
         m_model_spec.observations = m_model_spec.observations*m_scale;
         _scalechanged();
         _adjust_sample(first,last);
         sv = starting_values();
         res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
        }

      //--- 4. Configure Pre-sample Backcast Values
      vector bcast;
      if(!backcast.Size())
         bcast = m_vp.backCast(res);
      else
         bcast = m_vp.backCastTransform(backcast);

      // Direct routing failures for exceptional model layouts
      if(hascloseform)
         return false;
      if(total_params == 0.0)
         return false;

      // Initialize conditional variance pathways and bounds
      vector sigma2 = vector::Zeros(res.Size());
      m_backcast = bcast;
      vector sv_v = m_vp.startingValues(res);
      m_var_bounds = m_vp.varianceBounds(res);
      matrix vb = m_var_bounds;
      m_vp.computeVariance(sv_v,res,sigma2,bcast,vb);
      vector std_res = res/sqrt(sigma2);

      //--- 5. Collect and Format Inequality Constraint Systems ($A \cdot x \ge b$)
      Constraints cst[3];
      cst[0] = constraints();
      cst[1] = m_vp.constraints();
      cst[2] = m_distribution.constraints();

      vector num_cons(3);
      for(uint i = 0; i<cst.Size(); ++i)
         num_cons[i] = (double)cst[i]._one.Rows();

      // Allocate matrix rows matching stacked components, including linear boundary offsets
      matrix a = matrix::Zeros(ulong(num_cons.Sum()),ulong(total_params)+1);
      vector b = vector::Zeros(ulong(num_cons.Sum()));
      long ren, cen, rst, ct;
      vector nc, of;

      for(uint i = 0; i<3; ++i)
        {
         nc = np::sliceVector(num_cons,0,long(i+1));
         of = np::sliceVector(offsets,0,long(i+1));
         ren = long(nc.Sum());
         cen = long(of.Sum());
         rst = ren - long(num_cons[i]);
         ct = cen - long(offsets[i]);
         if(ren-rst>0)
           {
            np::matrixCopy(a,cst[i]._one,rst,ren,STEP,ct,cen);
            np::vectorCopy(b,cst[i]._two,rst,ren);
           }
        }

      // Convert structural constraint records to standard linear offset profiles
      a.Col(b,ulong(total_params));
      b = vector::Ones(a.Rows());

      //--- 6. Gather Boundary Matrices [Min / Max Box Constraints]
      matrix bound[3];
      bound[0] = bounds();
      bound[1] = m_vp.bounds(res);
      bound[2] = m_distribution.bounds(res);
      matrix all_bounds(2,bound[0].Cols()+bound[1].Cols()+bound[2].Cols());

      long cols = 0;
      for(uint i = 0; i<bound.Size(); ++i)
        {
         if(bound[i].Cols())
           {
            cols+=(i)?long(bound[i-1].Cols()):0;
            np::matrixCopyCols(all_bounds,bound[i],cols,long(cols+bound[i].Cols()));
           }
        }

      //--- 7. Validate Provided Starting Vector Against Constraints
      if(startingvalues.Size())
        {
         bool valid = (startingvalues.Size()==ulong(total_params));
         matrix aa = np::sliceMatrixCols(a,0,long(a.Cols()-1));
         vector bb = a.Col(a.Cols()-1);
         if(aa.Rows()&&valid)
           {
            vector temp = aa.MatMul(startingvalues) - bb;
            valid = (valid && temp.Min()>0.0);
           }
         for(uint i = 0; i<startingvalues.Size(); ++i)
            valid = (valid && all_bounds[0,i]<=startingvalues[i] && startingvalues[i]<=all_bounds[1,i]);
         if(!valid)
            startingvalues = EMPTY_VECTOR;
        }

      // Setup initial fallback structural seeds if vectors are unassigned
      if(!startingvalues.Size())
        {
         sv = starting_values();
         ulong ss = sv.Size();
         vector temp = m_distribution.startingValues(std_res);
         sv.Resize(ss+sv_v.Size()+temp.Size());
         if(sv_v.Size())
            np::vectorCopy(sv,sv_v,long(ss),long(ss+sv_v.Size()));
         if(temp.Size())
            np::vectorCopy(sv,temp,long(ss+sv_v.Size()));
        }
      else
         sv = startingvalues;

      //--- 8. Execute Target Core Log-Likelihood Extraction Engine
#ifdef __SLSQP__
      // Route evaluation through the internal SLSQP functional structure
      CObjectiveF fvec;
      fvec.setFuncParams(sigma2,bcast,vb,this);
      out_loglikelihood = fvec.orig_fun(params);
#else
      // Route evaluation through the standard ALGLIB functional object layout
      CObjectiveFunc fvec;
      fvec.SetParams(sigma2,bcast,vb,this);
      CObject ob;
      CRowDouble in_x = params;
      fvec.Func(in_x,out_loglikelihood,ob);
#endif

      return true;
     }
  };
//+------------------------------------------------------------------+
//| HAR model (Heterogeneous Autoregressive)                         |
//| Models the conditional mean by decomposing volatility into       |
//| daily, weekly, and monthly components to capture long memory.    |
//+------------------------------------------------------------------+
class HAR: public HARX
  {
public:
   // Default constructor: sets the model specification type to HAR
                     HAR(void)
     {
      m_model_spec.mean_model_type = MEAN_HAR;
     }

   // Overloaded constructor: initializes with a specific ArchParameters configuration
                     HAR(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_HAR;
      initialize(model_specification);
     }

   // Destructor: ensures proper cleanup of inherited resources
                    ~HAR(void)
     {
      deinitialize();
     }
  };

//+------------------------------------------------------------------+
//| The constant mean model                                          |
//| Represents a simple process: y_t = mu + epsilon_t                |
//+------------------------------------------------------------------+
class ConstantMean: public HARX
  {
public:
   // Default constructor: assigns constant mean model type
                     ConstantMean(void)
     {
      m_model_spec.mean_model_type = MEAN_CONSTANT;
     }

   // Overloaded constructor: initializes and sets model type
                     ConstantMean(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_CONSTANT;
      initialize(model_specification);
     }

   // Cleanup
                    ~ConstantMean(void)
     {
      deinitialize();
     }

   // Returns the count of parameters (only 1: the mean mu)
   virtual ulong     num_params(void) override
     {
      return 1;
     }

   // Computes model residuals (epsilon_t = y_t - mu)
   virtual vector    resids(vector& params, vector& y, matrix& regressors) override
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return EMPTY_VECTOR;
        }
      // Use internal fit data if y is empty, otherwise use provided vector
      if(!y.Size())
         return np::minus(m_fit_y,params);
      return np::minus(y,params);
     }

   // Generates synthetic data based on the constant mean and volatility process
   virtual matrix    simulate(vector &params, ulong nobs, ulong burn, vector &initial_vals, matrix &x, vector &initial_vals_vol)
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return EMPTY_MATRIX;
        }

      // Ensure no external regressors are passed for a simple constant mean simulation
      if(initial_vals.Size() || x.Cols())
        {
         Print(__FUNCTION__, " Both initial value and x must be none when simulating a constant mean process.");
         return matrix::Zeros(0,0);
        }

      // Extract segmented parameters (mean, volatility process, distribution)
      vector allparams[];
      if(!_parse_parameters(params,allparams))
        {
         Print(__FUNCTION__, " failed to parse parameters ");
         return matrix::Zeros(0,0);
        }

      // Initialize random number generator for innovations
      BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());

      // Simulate underlying volatility path
      matrix sim_values = m_vp.simulate(allparams[1],nobs+burn,Rng,burn,initial_vals_vol.Size()?initial_vals_vol[0]:0.0);

      // Construct return series by adding mean to simulated residuals (Row 0 = errors, Row 1 = variance)
      vector ers = sim_values.Row(0);
      vector y = ers + allparams[0];
      vector vol = sqrt(sim_values.Row(1));

      // Remove burn-in period samples
      y = np::sliceVector(y,long(burn));
      vol = np::sliceVector(vol,long(burn));
      ers = np::sliceVector(ers,long(burn));

      // Package results into matrix [y | volatility | residuals]
      matrix out(y.Size(),3);
      if(!out.Col(y,0) || !out.Col(vol,1) || !out.Col(ers,2))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return matrix::Zeros(0,0);
        }

      return out;
     }
  };
//+------------------------------------------------------------------+
//| The zero mean model                                              |
//| Assumes the return series has no deterministic mean (mu = 0).    |
//| Used primarily for high-frequency data or stationary returns.    |
//+------------------------------------------------------------------+
class ZeroMean: public HARX
  {
public:
   // Default constructor: assigns zero mean model type
                     ZeroMean(void)
     {
      m_model_spec.mean_model_type = MEAN_ZERO;
     }

   // Overloaded constructor: initializes with configuration
                     ZeroMean(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_ZERO;
      initialize(model_specification);
     }

   // Cleanup
                    ~ZeroMean(void)
     {
      deinitialize();
     }

   // Returns 0 parameters as there is no mean coefficient to estimate
   virtual ulong     num_params(void) override
     {
      return 0;
     }

   // Computes residuals (epsilon_t = y_t)
   virtual vector    resids(vector& params, vector& y, matrix& regressors) override
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return EMPTY_VECTOR;
        }
      // If y is not provided, return the internal model data directly
      if(!y.Size())
         return m_fit_y;
      return y;
     }

   // Generates synthetic data driven purely by the volatility process
   virtual matrix    simulate(vector &params, ulong nobs, ulong burn, vector &initial_vals, matrix &x, vector &initial_vals_vol)
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return EMPTY_MATRIX;
        }

      // Ensure no initial conditions or external regressors are used
      if(initial_vals.Size() || x.Cols())
        {
         Print(__FUNCTION__, " Both initial value and x must be none when simulating a zero mean process.");
         return matrix::Zeros(0,0);
        }

      vector allparams[];
      if(!_parse_parameters(params,allparams))
        {
         Print(__FUNCTION__, " failed to parse parameters ");
         return matrix::Zeros(0,0);
        }

      // Initialize generator for innovations
      BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());

      // Simulate volatility path (Row 0 = innovations, Row 1 = variance)
      matrix sim_values = m_vp.simulate(allparams[1],nobs+burn,Rng,burn,initial_vals_vol[0]);

      vector ers = sim_values.Row(0);
      vector y = ers; // y_t = epsilon_t
      vector vol = sqrt(sim_values.Row(1));

      // Remove burn-in samples
      y = np::sliceVector(y,long(burn));
      vol = np::sliceVector(vol,long(burn));
      ers = np::sliceVector(ers,long(burn));

      // Package output: [y | volatility | residuals]
      matrix out(y.Size(),3);
      if(!out.Col(y,0) || !out.Col(vol,1) || !out.Col(ers,2))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return matrix::Zeros(0,0);
        }

      return out;
     }
  };
//+------------------------------------------------------------------+
//| The AR model (Autoregressive)                                    |
//| Models y_t as a function of previous observations (lags).        |
//| Captures autocorrelation present in return series.               |
//+------------------------------------------------------------------+
class AR: public HARX
  {
public:
   // Default constructor: assigns AR model type
                     AR(void)
     {
      m_model_spec.mean_model_type = MEAN_AR;
     }

   // Overloaded constructor: initializes with configuration
                     AR(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_AR;
      initialize(model_specification);
     }

   // Clean up resources
                    ~AR(void)
     {
      deinitialize();
     }
  };

//+------------------------------------------------------------------+
//| The AR-X model (Autoregressive with Exogenous Variables)         |
//| Extends the AR model by incorporating external factors (X).      |
//| Useful for multifactor modeling and structural analysis.         |
//+------------------------------------------------------------------+
class ARX: public HARX
  {
public:
   // Default constructor: assigns ARX model type
                     ARX(void)
     {
      m_model_spec.mean_model_type = MEAN_ARX;
     }

   // Overloaded constructor: initializes with configuration
                     ARX(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_ARX;
      initialize(model_specification);
     }

   // Clean up resources
                    ~ARX(void)
     {
      deinitialize();
     }
  };
//+------------------------------------------------------------------