3524 lines
139 KiB
MQL5
3524 lines
139 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| mean.mqh |
|
|
//| Copyright 2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2025, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include"volatility.mqh"
|
|
#ifdef __SLSQP__
|
|
#include"..\..\slsqp.mqh"
|
|
#else
|
|
#include<Math\Alglib\delegatefunctions.mqh>
|
|
#include<Math\Alglib\optimization.mqh>
|
|
#endif
|
|
//+------------------------------------------------------------------+
|
|
//| Abstract base class for mean models |
|
|
//+------------------------------------------------------------------+
|
|
class CArchModel: public CObject
|
|
{
|
|
protected:
|
|
//--- Class Members
|
|
ArchParameters m_model_spec; // Structure containing execution specifications
|
|
vector m_params; // Optimized joint parameter values array
|
|
vector m_fit_indices; // Index flags mapping timestamps during optimization
|
|
vector m_backcast; // Initial pre-sample variance conditions vector
|
|
vector m_fit_y; // In-sample target timeline observation variables
|
|
matrix m_fit_x; // Independent exogenous variables design matrix
|
|
matrix m_fit_regressors; // Cached design matrix mapping lags/regressors
|
|
double m_scale; // Internal adjustment scaling multiplier
|
|
ulong m_nobs; // Number of total parsed dataset elements
|
|
ulong m_holdback; // Truncated initial sequence element offset
|
|
ulong m_maxlags; // Maximum depth threshold across mean/volatility processes
|
|
ulong m_num_params; // Total dedicated mean-model estimation parameter counts
|
|
bool m_converged; // Boolean indicating successful optimizer exit metrics
|
|
bool m_rescale; // Operational flag forcing internal auto-scaling routines
|
|
bool m_constant; // Structural toggle introducing intersect intercepts
|
|
bool m_rotated; // Indicator mapping variance targeting or HAR adjustments
|
|
string m_name; // Custom textual identifier for model variants
|
|
matrix m_regressors; // Dynamic matrix tracking user external series structures
|
|
matrix m_lags; // Historical lag tracking layout metrics
|
|
matrix m_var_bounds; // Linear constraint boundaries array for variance limits
|
|
bool m_delete_vp; // Lifecycle memory management tracking for Volatility object
|
|
bool m_delete_dist; // Lifecycle memory management tracking for Distribution object
|
|
|
|
CDistribution *m_distribution; // Object pointer managing target error innovations
|
|
CVolatilityProcess *m_vp; // Object pointer managing conditional variance steps
|
|
|
|
bool m_initialized; // System lock indicating model verification validity
|
|
|
|
//--- Calculate epsilon vector for numerical optimization step intervals
|
|
vector _get_epsilon(vector &x, double s, ulong n, double epsilon = EMPTY_VALUE)
|
|
{
|
|
vector out;
|
|
if(epsilon == EMPTY_VALUE)
|
|
// Machine epsilon scaling method based on parameter magnitude bounds
|
|
out = pow(DBL_EPSILON,1./s)*np::maximum(MathAbs(x),0.1);
|
|
else
|
|
{
|
|
out = vector::Zeros(n);
|
|
out.Fill(epsilon);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
//--- Approximates the partial derivatives (Jacobian/Gradient vector) of the objective function
|
|
matrix approx_fprime(vector &x,vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false,bool centered = false, double epsilon = EMPTY_VALUE)
|
|
{
|
|
ulong n = x.Size();
|
|
vector f0 = _loglikelihood(x,sigma2,backcast,varbounds,individual);
|
|
matrix grad = matrix::Zeros(n,f0.Size());
|
|
vector ei = vector::Zeros(n);
|
|
vector eps,row;
|
|
|
|
if(!centered) // Standard Forward Finite-Difference estimation
|
|
{
|
|
eps = _get_epsilon(x,2.,n,epsilon);
|
|
for(ulong i = 0; i<n; ++i)
|
|
{
|
|
ei[i] = eps[i];
|
|
row = (_loglikelihood(x+ei,sigma2,backcast,varbounds,individual) - f0)/eps[i];
|
|
ei[i] = 0.0;
|
|
grad.Row(row,i);
|
|
}
|
|
}
|
|
else // Higher accuracy O(h^2) Central Finite-Difference estimation
|
|
{
|
|
eps = _get_epsilon(x,3,n,epsilon)/2.;
|
|
for(ulong i = 0; i<n; ++i)
|
|
{
|
|
ei[i] = eps[i];
|
|
row = (_loglikelihood(x+ei,sigma2,backcast,varbounds,individual) - _loglikelihood(x-ei,sigma2,backcast,varbounds,individual))/(2. * eps[i]);
|
|
ei[i] = 0.0;
|
|
grad.Row(row,i);
|
|
}
|
|
}
|
|
return grad.Transpose();
|
|
}
|
|
|
|
//--- Approximates the structural Hessian matrix via localized finite-difference increments
|
|
matrix approx_hess(vector &x,vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false,double epsilon = EMPTY_VALUE)
|
|
{
|
|
ulong n = x.Size();
|
|
vector h = _get_epsilon(x,4,n,epsilon);
|
|
matrix ee;
|
|
ee.Diag(h); // Diagonal perturbation step sizing arrangement
|
|
matrix hess = h.Outer(h);
|
|
|
|
// Compute symmetric multi-directional cross-partial evaluations
|
|
for(ulong i = 0; i<n; ++i)
|
|
{
|
|
for(ulong j = i; j<n; ++j)
|
|
{
|
|
vector temp = (_loglikelihood(x+ee.Row(i)+ee.Row(j),sigma2,backcast,varbounds,individual) -
|
|
_loglikelihood(x+ee.Row(i)-ee.Row(j),sigma2,backcast,varbounds,individual) -
|
|
(_loglikelihood(x-ee.Row(i)+ee.Row(j),sigma2,backcast,varbounds,individual) -
|
|
_loglikelihood(x-ee.Row(i)-ee.Row(j),sigma2,backcast,varbounds,individual)))/(4.*hess[i,j]);
|
|
hess[i,j] = temp[0];
|
|
hess[j,i] = hess[i,j]; // Enforce mathematical matrix symmetry
|
|
}
|
|
}
|
|
return hess;
|
|
}
|
|
|
|
//--- Rescales the residuals vector to defend optimizer from underflow/overflow bounds
|
|
void _checkscale(vector& resids)
|
|
{
|
|
double ogscale = resids.Var();
|
|
double scale = ogscale;
|
|
double rescale = 1.;
|
|
|
|
// Loop iteratively to determine decimal step shifting scale bounds
|
|
while((scale>=1000.|| scale<0.1) && scale>0.0)
|
|
{
|
|
if(scale<1.0)
|
|
rescale*=10.0;
|
|
else
|
|
rescale/=10.0;
|
|
scale = ogscale*pow(rescale,2.0);
|
|
}
|
|
if(rescale == 1.0)
|
|
return;
|
|
if(!m_rescale)
|
|
{
|
|
Print(__FUNCTION__, " y is poorly scaled, which may affect convergence of the optimizer when estimating the model parameters.\n","Calculated scale is ",ogscale," . Parameter estimation works better when this value is between 1 and 1000. The recommended rescaling is ",rescale," * y ");
|
|
return;
|
|
}
|
|
m_scale = rescale;
|
|
}
|
|
|
|
//--- Virtual hook: triggered natively when localized array tracking scales change
|
|
virtual void _scalechanged(void)
|
|
{
|
|
return;
|
|
}
|
|
|
|
//--- Native placeholder computing the structural Coefficient of Determination ($R^2$)
|
|
virtual double _r2(vector& params)
|
|
{
|
|
return EMPTY_VALUE;
|
|
}
|
|
|
|
//--- Virtual hook: calculates structural analytical start guesses bypassing ARCH processes
|
|
virtual vector _fit_no_arch_normal_errors_params(void)
|
|
{
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
//--- Evaluates a baseline static Gaussian Log-Likelihood sequence
|
|
double _static_gaussian_loglikelihood(vector &resids_)
|
|
{
|
|
ulong nobs_ = resids_.Size();
|
|
double sigma2 = resids_.Dot(resids_)/double(nobs_);
|
|
double ll = -0.5 * double(nobs_) * log(2.0 * M_PI);
|
|
ll -= 0.5*double(nobs_)*log(sigma2);
|
|
ll -= 0.5*double(nobs_);
|
|
return ll;
|
|
}
|
|
|
|
//--- Unpacks the unified parameter array into mean, volatility, and error component slices
|
|
bool _parse_parameters(vector &x, vector& out[])
|
|
{
|
|
ulong km = m_num_params;
|
|
ulong kv = m_vp.numParams();
|
|
ArrayResize(out,3);
|
|
|
|
// Slice array segments explicitly based on model configuration widths
|
|
out[0] = (km)?np::sliceVector(x,0,long(km)):vector::Zeros(0);
|
|
out[1] = (kv)?np::sliceVector(x,long(km),long(km+kv)):vector::Zeros(0);
|
|
out[2] = ((km+kv)<x.Size())?np::sliceVector(x,long(km+kv)):vector::Zeros(0);
|
|
return true;
|
|
}
|
|
|
|
//--- Computes Negative Log-Likelihood aggregates across all linked subsystems
|
|
vector _loglikelihood(vector& parameters, vector& sigma2, vector& backcast, matrix& varbounds, bool individual=false)
|
|
{
|
|
vector prs[];
|
|
_parse_parameters(parameters,prs);
|
|
|
|
// Extract sequential system steps
|
|
vector resids = resids(prs[0],EMPTY_VECTOR,EMPTY_MATRIX);
|
|
sigma2 = m_vp.computeVariance(prs[1],resids,sigma2,backcast,varbounds);
|
|
vector llf = m_distribution.loglikelihood(prs[2],resids,sigma2,individual);
|
|
|
|
return (-1.*llf); // Inverse returns optimization-ready minimized objective shape
|
|
}
|
|
|
|
//--- Performs localized dataset boundary shifting adjustments
|
|
virtual void _adjust_sample(long& first, long& last)
|
|
{
|
|
return;
|
|
}
|
|
|
|
//--- Safely flushes dynamic sub-allocations on execution termination sequences
|
|
virtual void deinitialize(void)
|
|
{
|
|
if(m_delete_vp && CheckPointer(m_vp) == POINTER_DYNAMIC)
|
|
delete m_vp;
|
|
m_vp = NULL;
|
|
if(m_delete_dist && CheckPointer(m_distribution) == POINTER_DYNAMIC)
|
|
delete m_distribution;
|
|
m_distribution = NULL;
|
|
}
|
|
|
|
public:
|
|
//+------------------------------------------------------------------+
|
|
//| Constructor & Getters |
|
|
//+------------------------------------------------------------------+
|
|
CArchModel(void):m_vp(NULL),
|
|
m_initialized(false),
|
|
m_distribution(NULL),
|
|
m_nobs(0),
|
|
m_holdback(0),
|
|
m_converged(false),
|
|
m_rescale(false),
|
|
m_constant(true),
|
|
m_delete_vp(true),
|
|
m_delete_dist(true),
|
|
m_scale(1.0),
|
|
m_name(NULL),
|
|
m_fit_y(vector::Zeros(0)),
|
|
m_params(vector::Zeros(0)),
|
|
m_fit_indices(vector::Zeros(0)),
|
|
m_fit_regressors(matrix::Zeros(0,0)),
|
|
m_fit_x(matrix::Zeros(0,0)),
|
|
m_regressors(matrix::Zeros(0,0))
|
|
{
|
|
}
|
|
|
|
ulong nobs(void) { return m_nobs; }
|
|
ulong hold_back(void) { return m_holdback; }
|
|
bool converged(void) { return m_converged;}
|
|
bool rescaled(void) { return m_rescale; }
|
|
double get_scale(void) { return m_scale; }
|
|
vector get_y(void) { return m_model_spec.observations; }
|
|
vector get_params(void) { return m_params; }
|
|
ArchParameters get_specification(void) { return m_model_spec; }
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| File Input/Output Binary Serialization Handlers |
|
|
//+------------------------------------------------------------------+
|
|
virtual bool Save(const int file_handle) override
|
|
{
|
|
if(file_handle!=INVALID_HANDLE)
|
|
{
|
|
// Write validation marker block
|
|
if(FileWriteLong(file_handle,-1)==sizeof(long))
|
|
{
|
|
// Serialize Parameter structures sizes and contents
|
|
if(FileWriteInteger(file_handle,int(m_params.Size()),INT_VALUE)!=INT_VALUE)
|
|
return(false);
|
|
for(ulong i = 0; i<m_params.Size(); ++i)
|
|
{
|
|
if(FileWriteDouble(file_handle,m_params[i])!=sizeof(double))
|
|
{
|
|
Print(__FUNCTION__, " file write error ", GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// Serialize Mean Model configurations
|
|
if(FileWriteInteger(file_handle,int(m_model_spec.mean_model_type),INT_VALUE)!=INT_VALUE ||
|
|
FileWriteLong(file_handle,long(m_model_spec.holdout_size))!=sizeof(long) ||
|
|
FileWriteInteger(file_handle,int(m_model_spec.is_rescale_enabled),INT_VALUE)!=INT_VALUE ||
|
|
FileWriteInteger(file_handle,int(m_model_spec.include_constant),INT_VALUE)!=INT_VALUE ||
|
|
FileWriteInteger(file_handle,int(m_model_spec.use_har_rotation),INT_VALUE)!=INT_VALUE)
|
|
return(false);
|
|
|
|
if(FileWriteDouble(file_handle,m_model_spec.scaling_factor)!=sizeof(double))
|
|
return(false);
|
|
if(FileWriteLong(file_handle,long(m_model_spec.mean_lags.Size()))!=sizeof(long))
|
|
return(false);
|
|
|
|
for(ulong i = 0; i<m_model_spec.mean_lags.Size(); ++i)
|
|
{
|
|
if(FileWriteDouble(file_handle,m_model_spec.mean_lags[i])!=sizeof(double))
|
|
{
|
|
Print(__FUNCTION__, " file write error ", GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// Serialize Volatility specification tracking indices
|
|
if(FileWriteInteger(file_handle,int(m_model_spec.vol_model_type),INT_VALUE)!=INT_VALUE)
|
|
return(false);
|
|
if(FileWriteInteger(file_handle,int(m_model_spec.vol_rng_seed),INT_VALUE)!=INT_VALUE ||
|
|
FileWriteLong(file_handle,m_model_spec.sample_start_idx)!=sizeof(long) ||
|
|
FileWriteLong(file_handle,m_model_spec.sample_end_idx)!=sizeof(long) ||
|
|
FileWriteLong(file_handle,long(m_model_spec.min_bootstrap_sims))!=INT_VALUE ||
|
|
FileWriteLong(file_handle,long(m_model_spec.garch_p))!=sizeof(long) ||
|
|
FileWriteLong(file_handle,long(m_model_spec.garch_o))!=sizeof(long) ||
|
|
FileWriteLong(file_handle,long(m_model_spec.garch_q))!=sizeof(long))
|
|
return(false);
|
|
|
|
if(FileWriteDouble(file_handle,m_model_spec.vol_power)!=sizeof(double))
|
|
return false;
|
|
if(FileWriteInteger(file_handle,int(m_model_spec.dist_type),INT_VALUE)!=INT_VALUE)
|
|
return(false);
|
|
if(FileWriteInteger(file_handle,int(m_model_spec.dist_rng_seed),INT_VALUE)!=INT_VALUE)
|
|
return(false);
|
|
if(FileWriteLong(file_handle,long(m_model_spec.observations.Size()))!=sizeof(long))
|
|
return(false);
|
|
|
|
for(ulong i = 0; i<m_model_spec.observations.Size(); ++i)
|
|
{
|
|
if(FileWriteDouble(file_handle,m_model_spec.observations[i])!=sizeof(double))
|
|
{
|
|
Print(__FUNCTION__, " file write error ", GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// Serialize multi-dimensional Exogenous Matrix arrays
|
|
if(FileWriteLong(file_handle,long(m_model_spec.exog_data.Rows()))!=sizeof(long))
|
|
return(false);
|
|
if(FileWriteLong(file_handle,long(m_model_spec.exog_data.Cols()))!=sizeof(long))
|
|
return(false);
|
|
for(ulong i = 0; i<m_model_spec.exog_data.Rows(); ++i)
|
|
{
|
|
for(ulong j = 0; j<m_model_spec.exog_data.Cols(); ++j)
|
|
{
|
|
if(FileWriteDouble(file_handle,m_model_spec.exog_data[i,j])!=sizeof(double))
|
|
{
|
|
Print(__FUNCTION__, " file write error ", GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
virtual bool Load(const int file_handle) override
|
|
{
|
|
if(file_handle!=INVALID_HANDLE)
|
|
{
|
|
long ff = FileReadLong(file_handle);
|
|
if(ff == -1)
|
|
{
|
|
ResetLastError();
|
|
ulong size = ulong(FileReadLong(file_handle));
|
|
if(size)
|
|
{
|
|
m_params.Resize(size);
|
|
m_converged = true;
|
|
for(ulong i = 0; i<m_params.Size(); ++i)
|
|
m_params[i]=FileReadDouble(file_handle);
|
|
if(GetLastError())
|
|
{
|
|
Print(__FUNCTION__, " possible read error ", GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
m_converged = false;
|
|
m_params = vector::Zeros(0);
|
|
}
|
|
|
|
// Unpack deserialized variables back into configuration structure
|
|
m_model_spec.mean_model_type = ENUM_MEAN_MODEL(FileReadInteger(file_handle,INT_VALUE));
|
|
m_model_spec.holdout_size = ulong(FileReadLong(file_handle));
|
|
m_model_spec.is_rescale_enabled = bool(FileReadInteger(file_handle,INT_VALUE));
|
|
m_model_spec.include_constant = bool(FileReadInteger(file_handle,INT_VALUE));
|
|
m_model_spec.use_har_rotation = bool(FileReadInteger(file_handle,INT_VALUE));
|
|
|
|
if(GetLastError())
|
|
{
|
|
Print(__FUNCTION__, " possible read error ", GetLastError());
|
|
return false;
|
|
}
|
|
|
|
m_model_spec.scaling_factor=FileReadDouble(file_handle);
|
|
size = ulong(FileReadLong(file_handle));
|
|
if(size)
|
|
{
|
|
m_model_spec.mean_lags.Resize(size);
|
|
for(ulong i = 0; i<m_model_spec.mean_lags.Size(); ++i)
|
|
m_model_spec.mean_lags[i]=FileReadDouble(file_handle);
|
|
}
|
|
else
|
|
{
|
|
m_model_spec.mean_lags = vector::Zeros(0);
|
|
}
|
|
|
|
if(GetLastError())
|
|
{
|
|
Print(__FUNCTION__, " possible read error ", GetLastError());
|
|
return false;
|
|
}
|
|
|
|
m_model_spec.vol_model_type = ENUM_VOLATILITY_MODEL(FileReadInteger(file_handle,INT_VALUE));
|
|
m_model_spec.vol_rng_seed = (FileReadInteger(file_handle,INT_VALUE));
|
|
m_model_spec.sample_start_idx = long(FileReadLong(file_handle));
|
|
m_model_spec.sample_end_idx = long(FileReadLong(file_handle));
|
|
m_model_spec.min_bootstrap_sims = ulong(FileReadLong(file_handle));
|
|
m_model_spec.garch_p = ulong(FileReadLong(file_handle));
|
|
m_model_spec.garch_o = ulong(FileReadLong(file_handle));
|
|
m_model_spec.garch_q = ulong(FileReadLong(file_handle));
|
|
|
|
m_model_spec.vol_power = FileReadDouble(file_handle);
|
|
m_model_spec.dist_type = ENUM_DISTRIBUTION_MODEL(FileReadInteger(file_handle,INT_VALUE));
|
|
m_model_spec.dist_rng_seed = FileReadInteger(file_handle,INT_VALUE);
|
|
|
|
if(GetLastError())
|
|
{
|
|
Print(__FUNCTION__, " possible read error ", GetLastError());
|
|
return false;
|
|
}
|
|
|
|
size = ulong(FileReadLong(file_handle));
|
|
if(size)
|
|
{
|
|
m_model_spec.observations.Resize(size);
|
|
for(ulong i = 0; i<m_model_spec.observations.Size(); ++i)
|
|
m_model_spec.observations[i] = FileReadDouble(file_handle);
|
|
}
|
|
else
|
|
m_model_spec.observations = vector::Zeros(0);
|
|
|
|
long rws = FileReadLong(file_handle);
|
|
long cls = FileReadLong(file_handle);
|
|
m_model_spec.exog_data = matrix::Zeros(rws,cls);
|
|
if(m_model_spec.exog_data.Rows())
|
|
{
|
|
for(ulong i = 0; i<m_model_spec.exog_data.Rows(); ++i)
|
|
{
|
|
for(ulong j = 0; j<m_model_spec.exog_data.Cols(); ++j)
|
|
m_model_spec.exog_data[i,j] = FileReadDouble(file_handle);
|
|
}
|
|
}
|
|
|
|
if(GetLastError())
|
|
{
|
|
Print(__FUNCTION__, " possible read error ", GetLastError());
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
virtual bool is_initialized(void) { return m_initialized; }
|
|
~CArchModel(void) { deinitialize(); }
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Core Volatility and Distribution Access Pointers |
|
|
//+------------------------------------------------------------------+
|
|
CVolatilityProcess *volatility(void) { return m_vp; }
|
|
CDistribution *distribution(void) { return m_distribution; }
|
|
ENUM_MEAN_MODEL mean_model_type(void) { return m_model_spec.mean_model_type; }
|
|
ENUM_VOLATILITY_MODEL volatility_process_type(void) { return m_vp.volatilityprocess(); }
|
|
ENUM_DISTRIBUTION_MODEL distribution_type(void) { return m_distribution.distribution_type(); }
|
|
|
|
//--- Wrapper calling private log-likelihood objective processes
|
|
vector objective(vector& parameters, vector& sigma2, vector &backcast, matrix& varbounds, bool individual=false)
|
|
{
|
|
return _loglikelihood(parameters,sigma2,backcast,varbounds,individual);
|
|
}
|
|
|
|
//--- Wrapper tracking optimization Jacobian matrices
|
|
matrix jacobian(vector& parameters, vector& sigma2, vector &backcast, matrix& varbounds, bool individual=false,bool centered=false, double eps=EMPTY_VALUE)
|
|
{
|
|
return approx_fprime(parameters,sigma2,backcast,varbounds,individual,centered,eps);
|
|
}
|
|
|
|
virtual vector starting_values(void)
|
|
{
|
|
return _fit_no_arch_normal_errors_params();
|
|
}
|
|
|
|
virtual vector resids(vector& params, vector& y, matrix& regressors)
|
|
{
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
//--- Calculates the Parameter Covariance Matrix using robust Sandwich Estimation ($H^{-1} J H^{-1}$)
|
|
matrix compute_param_cov(vector& params, vector& backcast, bool robust = true)
|
|
{
|
|
vector sv_ = starting_values();
|
|
vector fy = vector::Zeros(0);
|
|
matrix fr = matrix::Zeros(0,0);
|
|
vector resids_ = resids(sv_,fy,fr);
|
|
matrix vb = m_vp.varianceBounds(resids_);
|
|
ulong nobs = resids_.Size();
|
|
|
|
if(!backcast.Size() && !m_backcast.Size())
|
|
{
|
|
backcast = m_vp.backCast(resids_);
|
|
m_backcast = backcast;
|
|
}
|
|
else
|
|
if(!backcast.Size())
|
|
backcast = m_backcast;
|
|
|
|
vector sgma2=vector::Zeros(nobs);
|
|
matrix hess = approx_hess(params,sgma2,backcast,vb,false);
|
|
hess/=double(nobs);
|
|
matrix inv_hess = hess.Inv(); // Invert the Hessian ($H^{-1}$)
|
|
|
|
if(robust) // Quasi-Maximum Likelihood (QMLE) Robust Sandwich Variance Covariance calculation
|
|
{
|
|
matrix scores = approx_fprime(params,sgma2,backcast,vb,true);
|
|
matrix score_cov = scores.Transpose().Cov(); // Inner Outer-Product Jacobian ($J$)
|
|
return (inv_hess.MatMul(score_cov).MatMul(inv_hess))/double(nobs);
|
|
}
|
|
else
|
|
return inv_hess/double(nobs);
|
|
}
|
|
|
|
virtual Constraints constraints(void)
|
|
{
|
|
Constraints out;
|
|
return out;
|
|
}
|
|
|
|
virtual matrix bounds(void)
|
|
{
|
|
vector temp = {AL_NEGINF,AL_POSINF};
|
|
return np::repeat_vector_as_rows_cols(temp,m_num_params,false);
|
|
}
|
|
|
|
virtual string name(void)
|
|
{
|
|
return m_name;
|
|
}
|
|
|
|
virtual ArchForecast forecast(vector& params,matrix& x[],ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
|
|
{
|
|
ArchForecast out;
|
|
return out;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Arch model fixed result struct |
|
|
//| Holds the optimization outputs and diagnostic metrics for a |
|
|
//| finalized ARCH/GARCH execution state. |
|
|
//+------------------------------------------------------------------+
|
|
struct ArchModelFixedResult
|
|
{
|
|
public:
|
|
double loglikelihood; // Maximized Joint Log-Likelihood objective value
|
|
vector params; // Final estimated parameter values vector
|
|
vector conditional_volatility; // Modeled time-varying conditional standard deviation (\sigma_t)
|
|
ulong nobs; // Total sample size count used in fitting
|
|
vector resid; // Extracted mean model residual errors (\epsilon_t)
|
|
|
|
//--- Default Constructor (Initializes to empty state)
|
|
ArchModelFixedResult(void)
|
|
{
|
|
loglikelihood = EMPTY_VALUE;
|
|
nobs = 0;
|
|
params = resid = conditional_volatility = vector::Zeros(0);
|
|
}
|
|
|
|
//--- Copy Constructor
|
|
ArchModelFixedResult(ArchModelFixedResult &other)
|
|
{
|
|
loglikelihood = other.loglikelihood;
|
|
nobs = other.nobs;
|
|
params = other.params;
|
|
resid = other.resid;
|
|
conditional_volatility = other.conditional_volatility;
|
|
}
|
|
|
|
//--- Parameterized Constructor
|
|
ArchModelFixedResult(vector &_params, vector &_resid, vector &_volatility, vector &_depvar, double _loglikelihood)
|
|
{
|
|
loglikelihood = _loglikelihood;
|
|
params = _params;
|
|
resid = _resid;
|
|
conditional_volatility = _volatility;
|
|
nobs = _resid.Size(); // Dynamically infers observation size from residual length
|
|
}
|
|
|
|
//--- Destructor
|
|
~ArchModelFixedResult(void)
|
|
{
|
|
}
|
|
|
|
//--- Overloaded Assignment Operator
|
|
void operator=(ArchModelFixedResult &other)
|
|
{
|
|
loglikelihood = other.loglikelihood;
|
|
nobs = other.nobs;
|
|
params = other.params;
|
|
resid = other.resid;
|
|
conditional_volatility = other.conditional_volatility;
|
|
}
|
|
|
|
//--- Akaike Information Criterion (AIC): Penalizes model complexity linearly
|
|
double aic(void)
|
|
{
|
|
return -2 * loglikelihood + 2 * double(num_params());
|
|
}
|
|
|
|
//--- Returns the total parameter count estimated across the complete system
|
|
ulong num_params(void)
|
|
{
|
|
return params.Size();
|
|
}
|
|
|
|
//--- Schwarz/Bayesian Information Criterion (BIC): Penalizes complexity logarithmic to sample size
|
|
double bic(void)
|
|
{
|
|
return -2.*loglikelihood+log(nobs)*double(num_params());
|
|
}
|
|
|
|
//--- Returns standardized residuals (z_t = \epsilon_t / \sigma_t) for distribution testing
|
|
vector std_resid(void)
|
|
{
|
|
return resid/conditional_volatility;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Arch model result struct |
|
|
//| Inherits base properties from ArchModelFixedResult to add |
|
|
//| statistical inference methods and parameter covariance parsing. |
|
|
//+------------------------------------------------------------------+
|
|
struct ArchModelResult: public ArchModelFixedResult
|
|
{
|
|
public:
|
|
long fit_indices[2]; // [0] = Start index, [1] = Stop index mapping the target data slice
|
|
matrix param_cov; // Calculated parameter Variance-Covariance matrix
|
|
double r2; // Unadjusted Coefficient of Determination (R-squared)
|
|
ENUM_COVAR_TYPE cov_type; // Strategy flag applied during variance calculations (Standard vs Robust)
|
|
|
|
//--- Default Constructor (Sets defaults for invalid metrics)
|
|
ArchModelResult(void)
|
|
{
|
|
fit_indices[0] = fit_indices[1] = 0;
|
|
param_cov = matrix::Zeros(0,0);
|
|
r2 = EMPTY_VALUE;
|
|
cov_type = WRONG_VALUE;
|
|
}
|
|
|
|
//--- Parameterized Constructor
|
|
ArchModelResult(vector& _params, matrix &_param_cov, double _r2, vector& _resid, vector & _volatility, ENUM_COVAR_TYPE _cov_type, double _loglikelihood, long fit_start, long fit_stop)
|
|
{
|
|
loglikelihood = _loglikelihood;
|
|
params = _params;
|
|
param_cov = _param_cov;
|
|
resid = _resid;
|
|
conditional_volatility = _volatility;
|
|
nobs = _resid.Size(); // Tracks observation sample depth internally
|
|
|
|
fit_indices[0] = fit_start;
|
|
fit_indices[1] = fit_stop;
|
|
r2 = _r2;
|
|
cov_type = _cov_type;
|
|
}
|
|
|
|
//--- Destructor
|
|
~ArchModelResult(void)
|
|
{
|
|
}
|
|
|
|
//--- Copy Constructor
|
|
ArchModelResult(ArchModelResult &other)
|
|
{
|
|
loglikelihood = other.loglikelihood;
|
|
nobs = other.nobs;
|
|
params = other.params;
|
|
resid = other.resid;
|
|
conditional_volatility = other.conditional_volatility;
|
|
|
|
fit_indices[0] = other.fit_indices[0];
|
|
fit_indices[1] = other.fit_indices[1];
|
|
r2 = other.r2;
|
|
cov_type = other.cov_type;
|
|
param_cov = other.param_cov;
|
|
}
|
|
|
|
//--- Overloaded Assignment Operator
|
|
void operator=(ArchModelResult &other)
|
|
{
|
|
loglikelihood = other.loglikelihood;
|
|
nobs = other.nobs;
|
|
params = other.params;
|
|
resid = other.resid;
|
|
conditional_volatility = other.conditional_volatility;
|
|
|
|
fit_indices[0] = other.fit_indices[0];
|
|
fit_indices[1] = other.fit_indices[1];
|
|
r2 = other.r2;
|
|
cov_type = other.cov_type;
|
|
param_cov = other.param_cov;
|
|
}
|
|
|
|
//--- Dynamically evaluates or updates parameter covariance states using standard or robust methods
|
|
matrix parameter_covariance(CArchModel &archmodel)
|
|
{
|
|
vector fitted_params = archmodel.get_params();
|
|
if(cov_type == COVAR_ROBUST)
|
|
param_cov = archmodel.compute_param_cov(fitted_params,EMPTY_VECTOR); // Triggers QMLE Sandwich Estimation
|
|
else
|
|
param_cov = archmodel.compute_param_cov(fitted_params,EMPTY_VECTOR,false); // Standard Hessian Inverse
|
|
return param_cov;
|
|
}
|
|
|
|
//--- Computes asymptotic symmetric confidence intervals based on Inverse Normal critical thresholds
|
|
matrix conf_int(double alpha = 0.05)
|
|
{
|
|
// Calculate two-tailed critical value Z_{1 - \alpha/2}
|
|
double cv = CNormalDistr::InvNormalCDF(1.0-alpha/2.0);
|
|
vector se = std_err();
|
|
matrix out = matrix::Zeros(params.Size(),2);
|
|
|
|
out.Col(params - cv * se, 0); // Lower Bound: Estimate - (Z * SE)
|
|
out.Col(params + cv * se, 1); // Upper Bound: Estimate + (Z * SE)
|
|
|
|
return out;
|
|
}
|
|
|
|
//--- Adjusted R-squared: factors parameter-to-sample loss penalty adjustments
|
|
double rsquared_adj(void)
|
|
{
|
|
return 1.0 - ((1.0-r2)*double(nobs-1)/double(nobs-num_params()));
|
|
}
|
|
|
|
//--- Isolates individual parameter standard errors by extracting the square root of the covariance diagonal
|
|
vector std_err(void)
|
|
{
|
|
return sqrt(param_cov.Diag());
|
|
}
|
|
|
|
//--- Computes two-tailed survival probability p-values under the Standard Normal curve
|
|
vector pvalues(void)
|
|
{
|
|
vector pvals = tvalues();
|
|
for(ulong i = 0; i<pvals.Size(); ++i)
|
|
// Mirror distributions for lower tail and double for a two-tailed probability assessment
|
|
pvals[i] = CNormalDistr::NormalCDF(-1.0*MathAbs(pvals[i]))*2.0;
|
|
|
|
return pvals;
|
|
}
|
|
|
|
//--- Calculates asymptotic t-statistics (Wald ratios) testing the null hypothesis H_0: \beta_i = 0
|
|
vector tvalues(void)
|
|
{
|
|
return params/std_err();
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| The Lagrange multiplier (LM) test for ARCH effects |
|
|
//| Computes Engle's LM test statistic ($T \times R^2$) to test for |
|
|
//| conditional heteroskedasticity in series observations/residuals. |
|
|
//+------------------------------------------------------------------+
|
|
WaldTestStatistic archlmtest(vector& residuals, vector& conditional_volatility, ulong lags, bool standardized = false)
|
|
{
|
|
WaldTestStatistic out;
|
|
vector resids = residuals;
|
|
|
|
//--- Standardize residuals if requested ($z_t = \epsilon_t / \sigma_t$)
|
|
if(standardized)
|
|
resids = resids/conditional_volatility;
|
|
|
|
//--- Data Integrity: Scrub and filter out missing values (NaNs)
|
|
if(resids.HasNan())
|
|
{
|
|
vector nresids = vector::Zeros(resids.Size() - resids.HasNan());
|
|
for(ulong i = 0, k = 0; i<resids.Size(); ++i)
|
|
if(MathClassify(resids[i]) == FP_NAN)
|
|
continue;
|
|
else
|
|
nresids[k++] = resids[i];
|
|
resids = nresids;
|
|
}
|
|
|
|
int nobs = (int)resids.Size();
|
|
vector resid2 = MathPow(resids,2.0); // Convert series into squared innovations ($\epsilon_t^2$)
|
|
|
|
//--- Dynamic Lag Selection: Default to Schwert rule-of-thumb if lags = 0
|
|
if(!lags)
|
|
lags = ulong(ceil(12.0*pow(nobs/100.0,1.0/4.0)));
|
|
|
|
//--- Ensure lags fit within mathematical bounds relative to sample size
|
|
lags = MathMax(MathMin(resids.Size()/2 - 1, lags), 1);
|
|
|
|
//--- Validation: Ensure adequate sample size exists for secondary regression
|
|
if(resid2.Size()<3)
|
|
{
|
|
Print(__FUNCTION__, " Test requires at least 3 non-nan observations ");
|
|
return out;
|
|
}
|
|
|
|
//--- Reshape squared series into a 1-column matrix for matrix lag operations
|
|
matrix matres = matrix::Zeros(resid2.Size(),1);
|
|
matres.Col(resid2,0);
|
|
matrix lag[];
|
|
|
|
//--- Construct the design lag matrices via lagmat (TRIM_BOTH drops pre-sample rows)
|
|
if(!lagmat(matres,lag,lags,TRIM_BOTH,ORIGINAL_SEP))
|
|
{
|
|
Print(__FUNCTION__, " lagmat failed ");
|
|
return out;
|
|
}
|
|
|
|
//--- Append an intercept/constant column to the independent regressor matrix (lag[0])
|
|
matrix lagout;
|
|
if(!addtrend(lag[0],lagout,TREND_CONST_ONLY,true,HAS_CONST_SKIP))
|
|
{
|
|
Print(__FUNCTION__, " addtrend failed ");
|
|
return out;
|
|
}
|
|
lag[0] = lagout;
|
|
|
|
//--- Fit the auxiliary OLS regression: $\epsilon_t^2 = \alpha_0 + \alpha_1 \epsilon_{t-1}^2 + ... + \alpha_p \epsilon_{t-p}^2 + v_t$
|
|
OLS ols;
|
|
if(!ols.Fit(lag[1].Col(0),lag[0]))
|
|
{
|
|
Print(__FUNCTION__, " OLS fitting failed ");
|
|
return out;
|
|
}
|
|
|
|
//--- Calculate Engle's LM test statistic: Sample Size ($T$) $\times$ R-squared ($R^2$)
|
|
out.stat = nobs*ols.Rsqe();
|
|
|
|
//--- Document testing hypothesis strings
|
|
if(standardized)
|
|
{
|
|
out.null = "Standardized residuals are homoskedastic";
|
|
out.alternative = "Standardized residuals are conditionally heteroskedastic";
|
|
}
|
|
else
|
|
{
|
|
out.null = "Residuals are homoskedastic";
|
|
out.alternative = "Residuals are conditionally heteroskedastic";
|
|
}
|
|
|
|
//--- Save degrees of freedom (corresponds asymptotically to a $\chi^2$ distribution with $p$ lags)
|
|
out.df = lags;
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Forecast helper function |
|
|
//| Pads the front of a forecast matrix array with NaN entries to |
|
|
//| align historical timeline intervals or pre-sample periods. |
|
|
//+------------------------------------------------------------------+
|
|
bool _forecast_pad(int count, matrix &forecasts[], matrix &out[])
|
|
{
|
|
//--- Expand target output array to hold both padding frames and raw forecasts
|
|
ArrayResize(out, count + (int)forecasts.Size());
|
|
|
|
//--- Construct a structural matrix populated cleanly with NaN values
|
|
matrix fill = matrix::Zeros(forecasts[0].Rows(), forecasts[0].Cols());
|
|
fill.Fill(double("nan"));
|
|
|
|
//--- Phase 1: Inject empty padding matrices up to the specified boundary offset
|
|
for(int i = 0; i < count; ++i)
|
|
out[i] = fill;
|
|
|
|
//--- Phase 2: Copy raw forecast entries into the remaining indexed positions
|
|
for(int i = count; i < int(out.Size()); ++i)
|
|
out[i] = forecasts[i - count];
|
|
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Forecasting helper function (Matrix Overload) |
|
|
//| Prepend empty rows (NaN) to a matrix to align structural lengths |
|
|
//| across historical timeline matrices. |
|
|
//+------------------------------------------------------------------+
|
|
matrix _forecast_pad(int count, matrix &forecasts)
|
|
{
|
|
//--- Initialize a destination matrix expanded by the offset row count
|
|
matrix fill = matrix::Zeros(count+forecasts.Rows(), forecasts.Cols());
|
|
fill.Fill(double("nan")); // Default fill everything with NaN
|
|
|
|
//--- Shift and copy original matrix rows sequentially into the padded matrix
|
|
for(ulong i = (ulong)count; i<fill.Rows(); ++i)
|
|
fill.Row(forecasts.Row(ulong(i-count)),i);
|
|
|
|
return fill;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Forecasting helper function (Vector Overload) |
|
|
//| Prepend empty values (NaN) to a vector to match data streams |
|
|
//| and account for pre-sample estimation gaps. |
|
|
//+------------------------------------------------------------------+
|
|
vector _forecast_pad(int count, vector &forecasts)
|
|
{
|
|
//--- Initialize a destination vector expanded by the offset size count
|
|
vector fill = vector::Zeros(count+forecasts.Size());
|
|
fill.Fill(double("nan")); // Default fill everything with NaN
|
|
|
|
//--- Shift and copy original vector elements into the padded vector positions
|
|
for(ulong i = (ulong)count; i<fill.Size(); ++i)
|
|
fill[i] = forecasts[ulong(i-count)];
|
|
|
|
return fill;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Generate mean forecasts from an AR-X model |
|
|
//| Operates recursively over an out-of-sample forecast horizon by |
|
|
//| rolling autoregressive lags forward and adding exogenous impacts. |
|
|
//+------------------------------------------------------------------+
|
|
matrix _ar_forecast(vector& y, ulong horizon, ulong start_index, double constant, vector & arp, matrix &x[], vector &exogp)
|
|
{
|
|
ulong t = y.Size();
|
|
ulong p = arp.Size();
|
|
ulong first,last;
|
|
vector col;
|
|
|
|
//--- Initialize destination forecast canvas matrix
|
|
//--- Columns span the historic autoregressive lags [0...p-1] plus the forward horizon [p...p+horizon-1]
|
|
matrix fcasts = matrix::Zeros(t - start_index,p+horizon);
|
|
|
|
//--- Seed the first 'p' columns of the matrix with actual historical lag observations
|
|
for(ulong i = 0; i<p; ++i)
|
|
{
|
|
first = start_index - p + i + 1;
|
|
last = t - p + i + 1;
|
|
col = np::sliceVector(y,long(first),long(last));
|
|
fcasts.Col(col,i);
|
|
}
|
|
|
|
//--- Reverse AR parameters to naturally match chronological history via dot-product matrix multiplications
|
|
vector arp_rev = arp;
|
|
np::reverseVector(arp_rev);
|
|
matrix mlt;
|
|
vector mat;
|
|
|
|
//--- Main loop tracking forward time-steps sequentially over the forecast horizon
|
|
for(ulong i = p; i<horizon+p; ++i)
|
|
{
|
|
//--- Step 1: Extract the rolling sliding window containing the preceding 'p' lags (mix of historical and generated forecasts)
|
|
mlt = ((i-p)!=i)?np::sliceMatrixCols(fcasts,long(i-p),long(i)):matrix::Zeros(0,0);
|
|
|
|
//--- Step 2: Calculate the standard Autoregressive linear combination ($AR(p) = \phi_1 y_{t-1} + ... + \phi_p y_{t-p}$)
|
|
mat = (arp_rev.Size()&&mlt.Cols())?mlt.MatMul(arp_rev):vector::Zeros(fcasts.Rows());
|
|
|
|
//--- Step 3: Accumulate structural intercept constants into the conditional mean column path
|
|
fcasts.Col(constant+mat,i);
|
|
|
|
//--- Step 4: Map exogenous features ($\beta \cdot X_t$) if external matrix trackers are explicitly provided
|
|
if(x.Size())
|
|
{
|
|
mlt = matrix::Zeros(x[0].Rows(),x.Size());
|
|
for(uint k = 0; k<x.Size(); ++k)
|
|
mlt.Col(x[k].Col(i-p),k); // Col(i-p) isolates the precise forward index for feature step $h$
|
|
|
|
//--- Append the combined dot-product exogenous additions back into current forecast matrix coordinate arrays
|
|
col = fcasts.Col(i);
|
|
col += mlt.MatMul(exogp);
|
|
fcasts.Col(col,i);
|
|
}
|
|
}
|
|
|
|
//--- Clean-up and trim away historical pre-sample columns, leaving only pure calculated out-of-sample forecasts
|
|
fcasts = np::sliceMatrixCols(fcasts,long(p));
|
|
return fcasts;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Convert AR coefficients to an Impulse Response Function (IRF) |
|
|
//| Recursively traces the long-term impact of a single structural |
|
|
//| shock unit ($\epsilon_0 = 1$) forward through time step horizons. |
|
|
//+------------------------------------------------------------------+
|
|
vector _ar_to_impulse(ulong steps, vector& params)
|
|
{
|
|
ulong p = params.Size();
|
|
vector impulse = vector::Zeros(steps);
|
|
|
|
//--- Seed the shock at time zero ($t_0 = 1$)
|
|
impulse[0] = 1.0;
|
|
|
|
//--- If there are no AR parameters, the impact is purely transitory (white noise)
|
|
if(p == 0)
|
|
return impulse;
|
|
|
|
long k, st;
|
|
vector rparams, imp;
|
|
|
|
//--- Recursively project the shock propagation across subsequent time intervals
|
|
for(long i = 1; i<long(steps); ++i)
|
|
{
|
|
//--- Dynamically bound indices to avoid reading beyond available AR parameters or past history
|
|
k = (long)fmin(long(p) - 1, i - 1);
|
|
st = (long)fmax(i - long(p), 0);
|
|
|
|
//--- Isolate historical impulse factors ($y_{t-1}, y_{t-2}, ...$)
|
|
imp = np::sliceVector(impulse, st, i);
|
|
|
|
//--- Reverse the applicable slice of AR parameters to align correctly chronologically
|
|
rparams = np::sliceVector(params, k, END_REVERSE, STEP_REVERSE);
|
|
|
|
//--- Compute the dot product to dynamically determine the current step's decay or amplification
|
|
if(imp.Size()==rparams.Size() && imp.Size())
|
|
impulse[i] = imp.Dot(rparams);
|
|
}
|
|
|
|
return impulse;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Heterogeneous Autoregression (HAR) |
|
|
//+------------------------------------------------------------------+
|
|
class HARX: public CArchModel
|
|
{
|
|
protected:
|
|
#ifdef __SLSQP__
|
|
//+------------------------------------------------------------------+
|
|
//| Objective function for SLSQP optimizer |
|
|
//| Extends CFunctor to provide a standard interface for evaluating |
|
|
//| the log-likelihood objective value and its respective gradients. |
|
|
//+------------------------------------------------------------------+
|
|
class CObjectiveF: public CFunctor
|
|
{
|
|
private:
|
|
vector m_sigma, m_bc; // Variance components and parameter constraint vectors
|
|
matrix m_vb; // Basis matrix/regressor blocks for variance computations
|
|
HARX* m_obj; // Pointer to the underlying HARX model instance
|
|
bool m_centered, m_individual; // Flags controlling structural type of gradient calculation
|
|
double m_epsilon; // Step size (perturbation value) used for numerical differentiation
|
|
vector obj_result; // Cached tracking vector containing objective evaluation scores
|
|
matrix obj_gradient; // Cached tracking matrix containing the calculated Jacobian/Gradient row blocks
|
|
public:
|
|
//--- Constructor
|
|
CObjectiveF(void) {}
|
|
|
|
//--- Destructor (Handles explicit resource resets and safe pointer detachment)
|
|
~CObjectiveF(void)
|
|
{
|
|
m_obj = NULL;
|
|
obj_result = vector::Zeros(1);
|
|
obj_gradient = matrix::Zeros(1,2);
|
|
}
|
|
|
|
//--- Sets internal parameter constants of the objective function
|
|
void setFuncParams(vector& s, vector& bc, matrix &vb, HARX &obj, bool individual=false, bool center = false, double ep=EMPTY_VALUE)
|
|
{
|
|
m_sigma = s;
|
|
m_bc = bc;
|
|
m_vb = vb;
|
|
m_obj = GetPointer(obj); // Safe encapsulation using structural runtime pointers
|
|
m_individual = individual;
|
|
m_centered = center;
|
|
m_epsilon = ep;
|
|
m_clip = true; // Activates parameter clipping boundaries to prevent parameter explosion
|
|
m_grad_options.method = GRAD_POINT_CALLABLE; // Tells the solver to utilize our explicit analytical/numerical override
|
|
}
|
|
|
|
//--- Objective function wrapper: Returns scalar evaluation ($f(x)$) for optimizer line search steps
|
|
virtual double orig_fun(vector& x) override
|
|
{
|
|
obj_result = m_obj.objective(x,m_sigma,m_bc,m_vb);
|
|
return obj_result[0];
|
|
}
|
|
|
|
//--- Gradient of objective function wrapper: Extracts and returns the local gradient vector ($\nabla f(x)$)
|
|
virtual vector grad_fun(vector& x) override
|
|
{
|
|
obj_gradient = m_obj.jacobian(x,m_sigma,m_bc,m_vb,m_individual,m_centered,m_epsilon);
|
|
return obj_gradient.Row(0); // Isolates the active row mapping the target scalar parameter bounds
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Inequality constraints class |
|
|
//| Formalizes boundary conditions structured as: $b - A \cdot x \ge 0$|
|
|
//+------------------------------------------------------------------+
|
|
class CIneqConstraints: public CConstraints
|
|
{
|
|
private:
|
|
matrix a; // Linear constraint coefficient matrix ($A$)
|
|
vector b; // Boundary threshold vector constants ($b$)
|
|
vector obj_result; // Cached tracking array mapping constraint slack values
|
|
public:
|
|
//--- Constructor
|
|
CIneqConstraints(void)
|
|
{
|
|
}
|
|
|
|
//--- Destructor
|
|
~CIneqConstraints(void)
|
|
{
|
|
}
|
|
|
|
//--- Registers the linear parameters defining the system boundaries
|
|
void setConstraints(matrix& a_, vector& b_)
|
|
{
|
|
a = a_;
|
|
b = b_;
|
|
m_m = int(a.Rows()); // Tracks total count of structural constraint equations
|
|
m_n = int(a.Cols()); // Tracks dimensional parameter inputs expected
|
|
obj_result = vector::Zeros(m_n);
|
|
m_options.method = GRAD_POINT_2; // Configures numerical finite-differencing strategy for boundaries
|
|
}
|
|
|
|
//--- Evaluates constraint spaces ($c(x) = b - A \cdot x$); positive values represent valid spaces
|
|
virtual vector orig_fun(vector& x) override
|
|
{
|
|
if(!a.Rows())
|
|
obj_result = vector::Zeros(0); // Handles unconstrained or open boundary executions safely
|
|
else
|
|
obj_result = b - a.MatMul(x); // Computes the linear coordinate residuals against boundary thresholds
|
|
return obj_result;
|
|
}
|
|
};
|
|
#else
|
|
//+------------------------------------------------------------------+
|
|
//| Objective function adapter for Alglib Optimizer |
|
|
//| Inherits from CNDimensional_Func to map multi-dimensional double |
|
|
//| arrays/rows to the underlying HARX model's math vector objective. |
|
|
//+------------------------------------------------------------------+
|
|
class CObjectiveFunc: public CNDimensional_Func
|
|
{
|
|
private:
|
|
vector m_sigma, m_bc; // Variance properties and optimization boundary limits
|
|
matrix m_vb; // Linear basis data partitions/regressors
|
|
HARX* m_obj; // Structural tracking link pointing back to the core HARX model
|
|
bool m_centered, m_individual; // Spatial properties configuration flags for differentiation
|
|
double m_epsilon; // Perturbation constraint index used during gradient calculations
|
|
public:
|
|
//--- Constructor
|
|
CObjectiveFunc(void) {}
|
|
|
|
//--- Destructor (Safely detaches the internal model pointer reference)
|
|
~CObjectiveFunc(void)
|
|
{
|
|
m_obj = NULL;
|
|
}
|
|
|
|
//--- Configures parameter constants required by the target objective calculation system
|
|
void SetParams(vector& s, vector& bc, matrix &vb, HARX &obj, bool individual=false, bool center = false, double ep=EMPTY_VALUE)
|
|
{
|
|
m_sigma = s;
|
|
m_bc = bc;
|
|
m_vb = vb;
|
|
m_obj = GetPointer(obj);
|
|
m_individual = individual;
|
|
m_centered = center;
|
|
m_epsilon = ep;
|
|
}
|
|
|
|
//--- Evaluates objective function via raw double array conversion interfaces
|
|
virtual void Func(double &x[], double &func, CObject &obj) override
|
|
{
|
|
vector temp;
|
|
temp.Assign(x); // Convert native Alglib parameter array layout into a native math vector
|
|
vector r = m_obj.objective(temp, m_sigma, m_bc, m_vb);
|
|
func = r[0]; // Isolates the primary log-likelihood target scalar metric
|
|
}
|
|
|
|
//--- Overloaded alternative evaluating objective space via Alglib's CRowDouble format
|
|
virtual void Func(CRowDouble &x, double &func, CObject &obj) override
|
|
{
|
|
vector temp = x.ToVector(); // Safe translation layer from specialized Alglib types to math vector profiles
|
|
vector r = m_obj.objective(temp, m_sigma, m_bc, m_vb);
|
|
func = r[0];
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Jacobian fvec vector implementation for Alglib |
|
|
//| Inherits from CNDimensional_Jac to handle simultaneous evaluation |
|
|
//| of the residual vector ($f_{vec}$) and its first-order derivative matrix. |
|
|
//+------------------------------------------------------------------+
|
|
class CObjectiveJacFVec: public CNDimensional_Jac
|
|
{
|
|
private:
|
|
vector m_sigma, m_bc; // Parameter limits and target variance matrices
|
|
matrix m_vb; // Basis matrix array partitions
|
|
HARX* m_obj; // Structural pointer pointing to the active execution engine
|
|
bool m_centered, m_individual; // Gradient computation alignment strategies
|
|
double m_epsilon; // Scaling value used for numeric parameter mutations
|
|
public:
|
|
//--- Constructor
|
|
CObjectiveJacFVec(void) {}
|
|
|
|
//--- Destructor
|
|
~CObjectiveJacFVec(void)
|
|
{
|
|
m_obj = NULL;
|
|
}
|
|
|
|
//--- Binds optimization parameters to internal tracker fields
|
|
void SetParams(vector& s, vector& bc, matrix &vb, HARX &obj, bool individual=false, bool center = false, double ep=EMPTY_VALUE)
|
|
{
|
|
m_sigma = s;
|
|
m_bc = bc;
|
|
m_vb = vb;
|
|
m_obj = GetPointer(obj);
|
|
m_individual = individual;
|
|
m_centered = center;
|
|
m_epsilon = ep;
|
|
}
|
|
|
|
//--- Evaluates both individual residual rows ($f_i$) and the Jacobian matrix simultaneously using native arrays
|
|
virtual void Jac(double &x[], double &fi[], CMatrixDouble &jac, CObject &obj) override
|
|
{
|
|
vector temp;
|
|
temp.Assign(x);
|
|
|
|
//--- Evaluate objective baseline at the current spatial parameters coordinate
|
|
vector r = m_obj.objective(temp, m_sigma, m_bc, m_vb);
|
|
fi[0] = r[0];
|
|
|
|
//--- Calculate the system Jacobian matrix ($\mathbf{J}_{i,j} = \frac{\partial f_i}{\partial x_j}$)
|
|
matrix mjac = m_obj.jacobian(temp, m_sigma, m_bc, m_vb, m_individual, m_centered, m_epsilon);
|
|
jac = mjac; // Implicit casting assignment to Alglib's internal CMatrixDouble format
|
|
}
|
|
|
|
//--- Overloaded variant mapping spatial dimensions directly from Alglib CRowDouble components
|
|
virtual void Jac(CRowDouble &x, CRowDouble &fi, CMatrixDouble &jac, CObject &obj) override
|
|
{
|
|
vector xv = x.ToVector();
|
|
|
|
//--- Evaluate the multidimensional residual mapping vector block
|
|
vector temp = m_obj.objective(xv, m_sigma, m_bc, m_vb);
|
|
fi = temp; // Directly overwrite the multi-row output reference structure
|
|
|
|
//--- Compute the matrix partial derivatives layout block
|
|
matrix mjac = m_obj.jacobian(xv, m_sigma, m_bc, m_vb, m_individual, m_centered, m_epsilon);
|
|
jac = mjac;
|
|
}
|
|
};
|
|
#endif
|
|
int m_extra_simulation_params;
|
|
//+------------------------------------------------------------------+
|
|
//| Calculate the Coefficient of Determination (R-squared) |
|
|
//| Computes $R^2 = 1 - \frac{SSR}{TSS}$ to determine the percentage |
|
|
//| of explained variance captured by the structural model. |
|
|
//+------------------------------------------------------------------+
|
|
virtual double _r2(vector ¶ms) override
|
|
{
|
|
vector y = m_fit_y;
|
|
bool constant = false;
|
|
matrix x = m_fit_regressors;
|
|
|
|
//--- Check for an explicit or implicit structural model intercept constant
|
|
if(x.Cols())
|
|
constant = (m_constant||implicit_constant(x));
|
|
|
|
//--- Center the dependent vector around its mean if a constant exists
|
|
if(constant)
|
|
{
|
|
//--- If the matrix is intercept-only, structural variance cannot be explained ($R^2 = 0$)
|
|
if(x.Cols() ==1)
|
|
return 0.0;
|
|
y = y - y.Mean(); // Transforms the vector to compute centered Sum of Squares
|
|
}
|
|
|
|
//--- Calculate Total Sum of Squares ($TSS = \sum (y_i - \bar{y})^2$) via vector dot-product
|
|
double tss = y.Dot(y);
|
|
if(tss<=0.)
|
|
return AL_NaN; // Guard against division by zero errors on flat data streams
|
|
|
|
//--- Extract model tracking errors / residuals ($\epsilon$)
|
|
vector e = resids(params,EMPTY_VECTOR,EMPTY_MATRIX);
|
|
|
|
//--- Return final metric: $R^2 = 1 - \frac{Residual\ Sum\ of\ Squares\ (e \cdot e)}{Total\ Sum\ of\ Squares\ (TSS)}$
|
|
return 1.-(e.Dot(e))/tss;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Dynamic sample framing array adjuster |
|
|
//| Re-calibrates active training timeline windows, accounting for |
|
|
//| designated lag pre-sample allocations or back-testing holdbacks. |
|
|
//+------------------------------------------------------------------+
|
|
bool _adjust_sample(long first_obs, long last_obs)
|
|
{
|
|
//--- Offset start boundary forward to shift past the warm-up/holdback lag horizon
|
|
first_obs+=long(m_holdback);
|
|
if(last_obs<= first_obs)
|
|
{
|
|
Print(__FUNCTION__, " invalid indices specified ");
|
|
return false;
|
|
}
|
|
|
|
//--- Cache target data window indices inside the fit index tracking vector
|
|
vector f = {double(first_obs),double(last_obs)};
|
|
m_fit_indices = f;
|
|
|
|
//--- Slice and align the dependent and independent series to the active estimation block
|
|
m_fit_y = adjust_y();
|
|
m_fit_x = adjust_x();
|
|
|
|
//--- Extract active explanatory feature rows (regressors) from the master collection matrix
|
|
m_fit_regressors = (m_regressors.Rows())?np::sliceMatrixRows(m_regressors,long(m_fit_indices[0]),long(m_fit_indices[1])):matrix::Zeros(0,0);
|
|
|
|
//--- Sync structural pointer/state indicators for the local model processing space
|
|
m_vp.start(long(first_obs));
|
|
m_vp.stop(long(last_obs));
|
|
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Adjust Dependent Variable Vector (y) |
|
|
//| Slices the master observation stream to isolate the active |
|
|
//| sample timeframe mapped by the current model estimation window. |
|
|
//+------------------------------------------------------------------+
|
|
vector adjust_y(void)
|
|
{
|
|
//--- Extract and slice data if the master observations vector contains entries
|
|
if(m_model_spec.observations.Size())
|
|
return np::sliceVector(m_model_spec.observations,long(m_fit_indices[0]),long(m_fit_indices[1]));
|
|
return vector::Zeros(0); // Return an empty vector frame if no observations exist
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Adjust Exogenous Regressor Matrix (X) |
|
|
//| Slices the master exogenous data matrix by either rows or |
|
|
//| columns to align tracking structures with the active sample window|
|
|
//+------------------------------------------------------------------+
|
|
matrix adjust_x(bool rows = true)
|
|
{
|
|
//--- Proceed only if the master matrix container contains data rows
|
|
if(m_model_spec.exog_data.Rows())
|
|
{
|
|
//--- Slice horizontally (default) to extract active temporal periods/rows
|
|
if(rows)
|
|
return np::sliceMatrixRows(m_model_spec.exog_data,long(m_fit_indices[0]),long(m_fit_indices[1]));
|
|
//--- Otherwise slice vertically to extract specific structural features/columns
|
|
else
|
|
return np::sliceMatrixCols(m_model_spec.exog_data,long(m_fit_indices[0]),long(m_fit_indices[1]));
|
|
}
|
|
else
|
|
return matrix::Zeros(0,0); // Return empty 0x0 matrix if no features exist
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Compute Baseline Parameters under Homoskedastic Normal Errors |
|
|
//| Generates starting parameters for the conditional mean equation |
|
|
//| via OLS using the Moore-Penrose Pseudo-Inverse ($\mathbf{X}^+$). |
|
|
//+------------------------------------------------------------------+
|
|
vector _fit_no_arch_normal_errors_params(void)
|
|
{
|
|
vector y = m_fit_y;
|
|
ulong nobs = y.Size();
|
|
|
|
//--- Data Check: Enforce degrees-of-freedom requirements ($N > K$)
|
|
if(nobs<m_num_params)
|
|
{
|
|
Print(__FUNCTION__, " insufficient data ");
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
matrix x = m_fit_regressors;
|
|
if(x.Cols())
|
|
{
|
|
ResetLastError();
|
|
|
|
//--- Compute Moore-Penrose Pseudo-Inverse: $\mathbf{X}^+ = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T$
|
|
matrix xpinv = np::pinv(x);
|
|
|
|
//--- Trapping singular matrices or calculation underflows
|
|
if(xpinv.HasNan())
|
|
{
|
|
Print(__FUNCTION__, " PInv() returned invalid number ", GetLastError());
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
//--- Calculate OLS coefficients: $\hat{\beta} = \mathbf{X}^+ \cdot \mathbf{y}$
|
|
return np::matmul(xpinv,y);
|
|
}
|
|
else
|
|
return vector::Zeros(0); // Return empty if there are no regressors in the system
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Estimates baseline model parameters with homoskedastic errors |
|
|
//| Fits OLS parameters, extracts residuals, and builds information |
|
|
//| matrices for either standard (Classic) or White's Robust errors. |
|
|
//+------------------------------------------------------------------+
|
|
ArchModelResult _fit_no_arch_normal_errors(ENUM_COVAR_TYPE cov_type = COVAR_ROBUST)
|
|
{
|
|
ArchModelResult out;
|
|
matrix x = m_fit_regressors;
|
|
vector y = m_fit_y;
|
|
vector fitted, rparams;
|
|
fitted = rparams = vector::Zeros(0);
|
|
matrix xpxi = matrix::Zeros(0,0);
|
|
|
|
//--- 1. Parameter Estimation via OLS
|
|
if(x.Cols())
|
|
{
|
|
matrix xpinv = np::pinv(x); // Moore-Penrose Pseudo-Inverse
|
|
rparams = np::matmul(xpinv,y); // $\hat{\beta} = \mathbf{X}^+ \cdot \mathbf{y}$
|
|
matrix xt = (x.Transpose().MatMul(x))/double(y.Size()); // Normalized $\frac{\mathbf{X}^T\mathbf{X}}{T}$
|
|
xpxi = xt.Inv(); // $(\mathbf{X}^T\mathbf{X})^{-1}$ tracking block
|
|
fitted = x.MatMul(rparams); // Predicted conditional mean outputs
|
|
}
|
|
|
|
//--- 2. Extract Residual Errors ($\epsilon$) and Variance Vector ($\sigma^2$)
|
|
vector e = np::minus(y,fitted);
|
|
double sigma2 = e.Dot(e)/double(y.Size()); // Homoskedastic variance estimate
|
|
vector params = rparams;
|
|
|
|
//--- Append the variance estimate ($\sigma^2$) as the final structural parameter row
|
|
rparams.Resize(params.Size()+1);
|
|
rparams[params.Size()] = sigma2;
|
|
|
|
//--- 3. Construct the Analytical Log-Likelihood Information Hessian Matrix
|
|
matrix hessian = matrix::Zeros(m_num_params+1,m_num_params+1);
|
|
np::matrixCopy(hessian,-1.*xpxi,0,long(m_num_params),1,0,long(m_num_params));
|
|
hessian[hessian.Rows()-1,hessian.Cols()-1] = -1.; // Map scale parameter derivative threshold
|
|
|
|
matrix param_cov = matrix::Zeros(0,0);
|
|
string nm;
|
|
|
|
//--- 4. Structural Covariance Type Resolution
|
|
switch(cov_type)
|
|
{
|
|
//--- Standard Asymptotic Covariance Matrix Evaluation
|
|
case COVAR_CLASSIC:
|
|
param_cov = sigma2 *(-1.*hessian);
|
|
param_cov[m_num_params,m_num_params] = 2.*pow(sigma2,2.); // Asymptotic variance of error variance
|
|
param_cov/=double(y.Size()); // Scale matrix by sample depth ($T$)
|
|
nm = "classic_ols";
|
|
break;
|
|
|
|
//--- White's Heteroskedasticity-Consistent (Sandwich) Covariance Estimator
|
|
case COVAR_ROBUST:
|
|
{
|
|
// Initialize score matrix to hold individual structural gradients/moment conditions
|
|
matrix scores = matrix::Zeros(y.Size(),m_num_params+1);
|
|
for(ulong i =0; i<m_num_params; ++i)
|
|
scores.Col(e*x.Col(i),i); // Mean equation parameter scores: $\epsilon_t \cdot \mathbf{X}_{t,i}$
|
|
|
|
// Variance parameter score column: $\epsilon_t^2 - \sigma^2$
|
|
scores.Col(pow(e,2.)-sigma2,scores.Cols()-1);
|
|
|
|
// Outer product of scores matrix ($\mathbf{B} = \frac{\mathbf{G}^T\mathbf{G}}{T}$)
|
|
matrix scorescov = scores.Transpose().MatMul(scores)/double(y.Size());
|
|
|
|
// Core Sandwich Assembly: $\mathbf{A}^{-1} \mathbf{B} \mathbf{A}^{-1}$
|
|
param_cov = hessian.MatMul(scorescov);
|
|
param_cov = param_cov.MatMul(hessian);
|
|
param_cov/=double(y.Size()); // Divide by overall sample length $T$
|
|
nm = "white";
|
|
}
|
|
break;
|
|
}
|
|
|
|
//--- 5. Diagnostics, Vector Alignment, and Output Structure Packaging
|
|
double r2 = _r2(params); // Compute R-squared metric
|
|
|
|
// Prepare full-sized array allocations aligned with the complete master timeline
|
|
vector resids_ = vector::Zeros(m_model_spec.observations.Size());
|
|
resids_.Fill(AL_NaN); // Pad unestimated rows with NaN
|
|
np::vectorCopy(resids_,e,long(m_fit_indices[0]),long(m_fit_indices[1])); // Splice in active window errors
|
|
|
|
vector vol = vector::Zeros(resids_.Size());
|
|
vol.Fill(AL_NaN);
|
|
double ss = sqrt(sigma2); // Calculate constant standard deviation ($\sigma$)
|
|
|
|
// Populate the volatility vector for the active estimation window slice
|
|
for(ulong i = ulong(m_fit_indices[0]); i<ulong(m_fit_indices[1]); vol[i] = ss, ++i);
|
|
|
|
// Compute final Gaussian Log-Likelihood objective score
|
|
double ll = _static_gaussian_loglikelihood(e);
|
|
|
|
// Save state variables to internal class properties
|
|
m_params = out.params = rparams;
|
|
m_converged = (m_params.Size()>0);
|
|
out.param_cov = param_cov;
|
|
out.r2 = r2;
|
|
out.resid = resids_;
|
|
out.conditional_volatility = vol;
|
|
out.loglikelihood = ll;
|
|
out.nobs = e.Size();
|
|
out.fit_indices[0] = long(m_fit_indices[0]);
|
|
out.fit_indices[1] = long(m_fit_indices[1]);
|
|
|
|
return out;
|
|
}
|
|
//--- Should be called whenever the model is initialized or changed
|
|
bool _init_model(void)
|
|
{
|
|
ResetLastError();
|
|
|
|
//--- 1. Process and Reformat Model Lag Specifications
|
|
if(m_model_spec.mean_lags.Size())
|
|
{
|
|
if(!_reformat_lags())
|
|
return false;
|
|
|
|
// Determine the largest autoregressive lag index to protect pre-sample sequences
|
|
m_maxlags = ulong(m_model_spec.mean_lags.Max());
|
|
|
|
// Enforce a minimum data holdback sample size to prevent index-out-of-bounds errors during lag construction
|
|
if(!m_holdback)
|
|
m_holdback = m_maxlags;
|
|
else
|
|
{
|
|
if(m_holdback<m_maxlags)
|
|
{
|
|
Print(__FUNCTION__, " hold_back is less then the minimum number given the lags selected");
|
|
m_holdback = m_maxlags;
|
|
}
|
|
}
|
|
}
|
|
|
|
//--- 2. Validate structural constraints across model options
|
|
if(!_check_specification())
|
|
return false;
|
|
|
|
ulong nobs_orig = m_model_spec.observations.Size();
|
|
vector reg_constant = vector::Zeros(0);
|
|
matrix reg_lags = matrix::Zeros(0,0);
|
|
|
|
//--- 3. Construct the Intercept Column vector (ones vector)
|
|
if(m_constant)
|
|
reg_constant = vector::Ones(nobs_orig);
|
|
|
|
//--- 4. Build Autoregressive Lag Arrays via lagmat
|
|
if(m_lags.Rows() && nobs_orig)
|
|
{
|
|
double maxlag = m_lags.Max();
|
|
matrix ymat(m_model_spec.observations.Size(),1);
|
|
ymat.Col(m_model_spec.observations,0);
|
|
matrix lagarray[];
|
|
|
|
// Build the pre-sample lag matrix pool (TRIM_FORWARD protects matrix dimensions)
|
|
if(!lagmat(ymat,lagarray,ulong(maxlag),TRIM_FORWARD,ORIGINAL_EX))
|
|
return false;
|
|
|
|
// Extract and aggregate specific lag column matrices (averaging or slicing subsets)
|
|
reg_lags = matrix::Zeros(nobs_orig,m_lags.Cols());
|
|
for(ulong i = 0; i<m_lags.Cols(); ++i)
|
|
{
|
|
matrix temp = np::sliceMatrixCols(lagarray[0],long(m_lags[0,i]),long(m_lags[1,i]));
|
|
vector vtemp = temp.Mean(1); // Compute mean across specified lag bounds if grouping is applied
|
|
reg_lags.Col(vtemp,i);
|
|
}
|
|
}
|
|
|
|
//--- 5. Initialize the Consolidated Design Matrix Container (m_regressors)
|
|
// Total column allocations: Constant (0 or 1) + AR Lag columns + Exogenous Feature columns
|
|
m_regressors = matrix::Zeros(nobs_orig,ulong(nobs_orig>0)*(ulong(m_constant)+reg_lags.Cols()+m_model_spec.exog_data.Cols()));
|
|
|
|
// Slot in the model intercept constant vector block at column index 0
|
|
if(m_constant && nobs_orig)
|
|
m_regressors.Col(reg_constant,0);
|
|
|
|
// Copy the constructed autoregressive series variables sequentially into adjacent channels
|
|
if(reg_lags.Cols())
|
|
for(ulong i = 0; i<reg_lags.Cols(); ++i)
|
|
m_regressors.Col(reg_lags.Col(i),i+ulong(m_constant));
|
|
|
|
// Append external exogenous design elements into the remaining trailing columns
|
|
if(m_model_spec.exog_data.Cols())
|
|
for(ulong i = 0; i<m_model_spec.exog_data.Cols(); ++i)
|
|
m_regressors.Col(m_model_spec.exog_data.Col(i),i+ulong(m_constant)+reg_lags.Cols());
|
|
|
|
//--- 6. Structural error boundary check before unlocking
|
|
if(GetLastError()!=0)
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Scale Change Event Listener Override |
|
|
//| Automatically re-initializes and updates the structural design |
|
|
//| matrix layouts if data scales or internal parameters are altered. |
|
|
//+------------------------------------------------------------------+
|
|
virtual void _scalechanged(void) override
|
|
{
|
|
//--- Trigger an inline model refresh to rebuild vectors and track matrix sizing
|
|
if(!_init_model())
|
|
Print(__FUNCTION__, " error ");
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Validate Model Structural Specifications |
|
|
//| Performs checking boundaries on data shapes to confirm exogenous |
|
|
//| matrix lengths perfectly align with the core observation series. |
|
|
//+------------------------------------------------------------------+
|
|
bool _check_specification(void)
|
|
{
|
|
//--- Evaluate matrix dimensional integrity if exogenous data blocks exist
|
|
if(m_model_spec.exog_data.Rows())
|
|
{
|
|
//--- Dimensional Gate: Enforce structural equality along the time axis (Rows == N)
|
|
if(m_model_spec.exog_data.Rows()!=m_model_spec.observations.Size())
|
|
{
|
|
Print(__FUNCTION__, " regressors shape does not match input data ");
|
|
return false; // Rejects initialization due to broken sample alignments
|
|
}
|
|
}
|
|
return true; // Pass specifications check
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Reformat input lags to simplify matrix operations |
|
|
//| Standardizes a 1-row lag vector or 2-row interval array into a |
|
|
//| normalized indexing matrix and validates entry uniqueness via rank.|
|
|
//+------------------------------------------------------------------+
|
|
bool _reformat_lags(void)
|
|
{
|
|
//--- If no lags are specified, bypass processing cleanly
|
|
if(!m_lags.Rows())
|
|
return true;
|
|
|
|
matrix lags = m_lags;
|
|
|
|
//--- Case 1: Input is a simple 1-Row array of specific lag milestones
|
|
if(lags.Rows()==1)
|
|
{
|
|
// Validate data inputs; lag values must be positive integers
|
|
if(lags.Min()<=0)
|
|
{
|
|
Print(__FUNCTION__, " invalid lags input ");
|
|
return false;
|
|
}
|
|
|
|
// Extract row, remove duplicate instances, and sort ascendingly
|
|
vector row = lags.Row(0);
|
|
row = np::unique(row);
|
|
|
|
// Convert to a 2-row matrix to represent continuous intervals [Start, Stop]
|
|
matrix temp = matrix::Zeros(2,row.Size());
|
|
if(!temp.Row(row,0) || !temp.Row(row,1))
|
|
return false;
|
|
|
|
//--- Process Rotated Lags: Maps cumulative/overlapping periods (e.g., HAR-Daily/Weekly/Monthly)
|
|
if(m_rotated)
|
|
{
|
|
vector v = np::sliceVector(row,0,long(row.Size()-1));
|
|
for(ulong i = 1; i<temp.Cols(); temp[0,i] = v[i-1], ++i);
|
|
temp[0,0] = 0.; // Sets base offset for the initial lag index
|
|
}
|
|
//--- Process Standard Lags: Disjoint individual time steps
|
|
else
|
|
temp.Row(vector::Zeros(temp.Cols()),0); // Base remains zero for direct offset parsing
|
|
|
|
m_lags = temp;
|
|
}
|
|
//--- Case 2: Input is a 2-Row matrix specifying range intervals [Start, Stop]
|
|
else
|
|
{
|
|
vector mins = lags.Min(1);
|
|
// Validate that all bounds are positive and upper boundaries exceed lower bounds
|
|
if(lags.Min()<=0 || mins[1]<mins[0])
|
|
{
|
|
Print(__FUNCTION__, " invalid lags input ");
|
|
return false;
|
|
}
|
|
|
|
// Normalize structure directions across rows and sort columns
|
|
if(!np::reverseMatrixRows(lags) || !np::sort(lags,true))
|
|
return false;
|
|
|
|
//--- Validate Multicollinearity / Redundant overlapping intervals
|
|
// Build a binary indicator profile matrix matching the max lag width
|
|
matrix testmat = matrix::Zeros(lags.Cols(),(ulong)lags.Max());
|
|
vector subtract = {1.,0.};
|
|
matrix sub = np::repeat_vector_as_rows_cols(subtract,lags.Cols(),false);
|
|
lags = lags-sub; // Zero-index adjustment subtraction layer
|
|
|
|
// Fill matrix coordinates representing the active ranges
|
|
for(ulong i = 0; i<lags.Cols(); ++i)
|
|
np::matrixFill(testmat,1.0,long(i),long(i+1),STEP,long(lags[0,i]),long(lags[1,i]));
|
|
|
|
// Check rank condition: Number of non-zero singular vectors must equal total column dimensions
|
|
ulong rank = np::matrix_rank(testmat);
|
|
if(rank != lags.Cols())
|
|
{
|
|
Print(__FUNCTION__," lags contains redundant entries ");
|
|
return false; // Rejects collinear/overlapping definitions
|
|
}
|
|
|
|
if(m_rotated)
|
|
Print(__FUNCTION__, " Rotation is not available for this model specification ");
|
|
|
|
m_lags = lags;
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Reformat and validate exogenous forecast parameters |
|
|
//| Enforces structural integrity on a 3-D container (array of |
|
|
//| matrices) containing multi-step exogenous values for forecasting. |
|
|
//+------------------------------------------------------------------+
|
|
bool _reformat_forecast_x(ulong start, ulong horizon, matrix& x[], matrix& out[])
|
|
{
|
|
//--- 1. Verification Gate: Confirm presence/absence of exogenous structures aligns with specification
|
|
if(!x.Size())
|
|
{
|
|
if(!m_model_spec.exog_data.Rows())
|
|
return true; // No exogenous data expected or provided; bypass cleanly
|
|
else
|
|
{
|
|
Print(__FUNCTION__, " error x function variable should not be empty since model specified with exogenous variables ");
|
|
return false;
|
|
}
|
|
}
|
|
else
|
|
if(!m_model_spec.exog_data.Rows())
|
|
{
|
|
Print(__FUNCTION__, " x is not None but the model does not contain any exogenous variables.");
|
|
return false;
|
|
}
|
|
|
|
//--- 2. Dimensional Gates: Check feature space and step horizon matches
|
|
ulong nx = m_model_spec.exog_data.Cols();
|
|
if(x.Size()!=uint(nx))
|
|
{
|
|
Print(__FUNCTION__, " the number of elements in x should match number of columns of m_model_spec.exog_data property ");
|
|
return false;
|
|
}
|
|
if(x[0].Cols()!=horizon)
|
|
{
|
|
Print(__FUNCTION__, " the number of values passed does not match the horizon of the forecasts ");
|
|
return false;
|
|
}
|
|
|
|
//--- 3. Validation: Verify active row depth conforms to the forecasting window index space
|
|
if(x[0].Rows()!=(m_model_spec.observations.Size()-start))
|
|
{
|
|
Print(__FUNCTION__," the shape of x does not have the proper dimensions ");
|
|
return false;
|
|
}
|
|
|
|
//--- 4. Memory Preparation and Row-Slicing Alignment
|
|
ArrayResize(out,x.Size());
|
|
|
|
// If the input matrices contain extra historical rows, slice starting from the active 'start' forecast index
|
|
if(x[0].Rows()>(m_model_spec.observations.Size()-start))
|
|
{
|
|
for(uint i = 0; i<x.Size(); ++i)
|
|
out[i] = np::sliceMatrixRows(x[i],long(start));
|
|
}
|
|
// Otherwise, the arrays are already perfectly trimmed; perform a direct copy swap
|
|
else
|
|
{
|
|
for(uint i = 0; i<x.Size(); ++i)
|
|
out[i] = x[i];
|
|
}
|
|
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Model transformation: HAR parameters to standard AR parameters |
|
|
//| Distributes cumulative HAR lag weights (e.g., weekly or monthly |
|
|
//| averages) uniformly across standard daily AR lag structures. |
|
|
//+------------------------------------------------------------------+
|
|
vector _har_to_ar(vector ¶ms)
|
|
{
|
|
//--- If no lags exist, isolate and return the intercept constant parameter block
|
|
if(!m_maxlags)
|
|
return np::sliceVector(params,0,long(m_constant));
|
|
|
|
//--- Strip away the intercept (if present) to isolate the raw HAR coefficients
|
|
vector har = np::sliceVector(params,long(m_constant));
|
|
vector ar = vector::Zeros(m_maxlags);
|
|
vector temp;
|
|
|
|
//--- 1. Linear Mapping: Disentangle intervals and distribute weights
|
|
for(ulong i = 0; i<m_lags.Cols(); ++i)
|
|
{
|
|
// Distribute the aggregated HAR weight uniformly across its interval span ($W_j = \phi_{HAR} / \Delta t$)
|
|
for(ulong j = ulong(m_lags[0,i]); j<ulong(m_lags[1,i]); ++j)
|
|
ar[j] += har[i]/(m_lags[1,i] - m_lags[0,i]);
|
|
}
|
|
|
|
//--- 2. Structure Reassembly: Re-attach the model intercept constant
|
|
if(m_constant)
|
|
{
|
|
vector a(ar.Size()+1);
|
|
a[0] = params[0]; // Set the intercept constant at position index 0
|
|
|
|
// Shift original standard AR lag values back into chronological vector channels
|
|
for(ulong i = 1; i<a.Size(); a[i] = ar[i-1], ++i);
|
|
ar = a;
|
|
}
|
|
|
|
return ar;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Runs a recursive simulation loop to calculate the model mean |
|
|
//| Generates a structural time-series sequence by combining path |
|
|
//| initializations, lag dependencies, exogenous variables, and errors.|
|
|
//+------------------------------------------------------------------+
|
|
vector _simulate_mean(vector ¶meters, matrix &x, vector &errs, vector& initial_values, vector& cond_var)
|
|
{
|
|
ulong max_lag = (!m_lags.Rows())?0:(ulong)m_lags.Max();
|
|
ulong nobs_and_burn = errs.Size();
|
|
vector y = vector::Zeros(nobs_and_burn);
|
|
vector iv;
|
|
ulong iv_size = initial_values.Size();
|
|
|
|
//--- 1. Condition and Seed Pre-Sample Initialization Arrays
|
|
if(iv_size && iv_size!=max_lag)
|
|
{
|
|
Print(__FUNCTION__, " initial_value has the wrong shape ");
|
|
return vector::Zeros(0);
|
|
}
|
|
else
|
|
if(!iv_size)
|
|
iv.Resize(max_lag); // If omitted, defaults to an array of zeros
|
|
else
|
|
iv = initial_values;
|
|
|
|
// Warm up the beginning of the time-series array with specified seed histories
|
|
for(ulong i = 0; i<max_lag; y[i] = iv[i],++i);
|
|
ulong kx = x.Cols();
|
|
|
|
ulong ind;
|
|
vector temp, col;
|
|
|
|
//--- 2. Main Chronological Simulation Iteration Loop
|
|
for(ulong t = max_lag; t<nobs_and_burn; ++t)
|
|
{
|
|
ind = 0; // Tracks parameter index mapping across structural boundaries
|
|
|
|
// Phase A: Apply Intercept Model Constants ($\mu$)
|
|
if(m_constant)
|
|
{
|
|
y[t] = parameters[ind];
|
|
ind+=1;
|
|
}
|
|
|
|
// Phase B: Compute Autoregressive / HAR Interval Lag Components
|
|
for(ulong lag = 0; lag<m_lags.Cols(); ++lag)
|
|
{
|
|
col = m_lags.Col(lag); // Extract the interval bounds [Start, Stop]
|
|
|
|
// Slice the active window tracking period chronologically relative to step 't'
|
|
temp = np::sliceVector(y,long(t - (ulong)col[1]), long(t - ulong(col[0])));
|
|
|
|
// Apply structural coefficient against the slice mean (supports standard and grouped/HAR components)
|
|
y[t] += parameters[ind] * temp.Mean();
|
|
ind+=1;
|
|
}
|
|
|
|
// Phase C: Incorporate Exogenous Regressor Linear Metrics ($\sum \beta_i X_{t,i}$)
|
|
for(ulong i = 0; i<kx; ++i)
|
|
{
|
|
y[t] += parameters[ind] * x[t,i];
|
|
ind += 1;
|
|
}
|
|
|
|
// Phase D: Inject Random/Conditional Error Innovations ($\epsilon_t$)
|
|
y[t] += errs[t];
|
|
}
|
|
|
|
return y;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Fit model without parameters |
|
|
//| Bypasses the optimizer when zero parameters require estimations. |
|
|
//| Directly extracts raw residuals and profiles the variance series.|
|
|
//+------------------------------------------------------------------+
|
|
ArchModelResult _fit_parameterless_model(ENUM_COVAR_TYPE cv, vector& backcast)
|
|
{
|
|
ArchModelResult out;
|
|
|
|
vector y = m_fit_y;
|
|
vector params = vector::Zeros(0); // No parameters to optimize
|
|
matrix param_cov = matrix::Zeros(0,0); // Covariance matrix is empty
|
|
ulong first_obs = (ulong)m_fit_indices[0];
|
|
ulong last_obs = (ulong)m_fit_indices[1];
|
|
|
|
//--- 1. Isolate and Extract Model Residuals
|
|
vector resids_final = vector::Zeros(y.Size());
|
|
if(!np::vectorCopy(resids_final, y, first_obs, last_obs))
|
|
{
|
|
Print(__FUNCTION__, " failed vector copy ");
|
|
return out;
|
|
}
|
|
|
|
//--- 2. Compute Target Variance Boundaries and Structural Conditions
|
|
matrix var_bounds = m_vp.varianceBounds(y);
|
|
vector vol = vector::Zeros(y.Size());
|
|
m_vp.computeVariance(params, y, vol, backcast, var_bounds);
|
|
// Transform raw variance values into volatility standard deviations ($\sigma_t = \sqrt{\sigma_t^2}$)
|
|
vol = sqrt(vol);
|
|
|
|
//--- 3. Slice and Align Volatility Elements inside the Timeline Window
|
|
vector vol_final = vector::Zeros(y.Size());
|
|
vol_final.Fill(AL_NaN); // Initialize array with NaN to pad unestimated spaces
|
|
if(!np::vectorCopy(vol_final, vol, first_obs, last_obs))
|
|
{
|
|
Print(__FUNCTION__, " failed vector copy ");
|
|
return out;
|
|
}
|
|
|
|
//--- 4. Calculate Summary Diagnostics and Core Model Log-Likelihood
|
|
double r2 = _r2(params);
|
|
vector ll = -1.0 * _loglikelihood(params, pow(vol, 2) * vector::Ones(ulong(m_fit_indices[1] - m_fit_indices[0])), backcast, var_bounds);
|
|
|
|
//--- 5. Package and Deliver Output Result Structures
|
|
out.params = params;
|
|
out.param_cov = param_cov;
|
|
out.r2 = r2;
|
|
out.resid = resids_final;
|
|
out.conditional_volatility = vol_final;
|
|
out.cov_type = cv;
|
|
out.loglikelihood = ll[0];
|
|
out.fit_indices[0] = long(m_fit_indices[0]);
|
|
out.fit_indices[1] = long(m_fit_indices[1]);
|
|
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Hidden class initialization routine |
|
|
//| Configures the complete execution framework by validating the |
|
|
//| mean specification, instantiating the chosen volatility process, |
|
|
//| and setting up the target innovation error distribution model. |
|
|
//+------------------------------------------------------------------+
|
|
bool _initialize(ArchParameters &vol_dist_params)
|
|
{
|
|
//--- 1. Verification Gate: Enforce structural class and specification alignment
|
|
if(m_model_spec.mean_model_type!=WRONG_VALUE && vol_dist_params.mean_model_type!=WRONG_VALUE && vol_dist_params.mean_model_type!=m_model_spec.mean_model_type)
|
|
{
|
|
m_initialized = false;
|
|
Print(__FUNCTION__" Incorrect initialization of mean model. \nMake sure input parameters correspond with the correct class ");
|
|
return false;
|
|
}
|
|
|
|
// Sync the specifications type if the target parameters structure is uninitialized
|
|
if(m_model_spec.mean_model_type!=WRONG_VALUE && vol_dist_params.mean_model_type==WRONG_VALUE)
|
|
vol_dist_params.mean_model_type = m_model_spec.mean_model_type;
|
|
|
|
// Bind the incoming parameter structure to the internal master specification state
|
|
m_model_spec = vol_dist_params;
|
|
|
|
//--- 2. Mean Model Configuration State Machine
|
|
switch(m_model_spec.mean_model_type)
|
|
{
|
|
case MEAN_CONSTANT:
|
|
m_model_spec.mean_lags = vector::Zeros(0);
|
|
m_name = "Constant Mean";
|
|
m_model_spec.include_constant = true;
|
|
m_model_spec.use_har_rotation = false;
|
|
break;
|
|
case MEAN_AR:
|
|
m_model_spec.use_har_rotation = false;
|
|
m_model_spec.exog_data = matrix::Zeros(0,0);
|
|
m_name = "AR";
|
|
break;
|
|
case MEAN_ARX:
|
|
m_model_spec.use_har_rotation = false;
|
|
m_name ="AR-X";
|
|
break;
|
|
case MEAN_HAR:
|
|
m_model_spec.exog_data = matrix::Zeros(0,0);
|
|
m_name = "HAR";
|
|
break;
|
|
case MEAN_HARX:
|
|
m_name = "HAR-X";
|
|
break;
|
|
case MEAN_ZERO:
|
|
m_model_spec.exog_data=matrix::Zeros(0,0);
|
|
m_model_spec.mean_lags = vector::Zeros(0);
|
|
m_model_spec.include_constant = false;
|
|
m_model_spec.use_har_rotation = false;
|
|
m_name = "Zero Mean";
|
|
break;
|
|
default:
|
|
m_name = "HAR-X";
|
|
break;
|
|
}
|
|
|
|
//--- 3. Set Timeline Estimation Boundaries and Format Structural Lags
|
|
m_fit_indices = vector::Zeros(2);
|
|
m_fit_indices[1] = (m_model_spec.observations.Size())?double(m_model_spec.observations.Size()):-1.;
|
|
|
|
if(m_model_spec.mean_lags.Size())
|
|
{
|
|
switch(m_model_spec.mean_model_type)
|
|
{
|
|
// Standard Autoregressive: Format structural arrays for point-to-point lag entries
|
|
case MEAN_AR:
|
|
case MEAN_ARX:
|
|
m_lags = matrix::Zeros(2,m_model_spec.mean_lags.Size());
|
|
m_lags.Row(m_model_spec.mean_lags,0);
|
|
m_lags.Row(m_model_spec.mean_lags,1);
|
|
break;
|
|
// Heterogeneous Autoregressive: Initialize arrays to track multi-day windows
|
|
case MEAN_HAR:
|
|
case MEAN_HARX:
|
|
m_lags = matrix::Zeros(1,m_model_spec.mean_lags.Size());
|
|
m_lags.Row(m_model_spec.mean_lags,0);
|
|
break;
|
|
}
|
|
}
|
|
else
|
|
m_lags = matrix::Zeros(0,0);
|
|
|
|
// Sync local runtime toggles with the verified specification snapshot
|
|
m_constant = m_model_spec.include_constant;
|
|
m_rescale = m_model_spec.is_rescale_enabled;
|
|
m_rotated = m_model_spec.use_har_rotation;
|
|
m_holdback = m_model_spec.holdout_size;
|
|
|
|
// Initialize core matrix calculations and allocate design variables
|
|
m_initialized = _init_model();
|
|
if(!m_initialized)
|
|
return false;
|
|
|
|
m_num_params = num_params();
|
|
|
|
//--- 4. Dynamic Polymorphic Volatility Process Object Allocation
|
|
if(CheckPointer(m_vp)==POINTER_DYNAMIC)
|
|
delete m_vp; // Prevent memory leaks by freeing old dynamic instances
|
|
m_vp = NULL;
|
|
|
|
switch(vol_dist_params.vol_model_type)
|
|
{
|
|
case VOL_CONST:
|
|
m_vp = new CConstantVariance(m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_ARCH:
|
|
m_vp = new CArchProcess(m_model_spec.garch_p,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_GARCH:
|
|
m_vp = new CGarchProcess(m_model_spec.garch_p,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_AVARCH:
|
|
m_vp = new CAvarchProcess(m_model_spec.garch_p,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_AVGARCH:
|
|
m_vp = new CAvgarchProcess(m_model_spec.garch_p,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_TARCH:
|
|
m_vp = new CTarchProcess(m_model_spec.garch_p,m_model_spec.garch_o,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_GJR_GARCH:
|
|
m_vp = new CGjrGarchProcess(m_model_spec.garch_p,m_model_spec.garch_o,m_model_spec.garch_q,m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_HARCH:
|
|
m_vp = new CHarchProcess(m_model_spec.harch_lags);
|
|
break;
|
|
case VOL_FIGARCH:
|
|
m_vp = new CFiGarchProcess(m_model_spec.garch_p,m_model_spec.garch_q,m_model_spec.vol_power,m_model_spec.figarch_truncation,m_model_spec.vol_rng_seed,m_model_spec.sample_start_idx,m_model_spec.sample_end_idx,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
case VOL_EGARCH:
|
|
m_vp = new CEgarchProcess(m_model_spec.garch_p,m_model_spec.garch_o,m_model_spec.garch_q,2.0,m_model_spec.vol_rng_seed,m_model_spec.sample_start_idx,m_model_spec.sample_end_idx,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
default:
|
|
m_vp = new CConstantVariance(m_model_spec.vol_rng_seed,m_model_spec.min_bootstrap_sims);
|
|
break;
|
|
}
|
|
|
|
//--- 5. Dynamic Innovation Error Distribution Object Allocation
|
|
if(CheckPointer(m_distribution)==POINTER_DYNAMIC)
|
|
delete m_distribution;
|
|
m_distribution = NULL;
|
|
|
|
switch(vol_dist_params.dist_type)
|
|
{
|
|
case DIST_NORMAL:
|
|
m_distribution = new CNormal();
|
|
break;
|
|
case DIST_STUDENT:
|
|
m_distribution = new CStudentsT();
|
|
break;
|
|
case DIST_SKEW_STUDENT:
|
|
m_distribution = new CSkewStudent();
|
|
break;
|
|
case DIST_GEN_ERROR:
|
|
m_distribution = new CGeneralizedError();
|
|
break;
|
|
default:
|
|
m_distribution = new CNormal();
|
|
break;
|
|
}
|
|
|
|
//--- 6. Combined Subsystem Initialization and Verification Audit
|
|
m_initialized = (
|
|
m_distribution!=NULL &&
|
|
m_vp!=NULL &&
|
|
m_vp.is_initialized() &&
|
|
m_distribution.initialize(m_model_spec.dist_init_params,
|
|
m_model_spec.dist_rng_seed)
|
|
);
|
|
|
|
//--- 7. Back-register Derived Subsystem Hyperparameters to Master Specifications
|
|
vector final_vol_specs = m_vp.garch_spec();
|
|
|
|
m_model_spec.garch_p = (ulong)final_vol_specs[0];
|
|
m_model_spec.garch_o = (ulong)final_vol_specs[1];
|
|
m_model_spec.garch_q = (ulong)final_vol_specs[2];
|
|
m_model_spec.vol_power = final_vol_specs[3];
|
|
|
|
return m_initialized;
|
|
}
|
|
public:
|
|
HARX(void):m_extra_simulation_params(0)
|
|
{
|
|
m_model_spec.mean_model_type = MEAN_HARX;
|
|
}
|
|
HARX(ArchParameters& model_specification)
|
|
{
|
|
model_specification.mean_model_type = MEAN_HARX;
|
|
initialize(model_specification);
|
|
}
|
|
|
|
~HARX(void)
|
|
{
|
|
deinitialize();
|
|
}
|
|
//--- Initialization method checks to ensure model correctly specified
|
|
bool initialize(ArchParameters& vol_dist_params)
|
|
{
|
|
return _initialize(vol_dist_params);
|
|
}
|
|
//--- Get exogenous variables if specified
|
|
matrix get_x(void) { return m_model_spec.exog_data; }
|
|
//--- Get regressors
|
|
matrix get_regressors(void) { return m_regressors; }
|
|
//--- Get observations
|
|
vector get_y(void) { return m_model_spec.observations; }
|
|
//--- Calculate residuals
|
|
virtual vector resids(vector ¶ms, vector &y, matrix ®ressors) override
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return EMPTY_VECTOR;
|
|
}
|
|
|
|
vector vec;
|
|
if(!y.Size())
|
|
{
|
|
vec = m_fit_regressors.MatMul(params);
|
|
return np::minus(m_fit_y,vec);
|
|
}
|
|
else
|
|
{
|
|
vec = regressors.MatMul(params);
|
|
return np::minus(y,vec);
|
|
}
|
|
}
|
|
//--- Get the number of model parameters
|
|
virtual ulong num_params(void)
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return 0;
|
|
}
|
|
return m_regressors.Cols();
|
|
}
|
|
//--- Manually set a previously initialized volatility process
|
|
bool set_volatility_process(CVolatilityProcess* &vp)
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return false;
|
|
}
|
|
if(CheckPointer(vp) == POINTER_INVALID)
|
|
{
|
|
Print(__FUNCTION__, " invalid pointer ");
|
|
return false;
|
|
}
|
|
|
|
if(!vp.is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " volatility process has not been initialized ");
|
|
return false;
|
|
}
|
|
|
|
if(CheckPointer(m_vp) == POINTER_DYNAMIC)
|
|
delete m_vp;
|
|
|
|
m_vp = NULL;
|
|
m_model_spec.vol_model_type = vp.volatilityprocess();
|
|
m_vp = vp;
|
|
m_delete_vp = false;
|
|
|
|
return true;
|
|
}
|
|
//--- Specify a previously initialized error distribution
|
|
bool set_distribution(CDistribution* &dist)
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return false;
|
|
}
|
|
if(CheckPointer(dist) == POINTER_INVALID)
|
|
{
|
|
Print(__FUNCTION__, " invalid pointer ");
|
|
return false;
|
|
}
|
|
|
|
if(!dist.is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " distribution has not been initialized ");
|
|
return false;
|
|
}
|
|
|
|
if(CheckPointer(m_distribution) == POINTER_DYNAMIC)
|
|
delete m_distribution;
|
|
|
|
m_distribution = NULL;
|
|
m_model_spec.dist_type = dist.distribution_type();
|
|
m_distribution = dist;
|
|
m_delete_dist = false;
|
|
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Simulates data from a linear regression, AR or HAR models |
|
|
//| Generates synthetic realization paths by combining a structural |
|
|
//| mean timeline with dynamic volatility pipelines and error shocks. |
|
|
//+------------------------------------------------------------------+
|
|
virtual matrix simulate(vector ¶ms, ulong nobs, ulong burn, vector &initial_vals, matrix &x, vector &initial_vals_vol)
|
|
{
|
|
//--- 1. Verification Gate: Check class instance state
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return EMPTY_MATRIX;
|
|
}
|
|
|
|
ulong kx = x.Cols();
|
|
ulong mc = 0;
|
|
|
|
//--- 2. Dimensional Gates: Verify exogenous matrix contains burn-in space allocations
|
|
if(kx && x.Rows()<nobs+burn)
|
|
{
|
|
Print(__FUNCTION__, " x must have nobs + burn rows");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
// Count parameters assigned exclusively to mean equation processes
|
|
if(m_model_spec.mean_lags.Size())
|
|
mc = ulong(m_constant)+m_lags.Cols()+kx+m_extra_simulation_params;
|
|
|
|
// Extract parameter allocation bounds across variance and error distributions
|
|
ulong vc = m_vp.numParams();
|
|
ulong dc = m_distribution.numParams();
|
|
ulong numparams = mc+vc+dc;
|
|
|
|
if(params.Size()!=numparams)
|
|
{
|
|
Print(__FUNCTION__, " params has the wrong number of elements.");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
//--- 3. Parameter Sub-vector Segmentation
|
|
vector allparams[];
|
|
if(!_parse_parameters(params,allparams))
|
|
return matrix::Zeros(0,0);
|
|
|
|
//--- 4. Structural Volatility and Innovation Error Modeling
|
|
// Initialize the specialized pseudorandom seed state engine based on target distribution rules
|
|
BootstrapRng rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
|
|
|
|
// Compute joint volatility pathways and return a multi-layered simulation snapshot matrix
|
|
matrix sim_data = m_vp.simulate(allparams[1],nobs+burn,rng,burn,initial_vals_vol.Size()?initial_vals_vol[0]:0.0);
|
|
|
|
// Extract unconditioned shock innovations and truncate the pre-sample burn-in allocation block
|
|
vector errors = sim_data.Row(0);
|
|
errors = np::sliceVector(errors,long(burn));
|
|
|
|
// Extract conditional variance, convert to volatility ($\sigma_t = \sqrt{\sigma_t^2}$), and drop burn-in entries
|
|
vector vol = sqrt(sim_data.Row(1));
|
|
vol = np::sliceVector(vol,long(burn));
|
|
|
|
//--- 5. Recursive Mean Synthesis Loop
|
|
vector y = _simulate_mean(allparams[0],x,errors,initial_vals,sim_data.Row(1));
|
|
y = np::sliceVector(y,long(burn)); // Trim transient burn-in rows from simulated dependent track
|
|
|
|
//--- 6. Construct Matrix Output Payload [Columns: Realized Series, Volatility, Error Shocks]
|
|
matrix out = matrix::Zeros(y.Size(),3);
|
|
if(!out.Col(y,0) || !out.Col(vol,1) || !out.Col(errors,2))
|
|
{
|
|
Print(__FUNCTION__, " ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
return out;
|
|
}
|
|
//--- Make forecasts from a model
|
|
ArchForecast forecast(ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
|
|
{
|
|
return forecast(EMPTY_VECTOR,EMPTY_MATRIX_ARRAY,horizon,start,method,simulations,seed);
|
|
}
|
|
//--- Forecast overload
|
|
ArchForecast forecast(matrix& x[],ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
|
|
{
|
|
return forecast(EMPTY_VECTOR,x,horizon,start,method,simulations,seed);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Multi-step conditional mean and volatility forecasting |
|
|
//| Generates multi-period analytic or simulation-based paths, |
|
|
//| handling rolling lag spaces and exogenous profile inputs ($X$). |
|
|
//+------------------------------------------------------------------+
|
|
virtual ArchForecast forecast(vector& params, matrix& x[], ulong horizon = 1, long start = -1, ENUM_FORECAST_METHOD method=FORECAST_ANALYTIC, ulong simulations=1000, uint seed=0)
|
|
{
|
|
ArchForecast out;
|
|
|
|
//--- 1. Verification Gate: Check object instance readiness
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return out;
|
|
}
|
|
|
|
// Automatically handle sample configuration boundaries if the active set is blank
|
|
if(!m_fit_y.Size())
|
|
{
|
|
long from = 0;
|
|
long to = long(m_model_spec.observations.Size());
|
|
if(!_adjust_sample(from,to))
|
|
return out;
|
|
}
|
|
|
|
//--- 2. Sample Alignment and Edge Case Boundary Assertions
|
|
long earliest = long(m_fit_indices[0]);
|
|
long default_start = long(m_fit_indices[1]);
|
|
default_start = MathMax(0,default_start-1);
|
|
long start_index = (start>-1 && start<=default_start)?start:default_start;
|
|
|
|
// Enforce the data boundary floor: Reject attempts to look behind the max lag memory window
|
|
if(start_index<(earliest-1))
|
|
{
|
|
Print(__FUNCTION__," Due to backcasting and/or data availability start cannot be less "\
|
|
"than the index of the largest value in the right-hand-side "\
|
|
"variables used to fit the first observation.");
|
|
return out;
|
|
}
|
|
|
|
//--- 3. Parameter Deconstruction and Signal Parsing
|
|
vector allparams[];
|
|
if(!CArchModel::_parse_parameters(params.Size()?params:m_params,allparams))
|
|
{
|
|
Print(__FUNCTION__, " failed to parse the parameters ");
|
|
return out;
|
|
}
|
|
|
|
// Extract structural residuals and calculate target variance backcast roots
|
|
vector res = resids(allparams[0],vector::Zeros(0),matrix::Zeros(0,0));
|
|
vector bc = m_vp.backCast(res);
|
|
vector ys = np::sliceVector(m_model_spec.observations,long(earliest));
|
|
matrix xs = np::sliceMatrixRows(m_regressors,long(earliest));
|
|
vector full_res = resids(allparams[0],ys,xs);
|
|
matrix vb = m_vp.varianceBounds(full_res,2.);
|
|
|
|
// Initialize random distribution engines
|
|
BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
|
|
long vstart = MathMax(0,start_index - earliest);
|
|
|
|
//--- 4. Generate Core Conditional Variance Forecast
|
|
VarianceForecast vfcast = m_vp.forecast(allparams[1],full_res,bc,vb,Rng,seed,vstart,horizon,method,simulations);
|
|
if(!vfcast.forecasts.Cols())
|
|
return out;
|
|
|
|
matrix var_fcasts = vfcast.forecasts;
|
|
if(start_index<earliest)
|
|
var_fcasts = _forecast_pad(int(earliest-start_index),var_fcasts);
|
|
|
|
// Convert HAR multi-day weights into clean point-to-point standard daily AR lag arrays
|
|
vector arp = _har_to_ar(allparams[0]);
|
|
ulong nexog = m_model_spec.exog_data.Cols();
|
|
vector exog_p = (nexog)?np::sliceVector(allparams[0],-1*long(nexog)):vector::Zeros(0);
|
|
double constant = m_constant?arp[0]:0.0;
|
|
vector dynp = (arp.Size()>ulong(m_constant))?np::sliceVector(arp,long(m_constant)):vector::Zeros(0);
|
|
|
|
// Align upcoming exogenous feature tracks ($X_{T+h}$)
|
|
matrix expected_x[];
|
|
if(!_reformat_forecast_x(ulong(start_index),horizon,x,expected_x))
|
|
{
|
|
Print(__FUNCTION__, " reformating error ");
|
|
return out;
|
|
}
|
|
|
|
//--- 5. Generate Conditional Mean Forecast Pathways
|
|
matrix mean_fcast = _ar_forecast(m_model_spec.observations,horizon,start_index,constant,dynp,expected_x,exog_p);
|
|
if(!mean_fcast.Cols())
|
|
return out;
|
|
|
|
//--- 6. Construct Long-Run Cumulative Volatility Profiles
|
|
// Maps localized shocks to future windows using an MA($\infty$) impulse map ($\psi_h$)
|
|
vector impulse = _ar_to_impulse(horizon,dynp);
|
|
matrix longrun_var_fcasts = var_fcasts;
|
|
matrix tempm;
|
|
vector tempv, lrf;
|
|
|
|
// Accumulate squared coefficients along the impulse timeline: $\sigma^2_{T+h} = \sum \psi_j^2 \cdot \sigma^2_{T+h-j}$
|
|
for(ulong i = 0; i<horizon; ++i)
|
|
{
|
|
tempm = np::sliceMatrixCols(var_fcasts,0,long(i+1));
|
|
tempv = np::sliceVector(impulse,long(i),END_REVERSE,STEP_REVERSE);
|
|
tempv = pow(tempv,2.);
|
|
lrf = tempm.MatMul(tempv);
|
|
longrun_var_fcasts.Col(lrf,i);
|
|
}
|
|
|
|
//--- 7. Path-Simulation Processing Layer [Bootstrap / Monte Carlo Simulation Paths]
|
|
matrix variance_paths[], mean_paths[];
|
|
matrix shocks[], long_run_variance_paths[];
|
|
|
|
if(method == FORECAST_BOOTSTRAP || method == FORECAST_SIMULATION)
|
|
{
|
|
// Copy simulated variance paths and raw standardized shock structures
|
|
ArrayResize(variance_paths,vfcast.forecastpaths.Size());
|
|
for(uint i = 0; i<variance_paths.Size(); ++i)
|
|
variance_paths[i] = vfcast.forecastpaths[i];
|
|
|
|
ArrayResize(shocks,vfcast.shocks.Size());
|
|
for(uint i = 0; i<shocks.Size(); ++i)
|
|
shocks[i] = vfcast.shocks[i];
|
|
|
|
if(start_index<earliest)
|
|
{
|
|
_forecast_pad(int(earliest-start_index),variance_paths,variance_paths);
|
|
_forecast_pad(int(earliest-start_index),shocks,shocks);
|
|
}
|
|
|
|
ArrayResize(long_run_variance_paths,variance_paths.Size());
|
|
for(uint i = 0; i<variance_paths.Size(); ++i)
|
|
long_run_variance_paths[i] = vfcast.forecastpaths[i];
|
|
|
|
matrix imps;
|
|
vector _impulses;
|
|
|
|
// Evaluate cumulative long-run variance profiles for simulated trials
|
|
for(ulong i = 0; i<horizon; ++i)
|
|
{
|
|
_impulses = np::sliceVector(impulse,long(i),END_REVERSE,STEP_REVERSE);
|
|
imps = matrix::Zeros(_impulses.Size(),1);
|
|
_impulses = pow(_impulses,2.);
|
|
imps.Col(_impulses,0);
|
|
for(uint k = 0; k<variance_paths.Size(); ++k)
|
|
{
|
|
tempm = np::sliceMatrixCols(variance_paths[k],BEGIN,long(i+1));
|
|
tempm = tempm.MatMul(imps);
|
|
long_run_variance_paths[k].Col(tempm.Col(0),i);
|
|
}
|
|
}
|
|
|
|
//--- 8. Recursive Mean Synthesis across Simulation Iterations
|
|
ulong t = m_model_spec.observations.Size();
|
|
ArrayResize(mean_paths,shocks.Size());
|
|
for(uint i =0; i<mean_paths.Size(); ++i)
|
|
mean_paths[i] = matrix::Zeros(shocks[0].Rows(),m_maxlags+horizon);
|
|
|
|
vector dynp_rev = dynp;
|
|
np::reverseVector(dynp_rev); // Match lag coefficients to reversed chronological rows
|
|
|
|
for(ulong i = start_index; i<t; ++i)
|
|
{
|
|
ulong ploc = i - start_index;
|
|
tempv = np::sliceVector(m_model_spec.observations,long(i-m_maxlags+1),long(i+1));
|
|
tempm = np::repeat_vector_as_rows_cols(tempv,mean_paths[ploc].Rows());
|
|
np::matrixCopy(mean_paths[ploc],tempm,BEGIN,END,STEP,BEGIN,long(m_maxlags));
|
|
|
|
// Loop step-by-step through the horizon window, updating mean pathways dynamically
|
|
for(ulong j = 0; j<horizon; ++j)
|
|
{
|
|
tempm = np::sliceMatrixCols(mean_paths[ploc],long(j),long(m_maxlags+j));
|
|
tempv = constant + tempm.MatMul(dynp_rev) + shocks[ploc].Col(j);
|
|
mean_paths[ploc].Col(tempv,m_maxlags+j);
|
|
|
|
// Inject multi-step exogenous model parameters if assigned
|
|
if(expected_x.Size())
|
|
{
|
|
tempv = vector::Zeros(mean_paths[ploc].Rows());
|
|
for(uint z = 0; z<expected_x.Size(); ++z)
|
|
tempv[z] = expected_x[z][ploc,j];
|
|
double dotproduct = tempv.Dot(exog_p);
|
|
for(ulong k = 0; k<mean_paths[ploc].Rows(); ++k)
|
|
mean_paths[ploc][k,m_maxlags+j] += dotproduct;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Strip off seed historical lag nodes to leave only pure out-of-sample forecast horizons
|
|
for(uint i = 0; i<mean_paths.Size(); ++i)
|
|
mean_paths[i] = np::sliceMatrixCols(mean_paths[i],long(m_maxlags));
|
|
}
|
|
|
|
//--- 9. Package final data arrays into the output container
|
|
return ArchForecast(start_index,mean_fcast,longrun_var_fcasts,var_fcasts,mean_paths,shocks,long_run_variance_paths,variance_paths);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Construct a fixed-parameter model execution result |
|
|
//| Bypasses estimation optimization to run evaluation diagnostics |
|
|
//| directly using a pre-determined user-defined parameter array. |
|
|
//+------------------------------------------------------------------+
|
|
ArchModelFixedResult fix(vector& params, long first_obs = 0, long last_obs = -1)
|
|
{
|
|
ArchModelFixedResult out;
|
|
|
|
//--- 1. Verification Gate: Check class instance state readiness
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return out;
|
|
}
|
|
|
|
// Synchronize data sequence indexes with given sample limits
|
|
if(!_adjust_sample(first_obs,last_obs))
|
|
return out;
|
|
|
|
//--- 2. Pre-Sample Backcast Conditioning and Target Variance Bounds
|
|
vector sv = starting_values();
|
|
vector res = resids(sv,EMPTY_VECTOR,EMPTY_MATRIX);
|
|
vector sigma2 = vector::Zeros(res.Size());
|
|
vector back_cast = m_vp.backCast(res);
|
|
m_backcast = back_cast; // Cache the historical conditioning backcast value
|
|
matrix vb = m_vp.varianceBounds(res);
|
|
|
|
// Compute the overall objective negative log-likelihood score vector
|
|
vector logl = _loglikelihood(params,sigma2,back_cast,vb)*(-1.0);
|
|
|
|
//--- 3. Parameter Deconstruction and Signal Slicing
|
|
vector parameters[];
|
|
if(!_parse_parameters(params,parameters))
|
|
return out;
|
|
|
|
//--- 4. Direct Volatility Evolution Trajectory Generation
|
|
vector vol = vector::Zeros(res.Size());
|
|
m_vp.computeVariance(parameters[1],res,vol,back_cast,vb);
|
|
|
|
// Convert raw variance states into volatility standard deviations ($\sigma_t = \sqrt{\sigma_t^2}$)
|
|
vol = sqrt(vol);
|
|
|
|
//--- 5. Synchronize Timeline Data Windows and Padding Assignments
|
|
first_obs = long(m_fit_indices[0]);
|
|
last_obs = long(m_fit_indices[1]);
|
|
|
|
// Copy core model residuals into a NaN-padded structural timeline vector
|
|
vector res_final = vector::Zeros(res.Size());
|
|
res_final.Fill(AL_NaN);
|
|
np::vectorCopy(res_final,res,long(first_obs),long(last_obs));
|
|
|
|
// Copy estimated conditional volatilities into a matching NaN-padded container
|
|
vector vol_final = vector::Zeros(vol.Size());
|
|
vol_final.Fill(AL_NaN);
|
|
np::vectorCopy(vol_final,vol,long(first_obs),long(last_obs));
|
|
|
|
//--- 6. Package final variables into the output object payload
|
|
out.params = params;
|
|
out.resid = res_final;
|
|
out.conditional_volatility = vol_final;
|
|
out.loglikelihood = logl[0];
|
|
|
|
return out;
|
|
}
|
|
#ifdef __SLSQP__
|
|
//--- Fit a model to data using SLSQP optimizer
|
|
ArchModelResult fit(double tol = 1e-6,uint maxits = 0,ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1)
|
|
{
|
|
vector starting_values,back_cast;
|
|
return fit(starting_values,back_cast,tol,maxits,cov_type,first,last);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Fit model using numerical optimization |
|
|
//| Orchestrates the optimization lifecycle by concatenating param |
|
|
//| segments, stacking dynamic boundaries, and passing everything |
|
|
//| into the SLSQP solver engine to extract optimal coefficients. |
|
|
//+------------------------------------------------------------------+
|
|
ArchModelResult fit(vector& startingvalues, vector& backcast, double tol = 1e-6, uint maxits = 0, ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1)
|
|
{
|
|
ArchModelResult out;
|
|
|
|
//--- 1. Verification Gates: Validate initialization and data state
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return out;
|
|
}
|
|
if(!m_model_spec.observations.Size())
|
|
{
|
|
Print(__FUNCTION__," Cannot estimate model without data.");
|
|
return out;
|
|
}
|
|
|
|
//--- 2. Parse Parameter Counts and Set Sample Boundaries
|
|
vector offsets(3);
|
|
offsets[0] = (double)m_num_params; // Mean equation parameters
|
|
offsets[1] = (double)m_vp.numParams(); // Volatility process parameters
|
|
offsets[2] = (double)m_distribution.numParams(); // Error distribution shape parameters
|
|
|
|
double total_params = offsets.Sum();
|
|
|
|
if(last<0 || last<first)
|
|
last = long(m_model_spec.observations.Size());
|
|
|
|
if(!_adjust_sample(first,last))
|
|
return out;
|
|
|
|
// Check for closed-form shortcuts (e.g., constant volatility, normal distribution)
|
|
bool hascloseform = (offsets[2] == 0.0 && m_vp.closedForm() && m_vp.volatilityprocess()==VOL_CONST);
|
|
|
|
//--- 3. Data Scaling Check and Adjustment
|
|
vector sv = starting_values();
|
|
vector res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
|
|
_checkscale(res);
|
|
|
|
if(m_scale!=1.)
|
|
{
|
|
m_model_spec.observations = m_model_spec.observations*m_scale;
|
|
_scalechanged();
|
|
_adjust_sample(first,last);
|
|
sv = starting_values();
|
|
res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
|
|
}
|
|
|
|
//--- 4. Configure Pre-sample Backcast Values
|
|
vector bcast;
|
|
if(!backcast.Size())
|
|
bcast = m_vp.backCast(res);
|
|
else
|
|
bcast = m_vp.backCastTransform(backcast);
|
|
|
|
// Direct routing shortcuts for exceptional model layouts
|
|
if(hascloseform)
|
|
return _fit_no_arch_normal_errors(cov_type);
|
|
if(total_params == 0.0)
|
|
return _fit_parameterless_model(cov_type,bcast);
|
|
// Initialize conditional variance pathways and bounds
|
|
vector sigma2 = vector::Zeros(res.Size());
|
|
m_backcast = bcast;
|
|
vector sv_v = m_vp.startingValues(res);
|
|
m_var_bounds = m_vp.varianceBounds(res);
|
|
matrix vb = m_var_bounds;
|
|
m_vp.computeVariance(sv_v,res,sigma2,bcast,vb);
|
|
vector std_res = res/sqrt(sigma2);
|
|
|
|
//--- 5. Collect and Vertically Stack Inequality Constraint Matrices
|
|
Constraints cst[3];
|
|
cst[0] = constraints();
|
|
cst[1] = m_vp.constraints();
|
|
cst[2] = m_distribution.constraints();
|
|
|
|
vector num_cons(3);
|
|
for(uint i = 0; i<cst.Size(); ++i)
|
|
num_cons[i] = (double)cst[i]._one.Rows();
|
|
|
|
matrix a = matrix::Zeros(ulong(num_cons.Sum()),ulong(total_params));
|
|
vector b = vector::Zeros(ulong(num_cons.Sum()));
|
|
long ren, cen, rst, ct;
|
|
vector nc, of;
|
|
|
|
for(uint i = 0; i<3; ++i)
|
|
{
|
|
nc = np::sliceVector(num_cons,0,long(i+1));
|
|
of = np::sliceVector(offsets,0,long(i+1));
|
|
ren = long(nc.Sum());
|
|
cen = long(of.Sum());
|
|
rst = ren - long(num_cons[i]);
|
|
ct = cen - long(offsets[i]);
|
|
if(ren-rst>0)
|
|
{
|
|
np::matrixCopy(a,cst[i]._one,rst,ren,STEP,ct,cen);
|
|
np::vectorCopy(b,cst[i]._two,rst,ren);
|
|
}
|
|
}
|
|
|
|
//--- 6. Gather and Horizontal Stitch Parameter Min/Max Bounds
|
|
matrix bound[3];
|
|
bound[0] = bounds();
|
|
bound[1] = m_vp.bounds(res);
|
|
bound[2] = m_distribution.bounds(res);
|
|
matrix all_bounds(2,bound[0].Cols()+bound[1].Cols()+bound[2].Cols());
|
|
|
|
long cols = 0;
|
|
for(uint i = 0; i<bound.Size(); ++i)
|
|
{
|
|
if(bound[i].Cols())
|
|
{
|
|
cols+=(i)?long(bound[i-1].Cols()):0;
|
|
np::matrixCopyCols(all_bounds,bound[i],cols,long(cols+bound[i].Cols()));
|
|
}
|
|
}
|
|
|
|
//--- 7. Validate or Automatically Seed Initial Starting Parameters
|
|
if(startingvalues.Size())
|
|
{
|
|
bool valid = (startingvalues.Size()==ulong(total_params));
|
|
if(a.Rows()&&valid)
|
|
{
|
|
vector temp = b - a.MatMul(startingvalues);
|
|
valid = (valid && temp.Max()<0.0);
|
|
}
|
|
for(uint i = 0; i<startingvalues.Size(); ++i)
|
|
valid = (valid && all_bounds[0,i]<=startingvalues[i] && startingvalues[i]<=all_bounds[1,i]);
|
|
if(!valid)
|
|
startingvalues = EMPTY_VECTOR;
|
|
}
|
|
|
|
if(!startingvalues.Size())
|
|
{
|
|
sv = starting_values();
|
|
ulong ss = sv.Size();
|
|
vector temp = m_distribution.startingValues(std_res);
|
|
sv.Resize(ss+sv_v.Size()+temp.Size());
|
|
if(sv_v.Size())
|
|
np::vectorCopy(sv,sv_v,long(ss),long(ss+sv_v.Size()));
|
|
if(temp.Size())
|
|
np::vectorCopy(sv,temp,long(ss+sv_v.Size()));
|
|
}
|
|
else
|
|
sv = startingvalues;
|
|
|
|
//--- 8. Prepare and Initialize Objective Functions and Constraints
|
|
CObject obj;
|
|
CObjectiveF obfunc;
|
|
CIneqConstraints obfunc_constraints;
|
|
|
|
obfunc.setFuncParams(sigma2,bcast,vb,this);
|
|
obfunc_constraints.setConstraints(a,b);
|
|
|
|
all_bounds = all_bounds.Transpose();
|
|
|
|
obfunc.setBounds(all_bounds);
|
|
obfunc_constraints.setBounds(all_bounds);
|
|
|
|
if(!obfunc.initialize(sv) || !obfunc_constraints.initialize(sv))
|
|
{
|
|
m_converged = false;
|
|
return out;
|
|
}
|
|
|
|
//--- 9. Execute Numerical Non-Linear Optimization Loop (SLSQP Solver)
|
|
CSlsqp slsqp_minim;
|
|
slsqp_minim.SetInequalityConstraints(obfunc_constraints);
|
|
|
|
if(maxits)
|
|
slsqp_minim.SetMaxEval(int(maxits));
|
|
if(fabs(tol))
|
|
{
|
|
slsqp_minim.SetAcc(tol);
|
|
slsqp_minim.SetXtolRel(tol);
|
|
}
|
|
|
|
OptimizeResult minim_result = slsqp_minim.Minimize(obfunc);
|
|
|
|
//--- 10. Process Optimization Solutions and Update Internal States
|
|
if(minim_result.return_code<0)
|
|
{
|
|
Print(__FUNCTION__, " termination reason ", minim_result.return_code);
|
|
m_converged = false;
|
|
return out;
|
|
}
|
|
else
|
|
{
|
|
m_params = minim_result.solution;
|
|
m_converged = true;
|
|
//Print(__FUNCTION__, " termination reason ", minim_result.return_code);
|
|
}
|
|
|
|
//--- 11. Extract Convergence Data Metrics and Pack Final Results Struct
|
|
vector allparams[];
|
|
vector llh = _loglikelihood(m_params,sigma2,bcast,vb);
|
|
|
|
_parse_parameters(m_params,allparams);
|
|
res = resids(allparams[0],EMPTY_VECTOR,EMPTY_MATRIX);
|
|
|
|
vector vol = vector::Zeros(res.Size());
|
|
m_vp.computeVariance(allparams[1],res,vol,bcast,vb);
|
|
vol = sqrt(vol);
|
|
|
|
double r2 = _r2(allparams[0]);
|
|
|
|
long f_obs = (long)m_fit_indices[0];
|
|
long l_obs = (long)m_fit_indices[1];
|
|
|
|
vector res_final, vol_final;
|
|
res_final = vol_final = m_model_spec.observations;
|
|
res_final.Fill(AL_NaN);
|
|
vol_final = res_final;
|
|
|
|
np::vectorCopy(res_final,res,f_obs,l_obs);
|
|
np::vectorCopy(vol_final,vol,f_obs,l_obs);
|
|
|
|
matrix pcov = compute_param_cov(m_params,bcast,(cov_type==COVAR_ROBUST));
|
|
|
|
return ArchModelResult(m_params,pcov,r2,res_final,vol_final,cov_type,-1.0*llh[0],f_obs,l_obs);
|
|
}
|
|
#else
|
|
//--- Fit model to data using Alglib optimizer
|
|
ArchModelResult fit(double scaling = 1.0, uint maxits = 0,ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1, double tol = 1e-9, bool guardsmoothness=false, double gradient_test_step = 0.0)
|
|
{
|
|
return fit(EMPTY_VECTOR,EMPTY_VECTOR,scaling,maxits,cov_type,first,last,tol,guardsmoothness,gradient_test_step);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Fit model using ALGLIB Nonlinear Constrained Optimization (AUL) |
|
|
//| Orchestrates estimation via Augmented Lagrangian Multipliers, |
|
|
//| integrating smoothness safeguards, gradient test steps, and |
|
|
//| custom interactive Jacobian evaluation loops. |
|
|
//+------------------------------------------------------------------+
|
|
ArchModelResult fit(vector& startingvalues, vector& backcast, double scaling = 1.0, uint maxits = 0, ENUM_COVAR_TYPE cov_type = COVAR_ROBUST, long first = 0, long last = -1, double tol = 1e-9, bool guardsmoothness=false, double gradient_test_step = 0.0)
|
|
{
|
|
ArchModelResult out;
|
|
|
|
//--- 1. Verification Gates: Validate initialization and data state
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return out;
|
|
}
|
|
if(!m_model_spec.observations.Size())
|
|
{
|
|
Print(__FUNCTION__," Cannot estimate model without data.");
|
|
return out;
|
|
}
|
|
|
|
//--- 2. Parse Parameter Counts and Set Sample Boundaries
|
|
vector offsets(3);
|
|
offsets[0] = (double)m_num_params; // Mean model parameter count
|
|
offsets[1] = (double)m_vp.numParams(); // Volatility process parameter count
|
|
offsets[2] = (double)m_distribution.numParams(); // Error distribution shape parameter count
|
|
|
|
double total_params = offsets.Sum();
|
|
|
|
if(last<0 || last<first)
|
|
last = long(m_model_spec.observations.Size());
|
|
|
|
if(!_adjust_sample(first,last))
|
|
return out;
|
|
|
|
// Check for closed-form shortcuts (e.g., constant volatility, normal distribution)
|
|
bool hascloseform = (offsets[2] == 0.0 && m_vp.closedForm() && m_vp.volatilityprocess()==VOL_CONST);
|
|
|
|
//--- 3. Data Scaling Check and Adjustment
|
|
vector sv = starting_values();
|
|
vector res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
|
|
_checkscale(res);
|
|
|
|
if(m_scale!=1.)
|
|
{
|
|
m_model_spec.observations = m_model_spec.observations*m_scale;
|
|
_scalechanged();
|
|
_adjust_sample(first,last);
|
|
sv = starting_values();
|
|
res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
|
|
}
|
|
|
|
//--- 4. Configure Pre-sample Backcast Values
|
|
vector bcast;
|
|
if(!backcast.Size())
|
|
bcast = m_vp.backCast(res);
|
|
else
|
|
bcast = m_vp.backCastTransform(backcast);
|
|
|
|
// Direct routing shortcuts for exceptional model layouts
|
|
if(hascloseform)
|
|
return _fit_no_arch_normal_errors(cov_type);
|
|
if(total_params == 0.0)
|
|
return _fit_parameterless_model(cov_type,bcast);
|
|
|
|
// Initialize conditional variance pathways and bounds
|
|
vector sigma2 = vector::Zeros(res.Size());
|
|
m_backcast = bcast;
|
|
vector sv_v = m_vp.startingValues(res);
|
|
m_var_bounds = m_vp.varianceBounds(res);
|
|
matrix vb = m_var_bounds;
|
|
m_vp.computeVariance(sv_v,res,sigma2,bcast,vb);
|
|
vector std_res = res/sqrt(sigma2);
|
|
|
|
//--- 5. Collect and Format Inequality Constraint Systems ($A \cdot x \ge b$)
|
|
Constraints cst[3];
|
|
cst[0] = constraints();
|
|
cst[1] = m_vp.constraints();
|
|
cst[2] = m_distribution.constraints();
|
|
|
|
vector num_cons(3);
|
|
for(uint i = 0; i<cst.Size(); ++i)
|
|
num_cons[i] = (double)cst[i]._one.Rows();
|
|
|
|
// Allocate space with an extra column to hold linear condition offsets natively
|
|
matrix a = matrix::Zeros(ulong(num_cons.Sum()),ulong(total_params)+1);
|
|
vector b = vector::Zeros(ulong(num_cons.Sum()));
|
|
long ren, cen, rst, ct;
|
|
vector nc, of;
|
|
|
|
for(uint i = 0; i<3; ++i)
|
|
{
|
|
nc = np::sliceVector(num_cons,0,long(i+1));
|
|
of = np::sliceVector(offsets,0,long(i+1));
|
|
ren = long(nc.Sum());
|
|
cen = long(of.Sum());
|
|
rst = ren - long(num_cons[i]);
|
|
ct = cen - long(offsets[i]);
|
|
if(ren-rst>0)
|
|
{
|
|
np::matrixCopy(a,cst[i]._one,rst,ren,STEP,ct,cen);
|
|
np::vectorCopy(b,cst[i]._two,rst,ren);
|
|
}
|
|
}
|
|
|
|
// Convert constraints into standard ALGLIB linear format: [A | -b] . x >= 0
|
|
a.Col(b,ulong(total_params));
|
|
b = vector::Ones(a.Rows());
|
|
|
|
//--- 6. Gather Boundary Matrices [Min / Max Box Constraints]
|
|
matrix bound[3];
|
|
bound[0] = bounds();
|
|
bound[1] = m_vp.bounds(res);
|
|
bound[2] = m_distribution.bounds(res);
|
|
matrix all_bounds(2,bound[0].Cols()+bound[1].Cols()+bound[2].Cols());
|
|
|
|
long cols = 0;
|
|
for(uint i = 0; i<bound.Size(); ++i)
|
|
{
|
|
if(bound[i].Cols())
|
|
{
|
|
cols+=(i)?long(bound[i-1].Cols()):0;
|
|
np::matrixCopyCols(all_bounds,bound[i],cols,long(cols+bound[i].Cols()));
|
|
}
|
|
}
|
|
|
|
//--- 7. Validate Provided Starting Vector Against Linear and Box Constraints
|
|
if(startingvalues.Size())
|
|
{
|
|
bool valid = (startingvalues.Size()==ulong(total_params));
|
|
matrix aa = np::sliceMatrixCols(a,0,long(a.Cols()-1));
|
|
vector bb = a.Col(a.Cols()-1);
|
|
if(aa.Rows()&&valid)
|
|
{
|
|
vector temp = aa.MatMul(startingvalues) - bb;
|
|
valid = (valid && temp.Min()>0.0);
|
|
}
|
|
for(uint i = 0; i<startingvalues.Size(); ++i)
|
|
valid = (valid && all_bounds[0,i]<=startingvalues[i] && startingvalues[i]<=all_bounds[1,i]);
|
|
if(!valid)
|
|
startingvalues = EMPTY_VECTOR;
|
|
}
|
|
|
|
// Compute structural initial guesses if the vector remains empty
|
|
if(!startingvalues.Size())
|
|
{
|
|
sv = starting_values();
|
|
ulong ss = sv.Size();
|
|
vector temp = m_distribution.startingValues(std_res);
|
|
sv.Resize(ss+sv_v.Size()+temp.Size());
|
|
if(sv_v.Size())
|
|
np::vectorCopy(sv,sv_v,long(ss),long(ss+sv_v.Size()));
|
|
if(temp.Size())
|
|
np::vectorCopy(sv,temp,long(ss+sv_v.Size()));
|
|
}
|
|
else
|
|
sv = startingvalues;
|
|
|
|
//--- 8. Initialize ALGLIB Nonlinear Constrained State Machine (MinNLC)
|
|
CMinNLCState optim_state;
|
|
CRowDouble _initial_values = sv;
|
|
|
|
CRowDouble _ones = vector::Ones(sv.Size());
|
|
if(scaling>0.0)
|
|
_ones*=scaling;
|
|
|
|
CRowDouble bl = all_bounds.Row(0);
|
|
CRowDouble bu = all_bounds.Row(1);
|
|
CMatrixDouble _constraints = a;
|
|
|
|
CMinNLC::MinNLCCreate(_initial_values.Size(),_initial_values,optim_state);
|
|
CMinNLC::MinNLCSetCond(optim_state,tol,int(maxits));
|
|
CMinNLC::MinNLCSetScale(optim_state,_ones);
|
|
CMinNLC::MinNLCSetSTPMax(optim_state,0.0);
|
|
CMinNLC::MinNLCSetAlgoAUL(optim_state,10,0); // Select Augmented Lagrangian Multipliers
|
|
CMinNLC::MinNLCSetBC(optim_state,bl,bu); // Bind simple box boundaries
|
|
|
|
CRowInt intones;
|
|
intones.Resize(_constraints.Rows());
|
|
intones.Fill(1); // Identify all rows as inequality expressions (>= 0)
|
|
CMinNLC::MinNLCSetLC(optim_state,_constraints,intones,intones.Size());
|
|
|
|
// Configure algorithmic safety features
|
|
if(guardsmoothness)
|
|
CMinNLC::MinNLCOptGuardSmoothness(optim_state,guardsmoothness);
|
|
if(gradient_test_step>0.0)
|
|
CMinNLC::MinNLCOptGuardGradient(optim_state,gradient_test_step);
|
|
|
|
CMinNLCReport rep;
|
|
CRowDouble solution;
|
|
|
|
// Set up objective function structural parameters
|
|
CObjectiveJacFVec fvec;
|
|
fvec.SetParams(sigma2,bcast,vb,this);
|
|
|
|
CNDimensional_Rep f_rep;
|
|
CObject object;
|
|
|
|
if(!CAp::Assert(GetPointer(fvec)!=NULL,"ALGLIB: error in 'minnlcoptimize()' (fvec is null)"))
|
|
{
|
|
m_converged = false;
|
|
return out;
|
|
}
|
|
|
|
//--- 9. Interactive Optimization Routine (Jacobian / Step-Evaluation Cycle)
|
|
while(CMinNLC::MinNLCIteration(optim_state))
|
|
{
|
|
// Handle request to re-evaluate function vector and full Jacobian matrix
|
|
if(optim_state.m_needfij)
|
|
{
|
|
fvec.Jac(optim_state.m_x,optim_state.m_fi,optim_state.m_j,object);
|
|
continue;
|
|
}
|
|
// Handle request to report update steps
|
|
if(optim_state.m_xupdated)
|
|
{
|
|
if(GetPointer(f_rep)!=NULL)
|
|
f_rep.Rep(optim_state.m_x,optim_state.m_f,object);
|
|
continue;
|
|
}
|
|
CAp::Assert(false,"ALGLIB: error in 'minnlcoptimize' (some derivatives were not provided?)");
|
|
break;
|
|
}
|
|
|
|
//--- 10. Extract Optimization Results and Assess Convergence Status
|
|
CMinNLC::MinNLCResults(optim_state,solution,rep);
|
|
if(rep.m_terminationtype<0)
|
|
{
|
|
Print(__FUNCTION__," failed convergence: ", rep.m_terminationtype);
|
|
m_converged = false;
|
|
return out;
|
|
}
|
|
else
|
|
m_converged = true;
|
|
|
|
//--- 11. Extract Convergence Data Metrics and Pack Final Results Struct
|
|
vector allparams[];
|
|
m_params = solution.ToVector();
|
|
|
|
// Calculate final target log-likelihood
|
|
vector llh = _loglikelihood(m_params,sigma2,bcast,vb);
|
|
|
|
_parse_parameters(m_params,allparams);
|
|
res = resids(allparams[0],EMPTY_VECTOR,EMPTY_MATRIX);
|
|
|
|
vector vol = vector::Zeros(res.Size());
|
|
m_vp.computeVariance(allparams[1],res,vol,bcast,vb);
|
|
vol = sqrt(vol);
|
|
|
|
double r2 = _r2(allparams[0]);
|
|
|
|
long f_obs = (long)m_fit_indices[0];
|
|
long l_obs = (long)m_fit_indices[1];
|
|
|
|
vector res_final, vol_final;
|
|
res_final = vol_final = m_model_spec.observations;
|
|
res_final.Fill(AL_NaN);
|
|
vol_final = res_final;
|
|
|
|
np::vectorCopy(res_final,res,f_obs,l_obs);
|
|
np::vectorCopy(vol_final,vol,f_obs,l_obs);
|
|
|
|
matrix pcov = compute_param_cov(m_params,bcast,(cov_type==COVAR_ROBUST));
|
|
|
|
return ArchModelResult(m_params,pcov,r2,res_final,vol_final,cov_type,-1.0*llh[0],f_obs,l_obs);
|
|
}
|
|
#endif
|
|
//+------------------------------------------------------------------+
|
|
//| Accessible loglikelihood evaluation routine |
|
|
//| Replicates the exact environment setup used during estimation |
|
|
//| to compute a clean out-of-sample or spot likelihood metric. |
|
|
//+------------------------------------------------------------------+
|
|
bool eval_loglikelihood(vector& params, double &out_loglikelihood)
|
|
{
|
|
long first = 0;
|
|
long last = -1;
|
|
|
|
vector startingvalues = EMPTY_VECTOR;
|
|
vector backcast = EMPTY_VECTOR;
|
|
|
|
//--- 1. Verification Gates: Validate initialization and data state
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return false;
|
|
}
|
|
if(!m_model_spec.observations.Size())
|
|
{
|
|
Print(__FUNCTION__," Cannot estimate model without data.");
|
|
return false;
|
|
}
|
|
|
|
//--- 2. Parse Parameter Counts and Set Sample Boundaries
|
|
vector offsets(3);
|
|
offsets[0] = (double)m_num_params; // Mean model parameter count
|
|
offsets[1] = (double)m_vp.numParams(); // Volatility process parameter count
|
|
offsets[2] = (double)m_distribution.numParams(); // Error distribution shape parameter count
|
|
|
|
double total_params = offsets.Sum();
|
|
|
|
if(last<0 || last<first)
|
|
last = long(m_model_spec.observations.Size());
|
|
|
|
if(!_adjust_sample(first,last))
|
|
return false;
|
|
|
|
// Check for closed-form shortcuts (e.g., constant volatility, normal distribution)
|
|
bool hascloseform = (offsets[2] == 0.0 && m_vp.closedForm() && m_vp.volatilityprocess()==VOL_CONST);
|
|
|
|
//--- 3. Data Scaling Check and Adjustment
|
|
vector sv = starting_values();
|
|
vector res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
|
|
_checkscale(res);
|
|
|
|
if(m_scale!=1.)
|
|
{
|
|
m_model_spec.observations = m_model_spec.observations*m_scale;
|
|
_scalechanged();
|
|
_adjust_sample(first,last);
|
|
sv = starting_values();
|
|
res = resids(sv,vector::Zeros(0),matrix::Zeros(0,0));
|
|
}
|
|
|
|
//--- 4. Configure Pre-sample Backcast Values
|
|
vector bcast;
|
|
if(!backcast.Size())
|
|
bcast = m_vp.backCast(res);
|
|
else
|
|
bcast = m_vp.backCastTransform(backcast);
|
|
|
|
// Direct routing failures for exceptional model layouts
|
|
if(hascloseform)
|
|
return false;
|
|
if(total_params == 0.0)
|
|
return false;
|
|
|
|
// Initialize conditional variance pathways and bounds
|
|
vector sigma2 = vector::Zeros(res.Size());
|
|
m_backcast = bcast;
|
|
vector sv_v = m_vp.startingValues(res);
|
|
m_var_bounds = m_vp.varianceBounds(res);
|
|
matrix vb = m_var_bounds;
|
|
m_vp.computeVariance(sv_v,res,sigma2,bcast,vb);
|
|
vector std_res = res/sqrt(sigma2);
|
|
|
|
//--- 5. Collect and Format Inequality Constraint Systems ($A \cdot x \ge b$)
|
|
Constraints cst[3];
|
|
cst[0] = constraints();
|
|
cst[1] = m_vp.constraints();
|
|
cst[2] = m_distribution.constraints();
|
|
|
|
vector num_cons(3);
|
|
for(uint i = 0; i<cst.Size(); ++i)
|
|
num_cons[i] = (double)cst[i]._one.Rows();
|
|
|
|
// Allocate matrix rows matching stacked components, including linear boundary offsets
|
|
matrix a = matrix::Zeros(ulong(num_cons.Sum()),ulong(total_params)+1);
|
|
vector b = vector::Zeros(ulong(num_cons.Sum()));
|
|
long ren, cen, rst, ct;
|
|
vector nc, of;
|
|
|
|
for(uint i = 0; i<3; ++i)
|
|
{
|
|
nc = np::sliceVector(num_cons,0,long(i+1));
|
|
of = np::sliceVector(offsets,0,long(i+1));
|
|
ren = long(nc.Sum());
|
|
cen = long(of.Sum());
|
|
rst = ren - long(num_cons[i]);
|
|
ct = cen - long(offsets[i]);
|
|
if(ren-rst>0)
|
|
{
|
|
np::matrixCopy(a,cst[i]._one,rst,ren,STEP,ct,cen);
|
|
np::vectorCopy(b,cst[i]._two,rst,ren);
|
|
}
|
|
}
|
|
|
|
// Convert structural constraint records to standard linear offset profiles
|
|
a.Col(b,ulong(total_params));
|
|
b = vector::Ones(a.Rows());
|
|
|
|
//--- 6. Gather Boundary Matrices [Min / Max Box Constraints]
|
|
matrix bound[3];
|
|
bound[0] = bounds();
|
|
bound[1] = m_vp.bounds(res);
|
|
bound[2] = m_distribution.bounds(res);
|
|
matrix all_bounds(2,bound[0].Cols()+bound[1].Cols()+bound[2].Cols());
|
|
|
|
long cols = 0;
|
|
for(uint i = 0; i<bound.Size(); ++i)
|
|
{
|
|
if(bound[i].Cols())
|
|
{
|
|
cols+=(i)?long(bound[i-1].Cols()):0;
|
|
np::matrixCopyCols(all_bounds,bound[i],cols,long(cols+bound[i].Cols()));
|
|
}
|
|
}
|
|
|
|
//--- 7. Validate Provided Starting Vector Against Constraints
|
|
if(startingvalues.Size())
|
|
{
|
|
bool valid = (startingvalues.Size()==ulong(total_params));
|
|
matrix aa = np::sliceMatrixCols(a,0,long(a.Cols()-1));
|
|
vector bb = a.Col(a.Cols()-1);
|
|
if(aa.Rows()&&valid)
|
|
{
|
|
vector temp = aa.MatMul(startingvalues) - bb;
|
|
valid = (valid && temp.Min()>0.0);
|
|
}
|
|
for(uint i = 0; i<startingvalues.Size(); ++i)
|
|
valid = (valid && all_bounds[0,i]<=startingvalues[i] && startingvalues[i]<=all_bounds[1,i]);
|
|
if(!valid)
|
|
startingvalues = EMPTY_VECTOR;
|
|
}
|
|
|
|
// Setup initial fallback structural seeds if vectors are unassigned
|
|
if(!startingvalues.Size())
|
|
{
|
|
sv = starting_values();
|
|
ulong ss = sv.Size();
|
|
vector temp = m_distribution.startingValues(std_res);
|
|
sv.Resize(ss+sv_v.Size()+temp.Size());
|
|
if(sv_v.Size())
|
|
np::vectorCopy(sv,sv_v,long(ss),long(ss+sv_v.Size()));
|
|
if(temp.Size())
|
|
np::vectorCopy(sv,temp,long(ss+sv_v.Size()));
|
|
}
|
|
else
|
|
sv = startingvalues;
|
|
|
|
//--- 8. Execute Target Core Log-Likelihood Extraction Engine
|
|
#ifdef __SLSQP__
|
|
// Route evaluation through the internal SLSQP functional structure
|
|
CObjectiveF fvec;
|
|
fvec.setFuncParams(sigma2,bcast,vb,this);
|
|
out_loglikelihood = fvec.orig_fun(params);
|
|
#else
|
|
// Route evaluation through the standard ALGLIB functional object layout
|
|
CObjectiveFunc fvec;
|
|
fvec.SetParams(sigma2,bcast,vb,this);
|
|
CObject ob;
|
|
CRowDouble in_x = params;
|
|
fvec.Func(in_x,out_loglikelihood,ob);
|
|
#endif
|
|
|
|
return true;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| HAR model (Heterogeneous Autoregressive) |
|
|
//| Models the conditional mean by decomposing volatility into |
|
|
//| daily, weekly, and monthly components to capture long memory. |
|
|
//+------------------------------------------------------------------+
|
|
class HAR: public HARX
|
|
{
|
|
public:
|
|
// Default constructor: sets the model specification type to HAR
|
|
HAR(void)
|
|
{
|
|
m_model_spec.mean_model_type = MEAN_HAR;
|
|
}
|
|
|
|
// Overloaded constructor: initializes with a specific ArchParameters configuration
|
|
HAR(ArchParameters& model_specification)
|
|
{
|
|
model_specification.mean_model_type = MEAN_HAR;
|
|
initialize(model_specification);
|
|
}
|
|
|
|
// Destructor: ensures proper cleanup of inherited resources
|
|
~HAR(void)
|
|
{
|
|
deinitialize();
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| The constant mean model |
|
|
//| Represents a simple process: y_t = mu + epsilon_t |
|
|
//+------------------------------------------------------------------+
|
|
class ConstantMean: public HARX
|
|
{
|
|
public:
|
|
// Default constructor: assigns constant mean model type
|
|
ConstantMean(void)
|
|
{
|
|
m_model_spec.mean_model_type = MEAN_CONSTANT;
|
|
}
|
|
|
|
// Overloaded constructor: initializes and sets model type
|
|
ConstantMean(ArchParameters& model_specification)
|
|
{
|
|
model_specification.mean_model_type = MEAN_CONSTANT;
|
|
initialize(model_specification);
|
|
}
|
|
|
|
// Cleanup
|
|
~ConstantMean(void)
|
|
{
|
|
deinitialize();
|
|
}
|
|
|
|
// Returns the count of parameters (only 1: the mean mu)
|
|
virtual ulong num_params(void) override
|
|
{
|
|
return 1;
|
|
}
|
|
|
|
// Computes model residuals (epsilon_t = y_t - mu)
|
|
virtual vector resids(vector& params, vector& y, matrix& regressors) override
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return EMPTY_VECTOR;
|
|
}
|
|
// Use internal fit data if y is empty, otherwise use provided vector
|
|
if(!y.Size())
|
|
return np::minus(m_fit_y,params);
|
|
return np::minus(y,params);
|
|
}
|
|
|
|
// Generates synthetic data based on the constant mean and volatility process
|
|
virtual matrix simulate(vector ¶ms, ulong nobs, ulong burn, vector &initial_vals, matrix &x, vector &initial_vals_vol)
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return EMPTY_MATRIX;
|
|
}
|
|
|
|
// Ensure no external regressors are passed for a simple constant mean simulation
|
|
if(initial_vals.Size() || x.Cols())
|
|
{
|
|
Print(__FUNCTION__, " Both initial value and x must be none when simulating a constant mean process.");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
// Extract segmented parameters (mean, volatility process, distribution)
|
|
vector allparams[];
|
|
if(!_parse_parameters(params,allparams))
|
|
{
|
|
Print(__FUNCTION__, " failed to parse parameters ");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
// Initialize random number generator for innovations
|
|
BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
|
|
|
|
// Simulate underlying volatility path
|
|
matrix sim_values = m_vp.simulate(allparams[1],nobs+burn,Rng,burn,initial_vals_vol.Size()?initial_vals_vol[0]:0.0);
|
|
|
|
// Construct return series by adding mean to simulated residuals (Row 0 = errors, Row 1 = variance)
|
|
vector ers = sim_values.Row(0);
|
|
vector y = ers + allparams[0];
|
|
vector vol = sqrt(sim_values.Row(1));
|
|
|
|
// Remove burn-in period samples
|
|
y = np::sliceVector(y,long(burn));
|
|
vol = np::sliceVector(vol,long(burn));
|
|
ers = np::sliceVector(ers,long(burn));
|
|
|
|
// Package results into matrix [y | volatility | residuals]
|
|
matrix out(y.Size(),3);
|
|
if(!out.Col(y,0) || !out.Col(vol,1) || !out.Col(ers,2))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
return out;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| The zero mean model |
|
|
//| Assumes the return series has no deterministic mean (mu = 0). |
|
|
//| Used primarily for high-frequency data or stationary returns. |
|
|
//+------------------------------------------------------------------+
|
|
class ZeroMean: public HARX
|
|
{
|
|
public:
|
|
// Default constructor: assigns zero mean model type
|
|
ZeroMean(void)
|
|
{
|
|
m_model_spec.mean_model_type = MEAN_ZERO;
|
|
}
|
|
|
|
// Overloaded constructor: initializes with configuration
|
|
ZeroMean(ArchParameters& model_specification)
|
|
{
|
|
model_specification.mean_model_type = MEAN_ZERO;
|
|
initialize(model_specification);
|
|
}
|
|
|
|
// Cleanup
|
|
~ZeroMean(void)
|
|
{
|
|
deinitialize();
|
|
}
|
|
|
|
// Returns 0 parameters as there is no mean coefficient to estimate
|
|
virtual ulong num_params(void) override
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
// Computes residuals (epsilon_t = y_t)
|
|
virtual vector resids(vector& params, vector& y, matrix& regressors) override
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return EMPTY_VECTOR;
|
|
}
|
|
// If y is not provided, return the internal model data directly
|
|
if(!y.Size())
|
|
return m_fit_y;
|
|
return y;
|
|
}
|
|
|
|
// Generates synthetic data driven purely by the volatility process
|
|
virtual matrix simulate(vector ¶ms, ulong nobs, ulong burn, vector &initial_vals, matrix &x, vector &initial_vals_vol)
|
|
{
|
|
if(!is_initialized())
|
|
{
|
|
Print(__FUNCTION__, " object instance was not successfully initialized ");
|
|
return EMPTY_MATRIX;
|
|
}
|
|
|
|
// Ensure no initial conditions or external regressors are used
|
|
if(initial_vals.Size() || x.Cols())
|
|
{
|
|
Print(__FUNCTION__, " Both initial value and x must be none when simulating a zero mean process.");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
vector allparams[];
|
|
if(!_parse_parameters(params,allparams))
|
|
{
|
|
Print(__FUNCTION__, " failed to parse parameters ");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
// Initialize generator for innovations
|
|
BootstrapRng Rng(m_distribution.distribution_type(),allparams[2],(int)m_distribution.randomState());
|
|
|
|
// Simulate volatility path (Row 0 = innovations, Row 1 = variance)
|
|
matrix sim_values = m_vp.simulate(allparams[1],nobs+burn,Rng,burn,initial_vals_vol[0]);
|
|
|
|
vector ers = sim_values.Row(0);
|
|
vector y = ers; // y_t = epsilon_t
|
|
vector vol = sqrt(sim_values.Row(1));
|
|
|
|
// Remove burn-in samples
|
|
y = np::sliceVector(y,long(burn));
|
|
vol = np::sliceVector(vol,long(burn));
|
|
ers = np::sliceVector(ers,long(burn));
|
|
|
|
// Package output: [y | volatility | residuals]
|
|
matrix out(y.Size(),3);
|
|
if(!out.Col(y,0) || !out.Col(vol,1) || !out.Col(ers,2))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
return out;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| The AR model (Autoregressive) |
|
|
//| Models y_t as a function of previous observations (lags). |
|
|
//| Captures autocorrelation present in return series. |
|
|
//+------------------------------------------------------------------+
|
|
class AR: public HARX
|
|
{
|
|
public:
|
|
// Default constructor: assigns AR model type
|
|
AR(void)
|
|
{
|
|
m_model_spec.mean_model_type = MEAN_AR;
|
|
}
|
|
|
|
// Overloaded constructor: initializes with configuration
|
|
AR(ArchParameters& model_specification)
|
|
{
|
|
model_specification.mean_model_type = MEAN_AR;
|
|
initialize(model_specification);
|
|
}
|
|
|
|
// Clean up resources
|
|
~AR(void)
|
|
{
|
|
deinitialize();
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| The AR-X model (Autoregressive with Exogenous Variables) |
|
|
//| Extends the AR model by incorporating external factors (X). |
|
|
//| Useful for multifactor modeling and structural analysis. |
|
|
//+------------------------------------------------------------------+
|
|
class ARX: public HARX
|
|
{
|
|
public:
|
|
// Default constructor: assigns ARX model type
|
|
ARX(void)
|
|
{
|
|
m_model_spec.mean_model_type = MEAN_ARX;
|
|
}
|
|
|
|
// Overloaded constructor: initializes with configuration
|
|
ARX(ArchParameters& model_specification)
|
|
{
|
|
model_specification.mean_model_type = MEAN_ARX;
|
|
initialize(model_specification);
|
|
}
|
|
|
|
// Clean up resources
|
|
~ARX(void)
|
|
{
|
|
deinitialize();
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------
|