ufranco/Article-22714-Volatility-Modeling-Part-3
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Article-22714-Volatility-Mo.../Include/slsqp_article/Arch/Univariate/volatility.mqh

1110 lines
41 KiB
MQL5
Raw Permalink Normal View History

first commit
2026-06-03 20:14:05 +02:00
//+------------------------------------------------------------------+
//| 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;
}
virtual vector garch_spec(void)
{
return vector::Zeros(4);
}
};
//+------------------------------------------------------------------+
//| 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;
}
virtual vector garch_spec(void) override
{
vector specs(4);
specs[0] = double(m_p);
specs[1] = double(m_o);
specs[2] = double(m_q);
specs[3] = double(m_power);
return specs;
}
};
//+------------------------------------------------------------------+
//| ARCH process |
//+------------------------------------------------------------------+
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(ulong(17),double(0.1),double(0.05));
vector svs[17];
vector bc = backCast(resids);
vector llfs = alphas;
matrix vb = varianceBounds(resids);
vector sv,temp;
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];
}
};
//+------------------------------------------------------------------+
//| GJR-GARCH process |
//+------------------------------------------------------------------+
class CGjrGarchProcess: public CGarchProcess
{
public:
CGjrGarchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_GJR_GARCH,true,false,0,"GJR-GARCH");
}
CGjrGarchProcess(ulong p, ulong o, ulong q, int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_GJR_GARCH,true,false,0,"GJR-GARCH",seed,p,o,q,2.0,0,-1,nboots);
}
};
//+------------------------------------------------------------------+
//| AVARCH process |
//+------------------------------------------------------------------+
class CAvarchProcess: public CGarchProcess
{
public:
CAvarchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_AVARCH,true,false,0,"AVARCH",0,1,0,0,1.0);
}
CAvarchProcess(ulong p,int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_AVARCH,true,false,0,"AVARCH",seed,p,0,0,1.0,0,-1,nboots);
}
};
//+------------------------------------------------------------------+
//| AVGARCH process |
//+------------------------------------------------------------------+
class CAvgarchProcess: public CGarchProcess
{
public:
CAvgarchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_AVGARCH,true,false,0,"AVGARCH",0,1,0,1,1.0);
}
CAvgarchProcess(ulong p,ulong q,int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_AVGARCH,true,false,0,"AVGARCH",seed,p,0,q,1.0,0,-1,nboots);
}
};
//+------------------------------------------------------------------+
//| TARCH process |
//+------------------------------------------------------------------+
class CTarchProcess: public CGarchProcess
{
public:
CTarchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_TARCH,true,false,0,"TARCH",0,1,1,1,1.0);
}
CTarchProcess(ulong p,ulong o,ulong q,int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_TARCH,true,false,0,"TARCH",seed,p,o,q,1.0,0,-1,nboots);
}
};
/*//+------------------------------------------------------------------+
//| power Arch process |
//+------------------------------------------------------------------+
class CPowerArchProcess: public CGarchProcess
{
public:
CPowerArchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_P_ARCH,true,false,0,"Power ARCH",0,1,0,0,1.5);
}
CPowerArchProcess(ulong p,double power,int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_P_ARCH,true,false,0,"Power ARCH",seed,p,0,0,power,0,-1,nboots);
}
};
//+------------------------------------------------------------------+
//| power Garch process |
//+------------------------------------------------------------------+
class CPowerGarchProcess: public CGarchProcess
{
public:
CPowerGarchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_P_GARCH,true,false,0,"Power GARCH",0,1,0,1,1.5);
}
CPowerGarchProcess(ulong p,ulong q,double power,int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_P_GARCH,true,false,0,"Power GARCH",seed,p,0,q,power,0,-1,nboots);
}
};
//+------------------------------------------------------------------+
//| asym power Garch process |
//+------------------------------------------------------------------+
class CAsymPowerGarchProcess: public CGarchProcess
{
public:
CAsymPowerGarchProcess(void)
{
m_initialized = CGarchProcess::_initialize(VOL_AP_GARCH,true,false,0,"Asym Power GARCH",0,1,0,1,1.5);
}
CAsymPowerGarchProcess(ulong p,ulong o,ulong q,double power,int seed = 0,ulong nboots=100)
{
m_initialized = CGarchProcess::_initialize(VOL_AP_GARCH,true,false,0,"Asym Power GARCH",seed,p,o,q,power,0,-1,nboots);
}
};*/
//+------------------------------------------------------------------+
Reference in a new issue Copy permalink