//+------------------------------------------------------------------+
//|                                                         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:

   ArchParameters    m_model_spec;

   vector            m_params,m_fit_indices,m_backcast;
   vector            m_fit_y;
   matrix            m_fit_x,m_fit_regressors;
   double            m_scale;
   ulong             m_nobs;
   ulong             m_holdback;
   ulong             m_maxlags;
   ulong             m_num_params;
   bool              m_converged,m_rescale,m_constant,m_rotated;
   string            m_name;
   matrix            m_regressors,m_lags,m_var_bounds;
   bool              m_delete_vp,m_delete_dist;

   CDistribution     *m_distribution;
   CVolatilityProcess *m_vp;

   bool               m_initialized;
   vector            _get_epsilon(vector &x, double s, ulong n, double epsilon = EMPTY_VALUE)
     {
      vector out;
      if(epsilon == EMPTY_VALUE)
         out = pow(DBL_EPSILON,1./s)*np::maximum(MathAbs(x),0.1);
      else
        {
         out = vector::Zeros(n);
         out.Fill(epsilon);
        }
      return out;
     }
   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)
        {
         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
        {
         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();
     }
   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);
      matrix hess = h.Outer(h);
      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];
           }
        }
      return hess;
     }
   void              _checkscale(vector& resids)
     {
      double ogscale = resids.Var();
      double scale = ogscale;
      double rescale = 1.;
      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 void      _scalechanged(void)
     {
      return;
     }
   virtual double    _r2(vector& params)
     {
      return EMPTY_VALUE;
     }
   virtual vector    _fit_no_arch_normal_errors_params(void)
     {
      return vector::Zeros(0);
     }
   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;
     }
   bool              _parse_parameters(vector &x, vector& out[])
     {
      ulong km = m_num_params;
      ulong kv = m_vp.numParams();
      ArrayResize(out,3);
      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;
     }

   vector            _loglikelihood(vector& parameters, vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false)
     {
      vector prs[];
      _parse_parameters(parameters,prs);
      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);
     }
   virtual void       _adjust_sample(long& first, long& last)
     {
      return;
     }
   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:
                     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; }
   virtual bool              Save(const int file_handle) override
     {
      if(file_handle!=INVALID_HANDLE)
        {
         //---
         if(FileWriteLong(file_handle,-1)==sizeof(long))
           {
            //---
            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;
                 }
              }
            //---
            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;
                 }
              }
            //---
            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;
                 }
              }
            //---
            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);
              }

            //---
            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();
     }

   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();
     }
   vector                   objective(vector& parameters, vector& sigma2, vector &backcast, matrix& varbounds, bool individual=false)
     {
      return _loglikelihood(parameters,sigma2,backcast,varbounds,individual);
     }
   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);
     }
   virtual 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();
      if(robust)
        {
         matrix scores = approx_fprime(params,sgma2,backcast,vb,true);
         matrix score_cov = scores.Transpose().Cov();
         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 ArchModelFixedResult
  {
public:
   double            loglikelihood;
   vector            params;
   vector            conditional_volatility;
   ulong             nobs;
   vector            resid;

                     ArchModelFixedResult(void)
     {
      loglikelihood = EMPTY_VALUE;
      nobs = 0;
      params = resid = conditional_volatility = vector::Zeros(0);
     }

                     ArchModelFixedResult(ArchModelFixedResult &other)
     {
      loglikelihood = other.loglikelihood;
      nobs = other.nobs;
      params = other.params;
      resid = other.resid;
      conditional_volatility = other.conditional_volatility;
     }
                     ArchModelFixedResult(vector &_params, vector &_resid, vector &_volatility, vector &_depvar, double _loglikelihood)
     {
      loglikelihood = _loglikelihood;
      params = _params;
      resid = _resid;
      conditional_volatility = _volatility;
     }
                    ~ArchModelFixedResult(void)
     {
     }
   void              operator=(ArchModelFixedResult &other)
     {
      loglikelihood = other.loglikelihood;
      nobs = other.nobs;
      params = other.params;
      resid = other.resid;
      conditional_volatility = other.conditional_volatility;
     }

   double            aic(void)
     {
      return -2 * loglikelihood + 2 * double(num_params());
     }
   ulong             num_params(void)
     {
      return params.Size();
     }
   double            bic(void)
     {
      return -2.*loglikelihood+log(nobs)*double(num_params());
     }
   vector            std_resid(void)
     {
      return resid/conditional_volatility;
     }
  };
//+------------------------------------------------------------------+
//| arch model result                                                |
//+------------------------------------------------------------------+
struct ArchModelResult: public ArchModelFixedResult
  {
public:
   long              fit_indices[2];
   matrix            param_cov;
   double            r2;
   ENUM_COVAR_TYPE   cov_type;

                     ArchModelResult(void)
     {
      fit_indices[0] = fit_indices[1] = 0;
      param_cov = matrix::Zeros(0,0);
      r2 = EMPTY_VALUE;
      cov_type = WRONG_VALUE;
     }
                     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;

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

                    ~ArchModelResult(void)
     {

     }

                     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;
     }

   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;
     }
   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);
      else
         param_cov =  archmodel.compute_param_cov(fitted_params,EMPTY_VECTOR,false);
      return param_cov;
     }
   matrix            conf_int(double alpha = 0.05)
     {
      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
      out.Col(params + cv *se,1);//upper

      return out;
     }
   double            rsquared_adj(void)
     {
      return 1.0 - ((1.0-r2)*double(nobs-1)/double(nobs-num_params()));
     }
   vector            std_err(void)
     {
      return sqrt(param_cov.Diag());
     }
   vector            pvalues(void)
     {
      vector pvals = tvalues();
      for(ulong i = 0; i<pvals.Size(); ++i)
         pvals[i] = CNormalDistr::NormalCDF(-1.0*MathAbs(pvals[i]))*2.0;

      return pvals;
     }
   vector            tvalues(void)
     {
      return params/std_err();
     }
  };
