//+------------------------------------------------------------------+
//|                                                   volatility.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "recursions.mqh"
#include "distribution.mqh"
//---
//+------------------------------------------------------------------+
//| Container for variance forecasts                                 |
//+------------------------------------------------------------------+
struct VarianceForecast
  {
   matrix            forecasts;
   matrix            forecastpaths[];
   matrix            shocks[];

                     VarianceForecast(void)
     {
      forecasts = matrix::Zeros(0,0);
     }
                     VarianceForecast(matrix &_forecasts, matrix &_forecastpaths[],matrix &_shocks[])
     {
      forecasts = _forecasts;
      ArrayResize(forecastpaths,_forecastpaths.Size());
      ArrayResize(shocks,_shocks.Size());
      for(uint i = 0 ; i<shocks.Size(); shocks[i] = _shocks[i], ++i);
      for(uint i = 0 ; i<forecastpaths.Size(); forecastpaths[i] = _forecastpaths[i], ++i);
     }
                     VarianceForecast(VarianceForecast &other)
     {
      forecasts = other.forecasts;
      ArrayResize(forecastpaths,other.forecastpaths.Size());
      ArrayResize(shocks,other.shocks.Size());
      for(uint i = 0 ; i<shocks.Size(); shocks[i] = other.shocks[i], ++i);
      for(uint i = 0 ; i<forecastpaths.Size(); forecastpaths[i] = other.forecastpaths[i], ++i);
     }
   void              operator=(VarianceForecast &other)
     {
      forecasts = other.forecasts;
      ArrayResize(forecastpaths,other.forecastpaths.Size());
      ArrayResize(shocks,other.shocks.Size());
      for(uint i = 0 ; i<shocks.Size(); shocks[i] = other.shocks[i], ++i);
      for(uint i = 0 ; i<forecastpaths.Size(); forecastpaths[i] = other.forecastpaths[i], ++i);
     }
  };
//+------------------------------------------------------------------+
//| base class volatility process                                    |
//+------------------------------------------------------------------+
class CVolatilityProcess
  {
protected:
   bool              m_updateable;
   bool              m_initialized;
   ulong             m_num_params;
   bool              m_closedform;
   ulong             m_bootstrap_obs;
   long              m_start;
   long              m_stop;
   CNormal           m_normal;
   VolatilityUpdater *m_volupdater;
   string            m_name;
   int               m_seed;
   ENUM_VOLATILITY_MODEL m_model;
   virtual double    _update(ulong index, vector& parameters, vector &resids, vector &sigma2, vector &backcast, vector &var_bounds)
     {
      return EMPTY_VALUE;
     }
   virtual bool      _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon)
     {
      return false;
     }
   virtual bool      _onestepforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, vector&out_sigma, matrix&out_forecasts)
     {
      ulong t = resids.Size();
      vector _resids=resids;
      _resids.Resize(t+1);
      _resids[t] = 0.;
      matrix _vb = varbounds;
      _vb.Resize(varbounds.Rows()+1,varbounds.Cols());
      vector temp = {0.,double("inf")};
      _vb.Row(temp,varbounds.Rows());
      vector sigma2 = vector::Zeros(t+1);
      computeVariance(parameters, _resids,sigma2,backcast,_vb);
      matrix forecasts = matrix::Zeros(t-_start,horizon);
      temp = np::sliceVector(sigma2,_start+1);
      forecasts.Col(temp,0);
      out_sigma = np::sliceVector(sigma2,0,-1);
      out_forecasts = forecasts;
      return true;
     }
   virtual double    _gaussianloglikelihood(vector &parameters,  vector& resids, vector&backcast, matrix& varbounds)
     {
      vector sigma2 = vector::Zeros(resids.Size());
      computeVariance(parameters,resids,sigma2,backcast,varbounds);
      vector empty  = vector::Zeros(0);
      vector v = m_normal.loglikelihood(empty,resids,sigma2);
      return v[0];
     }
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon)
     {
      return VarianceForecast();
     }
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      return VarianceForecast();
     }
   virtual VarianceForecast _bootstrapforecast(vector& parameters, vector& resids, vector& backcast, matrix &varbounds, ulong _start,ulong horizon, ulong simulations, int seed = 0)
     {
      VarianceForecast out;
      vector sigma2 = vector::Zeros(resids.Size());
      computeVariance(parameters,resids,sigma2,backcast,varbounds);
      vector std_resid = resids/sqrt(sigma2);
      if(_start<m_bootstrap_obs)
        {
         Print(__FUNCTION__, " start must include more than ",m_bootstrap_obs," observations");
         return out;
        }

      BootstrapRng Rng(std_resid,_start,seed);
      out = _simulationforecast(parameters,resids,backcast,varbounds,_start,horizon,simulations,Rng);
      return out;
     }
   virtual bool      _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE,bool updateable = true, bool closedform = false,ulong nparams = 0, string name = NULL,int seed = 0,ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0,long start=0, long stop=-1,ulong bootstrap_obs=100)
     {
      m_updateable = updateable;
      m_num_params = nparams;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;
      m_model = volmodel;
      return m_normal.initialize(vector::Zeros(0),seed);
     }
