//+------------------------------------------------------------------+
//|                                                        utils.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<Math\Alglib\solvers.mqh>
//+------------------------------------------------------------------+
//| Constant and trend used in regression model                      |
//+------------------------------------------------------------------+

enum ENUM_TREND
  {
   TREND_NONE=0,
   TREND_CONST_ONLY,
   TREND_LINEAR_ONLY,
   TREND_LINEAR_CONST,
   TREND_QUAD_LINEAR_CONST
  };
//+------------------------------------------------------------------+
//| Options for how to handle existing constants                     |
//+------------------------------------------------------------------+

enum ENUM_HAS_CONST
  {
   HAS_CONST_RAISE=0,
   HAS_CONST_SKIP,
   HAS_CONST_ADD
  };
//+------------------------------------------------------------------+
//| Options for  trimming invalid observations                       |
//+------------------------------------------------------------------+
enum ENUM_TRIM
  {
   TRIM_NONE=0,
   TRIM_FORWARD,
   TRIM_BACKWARD,
   TRIM_BOTH
  };
//+------------------------------------------------------------------+
//| options for how to handle original data set                      |
//+------------------------------------------------------------------+

enum ENUM_ORIGINAL
  {
   ORIGINAL_EX=0,
   ORIGINAL_IN,
   ORIGINAL_SEP
  };
//+------------------------------------------------------------------+
//| Yule walker calculation options                                  |
//+------------------------------------------------------------------+
enum ENUM_YW_METHOD
  {
   YW_ADJUSTED=0,//adjusted
   YW_MLE//mle
  };
//+------------------------------------------------------------------+
//|   YW result struct                                               |
//+------------------------------------------------------------------+
struct YWResult
  {
   double            sigma;
   vector            rho;
   matrix            inv;
                     YWResult(void)
     {
      sigma = EMPTY_VALUE;
      rho = vector::Zeros(0);
      inv = matrix::Zeros(0,0);
     }
                    ~YWResult(void)
     {
     }
                     YWResult(YWResult& other)
     {
      sigma = other.sigma;
      rho = other.rho;
      inv = other.inv;
     }
   void              operator=(YWResult& other)
     {
      sigma = other.sigma;
      rho = other.rho;
      inv = other.inv;
     }
  };
//+------------------------------------------------------------------+
//| sequence                                                         |
//+------------------------------------------------------------------+
//| The function generates a sequence of double numbers              |
//| with given starting and ending values and step.                  |
//|                                                                  |
//| Arguments:                                                       |
//| from        : Starting value of the sequence                     |
//| to          : Last value of the sequence                         |
//| step        : Step of the sequence                               |
//| result[]    : Array for the sequence values                      |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
bool Sequence(const double from,const double to,const double step,double &result[])
  {
//--- check NaN
   if(!MathIsValidNumber(from) || !MathIsValidNumber(to) || !MathIsValidNumber(step))
      return(false);

   if(to<from)
      return(false);
//--- calculate number of elements and prepare output array
   int count=1+int((to-from)/step);

   if(ArraySize(result)<count)
      if(ArrayResize(result,count)!=count)
         return(false);
//--- prepare sequence
   for(int i=0; i<count; i++)
      result[i]=from+i*step;
//---
   return(true);
  }
