//+------------------------------------------------------------------+
//|                                                Impact_Curves.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"
#resource "\\Files\\VolatilityModels\\sp500.csv" as string equity_data
#property script_show_inputs
#include<VolatilityModels\Arch\Univariate\mean.mqh>
//--- input parameters
input double ScaleFactor=100.;
input ulong _P_ = 1;
input ulong _O_ = 1;
input ulong _Q_ = 1;
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---
   matrix data = np::readcsv_from_string(equity_data,false,",",true,0);
   vector returns = log(data.Col(1));
   returns = np::diff(returns);
//---   
   returns*=ScaleFactor;
//---   
   ArchParameters gjr_spec,tarch_spec,garch_spec;
//---   
   gjr_spec.mean_model_type = MEAN_CONSTANT;
   gjr_spec.observations = returns;
   gjr_spec.garch_p = _P_;
   gjr_spec.garch_q = _Q_;
   gjr_spec.garch_o = _O_;
//---   
   garch_spec = gjr_spec;
   tarch_spec = gjr_spec;
//---   
   gjr_spec.vol_model_type = VOL_GJR_GARCH;
   garch_spec.vol_model_type = VOL_GARCH;
   tarch_spec.vol_model_type = VOL_TARCH;
//---
   ConstantMean gjr_model;
   ConstantMean garch_model;
   ConstantMean tarch_model;
//---
   if(!gjr_model.initialize(gjr_spec) || !garch_model.initialize(garch_spec) || !tarch_model.initialize(tarch_spec))
    return;
//---
   ArchModelResult gjr_params = gjr_model.fit(ScaleFactor);
//---
   ArchModelResult tarch_params = tarch_model.fit(ScaleFactor);
//---
   ArchModelResult garch_params = garch_model.fit(ScaleFactor);
//---
   vector shocks = np::linspace(-5.0,5.0,100);
//---
   vector gjr_response,tarch_response,garch_response;
//---
   gjr_response = gjr_impact(shocks,gjr_params.params);
   tarch_response = tarch_impact(shocks,tarch_params.params);
   garch_response = garch_impact(shocks,garch_params.params);
   
   matrix xaxis,yaxis;
   xaxis = np::vectorAsRowMatrix(shocks,3);
   yaxis = matrix::Zeros(3,shocks.Size());
   
   yaxis.Row(gjr_response,0);
   yaxis.Row(tarch_response,1);
   yaxis.Row(garch_response,2);
//---
   np::plotmatrices(xaxis,yaxis," Impact curves",false,"Shocks(%)","Predicted Volatility","GJR-GARCH(volatility),TARCH(volatility),GARCH(volatility)",CURVE_LINES,3,0,0,0,0,750,400,true,40);
  }
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
//| Function to calculate GJR-GARCH response                         |
//+------------------------------------------------------------------+
vector gjr_impact(vector& z, vector& p)
  {
   vector isless_than_zero = np::whereVectorIsLt(z,0.0);
   return p[1] + p[2] * pow(z,2.) + p[3] * (isless_than_zero) * pow(z,2.0);
  }
//+------------------------------------------------------------------+
//| Function to calculate GARCH response                             |
//+------------------------------------------------------------------+
vector garch_impact(vector& z, vector& p)
  {
   return p[1] + p[2] * pow(z,2.) + p[3] * pow(z,2.0);
  }
//+------------------------------------------------------------------+
//| Function to calculate TARCH response                             |
//+------------------------------------------------------------------+
vector tarch_impact(vector& z, vector& p)
  {
   vector isless_than_zero = np::whereVectorIsLt(z,0.0);
   return p[1] + p[2] * fabs(z) + p[3] * (isless_than_zero) * fabs(z);
  }