//+------------------------------------------------------------------+
//|                                                          acf.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<Math\Stat\ChiSquare.mqh>
//+------------------------------------------------------------------+
//| autocorrelation function result struct                           |
//+------------------------------------------------------------------+
struct ACFResult
  {
   vector            acf;
   vector            qstats;
   vector            pvalues;
   matrix            conf_intervals;
                     ACFResult(void)
     {
      acf = qstats = pvalues = vector::Zeros(0);
      conf_intervals = matrix::Zeros(0,0);
     }
                     ACFResult(vector& _acf, vector& _qstats, vector& _pvalues, matrix& _confintervals)
     {
      acf = _acf;
      qstats = _qstats;
      pvalues = _pvalues;
      conf_intervals = _confintervals;
     }
                     ACFResult(ACFResult& other)
     {
      acf = other.acf;
      qstats = other.qstats;
      pvalues = other.pvalues;
      conf_intervals = other.conf_intervals;
     }
   void              operator=(ACFResult& other)
     {
      acf = other.acf;
      qstats = other.qstats;
      pvalues = other.pvalues;
      conf_intervals = other.conf_intervals;
     }
  };
//+------------------------------------------------------------------+
//| ACF function                                                     |
//+------------------------------------------------------------------+
ACFResult acf(vector& x, ulong nlags = 0, double alpha = 0.05, bool bartlett_conf=false, bool demean= true, bool adjusted = false)
  {
   ACFResult out;
   ulong nobs = x.Size();

   if(!nobs)
     {
      Print(__FUNCTION__, " input container is empty ");
      return out;
     }

   if(!nlags)
      nlags = MathMin(ulong(10*log10(nobs)), nobs-1);

   vector avf = acovf(x,0,demean,adjusted);

   if(!avf.Size())
      return out;

   out.acf = np::sliceVector(avf,0,long(nlags+1));
   out.acf /=avf[0];
   int ec = 0;
   if(bartlett_conf)
     {
      vector varacf = vector::Ones(out.acf.Size())/double(nobs);
      varacf[0] = 0.;
      varacf[1] = 1./double(nobs);
      vector temp = np::sliceVector(out.acf,1,out.acf.Size()-1);
      vector tempsq = pow(temp,2.);
      temp = np::sliceVector(varacf,2);
      temp*=1.0+2.*tempsq.CumSum();
      np::vectorCopy(varacf,temp,2);

      vector interval = CNormalDistr::InvNormalCDF(1. - alpha/2.)*sqrt(varacf);
      out.conf_intervals = matrix::Zeros(2,varacf.Size());
      out.conf_intervals.Row(out.acf-interval,0);
      out.conf_intervals.Row(out.acf+interval,1);
     }
   else
     {
      double varacf = 1./double(x.Size());
      double nv = CNormalDistr::InvNormalCDF(1. - alpha/2.);
      double interval = nv*sqrt(varacf);
      out.conf_intervals = matrix::Zeros(2,out.acf.Size());
      out.conf_intervals.Row(out.acf-interval,0);
      out.conf_intervals.Row(out.acf+interval,1);
     }

   vector temp = np::sliceVector(out.acf,1);
   q_stat(temp,nobs,out.qstats,out.pvalues);

   return out;
  }
/*//+------------------------------------------------------------------+
//|Find the next regular number greater than or equal to target.     |
//+------------------------------------------------------------------+
long next_target(long target)
  {
   if(target<=6)
      return target;
   if(!(target&(target-1)))
      return target;
   double match = AL_POSINF;
   long p5 = 1;
   long p35,quotient,N,p2;
   while(p5<target)
     {
      p35 = p5;
      while(p35<target)
        {
         quotient =-(-target/p35);
         p2=(long)pow(2,(long)bit_length(quotient-1));
         N = p2*p35;
         if(N == target)
            return N;
         else
            if(double(N)<match)
               match = double(N);
         p35*=3;
         if(p35==target)
            return p35;
        }
      p5*=5;
      if(p5==target)
         return p5;
     }
   if(double(p5)<match)
      match = double(p5);
   return long(match);
  }
//+------------------------------------------------------------------+
//|  get the number bits for an integer                              |
//+------------------------------------------------------------------+
ulong bit_length(long x)
  {
   long var = MathAbs(x);
   ulong count = 1;
   while(var>=2)
     {
      var = var/2;
      ++count;
     }
   return count;
  }*/