//+------------------------------------------------------------------+
//|helper function : adds selected trend and\or constant             |
//+------------------------------------------------------------------+
bool addtrend(matrix &in,matrix &out,ENUM_TREND trend=TREND_CONST_ONLY, bool prepend=false, ENUM_HAS_CONST has_const=HAS_CONST_SKIP)
  {
//---
   ulong trendorder=0;
//---
   if(trend==TREND_NONE)
      return out.Copy(in);
//---
   if(trend==TREND_CONST_ONLY)
      trendorder = 0;
   else
      if(trend==TREND_LINEAR_CONST || trend==TREND_LINEAR_ONLY)
         trendorder = 1;
      else
         if(trend==TREND_QUAD_LINEAR_CONST)
            trendorder = 2;
//---
   ulong nobs = in.Rows();
//---
   vector tvector,temp;
//---
   if(!tvector.Resize(nobs) || !temp.Resize(nobs))
     {
      Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
      return false;
     }
//---
   double sequence[];
//---
   if(!Sequence(1,nobs,1,sequence) || !tvector.Assign(sequence))
     {
      Print(__FUNCTION__," ",__LINE__," Error ", ::GetLastError());
      return false;
     }
//---
   matrix trendmat;
//---
   if(!trendmat.Resize(tvector.Size(),(trend==TREND_LINEAR_ONLY)?1:trendorder+1))
     {
      Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
      return false;
     }
//---
   for(ulong i = 0; i<trendmat.Cols(); i++)
     {
      temp = MathPow(tvector,(trend==TREND_LINEAR_ONLY)?double(i+1):double(i));
      if(!trendmat.Col(temp,i))
        {
         Print(__FUNCTION__," ",__LINE__," Matrix Assign Error ", ::GetLastError());
         return false;
        }
     }
//---
   vector ptp = in.Ptp(0);
//---
   if(!ptp.Min())
     {
      if(has_const==HAS_CONST_RAISE)
        {
         Print("Input matrix contains one or more constant columns");
         return false;
        }
      if(has_const==HAS_CONST_SKIP)
        {
         matrix vsplitted[];

         ulong parts[] = {1,trendmat.Cols()-1};

         trendmat.Vsplit(parts,vsplitted);

         if(!trendmat.Copy(vsplitted[1]))
           {
            Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
            return false;
           }
        }
     }
//---
   int order = (prepend)?1:-1;
//---
   if(!out.Resize(trendmat.Rows(),trendmat.Cols()+in.Cols()))
     {
      Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
      return false;
     }
//---
   ulong j = 0;
   matrix mtemp;
//---
   if(prepend)
      mtemp.Copy(trendmat);
   else
      mtemp.Copy(in);
//---
   for(j = 0; j<mtemp.Cols(); j++)
     {
      vector col = mtemp.Col(j);

      if(!out.Col(col,j))
        {
         Print(__FUNCTION__," ",__LINE__," Matrix Assign Error ", ::GetLastError());
         return false;
        }

     }
//---
   ulong minus = j;
//---
   if(prepend)
      mtemp.Copy(in);
   else
      mtemp.Copy(trendmat);
//---
   for(; j<out.Cols(); j++)
     {
      vector col = mtemp.Col(j-minus);

      if(!out.Col(col,j))
        {
         Print(__FUNCTION__," ",__LINE__," Matrix Assign Error ", ::GetLastError());
         return false;
        }

     }
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//| helper function: transforms rhs matrix                           |
//+------------------------------------------------------------------+
bool lagmat(matrix &in,matrix &out[],ulong mlag,ENUM_TRIM trim=TRIM_BOTH,ENUM_ORIGINAL original=ORIGINAL_IN)
  {
//---
   ulong nobs=in.Rows();
   ulong nvars =in.Cols();
   ulong dropidx = 0;
//---
   if(mlag>=nobs)
     {
      Print("mlag should be < number of rows of the input matrix");
      return false;
     }
//---
   if(original==ORIGINAL_EX || original==ORIGINAL_SEP)
      dropidx = nvars;
//---
   matrix nn;
//---
   if(!nn.Resize(nobs + mlag,nvars * (mlag + 1)))
     {
      Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
      return false;
     }
//---
   nn.Fill(0.0);
//---
   ulong maxlag=mlag;
//---
   ulong row_end,row_start,col_end,col_start,j,z;
   row_end=row_start=col_end=col_start=j=z=0;
//---
   for(ulong k = 0; k<(maxlag+1); k++)
     {
      row_start = maxlag-k;
      row_end = nobs+maxlag-k;
      col_start = nvars*(maxlag-k);
      col_end = nvars*(maxlag-k+1);
      j = 0;
      for(ulong irow = row_start; irow < row_end; irow++, j++)
        {
         z = 0;
         for(ulong icol = col_start; icol < col_end; icol++, z++)
            nn[irow][icol]=in[j][z];
        }
     }
//---
   ulong startobs,stopobs;
   if(trim==TRIM_NONE || trim==TRIM_FORWARD)
      startobs=0;
   else
      startobs=maxlag;
//---
   if(trim==TRIM_NONE || trim==TRIM_BACKWARD)
      stopobs=nn.Rows();
   else
      stopobs=nobs;
//---
   if(dropidx && original==ORIGINAL_SEP)
      ArrayResize(out,2);
   else
      ArrayResize(out,1);
//---
   if(!out[0].Resize(stopobs-startobs,nn.Cols()-dropidx))
     {
      Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
      return false;
     }
//---
   for(ulong irow = startobs; irow<stopobs; irow++)
      for(ulong icol=dropidx; icol<nn.Cols(); icol++)
         out[0][irow-startobs][icol-dropidx]=nn[irow][icol];

//---
   if(out.Size()>1)
     {
      if(!out[1].Resize(stopobs-startobs,dropidx))
        {
         Print(__FUNCTION__," ",__LINE__," Resize Error ", GetLastError());
         return false;
        }
      //---
      for(ulong irow = startobs; irow<stopobs; irow++)
         for(ulong icol=0; icol<out[1].Cols(); icol++)
            out[1][irow-startobs][icol]=nn[irow][icol];
     }
//---
   return true;
  }
//+------------------------------------------------------------------+
//| toeplitz matrix construction                                     |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> toeplitz(vector<T>& x)
  {
   ulong n = x.Size();
   matrix<T> out(n,n);
   for(ulong i = 0; i<n; ++i)
      for(ulong j = 0; j<n; ++j)
         if(i<=j)
            out[i,j] = x[j-i];
         else
            out[i,j] = x[i-j];
   return out;
  }
//+------------------------------------------------------------------+
//| toeplitz matrix construction                                     |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> toeplitz(vector<T>& x, vector<T>& y)
  {
   ulong nx = x.Size();
   ulong ny = y.Size();
   vector<T> z = vector<T>::Zeros(nx+ny-1);
   ulong start = 1;
   for(ulong i = 0; i<z.Size(); ++i)
     {
      if(i<nx)
         z[i] = x[(nx-i)-1];
      else
         z[i] = y[start++];
     }

   matrix<T> out(nx,ny);
   for(ulong i = 0; i<nx; ++i)
     {
      for(ulong j = 0; j<ny; ++j)
        {
         if(i<j)
            out[i,j] = z[nx+(j-1)-i];
         else
            out[i,j] = z[nx-(i+1)+j];
        }
     }
   return out;
  }
//+---------------------------------------------------------------------------+
//| Estimate AR(p) parameters from a sequence using the Yule-Walker equations.|
//+---------------------------------------------------------------------------+
YWResult yule_walker(vector& x, long order=1, ENUM_YW_METHOD method =YW_ADJUSTED, long df = -1, bool inv = false, bool demean = true)
  {
   YWResult out;
   vector in = x;
   if(!in.Size())
     {
      Print(__FUNCTION__, " input sequence is empty ");
      return out;
     }

   if(demean)
      in = in - in.Mean();
   long n = df>0?df:long(in.Size());
   vector r  = vector::Zeros(order+1);
   r[0] = pow(in,2.0).Sum()/double(n);
   double adj_needed = (double)(method == YW_ADJUSTED);
   for(long k = 1; k<(order+1); ++k)
     {
      double sum = 0;
      for(long i = 0; i<(long(in.Size()) - k); ++i)
         sum+=(in[i]*in[k+i]);
      r[k] = sum/double(n - k*adj_needed);
     }
     
   vector rslice = vector::Zeros(r.Size()-1);
   for(ulong i = 0; i<rslice.Size(); rslice[i] = r[i], ++i);
   matrix R = toeplitz(rslice);
   for(ulong i = 0; i<rslice.Size(); rslice[i] = r[i+1], ++i);
   
   double b[],sol[];
   ArrayResize(b,(int)rslice.Size());
   for(uint i = 0; i<b.Size(); b[i] = rslice[i], ++i);
   
   CMatrixDouble Rr = R;
   CDenseSolverLSReport ls;
   int opinfo;
   CDenseSolver::RMatrixSolveLS(Rr,Rr.Rows(),Rr.Cols(),b,0,opinfo,ls,sol);
   if(opinfo<0)
     {
      Print(__FUNCTION__ " possible singular matrix. Using pinv ");
      out.rho = (R.PInv()).MatMul(rslice);
     }
   else
     out.rho.Assign(sol); 
   double sigmasq = r[0] - (rslice*out.rho).Sum();
   if(MathClassify(sigmasq)==FP_NORMAL && sigmasq>0)
      out.sigma = sqrt(sigmasq);
   else
      out.sigma = double("nan");
   if(inv)
      out.inv = R.Inv();
   return out;
  }

//+------------------------------------------------------------------+