public:
                     CVolatilityProcess(void):m_updateable(true),
                     m_num_params(0),
                     m_closedform(false),
                     m_bootstrap_obs(100),
                     m_start(0),
                     m_stop(-1),
                     m_initialized(false),
                     m_name(NULL),
                     m_volupdater(NULL),
                     m_model(WRONG_VALUE)
     {
     }
                     CVolatilityProcess(CVolatilityProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closeForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
     }

                    ~CVolatilityProcess(void)
     {
      if(!CheckPointer(m_volupdater)==POINTER_DYNAMIC)
         delete m_volupdater;
      m_volupdater = NULL;
     }
   void              operator=(CVolatilityProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closeForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
     }
   ENUM_VOLATILITY_MODEL volatilityprocess(void) { return m_model; }
   ulong             bootstraps(void) { return m_bootstrap_obs; }
   void              bootstraps(ulong bts)
     {
      if(bts>10)
         m_bootstrap_obs = bts;
      else
         return;
     }
   string            name(void) { return m_name; }
   bool              is_initialized(void) { return m_initialized; }
   virtual long      start(void) { return m_start; }
   virtual long      stop(void) { return m_stop; }
   virtual void      start(long _start) { m_start = _start; }
   virtual void      stop(long _stop) { m_stop = _stop; }
   virtual ulong     numParams(void) { return m_num_params; }
   virtual bool      upDateable(void) { return  m_updateable; }
   virtual bool      closeForm(void) { return m_closedform; }

   virtual VolatilityUpdater* volUpdater(void)
     {
      return m_volupdater;
     }
   virtual double    upDate(ulong _index, vector& parameters, vector& resids, vector& sigma2, vector& backcast, vector& varbounds)
     {
      return EMPTY_VALUE;
     }
   virtual matrix    varianceBounds(vector& resids, double power = 2.)
     {
      ulong nobs = resids.Size();
      ulong tau = MathMin(75,nobs);
      vector w = vector::Zeros(tau);
      for(ulong i = 0; i<tau; w[i] = pow(0.94,double(i)), ++i);
      w = w/w.Sum();
      double initial_value = w.Dot(pow(np::sliceVector(resids,0,long(tau)),2.));

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

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

      double var = resids.Var();
      double min_upper_bound = 1.+(pow(resids,2.)).Max();
      double lower_bound = var/1.e8;
      double upper_bound = 1.e7*(1.+(pow(resids,2.)).Max());
      vector col0 = varbounds.Col(0);
      vector col1 = varbounds.Col(1);

      if(!col0.Clip(lower_bound,DBL_MAX) || !col1.Clip(min_upper_bound,upper_bound) || !varbounds.Col(col0,0) || !varbounds.Col(col1,1))
        {
         Print(__FUNCTION__, " error ", __LINE__, " ", GetLastError());
         return matrix::Zeros(0,0);
        }

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

      return varbounds;
     }
   virtual vector    startingValues(vector& resids)
     {
      return vector::Zeros(0);
     }
   virtual vector    backCast(vector& resids)
     {
      ulong tau = MathMin(75,resids.Size());
      vector w = vector::Zeros(tau);
      for(ulong i = 0; i<tau; w[i] = pow(0.94,double(i)), ++i);
      w = w/w.Sum();
      vector out(1);
      out[0] = (pow(np::sliceVector(resids,0,long(tau)),2.0)*w).Sum();
      return out;
     }
   virtual vector    backCastTransform(vector& backcast)
     {
      if(backcast.Min()<0)
        {
         Print(__FUNCTION__, " User backcast value must be striclty positive ");
         return vector::Zeros(0);
        }
      return backcast;
     }
   virtual matrix    bounds(vector& resids)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds)
     {
      return vector::Zeros(0);
     }
   virtual Constraints constraints(void)
     {
      return Constraints();
     }
   VarianceForecast  forecast(vector &parameters, vector& resids, vector& _backcast, matrix& var_bounds, BootstrapRng &rng,int seed = 0, ulong _start=0, ulong horizon = 1, ENUM_FORECAST_METHOD ForecastingMethod = FORECAST_ANALYTIC, ulong simulations = 1000)
     {
      VarianceForecast out;
      if(!horizon)
        {
         Print(__FUNCTION__, " horizon must be >= 1");
         return out;
        }

      if(!_check_forecasting_method(ForecastingMethod,horizon))
        {
         Print(__FUNCTION__, " The specified forecasting method is not supported for this specific model ");
         return out;
        }

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

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

      return out;
     }

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

