//+------------------------------------------------------------------+
//|                                             APARCH_Benchmark.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"
#property script_show_inputs
//---
#include"Arch\univariate\mean.mqh"
//--- Script Inputs
input int      InpHoldoutSize   = 100;         // Forecast horizon
input int      InpInSampleSize  = 1000;        // In-Sample Training Bar Count
input ENUM_TIMEFRAMES InpTimeframe = PERIOD_D1;// Execution Timeframe
input double    InpScaleFactor  = 100.0;      // Scaling multiplier
input ENUM_FORECAST_METHOD InpForeCastMethod = FORECAST_ANALYTIC; // Forecasting method

// Struct to store final benchmarking metrics for each model
struct ModelPerformance
  {
   string   model_name;
   double   aic;
   double   bic;
   double   oos_mse;
   double   oos_mae;
   double   log_like;
  };

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   if(InpHoldoutSize<=0 || InpInSampleSize<50)
    {
     Alert("Invalid input parameters for InpHoldoutSize and/or InpInSampleSize\nBoth must be >0 and >=50 respectively.");
     return;
    }
   Print("==========================================================");
   Print(" Starting Volatility Model Evaluation: APARCH vs. Competitors");
   Print("==========================================================");

   // Prepare Returns Data
   ulong total_needed = (ulong)(InpInSampleSize + InpHoldoutSize + 1);
   vector close_prices;
   
   if(!close_prices.CopyRates(_Symbol, InpTimeframe, COPY_RATES_CLOSE, 0, (uint)total_needed))
     {
      Print("Error: Failed to fetch price data.");
      return;
     }

   // Compute Log Returns: r_t = ln(P_t / P_{t-1}) * scalefactor
   ulong n_returns = close_prices.Size() - 1;
   vector returns(n_returns);
   for(ulong i = 0; i < n_returns; ++i)
     {
      returns[i] = MathLog(close_prices[i + 1] / close_prices[i]) * InpScaleFactor;
     }

   // Split into In-Sample and Out-of-Sample arrays
   ulong in_sample_len = n_returns - (ulong)InpHoldoutSize;
   vector in_sample_returns = np::sliceVector(returns, 0, (long)in_sample_len);
   vector oos_returns       = np::sliceVector(returns, (long)in_sample_len, (long)n_returns);

   // Realized variance benchmark for OOS (using squared returns as proxy)
   vector oos_realized_variance = MathPow(oos_returns, 2.0);

   // Define Volatility Models to Compare
   ENUM_VOLATILITY_MODEL models[] =
     {
      VOL_GARCH,
      VOL_GJR_GARCH,
      VOL_EGARCH,
      VOL_TARCH,
      VOL_APARCH
     };

   string model_names[] =
     {
      "GARCH(1,1)",
      "GJR-GARCH(1,1,1)",
      "EGARCH(1,1,1)",
      "TARCH(1,1,1)",
      "APARCH(1,1,1)"
     };

   uint num_models = ArraySize(models);
   ModelPerformance results[];
   ArrayResize(results, num_models);

   // Loop Through and Fit Each Model
   for(uint i = 0; i < num_models; ++i)
     {
      PrintFormat("\n--- Fitting Model [%d/%d]: %s ---", i + 1, num_models, model_names[i]);

      // Populate base ARCH parameters structure
      ArchParameters params;
      params.observations       = in_sample_returns;
      params.mean_model_type    = MEAN_CONSTANT;
      params.vol_model_type     = models[i];
      params.dist_type          = DIST_NORMAL;
      params.garch_p            = 1;
      params.garch_o            = 1; // Asymmetry term enabled where applicable
      params.garch_q            = 1;

      // Special APARCH specific defaults if selected
      if(models[i] == VOL_APARCH)
        {
         params.aparch_delta       = EMPTY_VALUE; // Jointly estimate dynamic power parameter
         params.aparch_common_asym = false;
        }

      // Construct model instance (using CConstantMean model wrapper)
      ConstantMean model;
      if(!model.initialize(params))
        {
         PrintFormat("Failed to initialize model: %s", model_names[i]);
         continue;
        }

      // Optimize parameters over In-Sample dataset
      ArchModelResult fit_result = model.fit();

      if(fit_result.solver_return_code != 0)
        {
         PrintFormat("Warning: Solver ended with code %d for model %s", fit_result.solver_return_code, model_names[i]);
        }

      // Record In-Sample Information Criteria
      results[i].model_name = model_names[i];
      results[i].log_like   = fit_result.loglikelihood;
      results[i].aic        = fit_result.aic();
      results[i].bic        = fit_result.bic();

      // Evaluate Out-of-Sample (OOS) Volatility Forecasts
      matrix dummy_x[];
      vector fitted_params = model.get_params();
      
      // Perform multi-step analytical forecast across the holdout period
      ArchForecast fcast = model.forecast(fitted_params, dummy_x, (ulong)InpHoldoutSize, -1, InpForeCastMethod);

      if(fcast.variance.Rows() > 0)
        {
         // Extract 1-step to H-step variance forecast vector
         vector pred_variance = fcast.variance.Row(fcast.variance.Rows() - 1);

         // Calculate Forecast Error Metrics (Predicted Variance vs Realized Proxy)
         vector errors = oos_realized_variance - pred_variance;
         
         results[i].oos_mse = MathPow(errors, 2.0).Mean();
         results[i].oos_mae = MathAbs(errors).Mean();
        }
      else
        {
         results[i].oos_mse = DBL_MAX;
         results[i].oos_mae = DBL_MAX;
        }
     }

   // Print Comparison Summary Report Table
   Print("\n=========================================================================");
   Print("                        MODEL COMPARISON RESULTS                         ");
   Print("=========================================================================");
   PrintFormat("%-15s | %-12s | %-12s | %-12s | %-12s", "Model", "AIC", "BIC", "OOS MSE", "OOS MAE");
   Print("-------------------------------------------------------------------------");

   int best_aic_idx = 0;
   int best_mse_idx = 0;

   for(uint i = 0; i < num_models; ++i)
     {
      PrintFormat("%-15s | %-12.4f | %-12.4f | %-12.6f | %-12.6f",
                  results[i].model_name,
                  results[i].aic,
                  results[i].bic,
                  results[i].oos_mse,
                  results[i].oos_mae);

      if(results[i].aic < results[best_aic_idx].aic) best_aic_idx = (int)i;
      if(results[i].oos_mse < results[best_mse_idx].oos_mse) best_mse_idx = (int)i;
     }

   Print("-------------------------------------------------------------------------");
   PrintFormat("Optimal In-Sample Fit (Lowest AIC): %s", results[best_aic_idx].model_name);
   PrintFormat("Optimal Out-of-Sample Robustness (Lowest MSE): %s", results[best_mse_idx].model_name);
   Print("=========================================================================");
  }
//+------------------------------------------------------------------+