//+------------------------------------------------------------------+
//| The Lagrange multiplier test                                     |
//+------------------------------------------------------------------+
WaldTestStatistic archlmtest(vector& residuals,vector& conditional_volatility,ulong lags,bool standardized = false)
  {
   WaldTestStatistic out;
   vector resids = residuals;
   if(standardized)
      resids = resids/conditional_volatility;
   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);
   if(!lags)
      lags = ulong(ceil(12.0*pow(nobs/100.0,1.0/4.0)));
   lags = MathMax(MathMin(resids.Size()/2 - 1,lags),1);
   if(resid2.Size()<3)
     {
      Print(__FUNCTION__, " Test requires at least 3 non-nan observations ");
      return out;
     }
   matrix matres = matrix::Zeros(resid2.Size(),1);
   matres.Col(resid2,0);
   matrix lag[];
   if(!lagmat(matres,lag,lags,TRIM_BOTH,ORIGINAL_SEP))
     {
      Print(__FUNCTION__, " lagmat failed ");
      return out;
     }
   matrix lagout;
   if(!addtrend(lag[0],lagout,TREND_CONST_ONLY,true,HAS_CONST_SKIP))
     {
      Print(__FUNCTION__, " addtrend failed ");
      return out;
     }
   lag[0] = lagout;
   OLS ols;
   if(!ols.Fit(lag[1].Col(0),lag[0]))
     {
      Print(__FUNCTION__, " OLS fitting failed ");
      return out;
     }

   out.stat = nobs*ols.Rsqe();
   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";
     }

   out.df = lags;
   return out;
  }
//+------------------------------------------------------------------+
//| helper                                                           |
//+------------------------------------------------------------------+
bool _forecast_pad(int count, matrix &forecasts[], matrix &out[])
  {
   ArrayResize(out,count+(int)forecasts.Size());
   matrix fill = matrix::Zeros(forecasts[0].Rows(), forecasts[0].Cols());
   fill.Fill(double("nan"));

   for(int i = 0; i<count; ++i)
      out[i] = fill;
   for(int i = count; i<int(out.Size()); ++i)
      out[i] = forecasts[i-count];

   return true;
  }
//+------------------------------------------------------------------+
//| helper                                                           |
//+------------------------------------------------------------------+
matrix _forecast_pad(int count, matrix &forecasts)
  {
   matrix fill = matrix::Zeros(count+forecasts.Rows(), forecasts.Cols());
   fill.Fill(double("nan"));

   for(ulong i = (ulong)count; i<fill.Rows(); ++i)
      fill.Row(forecasts.Row(ulong(i-count)),i);

   return fill;
  }
//+------------------------------------------------------------------+
//| helper                                                           |
//+------------------------------------------------------------------+
vector _forecast_pad(int count, vector &forecasts)
  {
   vector fill = vector::Zeros(count+forecasts.Size());
   fill.Fill(double("nan"));

   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                      |
//+------------------------------------------------------------------+
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;
   matrix fcasts = matrix::Zeros(t - start_index,p+horizon);
   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);
     }
   vector arp_rev  = arp;
   np::reverseVector(arp_rev);
   matrix mlt;
   vector mat;
   for(ulong i = p; i<horizon+p; ++i)
     {
      mlt = ((i-p)!=i)?np::sliceMatrixCols(fcasts,long(i-p),long(i)):matrix::Zeros(0,0);
      mat = (arp_rev.Size()&&mlt.Cols())?mlt.MatMul(arp_rev):vector::Zeros(fcasts.Rows());
      fcasts.Col(constant+mat,i);
      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 = fcasts.Col(i);
         col += mlt.MatMul(exogp);
         fcasts.Col(col,i);
        }
     }
   fcasts = np::sliceMatrixCols(fcasts,long(p));
   return fcasts;
  }