public:
                     CConstantVariance(void)
     {
      m_initialized = CVolatilityProcess::_initialize(VOL_CONST,true,true,1,"Constant Variance");
     }
                     CConstantVariance(int seed=0,ulong nbootstraps=100)
     {
      m_initialized = CVolatilityProcess::_initialize(VOL_CONST,true,true,1,"Constant Variance",seed,0,0,0,0.0,0,-1,nbootstraps);
     }
                    ~CConstantVariance(void)
     {
     }
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) override
     {
      _sigma2.Fill(parameters[0]);
      return _sigma2;
     }
   virtual vector    startingValues(vector& resids) override
     {
      vector out(1);
      out[0] = resids.Var();
      return out;
     }
   virtual Constraints constraints(void) override
     {
      Constraints out;
      out._one = matrix::Ones(1,1);
      out._two = vector::Zeros(1);
      return out;
     }
   virtual vector    backCastTransform(vector& backcast) override
     {
      backcast=backCastTransform(backcast);
      return backcast;
     }
   virtual vector    backCast(vector& resids) override
     {
      vector out(1);
      out[0] = resids.Var();
      return out;
     }
   virtual matrix    bounds(vector& resids) override
     {
      double v = resids.Var();
      matrix out = matrix::Zeros(1,2);
      out[0,0] = v/100000.;
      out[0,1] = 10.0 * (v + pow(resids.Mean(), 2.0));
      return out;
     }
   virtual matrix       simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      vector errv = rng.rng(_nobs+burn);
      if(!errv.Size())
       {
        matrix container = rng.rng(_nobs+burn,1);
        errv = container.Col(0);
       }
      vector sigma2 = vector::Ones(_nobs+burn) * parameters[0];
      vector data =  sqrt(sigma2) * errv;
      sigma2 = np::sliceVector(sigma2,long(burn));
      data = np::sliceVector(data,long(burn));
      matrix out(data.Size(),2);
      out.Col(data,0);
      out.Col(sigma2,1);
      return out;
     }
   virtual string        parameterNames(void) override
     {
      return "sigma2";
     }
  };
