//+------------------------------------------------------------------+
//|                                                          OLS.mqh |
//|                                  Copyright 2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2023, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
//+------------------------------------------------------------------+
//| matrix multiplication                                            |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> matmul(const matrix<T>& matrix_a, const matrix<T>& matrix_b)
  {
   matrix<T> matrix_c = matrix<T>::Zeros(0,0);

   if(matrix_a.Cols()!=matrix_b.Rows())
      return(matrix_c);

   ulong M=matrix_a.Rows();
   ulong K=matrix_a.Cols();
   ulong N=matrix_b.Cols();
   matrix_c=matrix<T>::Zeros(M,N);

   for(ulong m=0; m<M; ++m)
      for(ulong k=0; k<K; ++k)
         for(ulong n=0; n<N; ++n)
            matrix_c[m][n]+=matrix_a[m][k]*matrix_b[k][n];

   return(matrix_c);
  }
//+------------------------------------------------------------------+
//| Matrix multiply with vector                                      |
//+------------------------------------------------------------------+
template<typename T>
vector<T> matmul(const matrix<T>& matrix_a, const vector<T>& vector_b)
  {
   matrix<T> matrix_b = matrix<T>::Zeros(vector_b.Size(),1);
   matrix_b.Col(vector_b,0);
   matrix<T>out = matmul(matrix_a,matrix_b);
   return out.Col(0);
  }
//+------------------------------------------------------------------+
//| Ordinary least squares class                                     |
//+------------------------------------------------------------------+
class OLS
  {
private:
   bool   m_const_prepended;    //check column for constant in design matrix
   matrix m_exog,               //design matrix
          m_pinv,               //pseudo-inverse of matrix
          m_cov_params,         //covariance of matrix
          m_m_error,            //error matrix
          m_norm_cov_params;    //normalized covariance matrix
   vector m_endog,              //dependent variables
          m_weights,            //weights
          m_singularvalues,     //singular values of solution
          m_params,             //coefficients of regression model(solution)
          m_tvalues,            //test statistics of model
          m_bse,                //standard errors of model
          m_const_cols,         //mark constant columns in design matrix
          m_resid;              //residuals of model
   ulong  m_obs,                //number of observations
          m_model_dof,          //degrees of freedom of model
          m_resid_dof,          //degrees of freedom of residuals
          m_kconstant,          //number of constants
          m_rank;               //rank of design matrix
   double m_aic,                //Akiake information criteria
          m_bic,                //Bayesian information criteria
          m_scale,              //scale of model
          m_llf,                //loglikelihood of model
          m_sse,                //sum of squared errors
          m_rsqe,               //r-squared of model
          m_centeredtss,        //centered sum of squares
          m_uncenteredtss;      //uncentered sum of squares
   uint              m_error;              //error flag
   // private methods
   ulong             countconstants(void);
   void              scale(void);
   void              sse(void);
   void              rsqe(void);
   void              centeredtss(void);
   void              uncenteredtss(void);
   void              aic(void);
   void              bic(void);
   void              bse(void);
   void              llf(void);
   void              tvalues(void);
   void              covariance_matrix(void);


public:
   //constructor
                     OLS(void);
   //destructor
                    ~OLS(void);
   //public methods
   bool              Fit(vector &y_vars,matrix &x_vars);
   double            Predict(vector &inputs);
   double            Predict(double _input);
   //get properties of OLS model
   ulong             ModelDOF(void) { if(m_error) return 0; else return m_model_dof;}
   ulong             ResidDOF(void) { if(m_error) return 0; else return m_resid_dof;}
   double            Scale(void)  { if(m_error) return EMPTY_VALUE; else return m_scale;    }
   double            Aic(void)    { if(m_error) return EMPTY_VALUE; else return m_aic;      }
   double            Bic(void)    { if(m_error) return EMPTY_VALUE; else return m_bic;    }
   double            Sse(void)    { if(m_error) return EMPTY_VALUE; else return m_sse;    }
   double            Rsqe(void)   { if(m_error) return EMPTY_VALUE; else return m_rsqe;   }
   double            C_tss(void)  { if(m_error) return EMPTY_VALUE; else return m_centeredtss;}
   double            Loglikelihood(void) { if(m_error) return EMPTY_VALUE; return m_llf; }
   vector            Tvalues(void) { if(m_error) return m_m_error.Col(0); return m_tvalues; }
   vector            Residuals(void) { if(m_error) return m_m_error.Col(0); return m_resid; }
   vector            ModelParameters(void) { if(m_error) return m_m_error.Col(0); return m_params; }
   vector            Bse(void) { if(m_error) return m_m_error.Col(0);  return m_bse; } 
   matrix            CovarianceMatrix(void) { if(m_error) return m_m_error; return m_cov_params; }
  };
