//+------------------------------------------------------------------+
//|                                                          OLS.mqh |
//|                                  Copyright 2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2023, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<Math\Stat\T.mqh>
#include<Math\Stat\Normal.mqh>
//+------------------------------------------------------------------+
//| Covariance type enumeration                                      |
//+------------------------------------------------------------------+
enum ENUM_COV_TYPE
  {
   COV_NON_ROBUST = 0,
   COV_HC0
  };
//+------------------------------------------------------------------+
//| 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 implementing classical Ordinary Least Squares (OLS) linear regression estimation.
// Uses Singular Value Decomposition (SVD) and Moore-Penrose pseudo-inverses for numerical stability.
class OLS
  {
private:
   bool              m_const_prepended;    // Flag indicating if a baseline constant column is tracked in the matrix
   matrix m_exog,               // Design Matrix ($X$ matrix of independent exogenous predictors)
          m_pinv,               // Calculated Moore-Penrose pseudo-inverse of the design matrix ($X^+$)
          m_cov_params,         // Scaled covariance matrix of the calculated parameters variance limits
          m_m_error,            // Fallback sentinel error matrix populated during failure states
          m_norm_cov_params,    // Normalized parameters covariance matrix base: $(X^T X)^{-1}$
          m_cov_params_default; // Default covariance matrix
   vector m_endog,              // Dependent response variables vector ($Y$ vector of endogenous data)
          m_weights,            // Weights profile configured for data vectors (initialized to 1s)
          m_singularvalues,     // Singular values diagonal array populated by the SVD decomposition pass
          m_params,             // Computed regression coefficients/slopes result vector ($\beta$ solution)
          m_tvalues,            // T-statistic values validating individual parameter statistical significance
          m_bse,                // Standard errors associated with calculated parameter estimates
          m_const_cols,         // Vector masking which specific columns serve as constant intercept vectors
          m_resid,              // Prediction residual error vector: $e = Y - \hat{Y}$
          m_pvalues;            // Parameter pvalues
   ulong  m_obs,                // Sample count: Number of observations rows in the matrix
          m_model_dof,          // Degrees of freedom of the model features (rank minus intercept tracking)
          m_resid_dof,          // Degrees of freedom of the calculated residuals ($N - \text{rank}$)
          m_kconstant,          // Number of constant intercept columns detected in the matrix model
          m_rank;               // Computed mathematical rank of the design matrix
   ENUM_COV_TYPE     m_cov_type;    // Covariance calculation method
   double m_aic,                // Akaike Information Criterion assessing comparative specification quality
          m_bic,                // Bayesian Information Criterion adding stricter parameters penalty layers
          m_scale,              // Variance scale parameter estimate of the model error component ($\sigma^2$)
          m_llf,                // Log-Likelihood function estimation value under normality assumptions
          m_sse,                // Sum of Squared Errors / Residual Sum of Squares (RSS)
          m_rsqe,               // Coefficient of Determination ($R^2$) indicating variance explanation rates
          m_centeredtss,        // Centered Total Sum of Squares tracking variance deviations about the mean
          m_uncenteredtss;      // Uncentered Total Sum of Squares used when no intercept constant exists
   uint              m_error;              // Boolean/integer indicator flag marking invalid calculation states

   // Private computation pipeline steps
   ulong              countconstants(void);   // Scans columns to locate and mask intercept structures
   void               scale(void);            // Computes residual variance scale factor
   void               sse(void);              // Accumulates squared error metric values
   void               rsqe(void);             // Computes the fit determination score ($R^2$)
   void               centeredtss(void);      // Computes variation tracking around response mean values
   void               uncenteredtss(void);    // Calculates total raw coordinate variations
   void               aic(void);              // Evaluates Akaike information scores
   void               bic(void);              // Evaluates Bayesian information scores
   void               bse(void);              // Extracts coefficient standard deviations profiles
   void               llf(void);              // Evaluates Gaussian log likelihood values
   void               tvalues(void);          // Scales parameter estimations by standard deviations
   void               covariance_matrix(void);// Scales normal distributions by the variance metric
   matrix             cov_params(void);       // Calculate covariance matrix
   matrix             hccm(vector &scale);    // Covariance helper
   matrix             cov_hco(void);          // Heteroscedasticity robust covariance matrix.
   vector             pvalues(bool use_t);          // Pvalue calculation

public:
   // Constructor & Destructor
                     OLS(void);
                    ~OLS(void);

   // Public Execution Controllers
   bool              Fit(vector &y_vars,matrix &x_vars,ENUM_COV_TYPE cov_type = COV_NON_ROBUST,bool use_t = false); // Core fitting routine executing regression logic
   double            Predict(vector &inputs);            // Generates predictions for multivariate input arrays
   double            Predict(double _input);             // Generates predictions for univariate input models

   // Public Property Accessors / Getters
   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; }
   vector            Pvalues(void) { if(m_error) return m_m_error.Col(0); return m_pvalues; }
  };