//+------------------------------------------------------------------+
//|  GARCH and related model estimation                              |
//+------------------------------------------------------------------+
class CGarchProcess: public CVolatilityProcess
  {
protected:
   ulong             m_p,m_o,m_q;
   double            m_power;
   matrix            _product(vector& one, vector& two, vector& three)
     {
      ulong size = one.Size();
      ulong rows = 3;
      matrix out(ulong(pow(size,rows)),rows);
      ulong shift = 0;
      for(ulong i = 0; i<size; ++i)
        {
         for(ulong j = 0; j<size; ++j)
           {
            for(ulong k = 0; k<size; ++k)
              {
               out[shift,0] = one[i];
               out[shift,1] = two[j];
               out[shift,2] = three[k];
               ++shift;
              }
           }
        }
      return out;
     }
   virtual bool      _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override
     {
      if(horizon == 1)
         return true;
      if(method == FORECAST_ANALYTIC && m_power!=2.)
        {
         Print(__FUNCTION__, " Analytic forecasts not available for horizon > 1 when power != 2");
         return false;
        }
      return true;
     }
   virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override
     {
      VarianceForecast out;
      ulong t = resids.Size();
      vector sigma2;
      matrix forecasts;
      _onestepforecast(parameters,resids,backcast,varbounds,_start,horizon,sigma2,forecasts);
      if(horizon == 1)
        {
         out.forecasts = forecasts;
         return out;
        }
      double omega = parameters[0];
      vector alpha = np::sliceVector(parameters,1,long(m_p+1));
      vector gamma = np::sliceVector(parameters,long(m_p+1),long(m_p+m_o+1));
      vector beta = np::sliceVector(parameters,long(m_p+m_o+1));
      ulong m = MathMax(MathMax(m_p,m_o),m_q);
      vector _resids,_asym_resids,_sigma2,temp;
      _resids = _asym_resids = _sigma2 = vector::Zeros(m+horizon);

      for(ulong i = ulong(_start); i<t; ++i)
        {
         if(i-m+1>=0)
           {
            temp = np::sliceVector(resids,long(i-m+1),long(i+1));
            np::vectorCopy(_resids,temp,0,long(m));
            temp = np::sliceVector(_resids,0,long(m));
            temp*=np::whereVectorIsLt(temp,0.0);
            np::vectorCopy(_asym_resids,temp,0,long(m));
            temp = np::sliceVector(sigma2,long(i-m+1),long(i+1));
            np::vectorCopy(_sigma2,temp,0,long(m));
           }
         else
           {
            np::fillVector(_resids,sqrt(backcast[0]),0,long(m-i-1));
            temp = np::sliceVector(resids,0,long(i+1));
            np::vectorCopy(_resids,temp,long(m-i-1),long(m));
            _asym_resids = _resids*(np::whereVectorIsLt(_resids,0.0));
            np::fillVector(_asym_resids,sqrt(0.5*backcast[0]),0,long(m-i-1));
            np::fillVector(_sigma2,backcast[0],0,long(m));
            temp = np::sliceVector(sigma2,0,long(i+1));
            np::vectorCopy(_sigma2,temp,long(m-i-10), long(m));
           }
         for(ulong h = 0; h<horizon; ++h)
           {
            ulong floc = i - ulong(_start);
            forecasts[floc,h] = omega;
            ulong sloc = h+m-1;
            for(ulong j = 0; j<m_p; ++j)
               forecasts[floc,h] += alpha[j]*pow(_resids[sloc-j],2.);
            for(ulong j = 0; j<m_o; ++j)
               forecasts[floc,h] += (gamma[j]*pow(_asym_resids[sloc-j],2.));
            for(ulong j = 0; j<m_q; ++j)
               forecasts[floc,h] += beta[j]*_sigma2[sloc-j];
            _resids[h+m] = sqrt(forecasts[floc,h]);
            _asym_resids[h+m] = sqrt(0.5*forecasts[floc,h]);
            _sigma2[h+m] = forecasts[floc,h];
           }
        }

      out.forecasts = forecasts;
      return out;
     }
   virtual bool      _simulatepaths(ulong m, vector& parameters, ulong horizon, matrix& std_shocks, matrix& scaled_forecast_paths, matrix& scaled_shock, matrix& asym_scaled_shock, vector& mean_fpaths, matrix& fpaths, matrix& shocks)
     {
      double omega = parameters[0];
      vector alpha = ((m_p+1)>1)?np::sliceVector(parameters,1,long(m_p+1)):EMPTY_VECTOR;
      vector gamma = ((m_p+1)<(m_p+m_o+1))?np::sliceVector(parameters,long(m_p+1),long(m_p+m_o+1)):EMPTY_VECTOR;
      vector beta = (parameters.Size()>(m_p+m_o+1))?np::sliceVector(parameters,long(m_p+m_o+1)):EMPTY_VECTOR;
      matrix shock = scaled_forecast_paths;
      vector temp;
      for(ulong h = 0; h<horizon; ++h)
        {
         ulong loc = h+m-1;
         temp = vector::Zeros(scaled_forecast_paths.Rows());
         temp.Fill(omega);
         scaled_forecast_paths.Col(temp,h+m);
         for(ulong j = 0; j<m_p; ++j)
           {
            temp = scaled_forecast_paths.Col(h+m);
            temp+=alpha[j]*scaled_shock.Col(loc-j);
            scaled_forecast_paths.Col(temp,h+m);
           }
         for(ulong j = 0; j<m_o; ++j)
           {
            temp = scaled_forecast_paths.Col(h+m);
            temp+=gamma[j]*asym_scaled_shock.Col(loc-j);
            scaled_forecast_paths.Col(temp,h+m);
           }
         for(ulong j = 0; j<m_q; ++j)
           {
            temp = scaled_forecast_paths.Col(h+m);
            temp+=beta[j]*scaled_forecast_paths.Col(loc-j);
            scaled_forecast_paths.Col(temp,h+m);
           }
         temp = std_shocks.Col(h)*pow(scaled_forecast_paths.Col(h+m),1./m_power);
         shock.Col(temp,h+m);
         vector lt_zero = np::whereVectorIsLt(shock.Col(h+m),0.);
         scaled_shock.Col(pow(fabs(shock.Col(h+m)),m_power),h+m);
         asym_scaled_shock.Col(scaled_shock.Col(h+m)*lt_zero,h+m);
        }

      fpaths = np::sliceMatrixCols(scaled_forecast_paths,long(m));
      fpaths = pow(fpaths,2./m_power);
      mean_fpaths = fpaths.Mean(0);
      shocks = np::sliceMatrixCols(shock,long(m));

      return true;
     }
   virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng)
     {
      VarianceForecast out;
      vector sigma2;
      matrix forecasts;
      _onestepforecast(parameters,resids,backcast,varbounds,_start,horizon,sigma2,forecasts);
      long t = (long)resids.Size();
      ArrayResize(out.forecastpaths, int(t-_start));
      ArrayResize(out.shocks,int(t-_start));
      ulong m = MathMax(MathMax(m_p,m_o),m_q);
      matrix scaled_forecast_paths = matrix::Zeros(simulations,m+horizon);
      matrix scaled_shock = matrix::Zeros(simulations, m+horizon);
      matrix asym_scaled_shock = matrix::Zeros(simulations,m+horizon);
      matrix temp_m;
      vector temp_v,mask;
      ulong count = 0;
      for(ulong i = ulong(_start); i<ulong(t); ++i)
        {
         matrix std_shocks = rng.rng(simulations,horizon);
         if(i-m<0)
           {
            np::fillMatrix(scaled_forecast_paths,pow(backcast[0],m_power/2.),BEGIN,END,STEP,BEGIN,long(m));
            np::fillMatrix(scaled_shock,pow(backcast[0],m_power/2.),BEGIN,END,STEP,BEGIN,long(m));
            np::fillMatrix(asym_scaled_shock,0.5*pow(backcast[0],m_power/2.),BEGIN,END,STEP,BEGIN,long(m));

            count = i+1;

            temp_v = np::sliceVector(sigma2,0,long(count));
            temp_v = pow(temp_v,m_power/2.);
            temp_m = np::repeat_vector_as_rows_cols(temp_v,scaled_forecast_paths.Rows());
            np::matrixCopy(scaled_forecast_paths,temp_m,BEGIN,END,STEP,long(m-count),long(m));

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

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

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

            temp_m = np::sliceMatrixCols(scaled_shock,0,long(m));
            temp_v = np::sliceVector(resids,long(i-m+1),long(i+1));
            mask = np::whereVectorIsLt(temp_v,0.0);
            temp_m = np::multiply(temp_m,mask);
            np::matrixCopy(asym_scaled_shock,temp_m,BEGIN,END,STEP,BEGIN,long(m));
           }
         vector f;
         matrix p, s;
         _simulatepaths(m,parameters,horizon,std_shocks,scaled_forecast_paths,scaled_shock,asym_scaled_shock,f,p,s);
         ulong loc = i - long(_start);
         out.forecasts.Row(f,loc);
         out.forecastpaths[loc] = p;
         out.shocks[loc] = s;
        }
      return out;
     }
   virtual bool      _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE,bool updateable = true, bool closedform = false,ulong nparams = 0,string name=NULL, int seed = 0,ulong p = 1, ulong o = 0, ulong q = 1, double power = 2.,long start=0, long stop=-1,ulong bootstrap_obs=100) override
     {
      m_model = volmodel;
      m_updateable = updateable;
      m_num_params = nparams;
      m_closedform = closedform;
      bootstraps(bootstrap_obs);
      m_start = start;
      m_stop = stop;
      m_name = name;
      m_seed = seed;

      if(m_normal.randomState()!=uint(m_seed))
         m_normal.initialize(vector::Zeros(0),seed);

      m_p = p;
      m_o = o;
      m_q = q;
      m_power = power;

      if(m_power<=0.0)
       {
        Print(__FUNCTION__, " invalid value for power variable ");
        return false;
       }
      switch(m_model)
        {
         case VOL_ARCH:
         case VOL_AVARCH:
            if(!m_p)
               m_p = 1;
            break;
         case VOL_GARCH:
         case VOL_AVGARCH:
            if(!m_p)
               m_p = 1;
            if(!m_q)
               m_q = 1;
            break;
         case VOL_GJR_GARCH:
         case VOL_TARCH:
            if(!m_p)
               m_p = 1;
            if(!m_q)
               m_q = 1;
            if(!m_o)
               m_o=1;
            break;
        }

      m_num_params = 1 + m_p + m_o + m_q;

      return (true);
     }
