//+------------------------------------------------------------------+
//|                                                   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 &parameters, 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 &parameters, vector& backcast, ulong _nobs);
   virtual void       update(ulong t, vector &parameters, 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 &parameters, 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 &parameters, 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 &parameters, 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 &parameters, vector& backcast, ulong _nobs) override
     {
      m_backcast = backcast[0];
      _resize(_nobs);
     }
   virtual void      update(ulong t, vector &parameters, 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]);
           }
     }
  };*/

//+------------------------------------------------------------------+