//+------------------------------------------------------------------+
//|  constructor                                                     |
//+------------------------------------------------------------------+
OLS::OLS(void)
  {
   m_kconstant = 0;
   m_obs = 0;
   m_model_dof = 0;
   m_resid_dof = 0;
   m_rank = 0;
   m_aic = 0;
   m_bic = 0;
   m_scale = 0;
   m_llf   = 0;
   m_error = 1; // Default to error state until data is successfully fitted

// Instantiate an fallback matrix populated with error placeholders
   m_m_error.Resize(2,2);
   m_m_error.Fill(EMPTY_VALUE);
  }

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

//+------------------------------------------------------------------+
//| Fit data                                                         |
//+------------------------------------------------------------------+
// Executes parameter estimation using linear algebraic decomposition techniques
bool OLS::Fit(vector &y_vars,matrix &x_vars,ENUM_COV_TYPE cov_type = COV_NON_ROBUST,bool use_t = false)
  {
// Synchronize raw input objects to internal structural data slots
   m_endog = y_vars;
   m_exog  = x_vars;
   m_cov_type = cov_type;
// Identify if an intercept constant is bundled into the independent matrix configurations
   m_kconstant = countconstants();
   m_error = 0;

// Reset structural state counters
   m_model_dof = m_resid_dof = m_rank = m_obs= 0;
   m_aic = m_bic= m_llf=0;

   m_obs=m_exog.Rows();          // Number of observation rows available for modeling
   m_weights = vector::Ones(m_obs); // Standard default OLS tracks even weighting profiles
   m_rank = m_exog.Rank();       // Extract initial matrix rank profile configuration counts
   m_model_dof = m_rank - m_kconstant;
   m_resid_dof = m_obs - m_rank;

// Enforce basic matrix dimension validation checks
   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();

// Execute SVD factor processing to assert geometric rank configurations stably
   if(!m_exog.SVD(U,V,m_singularvalues))
     {
      Print(__FUNCTION__," ",__LINE__," SVD operation failed : ",::GetLastError());
      m_error = 1;
      return false;
     }

// Extract the Moore-Penrose pseudo-inverse matrix: $X^+ = (X^T X)^{-1} X^T$
   m_pinv = m_exog.PInv();
   matrix pinv_t = m_pinv.Transpose();

// Compute normalized parameter covariances layout: $X^+ (X^+)^T = (X^T X)^{-1}$
   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;
     }

   m_rank = diag.Rank(); // Update matrix structural rank estimations via non-zero diagonal entries

// Compute parameter coefficients vector: $\hat{\beta} = X^+ Y$
   m_params = matmul(m_pinv,m_endog);

// Recompute final safe degrees of freedom allocations post-rank evaluation
   m_model_dof = m_rank - m_kconstant;
   m_resid_dof = m_obs - m_rank;

// Isolate model tracking residuals: $e = Y - X\hat{\beta}$
   m_resid = m_endog - matmul(m_exog,m_params);
// Calculate defaulte covariance
   switch(m_cov_type)
     {
      case COV_NON_ROBUST:
         m_cov_params_default = matrix::Zeros(0,0);
         break;
      case COV_HC0:
         m_cov_params_default = cov_hco();
         break;
     }