public:
                     CGarchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_GARCH,true,false,0,"GARCH",0,1,0,1,2.0);
     }

                     CGarchProcess(ulong p, ulong q, int seed = 0, ulong nbootstraps=100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_GARCH,true,false,0,"GARCH",seed,p,0,q,2.0,0,-1,nbootstraps);
     }

                     CGarchProcess(CGarchProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closeForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
      m_q = other.get_q();
      m_p = other.get_p();
      m_o = other.get_o();
      m_power = other.get_power();
     }

                    ~CGarchProcess(void)
     {
     }

   void              operator=(CGarchProcess& other)
     {
      m_updateable = other.upDateable();
      m_num_params = other.numParams();
      m_closedform = other.closeForm();
      m_bootstrap_obs = other.bootstraps();
      m_start = other.start();
      m_stop = other.stop();
      m_volupdater = other.volUpdater();
      m_name = other.name();
      m_model = other.volatilityprocess();
      m_q = other.get_q();
      m_p = other.get_p();
      m_o = other.get_o();
      m_power = other.get_power();
     }
   virtual ulong     get_p(void) { return m_p; }
   virtual ulong     get_o(void) { return m_o; }
   virtual ulong     get_q(void) { return m_q; }
   virtual double    get_power(void) { return m_power; }
   virtual matrix    varianceBounds(vector& resids, double power = 2.) override
     {
      return CVolatilityProcess::varianceBounds(resids,m_power);
     }
   virtual vector    computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) override
     {
      vector fresids = pow(MathAbs(resids),m_power);
      vector sresids(resids.Size());
      for(ulong i = 0; i<sresids.Size(); sresids[i] = np::sign(resids[i]), ++i);
      _sigma2 = garch_recursion(parameters,fresids,sresids,_sigma2,m_p,m_o,m_q,resids.Size(),_backcast[0],varbounds);
      double invp = 2./m_power;
      return pow(_sigma2,invp);
     }
   virtual vector    startingValues(vector& resids) override
     {
      vector alphas = {0.01, 0.05, 0.1, 0.2};
      vector gammas = alphas;
      vector abg = {0.5, 0.7, 0.9, 0.98};
      matrix abgs = _product(alphas,gammas,abg);
      double target = (pow(fabs(resids),m_power)).Mean();
      double scale = (pow(resids,2.)).Mean()/pow(target,2./m_power);
      target*=pow(scale,m_power/2.);
      vector svs[];
      ArrayResize(svs,int(abgs.Rows()));
      matrix vb = varianceBounds(resids);
      vector llfs = vector::Zeros(abgs.Rows());
      vector bc = backCast(resids);
      for(uint i = 0; i<svs.Size(); ++i)
        {
         vector row = abgs.Row(i);
         vector sv = (1.-row[2])*target*vector::Ones(m_p+m_o+m_q+1);
         if(m_p>0)
           {
            np::fillVector(sv,row[0]/double(m_p),long(1),long(1+m_p));
            row[2] -= row[0];
           }
         if(m_o>0)
           {
            np::fillVector(sv,row[1]/double(m_o),long(1+m_p),long(1+m_p+m_o));
            row[2] -= row[1]/2.;
           }
         if(m_q>0)
            np::fillVector(sv,row[2]/double(m_q),long(1+m_p+m_o),long(1+m_p+m_o+m_q));
         svs[i]=sv;
         llfs[i] = _gaussianloglikelihood(sv,resids,bc,vb);
        }
      ulong loc = llfs.ArgMax();
      return svs[loc];
     }
   virtual Constraints constraints(void) override
     {
      Constraints out;
      ulong k_arch = m_p+m_o+m_q;
      out._one = matrix::Zeros(k_arch+2,k_arch+1);
      for(ulong i = 0; i<(k_arch+1); ++i)
         out._one[i,i] = 1.;
      for(ulong i = 0; i<(m_o); ++i)
         if(i<m_p)
            out._one[i+m_p+1,i+1] = 1.;
      for(ulong i = 1; i<(k_arch+1); ++i)
         out._one[k_arch+1,i] = -1.;
      for(ulong i = m_p+1; i<(m_p+m_o+1); ++i)
         out._one[k_arch+1,i] = -0.5;
      out._two = vector::Zeros(k_arch+2);
      out._two[k_arch+1] = -1.;
      return out;
     }
   virtual vector    backCastTransform(vector& backcast) override
     {
      backcast= CVolatilityProcess::backCastTransform(backcast);
      if(backcast.Size())
         return pow(sqrt(backcast),m_power);
      else
         return vector::Zeros(0);
     }
   virtual vector    backCast(vector& resids) override
     {
      ulong tau = MathMin(75,resids.Size());
      vector out = np::arange(tau);
      for(ulong i = 0; i<out.Size(); out[i] = pow(0.94,out[i]), ++i);
      out = out/out.Sum();
      vector temp = np::sliceVector(resids,0,long(tau));
      out = pow(MathAbs(temp),m_power)*out;
      vector bc(1);
      bc[0] = out.Sum();
      return bc;
     }
   virtual matrix    bounds(vector& resids) override
     {
      double v = (pow(MathAbs(resids),m_power)).Mean();
      matrix out = matrix::Zeros(1+m_p+m_o+m_q,2);
      out[0,0] = 1e-8 * v;
      out[0,1] = 10.0 * (v);
      ulong from = 1;
      for(ulong i = 0; i<m_p; ++i)
        {
         out[from,0] = 0.0;
         out[from++,1] = 1.;
        }

      for(ulong i = 0; i<m_o; ++i)
        {
         if(i<m_p)
           {
            out[from,0] = 1.0;
            out[from++,1] = 2.;
           }
         else
           {
            out[from,0] = 0.0;
            out[from++,1] = 2.;
           }
        }

      for(ulong i = 0; i<m_q; ++i)
        {
         out[from,0] = 0.0;
         out[from++,1] = 1.;
        }

      return out.Transpose();
     }
   virtual string       parameterNames(void) override
     {
      return common_names(m_p,m_o,m_q);
     }
   virtual matrix       simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override
     {
      vector errors = rng.rng(_nobs+burn);

      if(initial_value == NULL)
        {
         vector scale = vector::Ones(parameters.Size());
         for(ulong i = m_p+1; i<(m_p+m_o+1); scale[i] = 0.5, ++i);
         vector temp = np::sliceVector(parameters,1) * np::sliceVector(scale,1);
         double persistence = temp.Sum();
         if((1.-persistence) >0)
            initial_value = parameters[0]/(1.-persistence);
         else
           {
            Print(__FUNCTION__, " InitialValueWarning ");
            initial_value = parameters[0];
           }
        }
      vector sigma2,data,fsigma,fdata;
      sigma2 = data = fsigma = fdata = vector::Zeros(_nobs+burn);
      ulong max_lag = MathMax(MathMax(m_p,m_o),m_q);
      np::fillVector(fsigma,initial_value,0,long(max_lag));
      double dv = pow(initial_value,2./m_power);
      np::fillVector(sigma2,dv,0,long(max_lag));
      for(ulong i = 0; i<max_lag; ++i)
        {
         data[i] = sqrt(sigma2[i])*errors[i];
         fdata[i] = pow(fabs(data[i]),m_power);
        }

      for(ulong t = max_lag; t<(_nobs+burn); ++t)
        {
         ulong loc = 0;
         fsigma[t] = parameters[loc];
         loc += 1;
         for(ulong j = 0; j<m_p; ++j)
           {
            fsigma[t] += parameters[loc]*fdata[t-1-j];
            loc+=1;
           }
         for(ulong j = 0; j<m_o; ++j)
           {
            fsigma[t] += parameters[loc]*fdata[t-1-j]*((data[t-1-j]<0.)?1.:0.0);
            loc+=1;
           }
         for(ulong j = 0; j<m_q; ++j)
           {
            fsigma[t] += parameters[loc]*fdata[t-1-j];
            loc+=1;
           }
         sigma2[t] = pow(fsigma[t],2./m_power);
         data[t] = errors[t]*sqrt(sigma2[t]);
         fdata[t] = pow(MathAbs(data[t]),m_power);
        }

      data = np::sliceVector(data,long(burn));
      sigma2 = np::sliceVector(sigma2,long(burn));

      matrix out(data.Size(),2);
      out.Col(data,0);
      out.Col(sigma2,1);

      return out;
     }
  };
