ufranco/Article-22714-SLSQP-MQL5-Volatility-Library-Integration
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Article-22714-SLSQP-MQL5-Vo.../Arch/Univariate/mean.mqh

2807 lines
94 KiB
MQL5
Raw Permalink Normal View History

modified code to enable easy compilation in shared projects folder
2026-06-05 14:42:49 +02:00
//+------------------------------------------------------------------+
//| 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();
}
};
//+------------------------------------------------------------------
Reference in a new issue Copy permalink