//+------------------------------------------------------------------+
//|                                           EGARCH_volatility.mq5  |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
//--- Preprocessor directives and indicator properties
#define __SLSQP__                                 // Enable Sequential Least Squares Programming optimization algorithm
#property indicator_separate_window               // Render indicator in an independent subwindow
#include "Arch\Univariate\mean.mqh"         // Include statistical mean model header dependencies
//--- Define buffer and plotting counts
#property indicator_buffers 7
#property indicator_plots   7
// --- Plot 1: Conditional Volatility (Standard Deviation)
#property indicator_label1  "ConditionalVolatility"
#property indicator_type1   DRAW_LINE
#property indicator_color1  clrBlue
#property indicator_style1  STYLE_SOLID
#property indicator_width1  1
// --- Plot 2: Conditional Variance
#property indicator_label2  "ConditionalVariance"
#property indicator_type2   DRAW_LINE
#property indicator_color2  clrGreen
#property indicator_style2  STYLE_SOLID
#property indicator_width2  1

// --- Plot 3: Standardized Residuals
#property indicator_label3  "StandardizedResiduals"
#property indicator_type3   DRAW_LINE
#property indicator_color3  clrRed
#property indicator_style3  STYLE_DASH
#property indicator_width3  1
// --- Plots 4-7: Estimated Model Parameters (Hidden from chart drawing, readable in Data Window)
#property indicator_label4  "Omega"
#property indicator_type4   DRAW_NONE

#property indicator_label5  "Alpha"
#property indicator_type5   DRAW_NONE

#property indicator_label6  "Gamma"
#property indicator_type6   DRAW_NONE

#property indicator_label7  "Beta"
#property indicator_type7   DRAW_NONE

//+------------------------------------------------------------------+
//| INPUT PARAMETERS                                                 |
//+------------------------------------------------------------------+
input int                   BarsToDraw        = 500;           // Number of historical bars to calculate and render
input ulong                 HistoryLen        = 500;           // Rolling sample size (lookback period) for EGARCH fitting
input double                ScaleFactor       = 100.;          // Multiplier applied to log returns (improves numerical optimizer stability)
input ENUM_MEAN_MODEL       MeanModel         = MEAN_CONSTANT; // Model type for time-series conditional mean
input bool                  MeanConstant      = true;          // Include intercept constant in mean specification
input string                MeanLags          = "";            // Comma-separated list of AR lag terms (e.g. "1,2")
ENUM_VOLATILITY_MODEL       VolatilityModel   = VOL_EGARCH;    // Volatility process specified as Exponential GARCH
ulong                 _P_               = 1;             // GARCH order (lagged conditional log-variances)
ulong                 _O_               = 1;             // Asymmetry order (leverage/asymmetric magnitude terms)
ulong                 _Q_               = 1;             // ARCH order (lagged innovations)
input int                   Volatility_Seed   = 0;             // Seed for volatility model random generator initialization
input ENUM_DISTRIBUTION_MODEL ErrorDistribution = DIST_NORMAL; // Innovation error distribution assumption
input int                   Distribution_Seed = 0;             // Seed for distribution model random generator initialization

//+------------------------------------------------------------------+
//| GLOBAL INDICATOR BUFFERS AND STATE VARIABLES                     |
//+------------------------------------------------------------------+
// Dynamic arrays mapped to indicator plot buffers
double         ConditionalVolatilityBuffer[];
double         ConditionalVarianceBuffer[];
double         StandardizedResidualsBuffer[];
double         OmegaBuffer[];
double         AlphaBuffer[];
double         GammaBuffer[];
double         BetaBuffer[];

// Computational vectors and model pointers
vector         returns = vector::Zeros(HistoryLen); // Rolling logarithmic returns vector
vector         stdresid = returns;                  // Model-standardized residuals vector
ArchParameters model_spec;                          // Struct holding model specifications and hyperparameters
HARX*          full_model;                          // Polymorphic pointer to mean model instance
vector         vol_params;                          // Vector slice containing extracted EGARCH parameters
long           volmodelparams, distmodelparams, allmodelparams;
ulong          output_size = 0;