//+------------------------------------------------------------------+
//| impulse                                                          |
//+------------------------------------------------------------------+
vector _ar_to_impulse(ulong steps,vector& params)
  {
   ulong p = params.Size();
   vector impulse = vector::Zeros(steps);
   impulse[0] = 1.;
   if(p == 0)
      return impulse;
   long k, st;
   vector rparams,imp;
   for(long i = 1; i<long(steps); ++i)
     {
      k =  fmin(long(p) -1, i - 1);
      st = fmax(i - long(p),0);
      imp  = np::sliceVector(impulse,st,i);
      rparams = np::sliceVector(params,k,END_REVERSE,STEP_REVERSE);
      if(imp.Size()==rparams.Size() && imp.Size())
         impulse[i] = imp.Dot(rparams);
     }
   return impulse;
  }
//+------------------------------------------------------------------+
//| Heterogeneous Autoregression (HAR)                               |
//+------------------------------------------------------------------+
class HARX: public CArchModel
  {
protected:
#ifdef __SLSQP__
   class CObjectiveF: public CFunctor
     {
   private:
      vector            m_sigma,m_bc;
      matrix            m_vb;
      HARX*             m_obj;
      bool              m_centered,m_individual;
      double            m_epsilon;
      vector            obj_result;
      matrix            obj_gradient;
   public:
      //--- constructor, destructor
                     CObjectiveF(void) {}
                    ~CObjectiveF(void)
        {
         m_obj = NULL;
         obj_result = vector::Zeros(1);
         obj_gradient = matrix::Zeros(1,2);
        }
      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);
         m_individual = individual;
         m_centered = center;
         m_epsilon = ep;
         m_clip = true;
         m_grad_options.method = GRAD_POINT_2;
        }
      virtual double    orig_fun(vector& x) override
        {
         obj_result = CheckPointer(m_obj)!=POINTER_INVALID?m_obj.objective(x,m_sigma,m_bc,m_vb):obj_result;
         return obj_result[0];
        }
      virtual vector grad_fun(vector& x) override
        {
         obj_gradient = CheckPointer(m_obj)!=POINTER_INVALID?m_obj.jacobian(x,m_sigma,m_bc,m_vb,m_individual,m_centered,m_epsilon):obj_gradient;
         return obj_gradient.Row(0);
        }
     };
   class CIneqConstraints: public CConstraints
     {
   private:
      matrix            a;
      vector            b;
      vector            obj_result;
   public:
                     CIneqConstraints(void)
        {
        }
                    ~CIneqConstraints(void)
        {
        }
      void              setConstraints(matrix& a_, vector& b_)
        {
         a = a_;
         b = b_;
         m_m = int(a.Rows());
         m_n = int(a.Cols());
         obj_result =  vector::Zeros(m_n);
         m_options.method = GRAD_POINT_2;
        }
      virtual vector    orig_fun(vector& x) override
        {
         if(!a.Rows())
            return obj_result;
         obj_result = b - a.MatMul(x);
         return obj_result;
        }
     };
