//+------------------------------------------------------------------+ //| volatility.mqh | //| Copyright 2025, MetaQuotes Ltd. | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Copyright 2025, MetaQuotes Ltd." #property link "https://www.mql5.com" #include "recursions.mqh" #include "distribution.mqh" //--- //+------------------------------------------------------------------+ //| Container for variance forecasts: Stores point-in-time predictions | //| and full stochastic simulation paths for GARCH-family models. | //+------------------------------------------------------------------+ struct VarianceForecast { matrix forecasts; // Holds point variance estimates [steps x horizons] matrix forecastpaths[]; // Array of simulation paths for each forecast origin matrix shocks[]; // Simulated residual shocks corresponding to paths // Default constructor VarianceForecast(void) { forecasts = matrix::Zeros(0, 0); } // Parameterized constructor to initialize from existing matrices VarianceForecast(matrix &_forecasts, matrix &_forecastpaths[], matrix &_shocks[]) { forecasts = _forecasts; // Copy simulation paths and associated shocks ArrayResize(forecastpaths, _forecastpaths.Size()); ArrayResize(shocks, _shocks.Size()); for(uint i = 0 ; i < shocks.Size(); shocks[i] = _shocks[i], ++i); for(uint i = 0 ; i < forecastpaths.Size(); forecastpaths[i] = _forecastpaths[i], ++i); } // Copy constructor VarianceForecast(VarianceForecast &other) { forecasts = other.forecasts; // Deep copy array members ArrayResize(forecastpaths, other.forecastpaths.Size()); ArrayResize(shocks, other.shocks.Size()); for(uint i = 0 ; i < shocks.Size(); shocks[i] = other.shocks[i], ++i); for(uint i = 0 ; i < forecastpaths.Size(); forecastpaths[i] = other.forecastpaths[i], ++i); } // Assignment operator for deep object copying void operator=(VarianceForecast &other) { forecasts = other.forecasts; // Resize and copy buffers from source object ArrayResize(forecastpaths, other.forecastpaths.Size()); ArrayResize(shocks, other.shocks.Size()); for(uint i = 0 ; i < shocks.Size(); shocks[i] = other.shocks[i], ++i); for(uint i = 0 ; i < forecastpaths.Size(); forecastpaths[i] = other.forecastpaths[i], ++i); } }; //+------------------------------------------------------------------+ //| Base class for all volatility processes (GARCH, EGARCH, etc.) | //| Enforces interface consistency for estimation and forecasting. | //+------------------------------------------------------------------+ class CVolatilityProcess { protected: bool m_updateable; // Can model parameters be updated? bool m_initialized; // Is the model ready for computation? ulong m_num_params; // Number of model parameters to optimize bool m_closedform; // Does model have an analytical forecast? ulong m_bootstrap_obs; // Minimum observations required for bootstrap long m_start, m_stop; // Data range indices CNormal m_normal; // Gaussian distribution reference VolatilityUpdater *m_volupdater; // Polymorphic helper for variance recursion string m_name; // Model identification string int m_seed; // Seed for stochastic simulations ENUM_VOLATILITY_MODEL m_model; // Model type identifier //--- Compute the variance for a single observation (Overridden by child) virtual double _update(ulong index, vector& parameters, vector &resids, vector &sigma2, vector &backcast, vector &var_bounds) { return EMPTY_VALUE; } //--- Verify if a forecasting method (e.g., Analytic vs Simulation) is valid virtual bool _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) { return false; } //--- One-step ahead forecast: Projects variance forward by one period virtual bool _onestepforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, vector&out_sigma, matrix&out_forecasts) { ulong t = resids.Size(); vector _resids = resids; _resids.Resize(t + 1); // Expand to include future residual _resids[t] = 0.; // Assume zero innovation for point forecast matrix _vb = varbounds; _vb.Resize(varbounds.Rows() + 1, varbounds.Cols()); vector temp = {0., double("inf")}; _vb.Row(temp, varbounds.Rows()); // Extend bounds vector sigma2 = vector::Zeros(t + 1); computeVariance(parameters, _resids, sigma2, backcast, _vb); matrix forecasts = matrix::Zeros(t - _start, horizon); temp = np::sliceVector(sigma2, _start + 1); forecasts.Col(temp, 0); out_sigma = np::sliceVector(sigma2, 0, -1); out_forecasts = forecasts; return true; } //--- Computes Gaussian Log-Likelihood, useful for initial parameter estimation virtual double _gaussianloglikelihood(vector ¶meters, vector& resids, vector&backcast, matrix& varbounds) { vector sigma2 = vector::Zeros(resids.Size()); computeVariance(parameters, resids, sigma2, backcast, varbounds); vector empty = vector::Zeros(0); vector v = m_normal.loglikelihood(empty, resids, sigma2); return v[0]; } //--- Analytic multi-step forecasts: Mathematical projection of variance virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) { return VarianceForecast(); // Default: must be implemented by subclasses } //--- Simulation-based forecasts: Generates Monte Carlo paths of future variance virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) { return VarianceForecast(); } //--- Bootstrap forecast: Uses historical residuals to simulate future volatility virtual VarianceForecast _bootstrapforecast(vector& parameters, vector& resids, vector& backcast, matrix &varbounds, ulong _start, ulong horizon, ulong simulations, int seed = 0) { VarianceForecast out; vector sigma2 = vector::Zeros(resids.Size()); computeVariance(parameters, resids, sigma2, backcast, varbounds); // Calculate standardized residuals to resample from vector std_resid = resids / sqrt(sigma2); // Safety check for data sufficiency if(_start < m_bootstrap_obs) { Print(__FUNCTION__, " start must include more than ", m_bootstrap_obs, " observations"); return out; } BootstrapRng Rng(std_resid, _start, seed); out = _simulationforecast(parameters, resids, backcast, varbounds, _start, horizon, simulations, Rng); return out; } //--- Internal initialization routine to set model hyper-parameters virtual bool _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100) { m_updateable = updateable; m_num_params = nparams; m_closedform = closedform; bootstraps(bootstrap_obs); m_start = start; m_stop = stop; m_name = name; m_seed = seed; m_model = volmodel; return m_normal.initialize(vector::Zeros(0), seed); } //+------------------------------------------------------------------+ //| Public Interface for CVolatilityProcess | //+------------------------------------------------------------------+ public: // Constructor: Initializes defaults for model configuration CVolatilityProcess(void): m_updateable(true), m_num_params(0), m_closedform(false), m_bootstrap_obs(100), m_start(0), m_stop(-1), m_initialized(false), m_name(NULL), m_volupdater(NULL), m_model(WRONG_VALUE) {} // Copy constructor: Performs shallow copy of configuration properties CVolatilityProcess(CVolatilityProcess& other) { m_updateable = other.upDateable(); m_num_params = other.numParams(); m_closedform = other.closedForm(); m_bootstrap_obs = other.bootstraps(); m_start = other.start(); m_stop = other.stop(); m_volupdater = other.volUpdater(); m_name = other.name(); m_model = other.volatilityprocess(); } // Destructor: Safely cleans up the dynamically allocated volupdater ~CVolatilityProcess(void) { if(CheckPointer(m_volupdater) == POINTER_DYNAMIC) delete m_volupdater; m_volupdater = NULL; } // Assignment operator: Updates object state from another process void operator=(CVolatilityProcess& other) { m_updateable = other.upDateable(); m_num_params = other.numParams(); m_closedform = other.closedForm(); m_bootstrap_obs = other.bootstraps(); m_start = other.start(); m_stop = other.stop(); m_volupdater = other.volUpdater(); m_name = other.name(); m_model = other.volatilityprocess(); } // --- Accessors --- ENUM_VOLATILITY_MODEL volatilityprocess(void) { return m_model; } ulong bootstraps(void) { return m_bootstrap_obs; } // Setters with validation void bootstraps(ulong bts) { if(bts > 10) m_bootstrap_obs = bts; } string name(void) { return m_name; } bool is_initialized(void) { return m_initialized; } // Index management virtual long start(void) { return m_start; } virtual long stop(void) { return m_stop; } virtual void start(long _start) { m_start = _start; } virtual void stop(long _stop) { m_stop = _stop; } // Model attributes virtual ulong numParams(void) { return m_num_params; } virtual bool upDateable(void) { return m_updateable; } virtual bool closedForm(void) { return m_closedform; } virtual VolatilityUpdater* volUpdater(void) { return m_volupdater; } // --- Core Methods --- // Compute variance for a specific index (Default interface) virtual double upDate(ulong _index, vector& parameters, vector& resids, vector& sigma2, vector& backcast, vector& varbounds) { return EMPTY_VALUE; } // Construct loose bounds to ensure variance remains positive/stable during optimization virtual matrix varianceBounds(vector& resids, double power = 2.) { ulong nobs = resids.Size(); ulong tau = MathMin(75, nobs); vector w = vector::Zeros(tau); for(ulong i = 0; i < tau; w[i] = pow(0.94, double(i)), ++i); w = w / w.Sum(); double initial_value = w.Dot(pow(np::sliceVector(resids, 0, long(tau)), 2.)); vector varbound = vector::Zeros(nobs); ewma_recursion(0.94, resids, varbound, resids.Size(), initial_value); matrix varbounds = matrix::Zeros(varbound.Size(), 2); varbounds.Col(varbound / 1.e6, 0); // Define lower bound varbounds.Col(varbound * 1.e6, 1); // Define upper bound double var = resids.Var(); double min_upper_bound = 1. + (pow(resids, 2.)).Max(); double lower_bound = var / 1.e8; double upper_bound = 1.e7 * (1. + (pow(resids, 2.)).Max()); vector col0 = varbounds.Col(0); vector col1 = varbounds.Col(1); // Clip values to prevent extreme numerical overflow/underflow if(!col0.Clip(lower_bound, DBL_MAX) || !col1.Clip(min_upper_bound, upper_bound) || !varbounds.Col(col0, 0) || !varbounds.Col(col1, 1)) { Print(__FUNCTION__, " error ", __LINE__, " ", GetLastError()); return matrix::Zeros(0, 0); } if(power != 2.) varbounds = pow(varbounds, power / 2.); return varbounds; } // --- Initialization and Estimation Helpers --- virtual vector startingValues(vector& resids) { return vector::Zeros(0); } // Backcast: Estimate initial variance before the start of the observed data virtual vector backCast(vector& resids) { ulong tau = MathMin(75, resids.Size()); vector w = vector::Zeros(tau); for(ulong i = 0; i < tau; w[i] = pow(0.94, double(i)), ++i); w = w / w.Sum(); vector out(1); out[0] = (pow(np::sliceVector(resids, 0, long(tau)), 2.0) * w).Sum(); return out; } virtual vector backCastTransform(vector& backcast) { if(backcast.Min() < 0) { Print(__FUNCTION__, " User backcast value must be strictly positive "); return vector::Zeros(0); } return backcast; } virtual matrix bounds(vector& resids) { return matrix::Zeros(0, 0); } virtual vector computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) { return vector::Zeros(0); } virtual Constraints constraints(void) { return Constraints(); } // --- Polymorphic Forecast Engine --- // Routes the forecast request to analytic, simulation, or bootstrap logic VarianceForecast forecast(vector ¶meters, vector& resids, vector& _backcast, matrix& var_bounds, BootstrapRng &rng, int seed = 0, ulong _start = 0, ulong horizon = 1, ENUM_FORECAST_METHOD ForecastingMethod = FORECAST_ANALYTIC, ulong simulations = 1000) { VarianceForecast out; if(!horizon) { Print(__FUNCTION__, " horizon must be >= 1"); return out; } if(!_check_forecasting_method(ForecastingMethod, horizon)) { Print(__FUNCTION__, " Method not supported for this model"); return out; } if(!_start) _start = resids.Size() - 1; switch(ForecastingMethod) { case FORECAST_ANALYTIC: out = _analyticforecast(parameters, resids, _backcast, var_bounds, _start, horizon); break; case FORECAST_SIMULATION: if(rng.is_initialized()) out = _simulationforecast(parameters, resids, _backcast, var_bounds, _start, horizon, simulations, rng); else Print(__FUNCTION__, " ERROR: BOOTSTRAP OBJECT IS NULL "); break; case FORECAST_BOOTSTRAP: if(_start < 10) Print(__FUNCTION__, " Bootstrap requires > 10 obs"); else if(double(horizon / _start) > 0.2) Print(__FUNCTION__, " Ratio of horizon-to-start < 20% required. "); else out = _bootstrapforecast(parameters, resids, _backcast, var_bounds, _start, horizon, simulations, seed); break; } return out; } virtual matrix simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) { return matrix::Zeros(0, 0); } virtual string parameterNames(void) { return NULL; } virtual vector garch_spec(void) { return vector::Zeros(4); } }; //+------------------------------------------------------------------+ //| Constant volatility process | //+------------------------------------------------------------------+ class CConstantVariance: public CVolatilityProcess { protected: //--- Verify the requested forecasting method as valid for the specification virtual bool _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override { return true; } //--- Analytic multi-step volatility forecasts from the model virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override { VarianceForecast out; long t = (long)resids.Size(); matrix forecasts = matrix::Zeros(t-_start,horizon); forecasts.Fill(parameters[0]); out.forecasts = forecasts; return out; } //--- Simulation-based volatility forecasts from the model virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) { VarianceForecast out; long t = (long)resids.Size(); matrix forecasts = matrix::Zeros(t-_start,horizon); ArrayResize(out.forecastpaths,int(t-_start)); ArrayResize(out.shocks,int(t-_start)); forecasts.Fill(parameters[0]); for(long i = 0; i 1 when power != 2"); return false; } return true; } //--- Analytic multi-step volatility forecasts from the model //--- Computes recursive variance projections using the GARCH/GJR-GARCH structure virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override { VarianceForecast out; ulong t = resids.Size(); vector sigma2; matrix forecasts; // Compute the immediate next-step forecast (h=1) _onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sigma2, forecasts); if(horizon == 1) { out.forecasts = forecasts; return out; } // Extract GARCH model components from parameter vector double omega = parameters[0]; vector alpha = np::sliceVector(parameters, 1, long(m_p + 1)); // ARCH terms vector gamma = np::sliceVector(parameters, long(m_p + 1), long(m_p + m_o + 1)); // Leverage/Asymmetry terms vector beta = np::sliceVector(parameters, long(m_p + m_o + 1)); // GARCH (persistence) terms // Determine maximum lag depth required for recursion ulong m = MathMax(MathMax(m_p, m_o), m_q); vector _resids, _asym_resids, _sigma2, temp; _resids = _asym_resids = _sigma2 = vector::Zeros(m + horizon); // Iterate through historical data to set up the recursive projection start point for(ulong i = ulong(_start); i < t; ++i) { // If historical data exists for the full lag depth 'm' if(i - m + 1 >= 0) { temp = np::sliceVector(resids, long(i - m + 1), long(i + 1)); np::vectorCopy(_resids, temp, 0, long(m)); temp = np::sliceVector(_resids, 0, long(m)); temp *= np::whereVectorIsLt(temp, 0.0); // Identify negative shocks for asymmetry np::vectorCopy(_asym_resids, temp, 0, long(m)); temp = np::sliceVector(sigma2, long(i - m + 1), long(i + 1)); np::vectorCopy(_sigma2, temp, 0, long(m)); } // If we are at the start of the series (padding with backcast values) else { np::fillVector(_resids, sqrt(backcast[0]), 0, long(m - i - 1)); temp = np::sliceVector(resids, 0, long(i + 1)); np::vectorCopy(_resids, temp, long(m - i - 1), long(m)); _asym_resids = _resids * (np::whereVectorIsLt(_resids, 0.0)); np::fillVector(_asym_resids, sqrt(0.5 * backcast[0]), 0, long(m - i - 1)); np::fillVector(_sigma2, backcast[0], 0, long(m)); temp = np::sliceVector(sigma2, 0, long(i + 1)); np::vectorCopy(_sigma2, temp, long(m - i - 10), long(m)); } // Perform recursive forecasting for the requested horizon for(ulong h = 0; h < horizon; ++h) { ulong floc = i - ulong(_start); forecasts[floc, h] = omega; ulong sloc = h + m - 1; // Add ARCH components: alpha * eps^2_{t-j} for(ulong j = 0; j < m_p; ++j) forecasts[floc, h] += alpha[j] * pow(_resids[sloc - j], 2.); // Add GJR/Asymmetry components: gamma * I_{t-j} * eps^2_{t-j} for(ulong j = 0; j < m_o; ++j) forecasts[floc, h] += (gamma[j] * pow(_asym_resids[sloc - j], 2.)); // Add GARCH persistence components: beta * sigma^2_{t-j} for(ulong j = 0; j < m_q; ++j) forecasts[floc, h] += beta[j] * _sigma2[sloc - j]; // Update buffers with the new forecast to use for next-step recursion _resids[h + m] = sqrt(forecasts[floc, h]); _asym_resids[h + m] = sqrt(0.5 * forecasts[floc, h]); _sigma2[h + m] = forecasts[floc, h]; } } out.forecasts = forecasts; return out; } //--- Simulates paths for forecasting operation using Monte Carlo methods //--- Projects future volatility by iteratively calculating shocks and updating variance virtual bool _simulatepaths(ulong m, vector& parameters, ulong horizon, matrix& std_shocks, matrix& scaled_forecast_paths, matrix& scaled_shock, matrix& asym_scaled_shock, vector& mean_fpaths, matrix& fpaths, matrix& shocks) { // Extract model parameters: Omega (constant), Alpha (ARCH), Gamma (Asymmetry/Leverage), Beta (GARCH) double omega = parameters[0]; vector alpha = ((m_p + 1) > 1) ? np::sliceVector(parameters, 1, long(m_p + 1)) : EMPTY_VECTOR; vector gamma = ((m_p + 1) < (m_p + m_o + 1)) ? np::sliceVector(parameters, long(m_p + 1), long(m_p + m_o + 1)) : EMPTY_VECTOR; vector beta = (parameters.Size() > (m_p + m_o + 1)) ? np::sliceVector(parameters, long(m_p + m_o + 1)) : EMPTY_VECTOR; matrix shock = scaled_forecast_paths; vector temp; // Iterate through the forecast horizon for(ulong h = 0; h < horizon; ++h) { ulong loc = h + m - 1; temp = vector::Zeros(scaled_forecast_paths.Rows()); temp.Fill(omega); scaled_forecast_paths.Col(temp, h + m); // Set baseline to constant omega // Add ARCH contributions: alpha * epsilon^2 for(ulong j = 0; j < m_p; ++j) { temp = scaled_forecast_paths.Col(h + m); temp += alpha[j] * scaled_shock.Col(loc - j); scaled_forecast_paths.Col(temp, h + m); } // Add Asymmetry/GJR-GARCH contributions: gamma * (indicator < 0) * epsilon^2 for(ulong j = 0; j < m_o; ++j) { temp = scaled_forecast_paths.Col(h + m); temp += gamma[j] * asym_scaled_shock.Col(loc - j); scaled_forecast_paths.Col(temp, h + m); } // Add GARCH persistence contributions: beta * sigma^2 for(ulong j = 0; j < m_q; ++j) { temp = scaled_forecast_paths.Col(h + m); temp += beta[j] * scaled_forecast_paths.Col(loc - j); scaled_forecast_paths.Col(temp, h + m); } // Generate the innovation (shock) for the current step based on projected volatility temp = std_shocks.Col(h) * pow(scaled_forecast_paths.Col(h + m), 1. / m_power); shock.Col(temp, h + m); // Calculate squared shock and flag asymmetry (negative shocks) for next recursion step vector lt_zero = np::whereVectorIsLt(shock.Col(h + m), 0.); scaled_shock.Col(pow(fabs(shock.Col(h + m)), m_power), h + m); asym_scaled_shock.Col(scaled_shock.Col(h + m) * lt_zero, h + m); } // Prepare final outputs: trim buffers, convert to variance scale, and calculate path means fpaths = np::sliceMatrixCols(scaled_forecast_paths, long(m)); fpaths = pow(fpaths, 2. / m_power); // Standardize back to variance scale mean_fpaths = fpaths.Mean(0); // Average of all simulation paths shocks = np::sliceMatrixCols(shock, long(m)); return true; } //--- Simulation-based volatility forecasts from the model //--- Performs Monte Carlo simulations to project future variance paths for GARCH/GJR-GARCH models virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) { VarianceForecast out; vector sigma2; matrix forecasts; // Compute base variance (one-step) to initialize simulation context _onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sigma2, forecasts); long t = (long)resids.Size(); ArrayResize(out.forecastpaths, int(t - _start)); ArrayResize(out.shocks, int(t - _start)); ulong m = MathMax(MathMax(m_p, m_o), m_q); // Determine maximum lag depth // Pre-allocate matrices for simulation state (paths, squared residuals, and asymmetry/leverage components) matrix scaled_forecast_paths = matrix::Zeros(simulations, m + horizon); matrix scaled_shock = matrix::Zeros(simulations, m + horizon); matrix asym_scaled_shock = matrix::Zeros(simulations, m + horizon); matrix temp_m; vector temp_v, mask; ulong count = 0; // Iterate over each historical observation to initialize simulation paths for(ulong i = ulong(_start); i < ulong(t); ++i) { matrix std_shocks = rng.rng(simulations, horizon); // Initialization for early periods (using backcast values if historical data is insufficient) if(i - m < 0) { np::fillMatrix(scaled_forecast_paths, pow(backcast[0], m_power / 2.), BEGIN, END, STEP, BEGIN, long(m)); np::fillMatrix(scaled_shock, pow(backcast[0], m_power / 2.), BEGIN, END, STEP, BEGIN, long(m)); np::fillMatrix(asym_scaled_shock, 0.5 * pow(backcast[0], m_power / 2.), BEGIN, END, STEP, BEGIN, long(m)); count = i + 1; // Populate initial buffers with historical sigma and residuals temp_v = np::sliceVector(sigma2, 0, long(count)); temp_v = pow(temp_v, m_power / 2.); temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_forecast_paths.Rows()); np::matrixCopy(scaled_forecast_paths, temp_m, BEGIN, END, STEP, long(m - count), long(m)); temp_v = np::sliceVector(resids, 0, long(count)); temp_v = pow(fabs(temp_v), m_power); temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_shock.Rows()); np::matrixCopy(scaled_shock, temp_m, BEGIN, END, STEP, long(m - count), long(m)); temp_v = np::sliceVector(resids, 0, long(count)); temp_v = pow(fabs(temp_v), m_power) * np::whereVectorIsLt(temp_v, 0.0); temp_m = np::repeat_vector_as_rows_cols(temp_v, asym_scaled_shock.Rows()); np::matrixCopy(asym_scaled_shock, temp_m, BEGIN, END, STEP, long(m - count), long(m)); } // Initialization for standard periods using historical data slices else { temp_v = np::sliceVector(sigma2, long(i - m + 1), long(i + 1)); temp_v = pow(temp_v, m_power / 2.); temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_forecast_paths.Rows()); np::matrixCopy(scaled_forecast_paths, temp_m, BEGIN, END, STEP, BEGIN, long(m)); temp_v = np::sliceVector(resids, long(i - m + 1), long(i + 1)); mask = np::whereVectorIsLt(temp_v, 0.0); temp_v = pow(fabs(temp_v), m_power); temp_m = np::repeat_vector_as_rows_cols(temp_v, scaled_shock.Rows()); np::matrixCopy(scaled_shock, temp_m, BEGIN, END, STEP, BEGIN, long(m)); temp_m = np::sliceMatrixCols(scaled_shock, 0, long(m)); temp_v = np::sliceVector(resids, long(i - m + 1), long(i + 1)); mask = np::whereVectorIsLt(temp_v, 0.0); temp_m = np::multiply(temp_m, mask); np::matrixCopy(asym_scaled_shock, temp_m, BEGIN, END, STEP, BEGIN, long(m)); } // Generate paths via internal simulation method vector f; matrix p, s; _simulatepaths(m, parameters, horizon, std_shocks, scaled_forecast_paths, scaled_shock, asym_scaled_shock, f, p, s); // Store results for the current path ulong loc = i - long(_start); out.forecasts.Row(f, loc); out.forecastpaths[loc] = p; out.shocks[loc] = s; } return out; } //--- Initialize class members and configure GARCH model structure virtual bool _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 1, ulong o = 0, ulong q = 1, double power = 2., long start = 0, long stop = -1, ulong bootstrap_obs = 100) override { // Assign basic model configuration m_model = volmodel; m_updateable = updateable; m_num_params = nparams; m_closedform = closedform; bootstraps(bootstrap_obs); m_start = start; m_stop = stop; m_name = name; m_seed = seed; // Synchronize random number generator seed if not already initialized if(m_normal.randomState() != uint(m_seed)) m_normal.initialize(vector::Zeros(0), seed); // Configure model lag orders and power transformation parameter m_p = p; // ARCH lags m_o = o; // Asymmetry/Leverage lags m_q = q; // GARCH (persistence) lags m_power = power; // Validate power parameter (must be positive, e.g., 2.0 for standard GARCH) if(m_power <= 0.0) { Print(__FUNCTION__, " invalid value for power variable "); return false; } // Set default lag values if not explicitly provided, based on the specific model type switch(m_model) { case VOL_ARCH: case VOL_AVARCH: if(!m_p) m_p = 1; break; case VOL_GARCH: case VOL_AVGARCH: if(!m_p) m_p = 1; if(!m_q) m_q = 1; break; case VOL_GJR_GARCH: case VOL_TARCH: if(!m_p) m_p = 1; if(!m_q) m_q = 1; if(!m_o) m_o = 1; break; } // Recalculate total parameter count: 1 (omega) + p + o + q m_num_params = 1 + m_p + m_o + m_q; return (true); } public: // Constructor: Default initialization for a standard GARCH(1,1) model CGarchProcess(void) { m_initialized = CGarchProcess::_initialize(VOL_GARCH, true, false, 0, "GARCH", 0, 1, 0, 1, 2.0); } // Constructor: Parameterized initialization for custom lag structures and seeding CGarchProcess(ulong p, ulong q, int seed = 0, ulong nbootstraps = 100) { m_initialized = CGarchProcess::_initialize(VOL_GARCH, true, false, 0, "GARCH", seed, p, 0, q, 2.0, 0, -1, nbootstraps); } // Copy Constructor: Creates a new instance by copying properties from an existing GARCH process CGarchProcess(CGarchProcess& other) { m_updateable = other.upDateable(); m_num_params = other.numParams(); m_closedform = other.closedForm(); m_bootstrap_obs = other.bootstraps(); m_start = other.start(); m_stop = other.stop(); m_volupdater = other.volUpdater(); m_name = other.name(); m_model = other.volatilityprocess(); m_q = other.get_q(); m_p = other.get_p(); m_o = other.get_o(); m_power = other.get_power(); } // Destructor ~CGarchProcess(void) { } // Assignment Operator: Copies properties from another GARCH process to this instance void operator=(CGarchProcess& other) { m_updateable = other.upDateable(); m_num_params = other.numParams(); m_closedform = other.closedForm(); m_bootstrap_obs = other.bootstraps(); m_start = other.start(); m_stop = other.stop(); m_volupdater = other.volUpdater(); m_name = other.name(); m_model = other.volatilityprocess(); m_q = other.get_q(); m_p = other.get_p(); m_o = other.get_o(); m_power = other.get_power(); } //--- Getters for model specifications (lags and power) virtual ulong get_p(void) { return m_p; } virtual ulong get_o(void) { return m_o; } virtual ulong get_q(void) { return m_q; } virtual double get_power(void) { return m_power; } //--- Construct boundaries for conditional variance estimation //--- Overrides base class to use the specific m_power exponent defined for this GARCH instance virtual matrix varianceBounds(vector& resids, double power = 2.) override { return CVolatilityProcess::varianceBounds(resids, m_power); } //--- Compute the variance for the GARCH model using recursive updates virtual vector computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) override { // Transform residuals based on the model power (e.g., standard GARCH uses power=2) vector fresids = pow(MathAbs(resids), m_power); // Capture sign of residuals to account for asymmetry (leverage effects) in GJR-GARCH models vector sresids(resids.Size()); for(ulong i = 0; i < sresids.Size(); sresids[i] = np::sign(resids[i]), ++i); // Run recursive GARCH filtering to estimate conditional variance series _sigma2 = garch_recursion(parameters, fresids, sresids, _sigma2, m_p, m_o, m_q, resids.Size(), _backcast[0], varbounds); // Scale back to original variance units (inverse power transformation) double invp = 2. / m_power; return pow(_sigma2, invp); } //--- Returns starting values for the GARCH model estimation //--- Performs a grid search over potential starting parameter combinations to improve solver convergence virtual vector startingValues(vector& resids) override { // Define candidate grid values for alphas, gammas, and persistence (betas) vector alphas = {0.01, 0.05, 0.1, 0.2}; vector gammas = alphas; vector abg = {0.5, 0.7, 0.9, 0.98}; // Generate Cartesian product of search candidates matrix abgs = _product(alphas, gammas, abg); // Scale target variance based on data properties to ensure reasonable starting scale double target = (pow(fabs(resids), m_power)).Mean(); double scale = (pow(resids, 2.)).Mean() / pow(target, 2. / m_power); target *= pow(scale, m_power / 2.); vector svs[]; ArrayResize(svs, int(abgs.Rows())); matrix vb = varianceBounds(resids); vector llfs = vector::Zeros(abgs.Rows()); vector bc = backCast(resids); // Evaluate Log-Likelihood for each combination in the grid for(uint i = 0; i < svs.Size(); ++i) { vector row = abgs.Row(i); // Construct starting vector: constant omega + lags vector sv = (1. - row[2]) * target * vector::Ones(m_p + m_o + m_q + 1); // Distribute coefficients across ARCH (p), Leverage (o), and GARCH (q) components if(m_p > 0) { np::fillVector(sv, row[0] / double(m_p), long(1), long(1 + m_p)); row[2] -= row[0]; } if(m_o > 0) { np::fillVector(sv, row[1] / double(m_o), long(1 + m_p), long(1 + m_p + m_o)); row[2] -= row[1] / 2.; // Adjust persistence for asymmetric impact } if(m_q > 0) np::fillVector(sv, row[2] / double(m_q), long(1 + m_p + m_o), long(1 + m_p + m_o + m_q)); svs[i] = sv; llfs[i] = _gaussianloglikelihood(sv, resids, bc, vb); } // Select the parameter set that yielded the highest Log-Likelihood ulong loc = llfs.ArgMax(); return svs[loc]; } //--- Construct parameter constraints arrays for parameter estimation //--- Enforces positivity constraints and stationarity: sum(alpha + gamma/2 + beta) < 1 virtual Constraints constraints(void) override { Constraints out; ulong k_arch = m_p + m_o + m_q; out._one = matrix::Zeros(k_arch + 2, k_arch + 1); // Enforce positivity: omega > 0 and all coefficients > 0 for(ulong i = 0; i < (k_arch + 1); ++i) out._one[i, i] = 1.; // Additional constraints for asymmetry terms (if applicable) for(ulong i = 0; i < (m_o); ++i) if(i < m_p) out._one[i + m_p + 1, i + 1] = 1.; // Enforce stationarity: sum(alpha_i) + 0.5 * sum(gamma_i) + sum(beta_i) < 1 // Implemented as linear inequality: - (sum(alpha) + 0.5*sum(gamma) + sum(beta)) >= -1 for(ulong i = 1; i < (k_arch + 1); ++i) out._one[k_arch + 1, i] = -1.; // Apply 0.5 multiplier for leverage (asymmetry) terms per GJR-GARCH definition for(ulong i = m_p + 1; i < (m_p + m_o + 1); ++i) out._one[k_arch + 1, i] = -0.5; out._two = vector::Zeros(k_arch + 2); out._two[k_arch + 1] = -1.; // Sets the stationarity limit return out; } //--- Transformation to apply to user-provided backcast values //--- Converts variance-scale backcast to the model's power-scale (e.g., standard deviation for power=2) virtual vector backCastTransform(vector& backcast) override { backcast = CVolatilityProcess::backCastTransform(backcast); if(backcast.Size()) return pow(sqrt(backcast), m_power); else return vector::Zeros(0); } //--- Construct values for backcasting to start the recursion //--- Uses an Exponentially Weighted Moving Average (EWMA) of past residuals for initialization virtual vector backCast(vector& resids) override { ulong tau = MathMin(75, resids.Size()); // Create decay weights for the EWMA vector out = np::arange(tau); for(ulong i = 0; i < out.Size(); out[i] = pow(0.94, out[i]), ++i); out = out / out.Sum(); // Calculate weighted average of squared (or power-transformed) residuals vector temp = np::sliceVector(resids, 0, long(tau)); out = pow(MathAbs(temp), m_power) * out; vector bc(1); bc[0] = out.Sum(); return bc; } //--- Returns parameter bounds for optimization //--- Defines the allowable range [min, max] for each GARCH parameter virtual matrix bounds(vector& resids) override { // Calculate the mean of transformed residuals as a scale factor for variance double v = (pow(MathAbs(resids), m_power)).Mean(); // Create a bounds matrix: rows = parameters (1 omega + p ARCH + o leverage + q GARCH), 2 columns (min, max) matrix out = matrix::Zeros(1 + m_p + m_o + m_q, 2); // Bounds for constant term (omega): near zero to 10x the mean variance out[0, 0] = 1e-8 * v; out[0, 1] = 10.0 * v; ulong from = 1; // Bounds for ARCH coefficients (alpha): range [0, 1] for(ulong i = 0; i < m_p; ++i) { out[from, 0] = 0.0; out[from++, 1] = 1.; } // Bounds for Leverage/Asymmetry coefficients (gamma): // For GJR-GARCH, if leverage term overlaps with ARCH lag, allow negative values (range [-1, 2]) // otherwise, standard positive constraint (range [0, 2]) for(ulong i = 0; i < m_o; ++i) { if(i < m_p) { out[from, 0] = -1.0; out[from++, 1] = 2.; } else { out[from, 0] = 0.0; out[from++, 1] = 2.; } } // Bounds for GARCH persistence coefficients (beta): range [0, 1] for(ulong i = 0; i < m_q; ++i) { out[from, 0] = 0.0; out[from++, 1] = 1.; } // Return transposed matrix to match the expected format [2 x N] for the optimizer return out.Transpose(); } //--- Names of model parameters virtual string parameterNames(void) override { return common_names(m_p,m_o,m_q); } //--- Simulate data from the GARCH/GJR-GARCH model virtual matrix simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override { // Generate random innovations (shocks) for the full period including burn-in vector errors = rng.rng(_nobs + burn); // Determine the initial long-run variance if not provided by the user if(initial_value == NULL) { // Calculate persistence = sum of ARCH + (0.5 * Leverage) + GARCH coefficients vector scale = vector::Ones(parameters.Size()); for(ulong i = m_p + 1; i < (m_p + m_o + 1); scale[i] = 0.5, ++i); vector temp = np::sliceVector(parameters, 1) * np::sliceVector(scale, 1); double persistence = temp.Sum(); // Unconditional Variance Formula: Omega / (1 - persistence) if((1. - persistence) > 0) initial_value = parameters[0] / (1. - persistence); else { Print(__FUNCTION__, " InitialValueWarning: Model parameters may imply non-stationarity"); initial_value = parameters[0]; } } vector sigma2, data, fsigma, fdata; sigma2 = data = fsigma = fdata = vector::Zeros(_nobs + burn); ulong max_lag = MathMax(MathMax(m_p, m_o), m_q); // Initialize pre-sample values to avoid cold-start bias np::fillVector(fsigma, initial_value, 0, long(max_lag)); double dv = pow(initial_value, 2. / m_power); np::fillVector(sigma2, dv, 0, long(max_lag)); for(ulong i = 0; i < max_lag; ++i) { data[i] = sqrt(sigma2[i]) * errors[i]; fdata[i] = pow(fabs(data[i]), m_power); } // Recursive simulation of the volatility process for(ulong t = max_lag; t < (_nobs + burn); ++t) { ulong loc = 0; fsigma[t] = parameters[loc]; // Intercept (Omega) loc += 1; // ARCH terms: Sum of lagged squared residuals for(ulong j = 0; j < m_p; ++j) { fsigma[t] += parameters[loc] * fdata[t - 1 - j]; loc += 1; } // Leverage (Asymmetry) terms: Gamma * I(res < 0) * res^2 for(ulong j = 0; j < m_o; ++j) { fsigma[t] += parameters[loc] * fdata[t - 1 - j] * ((data[t - 1 - j] < 0.) ? 1. : 0.0); loc += 1; } // GARCH (Persistence) terms: Beta * sigma^2 for(ulong j = 0; j < m_q; ++j) { fsigma[t] += parameters[loc] * fdata[t - 1 - j]; loc += 1; } // Transform internal power scale back to variance sigma2[t] = pow(fsigma[t], 2. / m_power); // Generate observation based on current volatility data[t] = errors[t] * sqrt(sigma2[t]); fdata[t] = pow(MathAbs(data[t]), m_power); } // Discard burn-in period to ensure the series has reached its stationary distribution data = np::sliceVector(data, long(burn)); sigma2 = np::sliceVector(sigma2, long(burn)); // Construct output matrix [data, variance] matrix out(data.Size(), 2); out.Col(data, 0); out.Col(sigma2, 1); return out; } //--- Get Garch specfication virtual vector garch_spec(void) override { vector specs(4); specs[0] = double(m_p); specs[1] = double(m_o); specs[2] = double(m_q); specs[3] = double(m_power); return specs; } }; //+------------------------------------------------------------------+ //| ARCH process | //| Implementation of the Autoregressive Conditional Heteroskedasticity | //| model, inheriting from the generalized GARCH framework. | //+------------------------------------------------------------------+ class CArchProcess: public CGarchProcess { public: // Constructor: Default initialization for a standard ARCH(1) model CArchProcess(void) { m_initialized = CGarchProcess::_initialize(VOL_ARCH, true, false, 0, "ARCH"); } // Constructor: Custom initialization for ARCH(p) model with seed and bootstrap settings CArchProcess(ulong p, int seed = 0, ulong nbootstraps = 100) { // ARCH models have no asymmetry (o=0) or GARCH (q=0) terms m_initialized = CGarchProcess::_initialize(VOL_ARCH, true, false, 0, "ARCH", seed, p, 0, 0, 2.0, 0, -1, nbootstraps); } // Destructor ~CArchProcess(void) { } //--- Returns starting values for the ARCH model optimization //--- Performs a targeted grid search to find the initial alpha parameters that maximize the log-likelihood virtual vector startingValues(vector& resids) override { // Create a grid of candidate ARCH persistence coefficients (alphas) to test vector alphas = np::arange(ulong(17), double(0.1), double(0.05)); vector svs[17]; vector bc = backCast(resids); vector llfs = alphas; matrix vb = varianceBounds(resids); vector sv, temp; // Iterate through candidates to find the best initial guess for(uint i = 0; i < svs.Size(); ++i) { // sv[0] is omega (constant), sv[1...p] are the ARCH coefficients // Initial omega is set based on the sample variance and assumed persistence sv = (1.0 - alphas[i]) * resids.Var() * vector::Ones(m_p + 1); // Evenly distribute the persistence coefficient (alpha) across all ARCH lags np::fillVector(sv, alphas[i] / double(m_p), 1); svs[i] = sv; // Evaluate log-likelihood for this parameter combination llfs[i] = _gaussianloglikelihood(sv, resids, bc, vb); } // Select the index that produced the highest log-likelihood ulong loc = llfs.ArgMax(); return svs[loc]; } }; //+------------------------------------------------------------------+ //| GJR-GARCH process | //| Models asymmetric volatility (leverage effect) using GARCH-style | //| recursion with quadratic variance (power=2.0). | //+------------------------------------------------------------------+ class CGjrGarchProcess: public CGarchProcess { public: // Default GJR-GARCH(1,1,1) constructor CGjrGarchProcess(void) { m_initialized = CGarchProcess::_initialize(VOL_GJR_GARCH, true, false, 0, "GJR-GARCH"); } // Custom GJR-GARCH(p,o,q) constructor CGjrGarchProcess(ulong p, ulong o, ulong q, int seed = 0, ulong nboots = 100) { m_initialized = CGarchProcess::_initialize(VOL_GJR_GARCH, true, false, 0, "GJR-GARCH", seed, p, o, q, 2.0, 0, -1, nboots); } }; //+------------------------------------------------------------------+ //| AVARCH process | //| Absolute Value ARCH process (power=1.0), modeling conditional | //| standard deviation instead of variance. | //+------------------------------------------------------------------+ class CAvarchProcess: public CGarchProcess { public: // Default AVARCH(1) constructor CAvarchProcess(void) { m_initialized = CGarchProcess::_initialize(VOL_AVARCH, true, false, 0, "AVARCH", 0, 1, 0, 0, 1.0); } // Custom AVARCH(p) constructor CAvarchProcess(ulong p, int seed = 0, ulong nboots = 100) { m_initialized = CGarchProcess::_initialize(VOL_AVARCH, true, false, 0, "AVARCH", seed, p, 0, 0, 1.0, 0, -1, nboots); } }; //+------------------------------------------------------------------+ //| AVGARCH process | //| Absolute Value GARCH process (power=1.0), extending AVARCH with | //| autoregressive variance components. | //+------------------------------------------------------------------+ class CAvgarchProcess: public CGarchProcess { public: // Default AVGARCH(1,1) constructor CAvgarchProcess(void) { m_initialized = CGarchProcess::_initialize(VOL_AVGARCH, true, false, 0, "AVGARCH", 0, 1, 0, 1, 1.0); } // Custom AVGARCH(p,q) constructor CAvgarchProcess(ulong p, ulong q, int seed = 0, ulong nboots = 100) { m_initialized = CGarchProcess::_initialize(VOL_AVGARCH, true, false, 0, "AVGARCH", seed, p, 0, q, 1.0, 0, -1, nboots); } }; //+------------------------------------------------------------------+ //| TARCH process | //| Threshold ARCH process (power=1.0), allowing for asymmetric | //| impacts on conditional standard deviation. | //+------------------------------------------------------------------+ class CTarchProcess: public CGarchProcess { public: // Default TARCH(1,1,1) constructor CTarchProcess(void) { m_initialized = CGarchProcess::_initialize(VOL_TARCH, true, false, 0, "TARCH", 0, 1, 1, 1, 1.0); } // Custom TARCH(p,o,q) constructor CTarchProcess(ulong p, ulong o, ulong q, int seed = 0, ulong nboots = 100) { m_initialized = CGarchProcess::_initialize(VOL_TARCH, true, false, 0, "TARCH", seed, p, o, q, 1.0, 0, -1, nboots); } }; //+------------------------------------------------------------------+ //| HARCH process | //| Heterogeneous Autoregressive Conditional Heteroskedasticity | //| Models volatility as a sum of lags over different frequencies | //+------------------------------------------------------------------+ class CHarchProcess: public CVolatilityProcess { protected: //--- Custom structure to hold Common Forecast Components (CFC) //--- This acts as a data container for efficient recursive volatility projection struct CFC { double Const; // Omega (long-term variance component) vector Arch; // Flattened ARCH coefficients vector matrix Resids2; // Matrix of squared residuals (historical + forecast buffer) // Default constructor initializing empty structures CFC(void) { Const = 0; Arch = vector::Zeros(0); Resids2 = matrix::Zeros(0, 0); } // Parameterized constructor for pre-populated components CFC(double c, vector& a, matrix& r) { Const = c; Arch = a; Resids2 = r; } // Copy constructor to handle structure duplication CFC(CFC& other) { Const = other.Const; Arch = other.Arch; Resids2 = other.Resids2; } // Assignment operator to facilitate component passing void operator=(CFC& other) { Const = other.Const; Arch = other.Arch; Resids2 = other.Resids2; } }; ulong m_num_lags; // Total number of frequency-based lag groups defined vector m_lags; // Vector containing the specific lag lengths for each group //--- Model transformation: Converts HARCH parameters (coefficients per frequency group) //--- into a standard ARCH(p) representation (one coefficient per individual lag) vector _harch_to_arch(vector& parameters) { // ARCH representation requires one parameter per lag up to the maximum lag defined vector arch_params = vector::Zeros(1 + int(m_lags.Max())); arch_params[0] = parameters[0]; // Copy the constant (omega) double lag = 0; double param = 0; // Distribute the HARCH parameter weight evenly across all lags within that group for(ulong i = 1; i < arch_params.Size(); ++i) { param = parameters[i]; lag = m_lags[i - 1]; // Uniform distribution of the impact: alpha_i = weight_group / length_group for(ulong j = 1; j < ulong(lag + 1.); arch_params[j] += param / lag, ++j); } return arch_params; } //--- Verify the requested forecasting method as valid for the specification //--- HARCH models support both analytic and simulation methods natively virtual bool _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override { return true; } //--- Intermediate forecast components: Prepares the history matrix for forecasting //--- Aligns squared residuals and backcast values to compute future conditional variance CFC _common_forecast_components(vector& parameters, vector& resids, vector& backcast, ulong horizon) { CFC out; // Transform HARCH parameters to an ARCH format for vectorization vector ap = _harch_to_arch(parameters); long t = long(resids.Size()); long m = long(m_lags.Max()); // Initialize the matrix with dimensions: Time x (MaxLag + Horizon) out.Resids2 = matrix::Zeros(t, m + horizon); // Seed the forecast horizon part of the matrix with initial backcast values for(long i = m; i < (long)out.Resids2.Rows(); ++i) for(long j = m; j < long(out.Resids2.Cols()); ++j) out.Resids2[i, j] = backcast[0]; // Populate historical squared residuals into the rolling buffer matrix vector sqresids = pow(resids, 2.); for(long i = 0; i < m; ++i) for(long j = m - i - 1, k = 0; j < long(out.Resids2.Rows()) && k < long(t - (m - i - 1)); ++k, ++j) out.Resids2[j, i] = sqresids[k]; // Store the constant and the ARCH terms in the CFC structure for use in _forecast out.Const = ap[0]; out.Arch = np::sliceVector(ap, 1); return out; } //--- Analytic multi-step volatility forecasts from the HARCH model //--- Computes expected future variance using the linear projection of squared residuals virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) override { // Retrieve the pre-calculated model constants and historical residual matrix CFC cfc = _common_forecast_components(parameters, resids, backcast, horizon); long m = long(m_lags.Max()); // Trim history to the requested start point cfc.Resids2 = np::sliceMatrixRows(cfc.Resids2, _start); // Reverse the ARCH coefficient vector to align with the rolling residual buffer (t-1, t-2, ...) vector arch_rev = cfc.Arch; if(!np::reverseVector(arch_rev)) return VarianceForecast(); // Iterative projection of future variance components for(long i = 0; i < long(horizon); ++i) { // Extract the relevant lag window (t-m to t-1) matrix cc = np::sliceMatrixCols(cfc.Resids2, i, m + i); // Forecast: Omega + Sum(Alpha_j * Resids^2_{t-j}) if(!cfc.Resids2.Col(cfc.Const + cc.MatMul(arch_rev), m + i)) { Print(__FUNCTION__, ": Failed Columns assignment cfc.Resids2 ", GetLastError()); return VarianceForecast(); } } // Extract the projected variance paths (excluding initial lags) matrix out = np::sliceMatrixCols(cfc.Resids2, m); matrix empty[]; return VarianceForecast(out, empty, empty); } //--- Simulation-based volatility forecasts from the HARCH model //--- Projects future variance paths using Monte Carlo integration virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) override { CFC cfc = _common_forecast_components(parameters, resids, backcast, horizon); long t = long(resids.Size()); long m = long(m_lags.Max()); matrix paths[], shocks[]; matrix temp_resids2 = matrix::Zeros(simulations, m + horizon); ArrayResize(paths, int(t - _start)); ArrayResize(shocks, int(t - _start)); for(uint i = 0; i < paths.Size(); paths[i] = shocks[i] = matrix::Zeros(simulations, horizon), ++i); // Reverse ARCH coefficients for matrix multiplication alignment vector arch_rev = np::sliceVector(cfc.Arch, BEGIN_REVERSE, END_REVERSE, -1); // Run simulations for each historical point for(long i = _start; i < t; ++i) { matrix stdshocks = rng.rng(simulations, horizon); matrix rd2 = np::sliceMatrixRows(cfc.Resids2, i, i + 1); temp_resids2 = rd2; long path_loc = i - _start; // Perform recursive step-ahead simulation for(long j = 0; j < long(horizon); ++j) { rd2 = np::sliceMatrixCols(temp_resids2, j, m + j); // Update variance: Forecast = Omega + Sum(Alpha * PastSqResids) if(!paths[path_loc].Col(cfc.Const + rd2.MatMul(arch_rev), j)) { Print(__FUNCTION__, ": Column assignment error paths", GetLastError()); return VarianceForecast(); } // Generate stochastic innovation based on the forecasted variance if(!shocks[path_loc].Col(stdshocks.Col(j) * sqrt(paths[path_loc].Col(j)), j)) { Print(__FUNCTION__, ": Column assignment error shocks", GetLastError()); return VarianceForecast(); } // Feed back the squared innovation into the residual buffer for next-step recursion if(!temp_resids2.Col(pow(shocks[path_loc].Col(j), 2.), m + j)) { Print(__FUNCTION__, ": Column assignment error temp_resids2", GetLastError()); return VarianceForecast(); } } } // Aggregate simulation results to find the mean volatility path matrix out(paths.Size(), paths[0].Cols()); for(uint i = 0; i < paths.Size(); ++i) if(!out.Row(paths[i].Mean(1), i)) { Print(__FUNCTION__, ": Row assignment error out", GetLastError()); return VarianceForecast(); } return VarianceForecast(out, paths, shocks); } //--- Initialize class members and configure HARCH model structure //--- HARCH models volatility as a weighted sum of returns at different frequency horizons virtual bool _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100) override { m_updateable = updateable; m_closedform = closedform; bootstraps(bootstrap_obs); m_start = start; m_stop = stop; m_name = name; m_seed = seed; m_model = volmodel; // Ensure the lag vector is populated; defaults to a single lag (standard ARCH(1)) if empty if(!m_lags.Size()) { m_lags = vector::Ones(1); for(ulong i = 0; i < m_lags.Size(); m_lags[i] = 1. + double(i), ++i); } // Validate that all defined lags represent positive integers if(int(m_lags.Max()) < 1) { Print(__FUNCTION__, ": Invalid input. All lags should resolve to > 0 integers "); return false; } // Set internal counters for lag depth and total parameters (Constant + 1 per lag group) m_num_lags = m_lags.Size(); m_num_params = m_num_lags + 1; // Initialize the random number generator for simulations return m_normal.initialize(vector::Zeros(0), seed); } public: // Constructor: Default initialization for a basic HARCH process CHarchProcess(void) : m_lags(vector::Ones(1)) { m_initialized = _initialize(VOL_HARCH, true, false, 0, "HARCH"); } // Constructor: Initializing with a custom vector of lag frequency groups CHarchProcess(vector& lags) { // If a single value is provided, treat it as the depth for a uniform lag range if(lags.Size() == 1) { m_lags = vector::Ones(ulong(lags[0])); for(ulong i = 0; i < m_lags.Size(); m_lags[i] = lags[0] + double(i), ++i); } else m_lags = lags; // Use the explicitly provided set of lags m_initialized = _initialize(VOL_HARCH, true, false, 0, "HARCH"); } //--- Returns bounds for parameters //--- Sets search intervals [min, max] for the HARCH parameters virtual matrix bounds(vector& resids) override { // Bounds matrix: rows = parameters (1 omega + number of lag groups), 2 columns (min, max) matrix out = matrix::Zeros(m_lags.Size() + 1, 2); // Estimate sample variance to provide a logical upper bound for the constant (omega) vector squared = pow(resids, 2.); out[0, 1] = 10. * squared.Mean(); // Set upper bounds for all ARCH coefficients (alpha_i) to 1.0 (enforcing local stability) for(ulong i = 1; i < out.Rows(); ++i) out[i, 1] = 1.0; // Transpose to match [2 x N] structure required by the optimization algorithm return out.Transpose(); } //--- Construct parameter constraints arrays for parameter estimation //--- Enforces positivity (omega, alpha > 0) and stationarity: sum(alpha_i) < 1 virtual Constraints constraints(void) override { // a[m_num_lags+1, m_num_lags+1] matrix for linear inequalities: A * params <= b matrix a = matrix::Zeros(m_num_lags + 2, m_num_lags + 1); // Enforce positivity: omega >= 0 and all alpha_i >= 0 // Implemented as -omega <= 0 and -alpha_i <= 0, or identity matrix for lower bounds for(ulong i = 0; i < m_num_lags + 1; a[i, i] = 1., ++i); // Enforce stationarity: sum(alpha_i) < 1 // Linear form: - (alpha_1 + alpha_2 + ... + alpha_k) >= -1 => (alpha_1 + ... + alpha_k) <= 1 for(ulong j = 1; j < a.Cols(); ++j) a[m_num_lags + 1, j] = -1.; vector b = vector::Zeros(m_num_lags + 2); b[m_num_lags + 1] = -1.; // Constraint threshold Constraints out; out._one = a; out._two = b; return out; } //--- Compute the conditional variance for the HARCH model //--- Filters the residuals through the heterogeneous lag structure to produce a variance series virtual vector computeVariance(vector& parameters, vector& resids, vector& sigma2, vector& backcast, matrix& varbounds) override { vector lags = m_lags; int nobs = int(resids.Size()); // Perform recursive filtering based on frequency-specific lag groups harch_recursion(parameters, resids, sigma2, lags, nobs, backcast[0], varbounds); return sigma2; } //--- Simulate synthetic data from the HARCH model //--- Generates a series of returns and their corresponding conditional variances virtual matrix simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override { // Generate random innovations (shocks) for the total period (including burn-in) vector errs = rng.rng(_nobs + burn); if(!errs.Size()) { matrix container = rng.rng(_nobs + burn, 1); errs = container.Col(0); } // Calculate unconditional long-run variance to seed the process if not provided if(!initial_value) { // Variance = Omega / (1 - sum of all ARCH parameters) vector vp = np::sliceVector(parameters, 1); double psum = vp.Sum(); if(1.0 - psum > 0) initial_value = parameters[0] / (1.0 - psum); else { Print(__FUNCTION__, ": Initial value WARNING: Parameters imply non-stationarity"); initial_value = parameters[0]; } } vector sigma2 = vector::Zeros(_nobs + burn); vector data = sigma2; long maxlag = long(m_lags.Max()); // Initialize pre-sample buffer to prevent cold-start effects for(long i = 0; i < maxlag; ++i) { sigma2[i] = initial_value; data[i] = sqrt(initial_value); } // Recursive simulation of the HARCH(p) process double param = 0; for(long t = maxlag; t < long(_nobs + burn); ++t) { sigma2[t] = parameters[0]; // Intercept term (Omega) // Aggregate volatility contributions from different frequency lag groups for(ulong i = 0; i < m_lags.Size(); ++i) { // Distribute group coefficient across the specified lag window length param = parameters[1 + i] / double(m_lags[i]); for(ulong j = 0; j < ulong(m_lags[i]); ++j) sigma2[t] += param * pow(data[t - 1 - j], 2.); } // Update observation based on conditional variance data[t] = errs[t] * sqrt(sigma2[t]); } // Remove burn-in period to ensure stationarity data = np::sliceVector(data, long(burn)); sigma2 = np::sliceVector(sigma2, long(burn)); // Return resulting data as a [returns, variance] matrix matrix out(data.Size(), 2); out.Col(data, 0); out.Col(sigma2, 1); return out; } //--- Returns starting values for the HARCH model estimation //--- Provides an initial guess based on assumed persistence to help the optimizer converge virtual vector startingValues(vector& resids) override { // Assume 90% persistence as a standard starting heuristic for long-memory models double alpha = 0.9; // Initialize omega (constant) based on sample variance, and split alpha persistence equally across all lag groups vector sv = (1. - alpha) * resids.Var() * vector::Ones(m_num_lags + 1); // Distribute the alpha weight evenly across the specified number of lag frequency groups for(ulong i = 1; i < sv.Size(); sv[i] = alpha / double(m_num_lags), ++i); return sv; } //--- Names of model parameters for logging and identification //--- Returns a comma-separated string listing omega followed by each frequency group label virtual string parameterNames(void) override { string names = "omega"; for(ulong i = 0; i < m_num_lags; ++i) // Format: alpha[lag_length] StringConcatenate(names, ",", StringFormat("%s%d%s", "alpha[", int(m_lags[i]), "]")); return names; } }; //+------------------------------------------------------------------+ //| FIGARCH process | //| Fractionally Integrated GARCH: Models long-memory volatility | //| behavior using fractional differencing of the variance process. | //+------------------------------------------------------------------+ class CFiGarchProcess: public CVolatilityProcess { protected: ulong m_p; // Order of the ARCH component ulong m_q; // Order of the GARCH component double m_power; // Power transformation (e.g., 2.0 for variance) ulong m_truncation; // Truncation lag for the infinite MA representation //--- Dynamically generate the model name string based on power and structure string _genetrateName(void) { // Standard variance-based models if(m_power == 2.) { if(m_q == 0) return "FIARCH"; else return "FIGARCH"; } // Absolute value-based models (e.g., power=1) else if(m_power == 1.) { if(m_q == 0) return "FIAVARCH"; else return "FIAVGARCH"; } // Generalized Power models else { if(m_q == 0) return StringFormat("Power FIARCH (power: {%.1f})", m_power); else return StringFormat("Power FIGARCH (power: {%.1f})", m_power); } } //--- Verify if the requested forecasting method is mathematically valid //--- Analytic methods require specific power=2.0 structures for multi-step persistence virtual bool _check_forecasting_method(ENUM_FORECAST_METHOD method, ulong horizon) override { // Single-step forecasts are always valid if(horizon == 1) return true; // Analytic multi-step forecasts are only supported for standard variance (power=2) if(method == FORECAST_ANALYTIC && m_power != 2.) { Print(__FUNCTION__, ": Analytic forecasts not available for horizon > 1 when power != 2"); return false; } return true; } //--- Analytic multi-step volatility forecasts from the model virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) { vector sg; matrix fc; //--- Establish the initial one-step look-ahead matrix boundaries to seed the multi-step recursion if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sg, fc)) return VarianceForecast(); //--- Short-circuit execution if only dealing with a single look-ahead window if(horizon == 1) return VarianceForecast(fc, EMPTY_MATRIX_ARRAY, EMPTY_MATRIX_ARRAY); //--- Extract model-specific parameters (skipping omega at parameters[0]) vector params = np::sliceVector(parameters, 1); //--- Pre-calculate the expansion weights up to the allocated truncation lag vector lam = figarch_weights(params, m_p, m_q, m_truncation); vector lam_rev = lam; //--- Reverse the weights vector to align chronological steps for a vector dot product (convolution) if(!np::reverseVector(lam)) return VarianceForecast(); long t = long(resids.Size()); double omega = parameters[0]; double beta = (m_q) ? parameters[parameters.Size() - 1] : 0.0; double omega_tilde = omega / (1. - beta); //--- Long-run intercept scaling //--- Allocate a unified workspace vector to hold past realized shocks and future forecast values sequentially vector temp_forecasts = vector::Zeros(m_truncation + horizon); vector resids2 = pow(resids, 2.0); //--- Pre-square raw structural innovations //--- Loop across the time-series sliding calculation window starting at '_start' for(long i = _start; i < t; ++i) { long fcast_loc = i - _start; long available = long(i + 1 - fmax(0, i - m_truncation + 1)); //--- Load historical squared innovations into the sliding workspace vector for(long j = long(m_truncation - available), k = (long)fmax(0, i - m_truncation + 1); j < long(m_truncation) && k < (i + 1); ++k, ++j) temp_forecasts[j] = resids2[k]; //--- Hand off backcast estimates if historical data depth does not span the truncation requirements if(available < long(m_truncation)) if(!np::fillVector(temp_forecasts, backcast[0], 0, long(m_truncation) - available)) return VarianceForecast(); //--- --- Multi-Step Analytic Loop --- //--- Step sequentially into the future horizon, resolving expectations via a convolution of reversed weights for(ulong h = 0; h < horizon; ++h) { vector lagged_forecasts = np::sliceVector(temp_forecasts, long(h), long(m_truncation + h)); temp_forecasts[m_truncation + h] = omega_tilde + lam_rev.Dot(lagged_forecasts); } //--- Transfer the computed future horizons into the destination forecast result matrix row if(!fc.Row(np::sliceVector(temp_forecasts, long(m_truncation)), fcast_loc)) return VarianceForecast(); } return VarianceForecast(fc, EMPTY_MATRIX_ARRAY, EMPTY_MATRIX_ARRAY); } //--- Simulation-based volatility forecasts from the model virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) { vector sig; matrix fc; //--- Baseline one-step layout verification if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sig, fc)) return VarianceForecast(); ulong t = resids.Size(); matrix paths[], shocks[]; ArrayResize(paths, int(t - _start)); ArrayResize(shocks, int(t - _start)); vector params = np::sliceVector(parameters, 1); vector lam = figarch_weights(params, m_p, m_q, m_truncation); vector lam_rev = lam; if(!np::reverseVector(lam_rev)) return VarianceForecast(); double omega = parameters[0]; double beta = (m_q) ? parameters[parameters.Size() - 1] : 0.0; double omega_tilde = omega / (1. - beta); //--- Allocate parallel simulation matrices [simulations rows x (truncation + horizon) columns] matrix fpath = matrix::Zeros(simulations, m_truncation + horizon); vector fresids = pow(fabs(resids), 2.); //--- Loop chronologically across the forecasting evaluation footprint for(long i = _start; i < long(t); ++i) { //--- Extract structured random standard normal variables from the random number generator matrix std_shocks = rng.rng(simulations, horizon); long available = i + 1 - fmax(0, i - long(m_truncation) + 1); vector frd = np::sliceVector(fresids, fmax(0, i + 1 - m_truncation), i + 1); //--- Inject identical historical realized shocks across all path simulations for(long ii = long(m_truncation - available); ii < long(m_truncation); ++ii) fpath.Col(frd, ii); //--- Pad empty history with backcast values where data boundary overflows if(available < long(m_truncation)) for(long ii = 0; ii < long(m_truncation - available); ++ii) fpath.Col(backcast, ii); //--- --- Horizon Path Integration Simulation --- for(ulong h = 0; h < horizon; ++h) { //--- Matrix column slice retrieves conditional tracking variables across all simulation paths simultaneously matrix lagged_forecasts = np::sliceMatrixCols(fpath, long(h), long(m_truncation + h)); vector temp = omega_tilde + lagged_forecasts.MatMul(lam_rev); //--- Invert the structural transformation power variable back into basic variance space vector sigma2 = pow(temp, 2. / m_power); long path_loc = i - _start; //--- Generate dynamic synthesized shocks and variance values for this horizon index shocks[path_loc].Col(std_shocks.Col(h) * sqrt(sigma2), h); paths[path_loc].Col(sigma2, h); //--- Expected forecast point equals the mathematical mean across all simulated iterations fc[path_loc, h] = sigma2.Mean(); //--- Feed back the generated structural shock into the path matrix for the next recursive horizon step fpath.Col(pow(fabs(shocks[path_loc].Col(h)), m_power), m_truncation + h); } } return VarianceForecast(fc, paths, shocks); } //--- Initialize class members virtual bool _initialize(ENUM_VOLATILITY_MODEL volmodel = WRONG_VALUE, bool updateable = true, bool closedform = false, ulong nparams = 0, string name = NULL, int seed = 0, ulong p = 0, ulong o = 0, ulong q = 0, double power = 0.0, long start = 0, long stop = -1, ulong bootstrap_obs = 100) override { //--- Synchronize configuration parameters to localized private class properties m_updateable = updateable; m_closedform = closedform; bootstraps(bootstrap_obs); m_start = start; m_stop = stop; m_name = name; m_seed = seed; m_model = volmodel; m_p = p; m_q = q; m_truncation = o; //--- Maps the truncation lag dimension m_power = power; //--- Declares structural model exponent (e.g. 2.0 for FIGARCH) m_num_params = 2 + m_p + m_q; //--- Core total parameter footprint (omega + phi + beta + d) //--- Initialize the distribution object wrapper return m_normal.initialize(vector::Zeros(0), seed); } public: // Constructor: Default initialization for a standard FIGARCH(1,d,1) model // Uses default truncation of 1000 lags to approximate fractional integration CFiGarchProcess(void) { m_initialized = _initialize(VOL_FIGARCH, true, false, 4, "FIGARCH", 0, 1, 1000, 1, 2.0); } // Constructor: Custom initialization for flexible FIGARCH(p,d,q) configurations // Allows user to define model orders, power transformation, and truncation depth CFiGarchProcess(ulong p, ulong q, double power, ulong truncation, int seed, long start, long stop, ulong boot) { // Passes parameters to the internal _initialize method // Note: 'truncation' maps to the 'o' (asymmetry/leverage-style) parameter slot in the base process m_initialized = _initialize(VOL_FIGARCH, true, false, 0, "FIGARCH", seed, p, truncation, q, power, start, stop, boot); } //--- Returns starting values for the ARCH model using a grid-search heuristic initialization virtual vector startingValues(vector& resids) override { //--- Establish parameter boundary array configurations for the grid search optimization initialization vector ds = {0.2, 0.5, 0.7}; vector phi_ratio, beta_ratio; if(m_p) { phi_ratio.Resize(3); phi_ratio[0] = 0.2; phi_ratio[1] = 0.5; phi_ratio[2] = 0.8; } else { phi_ratio.Resize(1); phi_ratio[0] = 0; } if(m_q) { beta_ratio.Resize(3); beta_ratio[0] = 0.1; beta_ratio[1] = 0.5; beta_ratio[2] = 0.9; } else { beta_ratio.Resize(1); beta_ratio[0] = 0; } //--- Establish baseline variance expectations derived directly from residuals footprint double target = pow(fabs(resids), m_power).Mean(); double scale = pow(resids, 2.).Mean() / (pow(target, 2.0 / m_power)); target *= pow(scale, m_power / 2.); //--- Map space to permute all parameter combinations matrix all_starting; all_starting.Resize(ds.Size() * phi_ratio.Size() * beta_ratio.Size(), 4); double beta, omega, d, phi, br, pr; vector temp = vector::Zeros(3); vector lam; ulong row = 0; //--- --- Grid Execution Matrix Search Loop --- for(ulong i = 0; i < ds.Size(); ++i) { d = ds[i]; for(ulong j = 0; j < phi_ratio.Size(); ++j) { pr = phi_ratio[j]; phi = (1. - d) / 2. * pr; //--- Force stationary constraints for(ulong k = 0; k < beta_ratio.Size(); ++k) { br = beta_ratio[k]; beta = (d + phi) * br; temp[0] = phi; temp[1] = d; temp[2] = beta; lam = figarch_weights(temp, 1, 1, int(m_truncation)); //--- Calculate corresponding consistent omega constant for the localized parameter combination omega = (1. - beta) * target * (1. - lam.Sum()); all_starting[row, 0] = omega; all_starting[row, 1] = phi; all_starting[row, 2] = d; all_starting[row++, 3] = beta; } } } //--- Extract unique configurations to minimize processing waste matrix sv = np::unique(all_starting, 1); //--- Strip unused parameter tracking columns if lag settings are set to zero (p=0 or q=0) if(!m_q) sv = np::sliceMatrixCols(sv, 0, long(sv.Cols() - 1)); if(!m_p) { long cols[]; ArrayResize(cols, int(sv.Cols() - 1)); cols[0] = 0; for(long i = 1; i < long(cols.Size()); cols[i] = i + 1, ++i); sv = np::selectMatrixCols(sv, cols); } matrix vb = varianceBounds(resids); vector bc = backCast(resids); vector llfs = vector::Zeros(sv.Rows()); //--- Score every starting candidate against the target log-likelihood function for(ulong i = 0; i < sv.Rows(); ++i) { vector row = sv.Row(i); llfs[i] = _gaussianloglikelihood(row, resids, bc, vb); } //--- Return the parameter configuration vector that optimized/maximized the log-likelihood (ArgMax) return sv.Row(llfs.ArgMax()); } //--- Construct values for backcasting to start the recursion virtual vector backCast(vector& resids) override { //--- Limit calculation depth up to a safe window threshold (e.g. 75 observations) ulong tau = MathMin(75, resids.Size()); vector w = vector::Zeros(tau); //--- Generate an exponentially decaying weights structure (RiskMetrics style, lambda = 0.94) for(ulong i = 0; i < w.Size(); w[i] = pow(0.94, double(i)), ++i); w = w / w.Sum(); //--- Normalize weight values to sum to 1.0 //--- Extract the start of the series, apply transformation power, and scale by decay weights vector sres = np::sliceVector(resids, 0, long(tau)); sres = pow(fabs(sres), m_power) * w; vector bc(1); bc[0] = sres.Sum(); //--- Aggregate into scalar backcast representation return bc; } //--- Transformation to apply to user-provided backcast values virtual vector backCastTransform(vector& backcast) override { //--- Defer to core base class method transformations first backcast = CVolatilityProcess::backCastTransform(backcast); //--- Apply specific fractional variance power adjustments to the processed vector values vector _backcast = pow(sqrt(backcast), m_power); return _backcast; } //--- Construct loose bounds for conditional variances virtual matrix varianceBounds(vector& resids, double power = 2.) override { //--- Pass computation parameters over to internal system tracker return CVolatilityProcess::varianceBounds(resids, m_power); } //--- Returns boundary matrix structures for parameters optimization constraints virtual matrix bounds(vector& resids) override { double eps_half = sqrt(2.220446049250313e-16); //--- System machine precision constant limit double v = (pow(fabs(resids), m_power)).Mean(); matrix bnds = matrix::Zeros(m_p + m_q + 2, 2); bnds[0, 1] = 10. * v; //--- Enforce an upper bound ceiling on baseline intercept omega constant for(ulong i = 1; i < m_p + 1; bnds[i, 1] = 0.5, ++i); //--- Limit short term memory parameters for(ulong i = m_p + 1; i < bnds.Rows(); bnds[i, 1] = 1. - eps_half, ++i); //--- Bound persistence variables below 1.0 return bnds.Transpose(); } //--- Compute the variance for the ARCH model virtual vector computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) { //--- Transform base residuals into process space based on model's power coefficient vector fresids = pow(fabs(resids), m_power); int nobs = (int)resids.Size(); //--- Execute the fractional integration recursion engine to write conditional variance calculations directly into _sigma2 figarch_recursion(parameters, fresids, _sigma2, int(m_p), int(m_q), nobs, m_truncation, _backcast[0], varbounds); //--- Apply inverse power mapping to re-align output to standard variance units double inv_power = 2. / m_power; _sigma2 = pow(_sigma2, inv_power); return _sigma2; } //--- Construct parameter constraints arrays for parameter estimation (Boundary bounding vectors) virtual Constraints constraints(void) override { //--- Initialize hard-coded inequality constraint matrices (A) for OLS/ML parameter bounding matrix a = { {1, 0, 0, 0}, {0, 1, 0, 0}, {0, -2, -1, 0}, {0, 0, 1, 0}, {0, 0, -1, 0}, {0, 0, 0, 1}, {0, 1, 1, -1}, }; matrix aa = { {1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 1, -1}, }; matrix aaa = { {1, 0, 0}, {0, 1, 0}, {0, -2, -1}, {0, 0, 1}, {0, 0, -1} }; //--- Associated boundary mapping verification vectors (B) where: A * parameters >= B vector b = {0, 0, -1, 0, -1, 0, 0}; vector bb = {0, 0, -1, 0, 0}; vector bbb = {0, 0, -1, 0, -1}; Constraints constr; constr._one = a; constr._two = b; //--- Adjust structural linear constraint layout mapping depending on specified lag dimension flags (p=0 or q=0) if(!m_q) { constr._one = aaa; constr._two = bbb; } if(!m_p) { constr._one = aa; constr._two = bb; } return constr; } //--- Simulate data paths from the model structure virtual matrix simulate(vector& parameters, ulong _nobs, BootstrapRng &rng, ulong burn = 500, double initial_value = NULL) override { vector params = np::sliceVector(parameters, 1); //--- Extract long memory hyperbolic expansion decay parameter array vector lam = figarch_weights(params, m_p, m_q, m_truncation); vector lam_rev = lam; if(!np::reverseVector(lam_rev)) return matrix::Zeros(0, 0); //--- Pre-allocate standardized random error distribution array including the requested burn-in footprint vector errs = rng.rng(ulong(m_truncation + _nobs + burn)); //--- Calculate stationary unconditional expectation values if no initial value parameter was provided if(initial_value == NULL) { double persistence = lam.Sum(); double beta = (m_q) ? parameters[parameters.Size() - 1] : 0.0; initial_value = parameters[0]; if(beta < 1.) initial_value /= 1. - beta; if(persistence < 1.) initial_value /= 1. - persistence; if(persistence >= 1. || beta >= 1.) Print(__FUNCTION__": WARNING: Initial value warning (non-stationary bounds detected)"); } if(initial_value == NULL) { Print(__FUNCTION__, ": Initial value is NULL"); return matrix::Zeros(0, 0); } //--- Allocate unified arrays tracking simulated states over time vector sigma2 = vector::Zeros(ulong(m_truncation + _nobs + burn)); vector data = sigma2; vector fsigma = data; vector fdata = fsigma; //--- Initialize index blocks spanning pre-sample truncation window constraints if(!np::vectorFill(fsigma, initial_value, 0, long(m_truncation)) || !np::vectorFill(sigma2, pow(initial_value, 2. / m_power), 0, long(m_truncation))) { Print(__FUNCTION__, ": vector fill error"); return matrix::Zeros(0, 0); } data = sqrt(sigma2) * errs; fdata = pow(fabs(data), m_power); double omega = parameters[0]; double beta = (m_q) ? parameters[parameters.Size() - 1] : 0; double omega_tilde = 0; if(beta < 1.) omega_tilde = omega / (1. - beta); else { Print(__FUNCTION__, ": Beta >= 1.0, using omega as intercept since long-run variance is ill-defined."); omega_tilde = omega; } //--- --- Simulation Time March Loop --- //--- Recursively project long-memory volatility evolution paths forward chronologically for(long t = long(m_truncation); t < long(m_truncation + _nobs + burn); ++t) { //--- Convolve the reversed weights against the sliding window vector slice of previous absolute values fsigma[t] = omega_tilde + lam_rev.Dot(np::sliceVector(fdata, long(t - m_truncation), t)); sigma2[t] = pow(fsigma[t], 2.0 / m_power); //--- Convert out of process shape parameters data[t] = errs[t] * sqrt(sigma2[t]); //--- Generate synthetic asset return realization fdata[t] = pow(fabs(data[t]), m_power); //--- Cache internal tracking absolute shape value } //--- Slice out and discard the initial initialization burn-in sequence block data = np::sliceVector(data, long(m_truncation + burn)); sigma2 = np::sliceVector(sigma2, long(m_truncation + burn)); //--- Package the output data streams into a 2-column matrix structure [Returns, Variance] matrix out = matrix::Zeros(data.Size(), 2); if(!out.Col(data, 0) || !out.Col(sigma2, 1)) { Print(__FUNCTION__, ": Column insertion error out ", GetLastError()); return matrix::Zeros(0, 0); } return out; } //--- Names of model parameters virtual string parameterNames(void) override { string names = "omega"; if(m_p) StringConcatenate(names,",","phi"); if(m_q) StringConcatenate(names,",","beta"); return names; } }; //+------------------------------------------------------------------+ //| EGARCH process | //| Models volatility as a sum of lags over different frequencies | //+------------------------------------------------------------------+ class CEgarchProcess: public CVolatilityProcess { protected: ulong m_p; // Order of the ARCH component ulong m_o; // Order of Assymetry component ulong m_q; // Order of the GARCH component double m_power; // Power transformation vector m_vectors[]; // Local container matrix _product(vector& one, vector& two, vector& three) { ulong rows = 3; matrix out(one.Size()*two.Size()*three.Size(),3); ulong shift = 0; for(ulong i = 0; i 1) { Print(__FUNCTION__, ": Analytic forecasts not available for horizon > 1 "); return false; } return true; } //--- Analytic multi-step volatility forecasts from the model virtual VarianceForecast _analyticforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon) { vector sg; matrix fc; //--- Establish the initial one-step look-ahead matrix boundaries to seed the multi-step recursion if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sg, fc)) return VarianceForecast(); return VarianceForecast(fc, EMPTY_MATRIX_ARRAY, EMPTY_MATRIX_ARRAY); } //--- Simulation-based volatility forecasts from the model virtual VarianceForecast _simulationforecast(vector& parameters, vector &resids, vector &backcast, matrix &varbounds, long _start, ulong horizon, ulong simulations, BootstrapRng &rng) { vector sig; matrix fc; //--- Baseline one-step layout verification if(!_onestepforecast(parameters, resids, backcast, varbounds, _start, horizon, sig, fc)) return VarianceForecast(); ulong t = resids.Size(); ulong m = fmax(fmax(m_p,m_o),m_q); vector Insigma2 = log(sig); vector e = resids/sqrt(sig); matrix Insigma2_mat = matrix::Zeros(t,m); if(backcast.Size()==m) { for(ulong i = 0; i0)?"EGARCH":"EARCH"; //--- Initialize the distribution object wrapper return m_normal.initialize(vector::Zeros(0), seed); } public: CEgarchProcess(void) { m_initialized = _initialize(VOL_EGARCH,true,false,0,NULL,0,1,0,1,2.0,0,-1,100); } CEgarchProcess(ulong p, ulong o, ulong q,double power, int seed, long start, long stop, ulong boot) { m_initialized = _initialize(VOL_EGARCH,true,false,0,NULL,seed,p,o,q,2.0,start,stop,boot); } //--- Returns starting values for the ARCH model using a grid-search heuristic initialization virtual vector startingValues(vector& resids) override { vector alphas = {0.01, 0.05, 0.1, 0.2}; vector gammas = {-0.1, 0.0, 0.1}; vector betas = {0.5, 0.7, 0.9, 0.98}; matrix agbs = _product(alphas,gammas,betas); double target = log(pow(resids,2).Mean()); vector svs[]; ArrayResize(svs, int(agbs.Rows())); matrix vb = varianceBounds(resids); vector llfs = vector::Zeros(agbs.Rows()); vector bc = backCast(resids); // Evaluate Log-Likelihood for each combination in the grid for(uint i = 0; i < svs.Size(); ++i) { vector row = agbs.Row(i); // Construct starting vector: constant omega + lags vector sv = (1. - row[2]) * target * vector::Ones(m_p + m_o + m_q + 1); // Distribute coefficients across ARCH (p), Leverage (o), and GARCH (q) components if(m_p > 0) np::fillVector(sv, row[0] / double(m_p), long(1), long(1 + m_p)); if(m_o > 0) np::fillVector(sv, row[1] / double(m_o), long(1 + m_p), long(1 + m_p + m_o)); if(m_q > 0) np::fillVector(sv, row[2] / double(m_q), long(1 + m_p + m_o), long(1 + m_p + m_o + m_q)); svs[i] = sv; llfs[i] = _gaussianloglikelihood(sv, resids, bc, vb); } // Select the parameter set that yielded the highest Log-Likelihood ulong loc = llfs.ArgMax(); return svs[loc]; } //--- Construct values for backcasting to start the recursion virtual vector backCast(vector& resids) override { return log(CVolatilityProcess::backCast(resids)); } //--- Transformation to apply to user-provided backcast values virtual vector backCastTransform(vector& backcast) override { //--- Defer to core base class method transformations first backcast = CVolatilityProcess::backCastTransform(backcast); return log(backcast); } //--- Construct loose bounds for conditional variances virtual matrix varianceBounds(vector& resids, double power = 2.) override { //--- Pass computation parameters over to internal system tracker return CVolatilityProcess::varianceBounds(resids, m_power); } //--- Returns boundary matrix structures for parameters optimization constraints virtual matrix bounds(vector& resids) override { double log_const = log(10000.0); double v = (pow(fabs(resids), m_power)).Mean(); matrix bnds = matrix::Zeros(1 + m_p + m_o + m_q, 2); bnds[0, 0] = log(v) - log_const; bnds[0, 1] = log(v) + log_const; for(ulong i = 1; i < m_p + m_o + 1; bnds[i, 0] = -1*double("inf"), bnds[i,1] = double("inf"), ++i); //--- Limit short term memory parameters for(ulong i = m_p + m_o + 1; i < bnds.Rows(); bnds[i, 1] = double(m_q), ++i); return bnds.Transpose(); } //--- Compute the variance for the ARCH model virtual vector computeVariance(vector& parameters, vector& resids, vector& _sigma2, vector& _backcast, matrix& varbounds) { //--- Transform base residuals into process space based on model's power coefficient //vector lnsigma2, std_resids, abs_std_resids; int nobs = (int)resids.Size(); if(!m_vectors.Size() || m_vectors[0].Size() != resids.Size()) { ArrayResize(m_vectors,3); m_vectors[0] = m_vectors[1] = m_vectors[2] = vector::Zeros(nobs); } egarch_recursion(parameters,resids,_sigma2,m_p,m_o,m_q,ulong(nobs),_backcast[0],varbounds,m_vectors[0],m_vectors[1],m_vectors[2]); return _sigma2; } //--- Construct parameter constraints arrays for parameter estimation (Boundary bounding vectors) virtual Constraints constraints(void) override { matrix a = matrix::Zeros(1,m_p+m_o+m_q+1); vector b = vector::Zeros(1); b[0] = -1.0; for(ulong i = m_p+m_o+1; i= 1.0, using omega as intercept since long-run variance is ill-defined."); omega_tilde = omega; } //--- --- Simulation Time March Loop --- ulong loc = 0; for(long t = long(max_lag); t < long(_nobs + burn); ++t) { loc = 0; lnsigma2[t] = parameters[loc]; loc+=1; for(ulong j = 0; j