Article-22714-Volatility-Modeling-Part-3/Include/slsqp_article/Arch/Univariate/ar.mqh

//+------------------------------------------------------------------+
//|                                                           ar.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"..\..\np.mqh"
#include"..\..\Regression\OLS.mqh"
#include"mean.mqh"
//+------------------------------------------------------------------+
//| info criterion                                                   |
//+------------------------------------------------------------------+
enum ENUM_INFO_CRITERION
 {
  AIC=0,//AIC
  BIC//BIC
 };

//+------------------------------------------------------------------+
//|Bits view                                                         |
//+------------------------------------------------------------------+
matrix bitview(vector& data, int basis, int view)
  {
   if(!data.Size())
      return matrix::Zeros(0,0);
   matrix out = matrix::Zeros(data.Size(),basis/view);
   int value = 0;
   int rem = 0;
   for(ulong i = 0; i<out.Rows(); ++i)
     {
      value = int(fabs(data[i]))/int(pow(2,8*view));
      rem = (int)MathMod(int(fabs(data[i])),int(pow(2,8*view)));
      out[i,0] = rem;
      for(ulong j = 1; j<out.Cols(); ++j)
        {
         if(value)
           {
            out[i,j] = (value<int(pow(2,8*view)))?value:MathMod(value,int(pow(2,8*view)));
           }
         else
            break;
         value/=(int)pow(2,8*view);
        }
     }
   return out;
  }
//+------------------------------------------------------------------+
//|unpack                                                            |
//+------------------------------------------------------------------+
matrix unpackbits(matrix& bview, ulong bits)
  {
   if(!bits || bits<bview.Cols() || !bview.Cols())
    return matrix::Zeros(0,0);
   matrix out = matrix::Zeros(bview.Rows(),bits);
   ulong factor = bits/bview.Cols();
   ulong k = 0;
   for(ulong i = 0; i<out.Rows(); ++i)
   {
    for(ulong j = 0; j<bview.Cols(); ++j)
     {
      ulong value = ulong(bview[i,j]);
      for(long k = 7; k>-1; --k)
        {
         if(value)
           out[i,j*8+k] = MathMod(value,2);
         else
           break;
         value/=2; 
        }      
     }
   }
   return out;
  }
//+------------------------------------------------------------------+
//|Autoregressive AR-X(p) model order selection.                     |
//+------------------------------------------------------------------+
vector ar_select_order(vector& endog, matrix& exog, ulong maxlag,bool use_constant=true, bool global_search = false, ENUM_INFO_CRITERION ic = BIC, ulong holdback=0, ulong period=0)
 {
  ulong nexog = exog.Cols();
  
  ArchParameters modelparams;
  modelparams.y = endog;
  modelparams.x = exog;
  modelparams.mean_constant = use_constant;
  
  ARX model;
  if(!model.initialize(modelparams))
    return vector::Zeros(0);
  
  vector y = endog;
  matrix x = model.get_regressors();
  
  ulong base_col = x.Cols() - nexog - maxlag;
  
  vector sel = vector::Ones(x.Cols());
  vector lags[];
  double ics[];
  int index[];
  matrix xx;
  OLS ols;
  if(!global_search)
   {
    np::fillVector(sel,0.0,long(base_col),long(base_col+maxlag));
    ArrayResize(lags,int(maxlag+1));
    for(ulong i = 0; i<maxlag+1; ++i)
      {
       index.Push(int(i));
       np::fillVector(sel,1.,long(base_col),long(base_col+i));
       if(!sel.Sum())
        continue;
       xx = matrix::Zeros(y.Size(),ulong(sel.Sum()));
       for(ulong k = 0,kk = 0; k<sel.Size(); ++k)
         {
          if(sel[k])
            xx.Col(x.Col(k),kk++);
         }
       lags[i] = np::range(int(i+1),1);
       if(!ols.Fit(y,xx))
         return vector::Zeros(0);
       ics.Push(ic==AIC?ols.Aic():ols.Bic());
      }
   }
  else
   {
    vector b = np::arange(ulong(pow(2,maxlag)));
    matrix bview = bitview(b,4,1);
    matrix bits = unpackbits(bview,32);
    matrix bb;
    for(ulong i = 0; i<4; ++i)
      {
       bb = np::sliceMatrixCols(bits,long(8*i),long(8*(i+1)));
       np::reverseMatrixCols(bb);
       np::matrixCopyCols(bits,bb,long(8*i),long(8*(i+1)));
      }
    matrix masks = np::sliceMatrixCols(bits,0,long(maxlag));
    ArrayResize(lags,int(masks.Rows()));
    for(ulong i = 0; i<masks.Rows(); ++i)
     {
      index.Push(int(i));
      np::vectorCopy(sel,masks.Row(i),long(base_col),long(base_col+maxlag));
      if(!sel.Sum())
        continue;
      xx = matrix::Zeros(y.Size(),ulong(sel.Sum()));
      
      for(ulong k = 0,kk = 0; k<sel.Size(); ++k)
         {
          if(sel[k])
            xx.Col(x.Col(k),kk++);
         }
         
      lags[i] = vector::Zeros(ulong(masks.Row(i).Sum()));
      
      for(ulong j = 0, k = 0; j<lags[i].Size(); ++j) 
           if(masks[i,j])
              lags[i][k++] = double(j+1);
              
      if(!ols.Fit(y,xx))
         return vector::Zeros(0);
      ics.Push(ic==AIC?ols.Aic():ols.Bic());
     } 
   }
  MathQuickSortAscending(ics,index,0,int(ics.Size()-1));
  return lags[index[0]];
 }
//+------------------------------------------------------------------+