#else
   class CObjectiveFunc: public CNDimensional_Func
     {
   private:
      vector            m_sigma,m_bc;
      matrix            m_vb;
      HARX*             m_obj;
      bool              m_centered,m_individual;
      double            m_epsilon;
   public:
      //--- constructor, destructor
                     CObjectiveFunc(void) {}
                    ~CObjectiveFunc(void)
        {
         m_obj = NULL;
        }
      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;
        }
      virtual void      Func(double &x[],double &func,CObject &obj) override
        {
         vector temp;
         temp.Assign(x);
         vector r = m_obj.objective(temp,m_sigma,m_bc,m_vb);
         func = r[0];
        }
      virtual void      Func(CRowDouble &x,double &func,CObject &obj) override
        {
         vector temp = x.ToVector();
         vector r = m_obj.objective(temp,m_sigma,m_bc,m_vb);
         func = r[0];
        }
     };

   //+------------------------------------------------------------------+
   //|jacobian fvec implementation                                      |
   //+------------------------------------------------------------------+
   class CObjectiveJacFVec: public CNDimensional_Jac
     {
   private:
      vector            m_sigma,m_bc;
      matrix            m_vb;
      HARX*             m_obj;
      bool              m_centered,m_individual;
      double            m_epsilon;
   public:
                     CObjectiveJacFVec(void) {}
                    ~CObjectiveJacFVec(void)
        {
         m_obj = NULL;
        }
      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;
        }
      virtual void      Jac(double &x[],double &fi[],CMatrixDouble &jac,CObject &obj) override
        {
         vector temp;
         temp.Assign(x);
         vector r = m_obj.objective(temp,m_sigma,m_bc,m_vb);
         fi[0] = r[0];

         matrix mjac = m_obj.jacobian(temp,m_sigma,m_bc,m_vb,m_individual,m_centered,m_epsilon);
         jac = mjac;
        }
      virtual void      Jac(CRowDouble &x,CRowDouble &fi,CMatrixDouble &jac,CObject &obj) override
        {
         vector xv = x.ToVector();
         vector temp = m_obj.objective(xv,m_sigma,m_bc,m_vb);
         fi = temp;

         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;
   virtual double    _r2(vector &params) override
     {
      vector y = m_fit_y;
      bool constant = false;
      matrix x = m_fit_regressors;
      if(x.Cols())
         constant = (m_constant||implicit_constant(x));
      if(constant)
        {
         if(x.Cols() ==1)
            return 0.0;
         y = y - y.Mean();
        }
      double tss = y.Dot(y);
      if(tss<=0.)
         return AL_NaN;
      vector e = resids(params,EMPTY_VECTOR,EMPTY_MATRIX);
      return 1.-(e.Dot(e))/tss;
     }
   bool              _adjust_sample(long first_obs, long last_obs)
     {
      first_obs+=long(m_holdback);
      if(last_obs<= first_obs)
        {
         Print(__FUNCTION__, " invalid indices specified ");
         return false;
        }

      vector f = {double(first_obs),double(last_obs)};
      m_fit_indices = f;
      m_fit_y = adjust_y();
      m_fit_x = adjust_x();
      m_fit_regressors = (m_regressors.Rows())?np::sliceMatrixRows(m_regressors,long(m_fit_indices[0]),long(m_fit_indices[1])):matrix::Zeros(0,0);
      m_vp.start(long(first_obs));
      m_vp.stop(long(last_obs));

      return true;
     }
   vector            adjust_y(void)
     {
      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);
     }
   matrix            adjust_x(bool rows = true)
     {
      if(m_model_spec.exog_data.Rows())
        {
         if(rows)
            return np::sliceMatrixRows(m_model_spec.exog_data,long(m_fit_indices[0]),long(m_fit_indices[1]));
         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);
     }
   vector            _fit_no_arch_normal_errors_params(void)
     {
      vector y = m_fit_y;
      ulong nobs = y.Size();
      if(nobs<m_num_params)
        {
         Print(__FUNCTION__, " insufficient data ");
         return vector::Zeros(0);
        }
      matrix x = m_fit_regressors;
      if(x.Cols())
        {
         ResetLastError();
         matrix xpinv = np::pinv(x);
         if(xpinv.HasNan())
           {
            Print(__FUNCTION__, " PInv() returned invalid number ", GetLastError());
            return vector::Zeros(0);
           }
         return np::matmul(xpinv,y);
        }
      else
         return vector::Zeros(0);
     }
   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);
      if(x.Cols())
        {
         matrix xpinv = np::pinv(x);
         rparams = np::matmul(xpinv,y);
         matrix xt = (x.Transpose().MatMul(x))/double(y.Size());
         xpxi = xt.Inv();
         fitted = x.MatMul(rparams);
        }
      vector e = np::minus(y,fitted);
      double sigma2 = e.Dot(e)/double(y.Size());
      vector params = rparams;
      rparams.Resize(params.Size()+1);
      rparams[params.Size()] = sigma2;
      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.;
      matrix param_cov = matrix::Zeros(0,0);
      string nm;
      switch(cov_type)
        {
         case COVAR_CLASSIC:
            param_cov = sigma2 *(-1.*hessian);
            param_cov[m_num_params,m_num_params] = 2.*pow(sigma2,2.);
            param_cov/=double(y.Size());
            nm = "classic_ols";
            break;
         case COVAR_ROBUST:
           {
            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);
            scores.Col(pow(e,2.)-sigma2,scores.Cols()-1);
            matrix scorescov = scores.Transpose().MatMul(scores)/double(y.Size());
            param_cov = hessian.MatMul(scorescov);
            param_cov = param_cov.MatMul(hessian);
            param_cov/=double(y.Size());
            nm = "white";
           }
         break;
        }
      double r2 = _r2(params);
      vector resids_ = vector::Zeros(m_model_spec.observations.Size());
      resids_.Fill(AL_NaN);
      np::vectorCopy(resids_,e,long(m_fit_indices[0]),long(m_fit_indices[1]));
      vector vol = vector::Zeros(resids_.Size());
      vol.Fill(AL_NaN);
      double ss = sqrt(sigma2);
      for(ulong i = ulong(m_fit_indices[0]); i<ulong(m_fit_indices[1]); vol[i] = ss, ++i);
      double ll = _static_gaussian_loglikelihood(e);

      m_params = out.params = rparams;
      m_converged = (m_params.Size()>0);
      out.param_cov = param_cov;
      out.r2 = r2;
      out.resid = resids_;
      out.loglikelihood = ll;
      out.fit_indices[0] = long(m_fit_indices[0]);
      out.fit_indices[1] = long(m_fit_indices[1]);

      return out;
     }
   bool              _init_model(void)
     {
      ResetLastError();
      if(m_model_spec.mean_lags.Size())
        {
         if(!_reformat_lags())
            return false;
         m_maxlags = ulong(m_model_spec.mean_lags.Max());
         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;
              }
           }
        }
      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);
      if(m_constant)
         reg_constant = vector::Ones(nobs_orig);
      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[];
         if(!lagmat(ymat,lagarray,ulong(maxlag),TRIM_FORWARD,ORIGINAL_EX))
            return false;
         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);
            reg_lags.Col(vtemp,i);
           }
        }
      m_regressors = matrix::Zeros(nobs_orig,ulong(nobs_orig>0)*(ulong(m_constant)+reg_lags.Cols()+m_model_spec.exog_data.Cols()));
      if(m_constant && nobs_orig)
         m_regressors.Col(reg_constant,0);
      if(reg_lags.Cols())
         for(ulong i = 0; i<reg_lags.Cols(); ++i)
            m_regressors.Col(reg_lags.Col(i),i+ulong(m_constant));
      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());
      if(GetLastError()!=0)
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }
      return true;
     }
   virtual void      _scalechanged(void) override
     {
      if(!_init_model())
         Print(__FUNCTION__, " error ");
     }
   bool              _check_specification(void)
     {
      if(m_model_spec.exog_data.Rows())
        {
         if(m_model_spec.exog_data.Rows()!=m_model_spec.observations.Size())
           {
            Print(__FUNCTION__, " regressors shape does not match input data ");
            return false;
           }
        }
      return true;
     }
   bool              _reformat_lags(void)
     {
      if(!m_lags.Rows())
         return true;
      matrix lags = m_lags;
      if(lags.Rows()==1)
        {
         if(lags.Min()<=0)
           {
            Print(__FUNCTION__, " invalid lags input ");
            return false;
           }
         vector row = lags.Row(0);
         row = np::unique(row);
         matrix temp = matrix::Zeros(2,row.Size());
         if(!temp.Row(row,0) || !temp.Row(row,1))
            return false;
         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.;
           }
         else
            temp.Row(vector::Zeros(temp.Cols()),0);
         m_lags = temp;
        }
      else
        {
         vector mins = lags.Min(1);
         if(lags.Min()<=0 || mins[1]<mins[0])
           {
            Print(__FUNCTION__, " invalid lags input ");
            return false;
           }

         if(!np::reverseMatrixRows(lags) || !np::sort(lags,true))
            return false;

         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;
         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]));
         ulong rank = np::matrix_rank(testmat);
         if(rank != lags.Cols())
           {
            Print(__FUNCTION__," lags contains redundant entries ");
            return false;
           }
         if(m_rotated)
            Print(__FUNCTION__, " Rotation is not available for this model specification ");
         m_lags = lags;
        }
      return true;
     }
   bool              _reformat_forecast_x(ulong start, ulong horizon, matrix& x[], matrix& out[])
     {
      if(!x.Size())
        {
         if(!m_model_spec.exog_data.Rows())
            return true;
         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;
           }
      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 forizon of the forecasts ");
         return false;
        }

      if(x[0].Rows()!=(m_model_spec.observations.Size()-start))
        {
         Print(__FUNCTION__," the shape of x does not have the proper dimensions ");
         return false;
        }

      ArrayResize(out,x.Size());
      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));
        }
      else
        {
         for(uint i = 0; i<x.Size(); ++i)
            out[i] = x[i];
        }

      return true;
     }
   vector            _har_to_ar(vector &params)
     {
      if(!m_maxlags)
         return np::sliceVector(params,0,long(m_constant));
      vector har = np::sliceVector(params,long(m_constant));
      vector ar = vector::Zeros(m_maxlags);
      vector temp;
      for(ulong i = 0; i<m_lags.Cols(); ++i)
        {
         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]);
        }
      if(m_constant)
        {
         vector a(ar.Size()+1);
         a[0] = params[0];
         for(ulong i = 1; i<a.Size(); a[i] = ar[i-1], ++i);
         ar = a;
        }
      return ar;
     }
   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();
      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);
         else
            iv = initial_values;
      for(ulong i = 0; i<max_lag; y[i] = iv[i],++i);
      ulong kx = x.Cols();

      ulong ind;
      vector temp,col;
      for(ulong t = max_lag; t<nobs_and_burn; ++t)
        {
         ind = 0;
         if(m_constant)
           {
            y[t] = parameters[ind];
            ind+=1;
           }
         for(ulong lag = 0; lag<m_lags.Cols(); ++lag)
           {
            col = m_lags.Col(lag);
            temp  = np::sliceVector(y,long(t - (ulong)col[1]), long(t - ulong(col[0])));
            y[t] += parameters[ind] * temp.Mean();
            ind+=1;
           }
         for(ulong i = 0; i<kx; ++i)
           {
            y[t] += parameters[ind] * x[t,i];
            ind += 1;
           }
         y[t] += errs[t];
        }
      return y;
     }
   ArchModelResult   _fit_parameterless_model(ENUM_COVAR_TYPE cv, vector& backcast)
     {
      ArchModelResult out;

      vector y = m_fit_y;
      vector params = vector::Zeros(0);
      matrix param_cov = matrix::Zeros(0,0);
      ulong first_obs = (ulong)m_fit_indices[0];
      ulong last_obs = (ulong)m_fit_indices[1];
      vector resids_final = vector::Zeros(y.Size());
      if(!np::vectorCopy(resids_final,y,first_obs,last_obs))
        {
         Print(__FUNCTION__,  " failed vector copy ");
         return out;
        }
      return out;
      matrix var_bounds = m_vp.varianceBounds(y);
      vector vol = vector::Zeros(y.Size());
      m_vp.computeVariance(params,y,vol,backcast,var_bounds);
      vol = sqrt(vol);
      vector vol_final = vector::Zeros(y.Size());
      vol_final.Fill(AL_NaN);
      if(!np::vectorCopy(vol_final,vol,first_obs,last_obs))
        {
         Print(__FUNCTION__,  " failed vector copy ");
         return out;
        }
      double r2 = _r2(params);
      vector ll = -1.*_loglikelihood(params,pow(vol,2)*vector::Ones(ulong(m_fit_indices[1] - m_fit_indices[0])),backcast,var_bounds);

      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;

     }
   bool              _initialize(ArchParameters &vol_dist_params)
     {
      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;
        }

      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_model_spec = vol_dist_params;

      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;

        }

      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)
           {
            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;
            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);
      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;
      m_initialized = _init_model();
      if(!m_initialized)
         return false;

      m_num_params = num_params();

      if(CheckPointer(m_vp)==POINTER_DYNAMIC)
         delete m_vp;
      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;
         default:
            m_vp = new CConstantVariance(m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
            break;
        }

      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;
        }

      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)
                      );

      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();
     }
   bool                      initialize(ArchParameters& vol_dist_params)
     {
      return _initialize(vol_dist_params);
     }
   matrix                    get_x(void) { return m_model_spec.exog_data; }
   matrix                    get_regressors(void) { return m_regressors; }
   vector                    get_y(void) { return m_model_spec.observations; }

   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);
        }
     }

   virtual ulong             num_params(void)
     {
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return 0;
        }
      return m_regressors.Cols();
     }
   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;
     }
   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;
     }
   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;
        }
      ulong kx = x.Cols();
      ulong mc = 0;
      if(kx && x.Rows()<nobs+burn)
        {
         Print(__FUNCTION__, " x must have nobs + burn rows");
         return matrix::Zeros(0,0);
        }
      if(m_model_spec.mean_lags.Size())
         mc = ulong(m_constant)+m_lags.Cols()+kx+m_extra_simulation_params;

      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);
        }
      vector allparams[];
      if(!_parse_parameters(params,allparams))
         return matrix::Zeros(0,0);
      BootstrapRng rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
      matrix sim_data = m_vp.simulate(allparams[1],nobs+burn,rng,burn,initial_vals_vol.Size()?initial_vals_vol[0]:0.0);
      vector errors = sim_data.Row(0);
      errors = np::sliceVector(errors,long(burn));
      vector vol = sqrt(sim_data.Row(1));
      vol = np::sliceVector(vol,long(burn));
      vector y = _simulate_mean(allparams[0],x,errors,initial_vals,sim_data.Row(1));
      y = np::sliceVector(y,long(burn));
      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;
     }
   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);
     }
   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);
     }
   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;
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return out;
        }

      if(!m_fit_y.Size())
        {
         long from = 0;
         long to = long(m_model_spec.observations.Size());
         if(!_adjust_sample(from,to))
            return out;
        }
      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;
      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;
        }
      vector allparams[];
      if(!CArchModel::_parse_parameters(params.Size()?params:m_params,allparams))
        {
         Print(__FUNCTION__, " failed to parse the parameters ");
         return out;
        }
      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.);
      BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
      long vstart = MathMax(0,start_index - earliest);
      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);
      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);
      matrix expected_x[];
      if(!_reformat_forecast_x(ulong(start_index),horizon,x,expected_x))
        {
         Print(__FUNCTION__, " reformating error ");
         return out;
        }

      matrix mean_fcast = _ar_forecast(m_model_spec.observations,horizon,start_index,constant,dynp,expected_x,exog_p);
      if(!mean_fcast.Cols())
         return out;
      vector impulse = _ar_to_impulse(horizon,dynp);
      matrix longrun_var_fcasts = var_fcasts;
      matrix tempm;
      vector tempv,lrf;
      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);
        }
      matrix variance_paths[],mean_paths[];
      matrix shocks[],long_run_variance_paths[];
      if(method == FORECAST_BOOTSTRAP || method == FORECAST_SIMULATION)
        {
         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;
         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);
              }
           }
         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);
         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));
            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);
               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;
                 }
              }
           }
         for(uint i = 0; i<mean_paths.Size(); ++i)
            mean_paths[i] = np::sliceMatrixCols(mean_paths[i],long(m_maxlags));
        }

      return ArchForecast(start_index,mean_fcast,longrun_var_fcasts,var_fcasts,mean_paths,shocks,long_run_variance_paths,variance_paths);
     }
   ArchModelFixedResult fix(vector& params,long first_obs = 0, long last_obs = -1)
     {
      ArchModelFixedResult out;
      if(!is_initialized())
        {
         Print(__FUNCTION__, " object instance was not successfully initialized ");
         return out;
        }
      if(!_adjust_sample(first_obs,last_obs))
         return out;

      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;
      matrix vb = m_vp.varianceBounds(res);
      vector logl= _loglikelihood(params,sigma2,back_cast,vb)*(-1.0);
      vector parameters[];
      if(!_parse_parameters(params,parameters))
         return out;
      vector vol = vector::Zeros(res.Size());
      m_vp.computeVariance(parameters[1],res,vol,back_cast,vb);
      vol = sqrt(vol);
      first_obs = long(m_fit_indices[0]);
      last_obs = long(m_fit_indices[1]);
      vector res_final = vector::Zeros(res.Size());
      res_final.Fill(AL_NaN);
      np::vectorCopy(res_final,res,long(first_obs),long(last_obs));
      vector vol_final = vector::Zeros(vol.Size());
      vol_final.Fill(AL_NaN);
      np::vectorCopy(vol_final,vol,long(first_obs),long(last_obs));

      out.params = params;
      out.resid = res_final;
      out.conditional_volatility = vol_final;
      out.loglikelihood = logl[0];

      return out;
     }