//+------------------------------------------------------------------+
//|  constructor                                                     |
//+------------------------------------------------------------------+
OLS::OLS(void)
  {
   m_kconstant = 0;
   m_obs = 0;
   m_model_dof = 0;
   m_resid_dof = 0;
   m_kconstant = 0;
   m_rank = 0;
   m_aic = 0;
   m_bic = 0;
   m_scale = 0;
   m_llf   = 0;
   m_error = 1;

   m_m_error.Resize(2,2);
   m_m_error.Fill(EMPTY_VALUE);

  }
//+------------------------------------------------------------------+
//|Destructor                                                        |
//+------------------------------------------------------------------+
OLS::~OLS(void)
  {

  }
//+------------------------------------------------------------------+
//| Fit data                                                         |
//+------------------------------------------------------------------+
bool OLS::Fit(vector &y_vars,matrix &x_vars)
  {
//---re-initialize internal properties
   m_endog = y_vars;
//---
   m_exog  = x_vars;
//---
   m_kconstant = countconstants();
//---
   m_error = 0;
//---
   m_model_dof = m_resid_dof = m_rank = m_obs= 0;
//---
   m_aic = m_bic= m_llf=0;
//---
   m_obs=m_exog.Rows();
//---
   m_weights = vector::Ones(m_obs);
//---
   m_rank = m_exog.Rank();
//---
   m_model_dof = m_rank - m_kconstant;
//---
   m_resid_dof = m_obs - m_rank;
//---
   if(y_vars.Size()!=x_vars.Rows())
     {
      Print(__FUNCTION__," ",__LINE__," Error Invalid inputs Rows in x_vars not equal to length of y_vars");
      m_error = 1;
      return false;
     }
//---
   matrix U,V;
//---
   ::ResetLastError();
//--- SVD operation
   if(!m_exog.SVD(U,V,m_singularvalues))
     {
      Print(__FUNCTION__," ",__LINE__," SVD operation failed : ",::GetLastError());
      m_error = 1;
      return false;
     }
//--- Compute the pseudo-inverse of a matrix by the Moore-Penrose method
   m_pinv = m_exog.PInv();
//--- transpose m_pinv
   matrix pinv_t = m_pinv.Transpose();
//--- normalize covariance matrix
   m_norm_cov_params = matmul(m_pinv,pinv_t);
//---
   matrix diag;
   if(!diag.Diag(m_singularvalues))
     {
      Print(__FUNCTION__," ",__LINE__," Diag operation failed : ",::GetLastError());
      m_error = 1;
      return false;
     }
//--- rank of diag
   m_rank = diag.Rank();
//--- estimate parameters of the model
   m_params = matmul(m_pinv,m_endog);
//--- caluclate its degrees of freedom
   m_model_dof = m_rank - m_kconstant;
//--- as well as those of the residuals
   m_resid_dof = m_obs - m_rank;
//---get the residuals
   m_resid = m_endog - matmul(m_exog,m_params);
//---compute scale of model
   scale();
//--- compute covariance matrix
   covariance_matrix();
//--- compute standard errors
   bse();
//--- compute sum of squared errors
   sse();
//--- compute loglikelihood function
   llf();
//--- compute centered sum of squares
   centeredtss();
//--- compute uncentered sum of squares
   uncenteredtss();
//--- compute r-squared
   rsqe();
//--- compute AIC
   aic();
//--- compute BIC
   bic();
//--- compute tvalues
   tvalues();
//---
   return true;
  }
//+------------------------------------------------------------------+
//| Predict next value based on model parameters                     |
//+------------------------------------------------------------------+
double OLS::Predict(vector &inputs)
  {
//---
   if(m_error)
     {
      Print("Invalid model");
      return EMPTY_VALUE;
     }
//---
   if(inputs.Size()!=(m_params.Size()-m_kconstant))
     {
      Print("invalid inputs: supplied vector does not match size of model parameters");
      return EMPTY_VALUE;
     }
//---
   double prediction = 0;
//---   
   for(ulong i = 0,k = 0;i<m_const_cols.Size(); i++)
      {
       if(!m_const_cols[i])
          prediction+=(m_params[i]*inputs[k++]); 
       else
          prediction+=(m_params[i]);   
      }           
//---
   return prediction;
  }