//+------------------------------------------------------------------+
//|acovf                                                             |
//+------------------------------------------------------------------+
vector acovf(vector &in, ulong lags = 0,bool demean= true, bool adjusted = false)
  {
   vector xo,acov;

   if(demean)
      xo = in - in.Mean();
   else
      xo = in;

   ulong n = in.Size();

   ulong laglen = lags;

   if(!laglen)
      laglen = n-1;
   else
      if(lags>n-1)
        {
         Print(__FUNCTION__, " nlag must be smaller than ",n - 1);
         return vector::Zeros(0);
        }

   acov = vector::Zeros(laglen+1);

   acov[0] = xo.Dot(xo);

   if(lags)
     {
      vector a,b;
      for(ulong i = 0; i<laglen; ++i)
        {
         a = np::sliceVector(xo,long(i+1));
         b = np::sliceVector(xo,0,long(xo.Size()-(i+1)));
         acov[i+1] = a.Dot(b);
        }
      if(adjusted)
         acov/=(double(n) - np::arange(laglen+1));
      else
         acov/=double(n);
      return acov;
     }

   vector xi,temp,d;

   if(adjusted)
     {
      xi = np::arange(n,1.);
      d = vector::Zeros((xi.Size()*2)-1);
      np::vectorCopy(d,xi,0,xi.Size());
      temp = np::sliceVector(xi,0,xi.Size()-1);
      np::reverseVector(temp);
      np::vectorCopy(d,temp,xi.Size());
     }
   else
      d = double(n)*vector::Ones(2*n-1);
   /*if(fft)
    {
     ulong nobs = xo.Size();
     n = ulong(next_target(long(2*nobs+1)));
     CRowComplex Frf;
     CRowDouble inx = xo;
     CFastFourierTransform::FFTR1D(inx,int(nobs),Frf);
     vectorc frf = Frf.ToVector();
     frf = frf*frf.Conjugate();
     Frf = frf;
     CFastFourierTransform::FFTR1DInv(Frf,Frf.Size(),inx);
     temp = inx.ToVector();
     temp = np::sliceVector(temp,0,nobs);
     acov = temp/np::sliceVector(d,long(nobs-1));
    }
   else
    {*/
   acov = xo.Correlate(xo,VECTOR_CONVOLVE_FULL);
   acov = np::sliceVector(acov,long(n-1))/np::sliceVector(d,long(n-1));
//}

   if(lags)
      return np::sliceVector(acov,0,long(laglen+1));

   return acov;
  }
//+------------------------------------------------------------------+
//|LjungBox q statistic                                              |
//+------------------------------------------------------------------+
bool q_stat(vector& x, ulong nvars, vector& out_stats, vector& out_pvals)
  {
   vector r = np::arange(x.Size(),1.) ;
   r = (1.0/(double(nvars) - r))*pow(x,2.);
   out_stats = double(nvars)*double(nvars+2)*r.CumSum();
   out_pvals = out_stats;
   r = np::arange(x.Size(),1.);
   int error_note = 0;
   for(ulong i = 0; i<out_pvals.Size(); ++i)
      out_pvals[i] =  1. - MathCumulativeDistributionChiSquare(out_stats[i],r[i],error_note);
   if(error_note)
      return false;
   return true;
  }
//+------------------------------------------------------------------+
//| plot the ACF                                                     |
//+------------------------------------------------------------------+

void plot_acf(ACFResult& acf,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(!acf.acf.Size())
    {
     Print(__FUNCTION__, " Results object is empty ");
     return;
    }
   
   if(StringLen(plot_name)<1)
    plot_name = "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)acf.acf.Size());
   ArrayResize(y,(int)acf.acf.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] = acf.acf[j];
            ctype = (i==0)?CURVE_HISTOGRAM:CURVE_POINTS;
           }
         else
           {
            ctype = CURVE_LINES;
            if(i == 2)
               y[j] = (acf.conf_intervals[1,j]-acf.conf_intervals[0,j])/2.;
            else
               y[j] = (acf.conf_intervals[0,j]-acf.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);
//---
}