#ifdef __SLSQP__
   ArchModelResult   fit(double tol = 1e-4,uint maxits = 0,ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1)
     {
      return fit(EMPTY_VECTOR,EMPTY_VECTOR,tol,maxits,cov_type,first,last);
     }
   ArchModelResult   fit(vector& startingvalues,vector& backcast,double tol = 1e-4, uint maxits = 0,ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1)
     {
      ArchModelResult out;
      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;
        }

      vector offsets(3);

      offsets[0] = (double)m_num_params;
      offsets[1] = (double)m_vp.numParams();
      offsets[2] = (double)m_distribution.numParams();

      double total_params = offsets.Sum();

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

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

      bool hascloseform = (offsets[2] == 0.0 && m_vp.closeForm() && m_vp.volatilityprocess()==VOL_CONST);

      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));
        }
      vector bcast;
      if(!backcast.Size())
         bcast = m_vp.backCast(res);
      else
         bcast = m_vp.backCastTransform(backcast);
      if(hascloseform)
         return _fit_no_arch_normal_errors(cov_type);
      if(total_params == 0.0)
         return _fit_parameterless_model(cov_type,bcast);
      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);

      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);
           }
        }

      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()));
           }
        }

      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;

      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;
        }

      CSlsqp slsqp_minim;
      slsqp_minim.SetInequalityConstraints(obfunc_constraints);

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

      OptimizeResult minim_result = slsqp_minim.Minimize(obfunc);

      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;
        }
      vector allparams[];
      //--- loglikelihood for estimated parameters
      vector llh = _loglikelihood(m_params,sigma2,bcast,vb);
      //Print(__FUNCTION__, " loglikelihood ", llh[0]);
      //---
      _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,l_obs;
      f_obs = (long)m_fit_indices[0];
      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
   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);
     }
   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;
      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;
        }

      vector offsets(3);

      offsets[0] = (double)m_num_params;
      offsets[1] = (double)m_vp.numParams();
      offsets[2] = (double)m_distribution.numParams();

      double total_params = offsets.Sum();

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

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

      bool hascloseform = (offsets[2] == 0.0 && m_vp.closeForm() && m_vp.volatilityprocess()==VOL_CONST);

      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));
        }
      vector bcast;
      if(!backcast.Size())
         bcast = m_vp.backCast(res);
      else
         bcast = m_vp.backCastTransform(backcast);
      if(hascloseform)
         return _fit_no_arch_normal_errors(cov_type);
      if(total_params == 0.0)
         return _fit_parameterless_model(cov_type,bcast);
      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);

      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)+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);
           }
        }

      a.Col(b,ulong(total_params));
      b = vector::Ones(a.Rows());

      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()));
           }
        }

      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;
        }

      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;


      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);

      CMinNLC::MinNLCSetBC(optim_state,bl,bu);

      CRowInt intones;
      intones.Resize(_constraints.Rows());
      intones.Fill(1);

      CMinNLC::MinNLCSetLC(optim_state,_constraints,intones,intones.Size());

      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 other objective function arguments
      CObjectiveJacFVec fvec;
      fvec.SetParams(sigma2,bcast,vb,this);
      //---
      CNDimensional_Rep f_rep;
      CObject object;
      //---optimize
      //CMinNLC::MinNLCOP(optim_state,fvec,f_rep,object);
      //--- check
      if(!CAp::Assert(GetPointer(fvec)!=NULL,"ALGLIB: error in 'minnlcoptimize()' (fvec is null)"))
        {
         m_converged = false;
         return out;
        }

      //--- cycle
      while(CMinNLC::MinNLCIteration(optim_state))
        {
         if(optim_state.m_needfij)
           {
            fvec.Jac(optim_state.m_x,optim_state.m_fi,optim_state.m_j,object);
            continue;
           }
         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;
        }
      //---get results of optimization
      CMinNLC::MinNLCResults(optim_state,solution,rep);
      //---check for failed convergence
      if(rep.m_terminationtype<0)
        {
         Print(__FUNCTION__," failed convergence ");
         m_converged = false;
         return out;
        }
      else
         m_converged = true;
      //---
      vector allparams[];
      //---
      m_params = solution.ToVector();
      //--- loglikelihood for estimated parameters
      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,l_obs;
      f_obs = (long)m_fit_indices[0];
      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
   bool              eval_loglikelihood(vector& params, double &out_loglikelihood)
     {

      long first = 0;
      long last = -1;

      vector startingvalues = EMPTY_VECTOR;
      vector backcast = EMPTY_VECTOR;

      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;
        }

      vector offsets(3);

      offsets[0] = (double)m_num_params;
      offsets[1] = (double)m_vp.numParams();
      offsets[2] = (double)m_distribution.numParams();

      double total_params = offsets.Sum();

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

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

      bool hascloseform = (offsets[2] == 0.0 && m_vp.closeForm() && m_vp.volatilityprocess()==VOL_CONST);

      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));
        }
      vector bcast;
      if(!backcast.Size())
         bcast = m_vp.backCast(res);
      else
         bcast = m_vp.backCastTransform(backcast);
      if(hascloseform)
         return false;
      if(total_params == 0.0)
         return false;
      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);

      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)+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);
           }
        }

      a.Col(b,ulong(total_params));
      b = vector::Ones(a.Rows());

      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()));
           }
        }

      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;
        }

      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;
      //--- set up other objective function arguments
