//+------------------------------------------------------------------+
//|                                                         pacf.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"..\..\np_graphics.mqh"
#include"..\..\Regression\utils.mqh"
//+------------------------------------------------------------------+
//| autocorrelation function result struct                           |
//+------------------------------------------------------------------+
struct PACFResult
  {
   vector            pacf;
   matrix            conf_intervals;
                     PACFResult(void)
     {
      pacf = vector::Zeros(0);
      conf_intervals = matrix::Zeros(0,0);
     }
                     PACFResult(vector& _pacf, matrix& _confintervals)
     {
      pacf = _pacf;
      conf_intervals = _confintervals;
     }
                     PACFResult(PACFResult& other)
     {
      pacf = other.pacf;
      conf_intervals = other.conf_intervals;
     }
   void              operator=(PACFResult& other)
     {
      pacf = other.pacf;
      conf_intervals = other.conf_intervals;
     }
  };
//+------------------------------------------------------------------+
//|  yule walker wrapper                                                      |
//+------------------------------------------------------------------+
vector pacf_yw(vector& x, ulong lags=0, ENUM_YW_METHOD method=YW_ADJUSTED)
  {
   ulong nobs = x.Size();
   if(!lags)
      lags = (ulong)MathMax(MathMin(ulong(10*log10(nobs)),nobs-1),1);
   vector pacf = vector::Zeros(lags+1);
   pacf[0] = 1.;
   vector inputx = x;
   for(ulong k = 1; k<lags+1; ++k)
     {
      YWResult yw = yule_walker(x,long(k),method);
      if(yw.rho.Size())
         pacf[k] = yw.rho[yw.rho.Size()-1];
     }
   return pacf;
  }
//+------------------------------------------------------------------+
//| Partial autocorrelation estimate.                                |
//+------------------------------------------------------------------+
PACFResult pacf(vector& x_data, ulong nlags = 0, double alpha = 0.05)
  {
   ulong nobs = x_data.Size();
   if(!nlags)
      nlags = (ulong)MathMin(ulong(10*log10(nobs)),nobs/2-1);
   nlags = (ulong)MathMax(nlags,1);

   if(nlags>(x_data.Size()/2))
     {
      Print(__FUNCTION__, " Can only compute partial correlations for lags up to 50% of the sample size.");
      return PACFResult();
     }
   vector ret = pacf_yw(x_data,nlags,YW_ADJUSTED);
   double varacf = 1.0/double(x_data.Size());
   double interval = CNormalDistr::InvNormalCDF(1.0 - alpha/2.0)*sqrt(varacf);
   matrix confint = matrix::Zeros(2,ret.Size());

   if(!confint.Row(ret-interval,0) ||
      !confint.Row(ret+interval,1))
     {
      Print(__FUNCTION__," error ", GetLastError());
      return PACFResult();
     }
     
   confint[0,0] = confint[1,0] = ret[0];

   return PACFResult(ret,confint);
  }

//+------------------------------------------------------------------+
//| plot the PACF                                                    |
//+------------------------------------------------------------------+

void plot_pacf(PACFResult& pacf,long chart_id = 0,int x_coordinate = 0, int y_coordinate = 0, int subwindow = 0,int width = 600, int height = 400,string plot_name = NULL, int fontsize = 15, long viewtimeseconds = 30, bool showchart=true )
{
   if(!pacf.pacf.Size())
    {
     Print(__FUNCTION__, " Results object is empty ");
     return;
    }
   
   if(StringLen(plot_name)<1)
    plot_name = "Partial Autocorrelation Function";
    
   CGraphic* graph = new CGraphic();
//---
   if(!chart_id)
     chart_id = ChartID();
//---
   ChartRedraw(chart_id);
//---
   ChartSetInteger(chart_id, CHART_SHOW, showchart);
//---
   if(!graph.Create(chart_id,plot_name,subwindow,x_coordinate,y_coordinate,width,height))
     {
      delete graph;
      Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
      return;;
     }
//---
   double x[], y[];
   ArrayResize(x,(int)pacf.pacf.Size());
   ArrayResize(y,(int)pacf.pacf.Size());
//---
   int ncurves = 4;
//---
   ENUM_CURVE_TYPE ctype = WRONG_VALUE;
//---
   for(int i=0; i<ncurves; ++i)
     {
      for(uint j = 0; j<x.Size(); ++j)
        {
         if(i<2)
           {
            x[j] = double(j);
            y[j] = pacf.pacf[j];
            ctype = (i==0)?CURVE_HISTOGRAM:CURVE_POINTS;
           }
         else
           {
            ctype = CURVE_LINES;
            if(i == 2)
               y[j] = (pacf.conf_intervals[1,j]-pacf.conf_intervals[0,j])/2.;
            else
               y[j] = (pacf.conf_intervals[0,j]-pacf.conf_intervals[1,j])/2.;
            
           }
        }
//---
      CCurve* newcurve = graph.CurveAdd(x, y,ctype,NULL);
//---
      switch(ctype)
       {
        case CURVE_HISTOGRAM:
          newcurve.Color(ColorToARGB(clrBlue));
          newcurve.HistogramWidth(2);
          break;
        case CURVE_POINTS:
          newcurve.PointsFill(true);
          newcurve.PointsColor(ColorToARGB(clrBlue));
          newcurve.Color(ColorToARGB(clrBlue));
          newcurve.PointsSize(8);
          break;
        default:
          newcurve.Color(ColorToARGB(clrRed));
          newcurve.LinesStyle(STYLE_DOT);
          break;
       }
     }
//---   
   graph.BackgroundMain(plot_name);
   graph.BackgroundMainSize(fontsize+5);
   graph.BackgroundMainColor(ColorToARGB(clrBlack));  
//---
   graph.HistoryNameWidth(0);
   graph.HistorySymbolSize(0);
   graph.HistoryNameSize(0);
//---
   graph.XAxis().Name("Lags");
   graph.XAxis().NameSize(fontsize);
//---
   graph.YAxis().Name("Correlation Coefficient");
   graph.YAxis().NameSize(fontsize);
//---
   graph.CurvePlotAll();
   graph.Update();
//---
   Sleep(int(viewtimeseconds)*1000);
   graph.Destroy();
   delete graph;
   ChartRedraw(chart_id);
//---
}