//+------------------------------------------------------------------+
//|                                                   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
const double LNSIGMA_MAX = log(DBL_MAX) - 0.1;
//+------------------------------------------------------------------+
//|  helper function                                                 |
//+------------------------------------------------------------------+
double bounds_check(double _sigma2, vector &varbounds)
  {
   if(_sigma2<varbounds[0])
      _sigma2 = varbounds[1];
   else
      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 utility                                            |
//+------------------------------------------------------------------+
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];
      for(ulong i = 0; i<p; ++i)
        {
         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;
  }
//+------------------------------------------------------------------+
//|Compute variance recursion for EGARCH models                      |
//+------------------------------------------------------------------+
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;
      for(ulong j = 0; j<p; ++j)
        {
         if((t - 1 - j) >= 0)
            lnsigma2[t] += parameters[loc] * (abs_std_resids[t - 1 - j] - SQRT2_OV_PI);
         loc += 1;
        }
      for(ulong j = 0; j<o; ++j)
        {
         if((t - 1 - j) >= 0)
            lnsigma2[t] += parameters[loc] * std_resids[t - 1 - j];
         loc += 1;
        }
      for(ulong j = 0; j<q; ++j)
        {
         if((t - 1 - j) < 0)
            lnsigma2[t] += parameters[loc] * backcast;
         else
            lnsigma2[t] += parameters[loc] * lnsigma2[t - 1 - j];
         loc += 1;
        }
      if(lnsigma2[t] > LNSIGMA_MAX)
         lnsigma2[t] = LNSIGMA_MAX;
      sigma2[t] = exp(lnsigma2[t]);
      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]);
           }
      std_resids[t] = resids[t] / sqrt(sigma2[t]);
      abs_std_resids[t] = MathAbs(std_resids[t]);
     }
   return sigma2;
  }

//+------------------------------------------------------------------+
//|Compute variance recursion for GARCH and related models           |
//+------------------------------------------------------------------+
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;
      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;
        }
      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] * (sresids[t - 1 - j] < 0));
         loc+=1;
        }
      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 EWMA/RiskMetrics Variance          |
//+------------------------------------------------------------------+
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);
   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 that all volatility updaters must inherit from.       |
//+------------------------------------------------------------------+
class VolatilityUpdater
  {
protected:
   bool m_initialized;
   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();
     }
                    ~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);
  };
//+------------------------------------------------------------------+
//| garch updater                                                    |
//+------------------------------------------------------------------+
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 = get_p();
      m_o = get_o();
      m_q = get_q();
      m_power = get_power();
      m_backcast = get_backcast();
      m_initialized = other.is_initialized();
     }
     
     void            operator=(GarchUpdater& other)
     {
      m_p = get_p();
      m_o = get_o();
      m_q = get_q();
      m_power = get_power();
      m_backcast = get_backcast();
      m_initialized = other.is_initialized();
     }
                    ~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];
     }
   virtual void      update(ulong t, vector &parameters, vector& resids, vector& sigma2, matrix& var_bounds) override
     {
      ulong loc = 0;
      sigma2[t] = parameters[loc];
      loc += 1;

      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;
        }
      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;
        }
      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]);
           }
     }
  };*/

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