#ifdef __SLSQP__
      CObjectiveF fvec;
      fvec.setFuncParams(sigma2,bcast,vb,this);
      out_loglikelihood = fvec.orig_fun(params);
#else
      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                                                        |
//+------------------------------------------------------------------+
class HAR: public HARX
  {
public:
                     HAR(void)
     {
      m_model_spec.mean_model_type = MEAN_HAR;
     }
                     HAR(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_HAR;
      initialize(model_specification);
     }
                    ~HAR(void)
     {
      deinitialize();
     }
  };
//+------------------------------------------------------------------+
//| The constant mean model                                          |
//+------------------------------------------------------------------+
class ConstantMean: public HARX
  {
public:
                     ConstantMean(void)
     {
      m_model_spec.mean_model_type = MEAN_CONSTANT;
     }
                     ConstantMean(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_CONSTANT;
      initialize(model_specification);
     }
                    ~ConstantMean(void)
     {
      deinitialize();
     }
   virtual ulong     num_params(void) override
     {
      return 1;
     }
   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.Size())
         return np::minus(m_fit_y,params);
      return np::minus(y,params);
     }
   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;
        }
      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);
        }

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

      BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
      matrix sim_values = m_vp.simulate(allparams[1],nobs+burn,Rng,burn,initial_vals_vol.Size()?initial_vals_vol[0]:0.0);

      vector ers = sim_values.Row(0);
      vector y = ers+allparams[0];
      vector vol = sqrt(sim_values.Row(1));
      y = np::sliceVector(y,long(burn));
      vol = np::sliceVector(vol,long(burn));
      ers = np::sliceVector(ers,long(burn));

      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                                               |