// Execute downstream statistical attribute calculation routines
   scale();
   covariance_matrix();
   bse();
   sse();
   llf();
   centeredtss();
   uncenteredtss();
   rsqe();
   aic();
   bic();
   tvalues();
   pvalues(use_t);

   return true;
  }

//+------------------------------------------------------------------+
//| Predict next value based on model parameters                     |
//+------------------------------------------------------------------+
// Evaluates dot product estimations while maintaining spatial transformations for masked constants
double OLS::Predict(vector &inputs)
  {
   if(m_error)
     {
      Print("Invalid model");
      return EMPTY_VALUE;
     }

// Ensure supplied non-constant features align with parameter dimensions
   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;

// Conditional dot-product tracking that maps slope parameters onto raw inputs or accumulates constant values
   for(ulong i = 0, k = 0; i<m_const_cols.Size(); i++)
     {
      if(!m_const_cols[i])
         prediction+=(m_params[i]*inputs[k++]); // Variable coordinate projection multiplication pass
      else
         prediction+=(m_params[i]);            // Intercept accumulation calculation pass
     }

   return prediction;
  }

//+------------------------------------------------------------------+
//| Predict next value based on model parameters                     |
//+------------------------------------------------------------------+
// Shortcut helper generating point calculations on single feature lines (Simple Linear Regression)
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)
      // Apply single feature multiplication adjusting indexes based on intercepted location patterns
      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                       |
//+------------------------------------------------------------------+
// Scans columns where the difference between max and min values is zero, indicating an invariant intercept vector
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);
      // Columns with zero variation and a maximum value of 1.0 are identified as constants
      if((temp.Max()-temp.Min()) == 0.0 && temp.Max() == 1.0)
        {
         count++;
         m_const_cols[i]=1.0; // Mask this column position index as an intercept
        }
     }

   return count;
  }

//+------------------------------------------------------------------+
//|scale factor for the covariance matrix                            |
//+------------------------------------------------------------------+
// Evaluates residual variance factor ($\sigma^2$): $\frac{\sum e_i^2}{\text{DOF}_{\text{residual}}}$
void OLS::scale(void)
  {
   m_scale =  m_resid.Dot(m_resid)/double(m_resid_dof);
  }

//+------------------------------------------------------------------+
//| the variance/covariance matrix                                   |
//+------------------------------------------------------------------+
// Scales the normalized parameter covariance matrix to capture the error variance structure
void OLS::covariance_matrix(void)
  {
   m_cov_params=m_norm_cov_params*m_scale;
  }

//+------------------------------------------------------------------+
//|Compute the variance/covariance matrix.                           |
//+------------------------------------------------------------------+
matrix OLS::cov_params(void)
  {
//if(!cov_p.Rows())
//  {
   matrix cov_p;
   if(m_cov_params_default.Rows())
      cov_p = m_cov_params_default;
   else
     {
      //if(scale == EMPTY_VALUE)
      //   scale = m_scale;
      cov_p = m_norm_cov_params * m_scale;
     }
   /*  }
   if(column.Size())
     {
      if(column.Size()==1)
        {
         matrix out = matrix::Zeros(1,1);
         ulong index = ulong(column[0]);
         out[0,0] = cov_p[index,index];
         return out;
        }
      else
        {
         matrix out = matrix::Zeros(column.Size(),column.Size());
         for(ulong i = 0; i<out.Rows(); ++i)
            for(ulong j = 0; j<out.Cols(); ++j)
               out[i,j] = cov_p[ulong(fabs(column[i])),ulong(fabs(column[j]))];
         return out;
        }
     }
   else
      if(r_matrix.Rows())
        {
         if(!other.Rows())
            other = r_matrix;
         return r_matrix.MatMul(cov_p.MatMul(other.Transpose()));
        }
      else*/
   return cov_p;
  }
//+------------------------------------------------------------------+
//|Compute Heteroscedasticity robust covariance matrix.              |
//+------------------------------------------------------------------+
matrix OLS::cov_hco(void)
  {
   vector het_scale = pow(m_resid,2.);
   return hccm(het_scale);
  }