//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int OnInit()
  {
   // Validate lookback window length
   if(HistoryLen < 30)
     {
      Print("Invalid input value for HistoryLen: Should be >= 30");
      return INIT_FAILED;
     }

   // --- Map dynamic arrays to indicator buffers
   SetIndexBuffer(0, ConditionalVolatilityBuffer, INDICATOR_DATA);
   SetIndexBuffer(1, ConditionalVarianceBuffer, INDICATOR_DATA);
   SetIndexBuffer(2, StandardizedResidualsBuffer, INDICATOR_DATA);
   SetIndexBuffer(3, OmegaBuffer, INDICATOR_DATA);
   SetIndexBuffer(4, AlphaBuffer, INDICATOR_DATA);
   SetIndexBuffer(5, GammaBuffer, INDICATOR_DATA);
   SetIndexBuffer(6, BetaBuffer, INDICATOR_DATA);

   // --- Set plot drawing offset to suppress initial uncalculated bars
   PlotIndexSetInteger(0, PLOT_DRAW_BEGIN, BarsToDraw);
   PlotIndexSetInteger(1, PLOT_DRAW_BEGIN, BarsToDraw);
   PlotIndexSetInteger(2, PLOT_DRAW_BEGIN, BarsToDraw);
   PlotIndexSetInteger(3, PLOT_DRAW_BEGIN, BarsToDraw);
   PlotIndexSetInteger(4, PLOT_DRAW_BEGIN, BarsToDraw);
   PlotIndexSetInteger(5, PLOT_DRAW_BEGIN, BarsToDraw);
   PlotIndexSetInteger(6, PLOT_DRAW_BEGIN, BarsToDraw);

   // --- Set standard empty value sentinel for plots
   PlotIndexSetDouble(0, PLOT_EMPTY_VALUE, EMPTY_VALUE);
   PlotIndexSetDouble(1, PLOT_EMPTY_VALUE, EMPTY_VALUE);
   PlotIndexSetDouble(2, PLOT_EMPTY_VALUE, EMPTY_VALUE);
   PlotIndexSetDouble(3, PLOT_EMPTY_VALUE, EMPTY_VALUE);
   PlotIndexSetDouble(4, PLOT_EMPTY_VALUE, EMPTY_VALUE);
   PlotIndexSetDouble(5, PLOT_EMPTY_VALUE, EMPTY_VALUE);
   PlotIndexSetDouble(6, PLOT_EMPTY_VALUE, EMPTY_VALUE);

   // --- Initialize buffers with empty markers
   ArrayInitialize(ConditionalVolatilityBuffer, EMPTY_VALUE);
   ArrayInitialize(ConditionalVarianceBuffer, EMPTY_VALUE);
   ArrayInitialize(StandardizedResidualsBuffer, EMPTY_VALUE);
   ArrayInitialize(OmegaBuffer, EMPTY_VALUE);
   ArrayInitialize(AlphaBuffer, EMPTY_VALUE);
   ArrayInitialize(GammaBuffer, EMPTY_VALUE);
   ArrayInitialize(BetaBuffer, EMPTY_VALUE);

   // --- Configure model specifications
   model_spec.mean_model_type = MeanModel;

   // Parse comma-separated mean lag configuration string if provided
   if(StringLen(MeanLags))
     {
      string lag_info[];
      int nlags = StringSplit(MeanLags, StringGetCharacter(",", 0), lag_info);
      if(nlags > 0)
        {
         for(uint i = 0; i < uint(nlags); ++i)
           {
            if(StringLen(lag_info[i]) > 0)
              {
               if(model_spec.mean_lags.Resize(model_spec.mean_lags.Size() + 1, 3))
                  model_spec.mean_lags[model_spec.mean_lags.Size() - 1] = StringToDouble(lag_info[i]);
               else
                 {
                  Print(" error ", GetLastError());
                  return INIT_FAILED;
                 }
              }
           }
        }
     }

   // Set EGARCH structural dynamics and seeds
   model_spec.vol_rng_seed = Volatility_Seed;
   model_spec.garch_o     = _O_;
   model_spec.garch_p     = _P_;
   model_spec.garch_q     = _Q_;
   model_spec.dist_type   = ErrorDistribution;
   model_spec.dist_rng_seed = Distribution_Seed;

   // Factory instantiation of the requested mean model class
   switch(MeanModel)
     {
      case MEAN_CONSTANT:
         full_model = new ConstantMean();
         break;
      case MEAN_ZERO:
         full_model = new ZeroMean();
         break;
      case MEAN_AR:
         full_model = new AR();
         break;
      default:
         full_model = new ConstantMean();
         break;
     }

   // Ensure memory allocation for mean model succeeded
   if(CheckPointer(full_model) == POINTER_INVALID)
      return INIT_FAILED;

   model_spec.vol_model_type = VolatilityModel;

   return(INIT_SUCCEEDED);
  }