//+------------------------------------------------------------------+
class ZeroMean: public HARX
  {
public:
                     ZeroMean(void)
     {
      m_model_spec.mean_model_type = MEAN_ZERO;
     }
                     ZeroMean(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_ZERO;
      initialize(model_specification);
     }
                    ~ZeroMean(void)
     {
      deinitialize();
     }
   virtual ulong     num_params(void) override
     {
      return 0;
     }
   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.Size())
         return m_fit_y;
      return y;
     }
   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;
        }
      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);
        }

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

      BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
      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;
      vector vol = sqrt(sim_values.Row(1));
      y = np::sliceVector(y,long(burn));
      vol = np::sliceVector(vol,long(burn));
      ers = np::sliceVector(ers,long(burn));

      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                                                     |
//+------------------------------------------------------------------+
class AR: public HARX
  {
public:
                     AR(void)
     {
      m_model_spec.mean_model_type = MEAN_AR;
     }
                     AR(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_AR;
      initialize(model_specification);
     }
                    ~AR(void)
     {
      deinitialize();
     }
  };
//+------------------------------------------------------------------+
//| The AR-X model                                                   |
//+------------------------------------------------------------------+
class ARX: public HARX
  {
public:
                     ARX(void)
     {
      m_model_spec.mean_model_type = MEAN_ARX;
     }
                     ARX(ArchParameters& model_specification)
     {
      model_specification.mean_model_type = MEAN_ARX;
      initialize(model_specification);
     }
                    ~ARX(void)
     {
      deinitialize();
     }
  };
//+------------------------------------------------------------------