509 lines
17 KiB
MQL5
509 lines
17 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| recursions.mqh |
|
|
//| Copyright 2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2025, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include"base.mqh"
|
|
//---
|
|
#define SQRT2_OV_PI 0.79788456080286541 // Constant used for Gaussian EGARCH scaling
|
|
const double LNSIGMA_MAX = log(DBL_MAX) - 0.1; // Overflow protection for log-variance
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Helper: Adjusts sigma2 to ensure it stays within defined bounds |
|
|
//| Prevents variance from collapsing to zero or exploding to inf. |
|
|
//+------------------------------------------------------------------+
|
|
double bounds_check(double _sigma2, vector &varbounds)
|
|
{
|
|
// Enforce floor boundary
|
|
if(_sigma2 < varbounds[0])
|
|
_sigma2 = varbounds[0];
|
|
else
|
|
// Enforce ceiling boundary with soft-clip for numerical stability
|
|
if(_sigma2 > varbounds[1])
|
|
{
|
|
if(MathClassify(_sigma2) == FP_NORMAL)
|
|
_sigma2 = varbounds[1] + log(_sigma2/varbounds[1]);
|
|
else
|
|
_sigma2 = varbounds[1] + 1000;
|
|
}
|
|
return _sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| ARCH recursion: Sigma^2 = omega + alpha_1 * res^2_{t-1} |
|
|
//+------------------------------------------------------------------+
|
|
vector arch_recursion(vector& parameters, vector& resids, vector& sigma2, ulong p, ulong nobs, double backcast, matrix &var_bounds)
|
|
{
|
|
for(ulong t = 0; t < nobs; ++t)
|
|
{
|
|
sigma2[t] = parameters[0]; // Constant term (omega)
|
|
for(ulong i = 0; i < p; ++i)
|
|
{
|
|
// Use backcast for pre-sample values
|
|
if((t - i - 1) < 0)
|
|
sigma2[t] += parameters[i + 1] * backcast;
|
|
else
|
|
sigma2[t] += parameters[i + 1] * pow(resids[t - i - 1], 2);
|
|
}
|
|
sigma2[t] = bounds_check(sigma2[t], var_bounds.Row(t));
|
|
}
|
|
return sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| EGARCH recursion: Models log(sigma^2) to ensure positivity |
|
|
//| Handles asymmetries using std_resids and abs_std_resids. |
|
|
//+------------------------------------------------------------------+
|
|
vector egarch_recursion(vector& parameters, vector& resids, vector& sigma2, ulong p, ulong o, ulong q, ulong nobs, double backcast, matrix& var_bounds, vector& lnsigma2, vector& std_resids, vector& abs_std_resids)
|
|
{
|
|
for(ulong t = 0; t < nobs; ++t)
|
|
{
|
|
ulong loc = 0;
|
|
lnsigma2[t] = parameters[loc];
|
|
loc += 1;
|
|
|
|
// Compute symmetric innovation component
|
|
for(ulong j = 0; j < p; ++j)
|
|
{
|
|
if((long(t) - 1 - long(j)) >= 0)
|
|
lnsigma2[t] += parameters[loc] * (abs_std_resids[t - 1 - j] - SQRT2_OV_PI);
|
|
loc += 1;
|
|
}
|
|
// Compute asymmetric (leverage) innovation component
|
|
for(ulong j = 0; j < o; ++j)
|
|
{
|
|
if((long(t) - 1 - long(j)) >= 0)
|
|
lnsigma2[t] += parameters[loc] * std_resids[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
// Compute persistence (AR) component
|
|
for(ulong j = 0; j < q; ++j)
|
|
{
|
|
if((long(t) - 1 - long(j)) < 0)
|
|
lnsigma2[t] += parameters[loc] * backcast;
|
|
else
|
|
lnsigma2[t] += parameters[loc] * lnsigma2[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
|
|
// Numerical clamping for log-variance
|
|
if(lnsigma2[t] > LNSIGMA_MAX)
|
|
lnsigma2[t] = LNSIGMA_MAX;
|
|
|
|
sigma2[t] = exp(lnsigma2[t]);
|
|
|
|
// Boundary enforcement in variance space
|
|
if(sigma2[t] < var_bounds[t, 0])
|
|
{
|
|
sigma2[t] = var_bounds[t, 0];
|
|
lnsigma2[t] = log(sigma2[t]);
|
|
}
|
|
else if(sigma2[t] > var_bounds[t, 1])
|
|
{
|
|
sigma2[t] = var_bounds[t, 1] + log(sigma2[t]) - log(var_bounds[t, 1]);
|
|
lnsigma2[t] = log(sigma2[t]);
|
|
}
|
|
|
|
// Update standardized residuals for next step
|
|
std_resids[t] = resids[t] / sqrt(sigma2[t]);
|
|
abs_std_resids[t] = MathAbs(std_resids[t]);
|
|
}
|
|
return sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| GARCH recursion: Standard variance update with GJR terms |
|
|
//+------------------------------------------------------------------+
|
|
vector garch_recursion(vector ¶meters, vector& fresids, vector& sresids, vector& sigma2, ulong p, ulong o, ulong q, ulong nobs, double backcast, matrix &var_bounds)
|
|
{
|
|
for(long t = 0; t < long(nobs); ++t)
|
|
{
|
|
long loc = 0;
|
|
sigma2[t] = parameters[loc];
|
|
loc += 1;
|
|
// ARCH terms (p)
|
|
for(long j = 0; j < long(p); ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
sigma2[t] += parameters[loc] * backcast;
|
|
else
|
|
sigma2[t] += parameters[loc] * fresids[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
// GJR terms (o) - Asymmetric leverage effect
|
|
for(long j = 0; j < long(o); ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
sigma2[t] += parameters[loc] * 0.5 * backcast;
|
|
else
|
|
sigma2[t] += (parameters[loc] * fresids[t - 1 - j] * double((sresids[t - 1 - j] < 0)));
|
|
loc += 1;
|
|
}
|
|
// GARCH terms (q) - Persistence
|
|
for(long j = 0; j < long(q); ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
sigma2[t] += parameters[loc] * backcast;
|
|
else
|
|
sigma2[t] += parameters[loc] * sigma2[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
sigma2[t] = bounds_check(sigma2[t], var_bounds.Row(t));
|
|
}
|
|
return sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Compute variance recursion for HARCH (Heterogeneous ARCH) |
|
|
//+------------------------------------------------------------------+
|
|
vector harch_recursion(vector& parameters, vector& resids, vector& sigma2, vector& lags, int nobs, double backcast, matrix& var_bounds)
|
|
{
|
|
for(int t = 0; t < nobs; ++t)
|
|
{
|
|
sigma2[t] = parameters[0];
|
|
double param;
|
|
for(int i = 0; i < int(lags.Size()); ++i)
|
|
{
|
|
param = parameters[i + 1] / lags[i]; // Normalize params by lag length
|
|
for(int j = 0; j < int(lags[i]); ++j)
|
|
{
|
|
if((t - j - 1) >= 0)
|
|
sigma2[t] += param * resids[t - j - 1] * resids[t - j - 1];
|
|
else
|
|
sigma2[t] += param * backcast;
|
|
}
|
|
}
|
|
sigma2[t] = bounds_check(sigma2[t], var_bounds.Row(t));
|
|
}
|
|
return sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| FIGARCH weight helper: Calculates fractional integration weights |
|
|
//+------------------------------------------------------------------+
|
|
vector figarch_weights(vector& _params, ulong p, ulong q, ulong trunc_lag)
|
|
{
|
|
double phi = p ? _params[0] : 0.0;
|
|
double d = p ? _params[1] : _params[0];
|
|
double beta = q ? _params[p + q] : 0.0;
|
|
|
|
vector lam = vector::Zeros(trunc_lag);
|
|
vector delta = vector::Zeros(trunc_lag);
|
|
|
|
lam[0] = phi - beta + d;
|
|
delta[0] = d;
|
|
// Iterative calculation of long-memory weights
|
|
for(ulong i = 1; i < trunc_lag; ++i)
|
|
{
|
|
delta[i] = double(i - d) / double(i + 1) * delta[i - 1];
|
|
lam[i] = beta * lam[i - 1] + (delta[i] - phi * delta[i - 1]);
|
|
}
|
|
return lam;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Compute variance recursion for FIGARCH Variance |
|
|
//+------------------------------------------------------------------+
|
|
vector figarch_recursion(vector& parameters, vector& fresids, vector& sigma2, ulong p, ulong q, ulong nobs, ulong trunc_lag, double backcast, matrix& var_bounds)
|
|
{
|
|
double omega = parameters[0];
|
|
double beta = q ? parameters[1 + p + q] : 0.0;
|
|
double omega_tilde = omega / (1. - beta);
|
|
vector sparams = vector::Zeros(parameters.Size() - 1);
|
|
for(ulong i = 0; i < sparams.Size(); sparams[i] = parameters[i + 1], ++i);
|
|
|
|
vector lam = figarch_weights(sparams, p, q, trunc_lag);
|
|
ulong stop = 0;
|
|
for(ulong t = 0; t < nobs; ++t)
|
|
{
|
|
// Calculate contribution of initial backcast
|
|
double bc_weight = 0.0;
|
|
for(ulong i = t; i < trunc_lag; ++i)
|
|
bc_weight += lam[i];
|
|
|
|
sigma2[t] = omega_tilde + bc_weight * backcast;
|
|
|
|
// Apply weighted residuals (memory)
|
|
stop = fmin(t, trunc_lag);
|
|
for(ulong i = 0; i < stop; ++i)
|
|
sigma2[t] += lam[i] * fresids[t - i - 1];
|
|
|
|
sigma2[t] = bounds_check(sigma2[t], var_bounds.Row(t));
|
|
}
|
|
return sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| EWMA recursion (RiskMetrics approach) |
|
|
//+------------------------------------------------------------------+
|
|
vector ewma_recursion(double lam, vector& resids, vector& sigma2, ulong nobs, double _backcast)
|
|
{
|
|
matrix varbs = matrix::Ones(nobs, 2);
|
|
vector temp = {-1., 1.7e308};
|
|
varbs = np::multiply(varbs, temp);
|
|
|
|
// Map EWMA parameters to GARCH(1,0,1) structure
|
|
temp.Resize(3, 1);
|
|
temp[0] = 0.0;
|
|
temp[1] = 1.0 - lam;
|
|
temp[2] = lam;
|
|
|
|
vector resids2 = pow(resids, 2.);
|
|
sigma2 = garch_recursion(temp, resids2, resids, sigma2, 1, 0, 1, nobs, _backcast, varbs);
|
|
return sigma2;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Base class for volatility updaters (polymorphism support) |
|
|
//+------------------------------------------------------------------+
|
|
class VolatilityUpdater
|
|
{
|
|
protected:
|
|
bool m_initialized;
|
|
// Internal test wrapper
|
|
virtual void update_tester(ulong t, vector& parameters, vector& resids, vector& _sigma2, matrix& varbounds)
|
|
{
|
|
update(t, parameters, resids, _sigma2, varbounds);
|
|
}
|
|
public:
|
|
VolatilityUpdater(void):m_initialized(false) {}
|
|
VolatilityUpdater(VolatilityUpdater& other) { m_initialized = other.is_initialized(); }
|
|
void operator=(VolatilityUpdater& other) { m_initialized = other.is_initialized(); }
|
|
virtual ~VolatilityUpdater(void) {}
|
|
|
|
virtual bool is_initialized(void) { return m_initialized; }
|
|
virtual void initialize_update(vector ¶meters, vector& backcast, ulong _nobs);
|
|
virtual void update(ulong t, vector ¶meters, vector& resids, vector& _sigma2, matrix& varbounds);
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| GarchUpdater class: Manages specific GARCH recursion settings |
|
|
//+------------------------------------------------------------------+
|
|
class GarchUpdater: public VolatilityUpdater
|
|
{
|
|
private:
|
|
ulong m_p, m_o, m_q;
|
|
double m_power;
|
|
double m_backcast;
|
|
public:
|
|
GarchUpdater(void) {}
|
|
GarchUpdater(ulong p, ulong o, ulong q, double power) { initialize(p, o, q, power); }
|
|
GarchUpdater(GarchUpdater& other)
|
|
{
|
|
m_p = other.get_p();
|
|
m_o = other.get_o();
|
|
m_q = other.get_q();
|
|
m_power = other.get_power();
|
|
m_backcast = other.get_backcast();
|
|
m_initialized = other.is_initialized();
|
|
}
|
|
|
|
void operator=(GarchUpdater& other)
|
|
{
|
|
m_p = other.get_p();
|
|
m_o = other.get_o();
|
|
m_q = other.get_q();
|
|
m_power = other.get_power();
|
|
m_backcast = other.get_backcast();
|
|
m_initialized = other.is_initialized();
|
|
}
|
|
|
|
virtual ~GarchUpdater(void) {}
|
|
|
|
ulong get_p(void) { return m_p; }
|
|
ulong get_o(void) { return m_o; }
|
|
ulong get_q(void) { return m_q; }
|
|
double get_power(void) { return m_power; }
|
|
double get_backcast(void) { return m_backcast; }
|
|
|
|
bool initialize(ulong p, ulong o, ulong q, double power)
|
|
{
|
|
m_p = p;
|
|
m_o = o;
|
|
m_backcast = -1.;
|
|
m_q = q;
|
|
m_power = power;
|
|
m_initialized = (m_p || m_o);
|
|
return m_initialized;
|
|
}
|
|
|
|
virtual void initialize_update(vector& parameters, vector& backcast, ulong nobs) override
|
|
{
|
|
m_backcast = backcast[0];
|
|
}
|
|
|
|
// Core update loop for GARCH volatility
|
|
virtual void update(ulong t, vector ¶meters, vector& resids, vector& sigma2, matrix& var_bounds) override
|
|
{
|
|
ulong loc = 0;
|
|
sigma2[t] = parameters[loc];
|
|
loc += 1;
|
|
|
|
// ARCH terms with power scaling
|
|
for(ulong j = 0; j < m_p; ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
sigma2[t] += parameters[loc] * m_backcast;
|
|
else
|
|
sigma2[t] += parameters[loc] * pow(MathAbs(resids[t - 1 - j]), m_power);
|
|
loc += 1;
|
|
}
|
|
// GJR terms with power scaling
|
|
for(ulong j = 0; j < m_o; ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
sigma2[t] += parameters[loc] * 0.5 * m_backcast;
|
|
else
|
|
sigma2[t] += (parameters[loc] * pow(MathAbs(resids[t - 1 - j]), m_power) * (resids[t - 1 - j] < 0));
|
|
loc += 1;
|
|
}
|
|
// GARCH persistence terms
|
|
for(ulong j = 0; j < m_q; ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
sigma2[t] += parameters[loc] * m_backcast;
|
|
else
|
|
sigma2[t] += parameters[loc] * sigma2[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
sigma2[t] = bounds_check(sigma2[t], var_bounds.Row(t));
|
|
}
|
|
};
|
|
/*+------------------------------------------------------------------+
|
|
//| EWMAUpdater |
|
|
//+------------------------------------------------------------------+
|
|
class EWMAUpdater: public VolatilityUpdater
|
|
{
|
|
private:
|
|
bool m_estimate_lam;
|
|
double m_backcast;
|
|
vector m_params;
|
|
EWMAUpdater(void)
|
|
{
|
|
|
|
}
|
|
~EWMAUpdater(void)
|
|
{
|
|
|
|
}
|
|
virtual bool initialize(ulong p,ulong o, ulong q,double power) override
|
|
{
|
|
Print(__FUNCTION__, " wrong function call ");
|
|
return false;
|
|
}
|
|
virtual bool initialize(double lam) override
|
|
{
|
|
m_estimate_lam = (lam == EMPTY_VALUE);
|
|
m_params = vector::Zeros(3);
|
|
if(!m_estimate_lam)
|
|
{
|
|
m_params[1] = 1.-lam;
|
|
m_params[2] = lam;
|
|
}
|
|
|
|
m_initialized = true;
|
|
return m_initialized;
|
|
}
|
|
virtual void initialize_update(vector ¶meters, vector& backcast, ulong _nobs) override
|
|
{
|
|
if(m_estimate_lam)
|
|
{
|
|
m_params[1] = 1.0 - parameters[0];
|
|
m_params[2] = parameters[0];
|
|
}
|
|
m_backcast = backcast[0];
|
|
}
|
|
virtual void update(ulong t, vector ¶meters, vector& resids, vector& sigma2, matrix& var_bounds) override
|
|
{
|
|
sigma2[t] = m_params[0];
|
|
if(t == 0)
|
|
sigma2[t] += m_backcast;
|
|
else
|
|
sigma2[t] += (
|
|
m_params[1] * resids[t - 1] * resids[t - 1]
|
|
+ m_params[2] * sigma2[t - 1]
|
|
);
|
|
sigma2[t] = bounds_check(sigma2[t], var_bounds.Row(t));
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//|EGARCHUpdater |
|
|
//+------------------------------------------------------------------+
|
|
class EGARCHUpdater: public VolatilityUpdater
|
|
{
|
|
private:
|
|
ulong m_p, m_o, m_q;
|
|
double m_backcast;
|
|
vector m_insigma2, m_std_resids, m_abs_std_resids;
|
|
void _resize(ulong nobs)
|
|
{
|
|
if(m_insigma2.Size()<nobs)
|
|
m_insigma2 = m_std_resids = m_abs_std_resids = vector::Zeros(nobs);
|
|
}
|
|
public:
|
|
EGARCHUpdater(ulong p, ulong o, ulong q)
|
|
{
|
|
m_p = p;
|
|
m_o = o;
|
|
m_q = q;
|
|
m_backcast = 9999.99;
|
|
m_insigma2 = m_std_resids = m_abs_std_resids = vector::Zeros(0);
|
|
}
|
|
~EGARCHUpdater(void)
|
|
{
|
|
}
|
|
virtual void initialize_update(vector ¶meters, vector& backcast, ulong _nobs) override
|
|
{
|
|
m_backcast = backcast[0];
|
|
_resize(_nobs);
|
|
}
|
|
virtual void update(ulong t, vector ¶meters, vector& resids, vector& sigma2, matrix& var_bounds) override
|
|
{
|
|
if(t)
|
|
{
|
|
m_std_resids[t - 1] = resids[t - 1] / sqrt(sigma2[t - 1]);
|
|
m_abs_std_resids[t - 1] = MathAbs(m_std_resids[t - 1]);
|
|
}
|
|
|
|
m_insigma2[t] = parameters[0];
|
|
ulong loc = 1;
|
|
for(ulong j = 0; j<m_p; ++j)
|
|
{
|
|
if((t - 1 - j) >= 0)
|
|
m_insigma2[t] += parameters[loc] * (m_abs_std_resids[t - 1 - j] - SQRT2_OV_PI);
|
|
loc += 1;
|
|
}
|
|
for(ulong j = 0; j<m_o; ++j)
|
|
{
|
|
if((t - 1 - j) >= 0)
|
|
m_insigma2[t] += parameters[loc] * m_std_resids[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
for(ulong j = 0; j<m_q; ++j)
|
|
{
|
|
if((t - 1 - j) < 0)
|
|
m_insigma2[t] += parameters[loc] * m_backcast;
|
|
else
|
|
m_insigma2[t] += parameters[loc] * m_insigma2[t - 1 - j];
|
|
loc += 1;
|
|
}
|
|
if(m_insigma2[t] > LNSIGMA_MAX)
|
|
m_insigma2[t] = LNSIGMA_MAX;
|
|
sigma2[t] = exp(m_insigma2[t]);
|
|
if(sigma2[t] < var_bounds[t, 0])
|
|
{
|
|
sigma2[t] = var_bounds[t, 0];
|
|
m_insigma2[t] = log(sigma2[t]);
|
|
}
|
|
else
|
|
if(sigma2[t] > var_bounds[t, 1])
|
|
{
|
|
sigma2[t] = var_bounds[t, 1] + log(sigma2[t]) - log(var_bounds[t, 1]);
|
|
m_insigma2[t] = log(sigma2[t]);
|
|
}
|
|
}
|
|
};*/
|
|
|
|
//+------------------------------------------------------------------+
|