//+------------------------------------------------------------------+
//|Covariance matrix calculation helper                              |
//+------------------------------------------------------------------+
matrix OLS::hccm(vector &scale)
  {
   matrix scale_mat = matrix::Zeros(m_pinv.Cols(),m_pinv.Rows());
   for(ulong i = 0; i<m_pinv.Rows(); ++i)
      scale_mat.Col(scale*m_pinv.Row(i),i);
   return m_pinv.MatMul(scale_mat);
  }
//+------------------------------------------------------------------+
//| standard errors of the parameter estimates                       |
//+------------------------------------------------------------------+
// Extracts standard errors from the square root of the main diagonal elements of the covariance matrix
void OLS::bse(void)
  {
   matrix cp = cov_params();
   m_bse=cp.Diag();
   m_bse=MathSqrt(m_bse);
  }

//+------------------------------------------------------------------+
//| sum of squared errors                                            |
//+------------------------------------------------------------------+
// Calculates total residual error metrics squared: $\sum e_i^2$
void OLS::sse(void)
  {
   vector sqresid=MathPow(m_resid,2);
   m_sse = sqresid.Sum();
  }

//+------------------------------------------------------------------+
//|likelihood function for the OLS model                             |
//+------------------------------------------------------------------+
// Computes log-likelihood under standard normal error assumptions:
// $\ln L = -\frac{N}{2}\ln(2\pi) - -\frac{N}{2}\ln(\frac{\text{SSE}}{N}) - \frac{N}{2}$
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;
  }

//+------------------------------------------------------------------+
//| Akaike's information criteria                                    |
//+------------------------------------------------------------------+
// Computes the information quality metric with parameters scale penalty: $\text{AIC} = -2\ln L + 2K$
void OLS::aic(void)
  {
   double params = double(m_model_dof)+double(m_kconstant);
   m_aic = -2*m_llf+2*params;
  }

//+------------------------------------------------------------------+
//|Bayes' information criteria                                       |
//+------------------------------------------------------------------+
// Computes the information quality metric with sample size weight scale penalty: $\text{BIC} = -2\ln L + K\ln N$
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                        |
//+------------------------------------------------------------------+
// Computes $t$-scores to check the null hypothesis ($H_0: \beta_i = 0$): $t = \frac{\hat{\beta_i}}{\text{SE}(\hat{\beta_i})}$
void OLS::tvalues(void)
  {
   m_tvalues = m_params/m_bse;
  }

//+------------------------------------------------------------------+
//|total  sum of squares centered about the mean                     |
//+------------------------------------------------------------------+
// Computes centered variation metrics for models containing an intercept constant: $\sum (Y_i - \bar{Y})^2$
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|
//+------------------------------------------------------------------+
// Computes uncentered variation metrics for regression models that pass strictly through the origin: $\sum Y_i^2$
void OLS::uncenteredtss(void)
  {
   m_uncenteredtss = m_endog.Dot(m_endog);
  }
//+------------------------------------------------------------------+
//| Pvalue calculation                                               |
//+------------------------------------------------------------------+
vector OLS::pvalues(bool use_t)
  {
   m_pvalues = vector::Zeros(m_tvalues.Size());
   int error = 0;
   if(use_t)
     {
      for(ulong i = 0; i<m_pvalues.Size(); ++i)
         m_pvalues[i] = (1.0 - (double)MathCumulativeDistributionT((double)fabs(m_tvalues[i]),double(m_resid_dof),error)) * 2.0;
     }
   else
     {
      for(ulong i = 0; i<m_pvalues.Size(); ++i)
         m_pvalues[i] = (1.0 - (double)MathCumulativeDistributionNormal((double)fabs(m_tvalues[i]),double(0.0),double(1.0),error)) * 2.0;
     }
   return m_pvalues;
  }
//+------------------------------------------------------------------+
//|R-squared of the model                                            |
//+------------------------------------------------------------------+
// Computes $R^2 = 1 - \frac{\text{SSE}}{\text{TSS}}$, selecting the appropriate TSS variant depending on intercept presence
void OLS::rsqe(void)
  {
   m_rsqe = (m_kconstant)? 1 - m_sse/m_centeredtss: 1 - m_sse/m_uncenteredtss;
  }
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+