//+------------------------------------------------------------------+
//|                                                        ljung.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"acf.mqh"
///+------------------------------------------------------------------+
//| Ljung-Box and Box-Pierce test of autocorrelation in residuals.   |
//| Returns a matrix: Col 0 = LB Stat, Col 1 = LB p-val              |
//| (Optional Box-Pierce) Col 2 = BP Stat, Col 3 = BP p-val          |
//+------------------------------------------------------------------+
matrix ljungboxtest(vector& x, ulong lags=0, ulong model_df=0, ulong period = 0, bool demean=false, bool boxpierce=false, bool autolag=false)
  {
   //--- Validation: Seasonal periods must encompass at least 2 observations
   if(period==1)
     {
      Print(__FUNCTION__, " period must be >= 2");
      return matrix::Zeros(0,0);
     }

   ulong nobs = x.Size();
   vector laggs;
   
   //+------------------------------------------------------------------+
   //| Automated Lag Selection (Information Criterion Style)            |
   //+------------------------------------------------------------------+
   if(autolag)
     {
      ulong maxlag = nobs-1;
      ACFResult sacf = acf(x,maxlag,0.05,false,demean,false);
      if(!sacf.acf.Size())
         return matrix::Zeros(0,0);
         
      vector ssacf(maxlag);
      vector r(maxlag);
      
      // Shift ACF results to isolate positive lags (drop lag 0)
      for(ulong i = 1; i<(maxlag+1); ssacf[i-1]=sacf.acf[i], ++i);
      for(ulong i = 0; i<(maxlag); r[i] = double(i+1), ++i);
      
      vector qsacf;
      // Calculate standard running test statistics for automated profiling
      if(!boxpierce)
        {
         vector n = pow(ssacf,2.)/(double(nobs)-r);
         qsacf = nobs*(nobs+2.)*n.CumSum();
        }
      else
        {
         vector n = pow(ssacf,2.);
         qsacf = nobs*n.CumSum();
        }
        
      // Dynamic penalty adjustments based on information thresholds
      double q = 2.4;
      double threshold = sqrt(q*log(nobs));
      double threshold_metric = MathAbs(sacf.acf).Max()*sqrt(nobs);
      
      r = vector::Zeros(nobs-1);
      for(ulong i = 0; i<r.Size(); r[i] = double(i+1), ++i);
      
      // Apply HQ / BIC style sequence selection penalties
      if(threshold_metric<=threshold)
         qsacf = qsacf - (r*log(nobs)); // Severe penalty (No strong autocorrelation detected)
      else
         qsacf = qsacf - (2.*r);        // Lighter AIC-style penalty
         
      // Select lag index maximizing the penalized metric sequence
      ulong m = qsacf.ArgMax();
      laggs=vector::Zeros(m);
      for(ulong i = 0; i<laggs.Size(); laggs[i] = double(i+1), ++i);
     }
   //+------------------------------------------------------------------+
   //| Manual or Rule-of-Thumb Lag Selection                            |
   //+------------------------------------------------------------------+
   else
     {
      if(period) // Contextual seasonal spacing rule
        {
         ulong m = MathMin(nobs/5,2*period);
         laggs = vector::Zeros(m);
         for(ulong i = 0; i<laggs.Size(); laggs[i] = double(i+1), ++i);
        }
      else
        {
         if(!lags) // Default rule-of-thumb: min(N/5, 10)
           {
            ulong m = MathMin(nobs/5,10);
            laggs = vector::Zeros(m);
            for(ulong i = 0; i<laggs.Size(); laggs[i] = double(i+1), ++i);
           }
         else // Explicitly defined user bounds
           {
            laggs = vector::Zeros(lags);
            for(ulong i = 0; i<laggs.Size(); laggs[i] = double(i+1), ++i);
           }
        }
     }
     
   //+------------------------------------------------------------------+
   //| Execution: Core Statistic Computations                           |
   //+------------------------------------------------------------------+
   ulong maxlag = (ulong)laggs.Max();
   ACFResult sacf = acf(x,maxlag,0.05,false,demean,false);
   if(!sacf.acf.Size())
      return matrix::Zeros(0,0);
      
   vector ssacf(maxlag);
   vector r(maxlag);
   
   // Parse execution arrays
   for(ulong i = 1; i<(maxlag+1); ssacf[i-1]=sacf.acf[i], ++i);
   for(ulong i = 0; i<(maxlag); r[i] = double(i+1), ++i);
   
   // 1. Compute Ljung-Box Q-Statistic cumulative terms
   vector ssacf2 = pow(ssacf,2.)/(nobs-r);
   vector temp = ssacf2.CumSum();
   vector qljungbox(temp.Size());
   for(ulong i = 0; i<laggs.Size(); qljungbox[i] = nobs*(nobs+2)*temp[ulong(laggs[i])-1], ++i);
   
   // 2. Adjust degrees of freedom for estimated model parameters (e.g., ARMA(p,q) coefficients)
   vector adj_lags = laggs - double(model_df);
   vector pval = qljungbox;
   pval.Fill(double("nan"));
   int e = 0;
   
   // 3. Compute Chi-Square survival probability (p-values) for Ljung-Box
   for(ulong i = 0; i<adj_lags.Size(); ++i)
      if(adj_lags[i]>0)
         pval[i] = 1. - MathCumulativeDistributionChiSquare(qljungbox[i],adj_lags[i],e);
         
   // Build base matrix packaging results
   matrix out;
   out = matrix::Zeros(adj_lags.Size(),2);
   out.Col(qljungbox,0);
   out.Col(pval,1);
   
   //+------------------------------------------------------------------+
   //| Optional: Compute Classic Box-Pierce Statistics                  |
   //+------------------------------------------------------------------+
   if(boxpierce)
     {
      out.Resize(out.Rows(),4); // Expand matrix allocation to hold 4 metrics
      temp = pow(ssacf,2.0);
      temp = temp.CumSum();
      vector qboxpierce = temp;
      
      // Calculate classic BP statistics (Unscaled by sample denominator increments)
      for(ulong i = 0; i<laggs.Size(); qboxpierce[i] = nobs*temp[ulong(laggs[i])-1], ++i);
      for(ulong i = 0; i<adj_lags.Size(); ++i)
         if(adj_lags[i]>0)
            pval[i] = 1. - MathCumulativeDistributionChiSquare(qboxpierce[i],adj_lags[i],e);

      // Append Box-Pierce results to final columns
      out.Col(qboxpierce,2);
      out.Col(pval,3);
     }
     
   return out;
  }
//+------------------------------------------------------------------+