//+------------------------------------------------------------------+
//| Predict next value based on model parameters                     |
//+------------------------------------------------------------------+
double OLS::Predict(double _input)
  {
//---
   if(m_error || (m_params.Size()-m_kconstant)>1)
     {
      if(m_error)
         Print("Invalid model");
      else
         Print("invalid inputs: insufficient number of predictors supplied");   
      return EMPTY_VALUE;
     }
//---
  double prediction = 0;
//---
  if(m_kconstant)
      prediction=(m_const_cols[0])?m_params[0]+(m_params[1]*_input):m_params[1]+(m_params[0]*_input);
  else
      prediction=m_params[0]*_input;    
//---
   return prediction;
  }
//+------------------------------------------------------------------+
//|Count the number of constants in RHS matrix                       |
//+------------------------------------------------------------------+
ulong OLS::countconstants(void)
  {
   ulong count = 0;
   vector temp;
   
   m_const_cols = vector::Zeros(m_exog.Cols());
   
   for(ulong i =0; i<m_exog.Cols(); i++)
     {
      temp = m_exog.Col(i);
      if((temp.Max()-temp.Min()) == 0.0 && temp.Max() == 1.0)
         {
          count++;
          m_const_cols[i]=1.0;
         }

     }
     
   return count;
  }
//+------------------------------------------------------------------+
//|scale factor for the covariance matrix                            |
//+------------------------------------------------------------------+
void OLS::scale(void)
  {
   m_scale =  m_resid.Dot(m_resid)/double(m_resid_dof);
  }
//+------------------------------------------------------------------+
//| the variance/covariance matrix                                   |
//+------------------------------------------------------------------+
void OLS::covariance_matrix(void)
  {
//---
   m_cov_params=m_norm_cov_params*m_scale;
  }
//+------------------------------------------------------------------+
//| standard errors of the parameter estimates                       |
//+------------------------------------------------------------------+
void OLS::bse(void)
  {
   m_bse=m_cov_params.Diag();
   m_bse=MathSqrt(m_bse);
  }
//+------------------------------------------------------------------+
//| sum of squared errors                                            |
//+------------------------------------------------------------------+
void OLS::sse(void)
  {
   vector sqresid=MathPow(m_resid,2);

   m_sse = sqresid.Sum();

  }
//+------------------------------------------------------------------+
//|likelihood function for the OLS model                             |
//+------------------------------------------------------------------+
void OLS::llf(void)
  {
   double obs2=double(m_obs)/2.0;
   m_llf = -obs2*log(2.0*M_PI) - obs2*log(m_sse / double(m_obs)) - obs2;
//m_llf = -1*obs2*MathLog(2.0*M_PI*m_scale) - m_sse/(2.0*m_scale);
  }
//+------------------------------------------------------------------+
//| Akaike's information criteria                                    |
//+------------------------------------------------------------------+
void OLS::aic(void)
  {
   double params = double(m_model_dof)+double(m_kconstant);

   m_aic = -2*m_llf+2*params;

  }
//+------------------------------------------------------------------+
//|Bayes' information criteria                                       |
//+------------------------------------------------------------------+
void OLS::bic(void)
  {
   double params = double(m_model_dof)+double(m_kconstant);

   m_bic = -2*m_llf+MathLog(m_obs)*params;

  }
//+------------------------------------------------------------------+
//|t-statistic for a given parameter estimate                        |
//+------------------------------------------------------------------+
void OLS::tvalues(void)
  {
   m_tvalues = m_params/m_bse;
  }
//+------------------------------------------------------------------+
//|total  sum of squares centered about the mean                     |
//+------------------------------------------------------------------+
void OLS::centeredtss(void)
  {
   vector centered_endog = m_endog - m_endog.Mean();
//---
   m_centeredtss = centered_endog.Dot(centered_endog);
  }
//+------------------------------------------------------------------+
//| The sum of the squared values of the endogenous response variable|
//+------------------------------------------------------------------+
void OLS::uncenteredtss(void)
  {
   m_uncenteredtss = m_endog.Dot(m_endog);
  }
//+------------------------------------------------------------------+
//|R-squared of the model                                            |
//+------------------------------------------------------------------+
void OLS::rsqe(void)
  {
   m_rsqe = (m_kconstant)? 1 - m_sse/m_centeredtss: 1 - m_sse/m_uncenteredtss;
  }
//+------------------------------------------------------------------+