//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
void OnDeinit(const int reason)
  {
   // Free dynamic memory allocated for the mean model instance
   if(CheckPointer(full_model) == POINTER_DYNAMIC)
      delete full_model;
  }

//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int OnCalculate(const int32_t rates_total,
                const int32_t prev_calculated,
                const datetime &time[],
                const double &open[],
                const double &high[],
                const double &low[],
                const double &close[],
                const long &tick_volume[],
                const long &volume[],
                const int32_t &spread[])
  {
   // Verify total bars available meet the required lookback + rendering window
   if(rates_total < int32_t(HistoryLen + BarsToDraw))
     {
      Print("Not enough bars for indicator calculation");
      return -1;
     }

   // Determine starting index for incremental bar calculation
   int32_t limit = 0;
   if(prev_calculated <= 0)
      limit = rates_total - int32_t(fabs(BarsToDraw)); // First run: calculate specified historical depth
   else
      limit = prev_calculated - 1;                     // Subsequent runs: update only latest bar(s)

   // Main calculation loop iterating through historical price bars
   for(int32_t shift = limit; shift < rates_total; ++shift)
     {
      // Reset current bar values to empty defaults
      ConditionalVarianceBuffer[shift] = ConditionalVolatilityBuffer[shift] = StandardizedResidualsBuffer[shift] = EMPTY_VALUE;
      
      int32_t from = (shift - int32_t(HistoryLen)) + 1;
      
      // Calculate logarithmic returns over the rolling lookback window
      for(int32_t i = from, k = 0; k < int32_t(HistoryLen); ++i, ++k)
         returns[k] = log(close[i] / close[i - 1]);
      
      // Scale returns to assist optimizer convergence
      returns *= fabs(ScaleFactor);
      
      // Load current sample window into specification structure
      model_spec.observations = returns;
      
      // Re-initialize model state with new sample window
      if(!full_model.initialize(model_spec))
        {
         Print(" initialization error ");
         continue;
        }
      
      // Fit EGARCH model via maximum likelihood optimization
      ArchModelResult result = full_model.fit();
      output_size = result.conditional_volatility.Size();
      if(!output_size)
        {
         Print(" model fit error ");
         continue;
        }

      // Extract standardized residuals and terminal fitted values
      stdresid = result.std_resid();
      
      // Convert fitted log-volatility back to original scale (exp) and compute variance/volatility
      ConditionalVarianceBuffer[shift]   = exp(pow(result.conditional_volatility[output_size - 1], 2));
      ConditionalVolatilityBuffer[shift] = exp(result.conditional_volatility[output_size - 1]);
      StandardizedResidualsBuffer[shift] = stdresid[output_size - 1];

      // Retrieve model parameter counts on first successful fit
      if(!vol_params.Size())
        {
         volmodelparams  = long(full_model.volatility().numParams());
         distmodelparams = long(full_model.distribution().numParams());
         allmodelparams  = long(result.params.Size());
        }

      // Slice out EGARCH specific parameters (Omega, Alpha, Gamma, Beta) from full parameter vector
      vol_params = np::sliceVector(result.params, allmodelparams - (volmodelparams + distmodelparams), allmodelparams - distmodelparams);
      // Store estimated parameters into indicator buffers
      OmegaBuffer[shift] = vol_params[0];
      AlphaBuffer[shift] = vol_params[1];
      GammaBuffer[shift] = vol_params[2];
      BetaBuffer[shift] =  vol_params[3];
     }

   // Return calculated count to optimize subsequent iteration calls
   return(rates_total);
  }
//+------------------------------------------------------------------+