//+------------------------------------------------------------------+
//| ARCH process                                                     |
//+------------------------------------------------------------------+
class CArchProcess: public CGarchProcess
  {
public:
                     CArchProcess(void)
     {
      m_initialized = CGarchProcess::_initialize(VOL_ARCH,true,false,0,"ARCH");
     }
                     CArchProcess(ulong p,int seed =0, ulong nbootstraps=100)
     {
      m_initialized = CGarchProcess::_initialize(VOL_ARCH,true,false,0,"ARCH",seed,p,0,0,2.0,0,-1,nbootstraps);
     }
                    ~CArchProcess(void)
     {
     }
   virtual vector    startingValues(vector& resids) override
     {
      vector alphas = np::arange(17,0.1,0.05);
      vector svs[17];
      vector bc = backCast(resids);
      vector llfs = alphas;
      matrix vb = varianceBounds(resids);
      vector sv,temp;
      for(uint i = 0; i<svs.Size(); ++i)
        {
         sv = (1.0 - alphas[i])*resids.Var() * vector::Ones(m_p+1);
         np::fillVector(sv,alphas[i]/double(m_p),1);
         svs[i] = sv;
         llfs[i] = _gaussianloglikelihood(sv,resids,bc,vb);
        }
      ulong loc = llfs.ArgMax();
      return svs[loc];
     }
  };

//+------------------------------------------------------------------+