//+------------------------------------------------------------------+
//|                                                EconometricsM.mqh |
//|                             Copyright 2000-2026, MetaQuotes Ltd. |
//|                                                     www.mql5.com |
//+------------------------------------------------------------------+
//---
#include "EconometricsA.mqh"
//---
//+------------------------------------------------------------------+
//| struct for coefficients testing                                  |
//+------------------------------------------------------------------+
struct CoefficientStats
  {
   double            estimate;            // point estimate
   double            std_error;           // standard error
   double            t_stat;              // t-statistic
   double            p_value;             // one-sided p-value
   double            conf_low;            // low bound of confidence interval
   double            conf_high;           // high bound of confidence interval
  };
//+------------------------------------------------------------------+
//| struct for prognosed data                                        |
//+------------------------------------------------------------------+
struct SPrognose
  {
   vector            xnew;                // regressors for prognose (xnew[0]==1!)
   double            point_progn;         // point prognose (mean)
   double            conf_low;            // low bound of confidence interval
   double            conf_high;           // high bound of confidence interval
   double            progn_low;           // low bound of prognose interval
   double            progn_high;          // high bound of prognose interval
   void              print()              // print struct
     {
      PrintFormat("point prognose: %.3f",point_progn);
      PrintFormat("confidence interval, low bound: %.3f, high bound: %.3f",conf_low,conf_high);
      PrintFormat("prognose interval, low bound: %.3f, high bound: %.3f",progn_low,progn_high);
     };
  };
//+------------------------------------------------------------------+
//| struct for prognosed data for TSLS                               |
//+------------------------------------------------------------------+
struct SPrognose_TSLS
  {
   vector            xnew_exo;            // regressors for prognose (xnew_exo[0]==1!)
   vector            xnew_endo;           // regressors for prognose
   vector            znew_iv;             // instruments for prognose
   double            point_progn;         // point prognose (mean)
   double            conf_low;            // low bound of confidence interval
   double            conf_high;           // high bound of confidence interval
   double            progn_low;           // low bound of prognose interval
   double            progn_high;          // high bound of prognose interval
   void              print()              // print struct
     {
      PrintFormat("point prognose: %.3f",point_progn);
      PrintFormat("confidence interval, low bound: %.3f, high bound: %.3f",conf_low,conf_high);
      PrintFormat("prognose interval, low bound: %.3f, high bound: %.3f",progn_low,progn_high);
     };
  };
//+------------------------------------------------------------------+
//| Computation of parameters and residuals for multiple regression  |
//| Inputs: y - target variable, X - regressors                      |
//| Outputs: b - parameters, e - residuals, c - residuals' SD        |
//+------------------------------------------------------------------+
void regression(vector& y, matrix& X, vector& b, vector& e, double& c)
  {
//--- Input data validation
   ulong k=X.Cols();
   if(k<1)
     {Print("error: empty X");  return;}
   ulong n=X.Rows();
   if(n<k)
     {Print("error: not enough samples"); return;}
   if(X.Rank()<k)
     {Print("error: low rank of X"); return;}
   if(y.Size()!=n)
     {Print("error: wrong size of y"); return;}
//--- Parameters and residuals computation
   b=X.Transpose().MatMul(X).Inv().MatMul(X.Transpose()).MatMul(y);
   e=y-X.MatMul(b);
   c=e.Std((int)k);
  }
//+------------------------------------------------------------------+
//| Computation of parameters and residuals for multiple regression  |
//| Inputs: y - target variable, X - regressors                      |
//| Outputs: b - parameters, e - residuals, c - residuals' SD        |
//|          XX - Inverse(Transpose(X)*X) (need for stat calculation)|
//+------------------------------------------------------------------+
void regression(vector& y, matrix& X, vector& b, vector& e, double& c, matrix& XX)
  {
//--- Input data validation
   ulong k=X.Cols();
   if(k<1)
     {Print("error: empty X");  return;}
   ulong n=X.Rows();
   if(n<k)
     {Print("error: not enough samples"); return;}
   if(X.Rank()<k)
     {Print("error: low rank of X"); return;}
   if(y.Size()!=n)
     {Print("error: wrong size of y"); return;}
//--- Parameters and residuals computation
   XX=X.Transpose().MatMul(X).Inv();
   b=XX.MatMul(X.Transpose()).MatMul(y);
   e=y-X.MatMul(b);
   c=e.Std((int)k);
  }
//+------------------------------------------------------------------+
//| Prognose computation                                             |
//| Inputs: XX - Inverse(Transpose(X)*X), b - parameters vector,     |
//|         c - unbiased residuals SD, n - sample size               |
//|         conf_level - confidence level for interval estimation    |
//|         prgns.xnew - point for prognose (prgns.xnew[0]==1!)      |
//| Outputs: prgns - point and interval prognoses for prgns.xnew     |
//+------------------------------------------------------------------+
void prognose(matrix& XX, vector& b, double c, ulong n, SPrognose& prgns, double conf_level=0.95)
  {
//--- Input data validation
   ulong k=XX.Cols();
//--- Input data validation
   if(k<1)
     {Print("error: empty X");  return;}
   if(b.Size()!=k)
     {Print("error: wrong b size");  return;}
   if(prgns.xnew.Size()!=k)
     {Print("error: wrong prgns.xnew size");  return;}
   prgns.point_progn=b.MatMul(prgns.xnew);
   double t,h=XX.MatMul(prgns.xnew).MatMul(prgns.xnew), d;
   int err=0;
//--- t = MathQuantileT((1.0 - conf_level) / 2.0, n - k, err);  // faulty library function
   t = MathQuantileT_TMP((1.0 - conf_level) / 2.0, n - k, err); // custom implementation
   if(err != 0)
     {Print("MathQuantileT() error ", err); return;}
   t=MathAbs(t);
   d=t*c*MathSqrt(h);
   prgns.conf_low=prgns.point_progn-d;
   prgns.conf_high=prgns.point_progn+d;
   d=t*c*MathSqrt(1.0+h);
   prgns.progn_low=prgns.point_progn-d;
   prgns.progn_high=prgns.point_progn+d;
  }
//+------------------------------------------------------------------+
//| Computation of parameters stat                                   |
//| Inputs: XX - Inverse(Transpose(X)*X), b - parameters vector,     |
//|         c - unbiased residuals SD, n - sample size               |
//|         conf_level - confidence level for interval estimation    |
//| Outputs: cs - parameters statistics                              |
//+------------------------------------------------------------------+
void parameter_stat(matrix& XX, vector& b, double c, ulong n, CoefficientStats& cs[], double conf_level=0.95)
  {
//--- Input data validation
   ulong k=XX.Cols();
   if(k<1)
     {Print("error: empty X");  return;}
   vector d=XX.Diag();
   ArrayResize(cs,(int)k);
   double t;
   int err=0;
//--- t = MathQuantileT((1.0 - conf_level) / 2.0, n - k, err);  // faulty library function
   t = MathQuantileT_TMP((1.0 - conf_level) / 2.0, n - k, err); // custom implementation
   if(err != 0)
     {Print("MathQuantileT() error ", err); return;}
   t=MathAbs(t);
   for(int i=0;i<(int)k;++i)
     {
      cs[i].estimate=b[i];
      cs[i].std_error=c*MathSqrt(d[i]);
      cs[i].t_stat=cs[i].estimate/cs[i].std_error;
      cs[i].p_value=MathCumulativeDistributionT(MathAbs(cs[i].t_stat),n-k,false,false,err);
      if(err!=0)
        {Print("error: MathCumulativeDistributionT() error ",err); return;}
      cs[i].conf_low=cs[i].estimate-t*cs[i].std_error;
      cs[i].conf_high=cs[i].estimate+t*cs[i].std_error;
     }
  }
//+------------------------------------------------------------------+
//| Computation of parameters and residuals for                      |
//| simple linear regression                                         |
//+------------------------------------------------------------------+
void regression1(vector& y, vector& x, double& a, double& b, double& c, vector& e)
  {
   double ya[],xa[],ea[];
   vector2array(y,ya);
   vector2array(x,xa);
   regression1(ya,xa,a,b,c,ea);
   e.Assign(ea);
  }
//+------------------------------------------------------------------+
//| Durbin-Wu-Hausman (DWH) Endogeneity Test (Regression-Based)      |
//| H0: Regressors are exogenous (OLS is consistent and efficient)   |
//| H1: Regressors are endogenous (OLS is biased, 2SLS is required)  |
//| Returns: p-value of the test (p-value < 0.05 means endogeneity)  |
//+------------------------------------------------------------------+
double DWH_test(const vector &y, const matrix &X_exo, const matrix &X_endo, const matrix &Z_iv)
  {
   ulong n = y.Size();
   ulong p = X_exo.Cols();
   ulong r = X_endo.Cols();
   ulong m = Z_iv.Cols();

//--- 1. Input data validation
   if(n < 1 || X_exo.Rows() != n || X_endo.Rows() != n || Z_iv.Rows() != n)
     {
      Print("DWH Test error: matrix row dimensions mismatch or empty data.");
      return 1.0;
     }
   if(m < r)
     {
      Print("DWH Test error: too few instruments (m < r). Test cannot be performed.");
      return 1.0;
     }

//--- 2. Construct the full instrument matrix Z = [X_exo | Z_iv]
   matrix Z = X_exo.Concat(Z_iv, 1);
   if(Z.Rank() < Z.Cols())
     {
      Print("DWH Test error: multicollinearity in the instrument matrix Z.");
      return 1.0;
     }

//--- 3. FIRST STAGE: Regress each endogenous variable on all instruments to get residuals
   matrix V_hat;
   V_hat.Resize(n, r);

   matrix Pl = Z.Transpose().MatMul(Z).Inv().MatMul(Z.Transpose());
   vector bl, el, ytmp;

   for(ulong i = 0; i < r; ++i)
     {
      ytmp = X_endo.Col(i);
      bl = Pl.MatMul(ytmp);       //--- Coefficients of the first stage
      el = ytmp - Z.MatMul(bl);    //--- Residuals containing the "toxic" endogenous part
      V_hat.Col(el, i);            //--- Save residuals as a column
     }

//--- 4. SECOND STAGE: Fit the Unrestricted (Augmented) Model via OLS: Y ~ X_exo + X_endo + V_hat
   matrix X_orig = X_exo.Concat(X_endo, 1);
   matrix X_augmented = X_orig.Concat(V_hat, 1);

   ulong k_aug = X_augmented.Cols(); // Total variables in the augmented model (p + r + r)

//--- Validate sample size after augmentation
   if(n <= k_aug)
     {
      Print("DWH Test error: too few observations for the augmented model.");
      return 1.0;
     }

   matrix X_aug_T = X_augmented.Transpose();
   matrix XX_aug_inv = X_aug_T.MatMul(X_augmented).Inv();
   vector b_aug = XX_aug_inv.MatMul(X_aug_T).MatMul(y);

//--- Calculate Sum of Squares for the Unrestricted Model (RSS_unrestricted)
   vector e_aug = y - X_augmented.MatMul(b_aug);
   double rss_unrestricted = e_aug.MatMul(e_aug);

//--- 5. Fit the Restricted Model (Standard OLS): Y ~ X_exo + X_endo
   matrix X_orig_T = X_orig.Transpose();
   vector b_ols = X_orig_T.MatMul(X_orig).Inv().MatMul(X_orig_T).MatMul(y);

//--- Calculate Sum of Squares for the Restricted Model (RSS_restricted)
   vector e_ols = y - X_orig.MatMul(b_ols);
   double rss_restricted = e_ols.MatMul(e_ols);

//--- 6. Compute the F-statistic for the joint significance of V_hat coefficients
//--- Degrees of freedom: numerator = r (number of restrictions), denominator = n - k_aug
   double h_numerator = (rss_restricted - rss_unrestricted) / (double)r;
   double h_denominator = rss_unrestricted / (double)(n - k_aug);

   if(h_denominator <= 0.0)
     {
      Print("DWH Test warning: perfect fit or zero residual variance in the augmented model.");
      return 1.0;
     }

   double hausman_F = h_numerator / h_denominator;
   int h_err = 0;

//--- Compute the right-tailed p-value from the Fisher F-distribution
   double f_cdf = MathCumulativeDistributionF(hausman_F, (double)r, (double)(n - k_aug), h_err);
   double p_value = 1.0 - f_cdf;

//--- 7. Print diagnostic results to the terminal log
   PrintFormat("DWH Test | F-stat: %.4f | p-value: %.6f", hausman_F, p_value);
   if(p_value < 0.05)
      Print("DWH Test Result: Reject H0. Endogeneity is significant. Use 2SLS.");
   else
      Print("DWH Test Result: Fail to reject H0. No significant endogeneity. OLS is preferred.");

   return p_value;
  }
//+------------------------------------------------------------------+
//| Two stage least square algorithm (TSLS)                          |
//+------------------------------------------------------------------+
void TSLS(vector& y, matrix& X_exo, matrix& X_endo, matrix& Z_iv,
          SPrognose_TSLS& prog, CoefficientStats& stats[], double confidence_level = 0.95)
  {
//--- Extract data dimensions: n - observations, p - exogenous, r - endogenous, m - clean instruments
   ulong n=y.Size(), q=X_exo.Cols(), r=X_endo.Cols(), m=Z_iv.Cols();

//--- Input data validation
   if(n<1)
     {Print("TSLS error: empty y");  return;}
   if(X_exo.Rows()!=n)
     {Print("TSLS error: wrong X_exo rows number");  return;}
   if(X_endo.Rows()!=n)
     {Print("TSLS error: wrong X_endo rows number");  return;}
   if(Z_iv.Rows()!=n)
     {Print("TSLS error: wrong Z_iv rows number");  return;}
   if(n <= (q + r) || n <= (q + m))
     {Print("TSLS error: too few observations (n) for the number of variables"); return;}
   if(m<r)
     {Print("TSLS error: too few instrumental variables");  return;}
   if(X_exo.Rank()<q)
     {Print("TSLS error: multicollinearity in the X_exo matrix");  return;}

//--- Validate dimensional consistency for the out-of-sample forecast fields inside the structure
   if(prog.xnew_exo.Size()!=q)
     {Print("TSLS error: wrong out-of-sample prog.xnew_exo vector size"); return;}
   if(prog.znew_iv.Size()!=m)
     {Print("TSLS error: wrong out-of-sample prog.znew_iv vector size"); return;}

//--- Z - regressors for the first TSLS stage (combined exogenous and clean instruments)
   matrix Z=X_exo.Concat(Z_iv,1);
   if(Z.Rank()<Z.Cols())
     {Print("TSLS error: multicollinearity in the Z matrix");  return;}

//--- X - regressors for the second TSLS stage (initialized with exogenous variables)
   matrix X=X_exo;
   X.Resize(n,q+r);

//--- Initialize the first-stage coefficient matrix locally to use later in forecasting
   matrix B_first;
   B_first.Resize(Z.Cols(), r);

//--- Precalculate projection matrices for the full and restricted models to optimize performance
   matrix Pl=Z.Transpose().MatMul(Z).Inv().MatMul(Z.Transpose());
   matrix Ps=X_exo.Transpose().MatMul(X_exo).Inv().MatMul(X_exo.Transpose());
   vector bl,bs,el,es,ytmp;
   double F=0.0,Fmin=10.0,rl,rs;

//--- First TSLS stage: instrumenting each endogenous variable and checking instrument relevance
   for(ulong i=0;i<r;++i)
     {
      //--- Extract the current endogenous regressor
      ytmp=X_endo.Col(i);

      //--- Fit the full model using all instruments (Z)
      bl=Pl.MatMul(ytmp);
      el=ytmp-Z.MatMul(bl);
      rl=el.MatMul(el);

      //--- Fit the restricted model using only exogenous variables (X_exo)
      bs=Ps.MatMul(ytmp);
      es=ytmp-X_exo.MatMul(bs);
      rs=es.MatMul(es);

      //--- Calculate the F-statistic to test the joint significance of clean instruments
      F=(rs/rl-1.0)*((double)(n-q-m)/m);
      if(F<Fmin)
        {PrintFormat("TSLS error: instruments are weak for endogenous regressor %d. F = %.2f",i+1,F); return;}

      //--- Store the first-stage coefficients for this endogenous variable as a column
      B_first.Col(bl, i);

      //--- Store the predicted (cleaned) endogenous variable into the second-stage regressor matrix
      X.Col(Z.MatMul(bl),q+i);
     }

//--- Second TSLS stage: estimating final coefficients and calculating statistics

//--- 1. Check the rank of the cleaned regressor matrix X before inversion to avoid singularity
   ulong total_cols = q + r;
   if(X.Rank() < total_cols)
     {
      Print("TSLS error: multicollinearity in the cleaned X matrix on the second stage");
      return;
     }

//--- 2. Compute the inverted matrix of the second stage (saved for standard error calculation)
   matrix XT = X.Transpose();
   matrix XX_inv = XT.MatMul(X).Inv();

//--- 3. Calculate the final point estimates (beta vector)
   vector b_vec = XX_inv.MatMul(XT).MatMul(y);

//--- 4. Calculate the true residual variance (sigma_sq) using the original uncleaned regressors
   matrix X_orig = X_exo.Concat(X_endo, 1);

//--- 5. Compute true model residuals based on actual historical data
   vector e = y - X_orig.MatMul(b_vec);

//--- 6. Calculate degrees of freedom (N - K) and residual variance
   long df = (long)n - (long)total_cols;
   double sigma_sq = (e.MatMul(e)) / (double)df;

//--- 7. Compute the coefficient covariance matrix V_beta (uses second-stage XX_inv matrix)
   matrix V_beta = XX_inv * sigma_sq;
   vector variances = V_beta.Diag();

//--- 8. Prepare the critical t-value for the specified confidence level (two-sided interval)
   double alpha = 1.0 - confidence_level;
   int err_code = 0;
   double t_crit = MathQuantileT(1.0 - (alpha / 2.0), (double)df, err_code);
   if(err_code!=0)
     {Print("MathQuantileT() error: ",err_code,", MathQuantileT_TMP() trying"); err_code=0; t_crit = MathQuantileT_TMP(1.0 - (alpha / 2.0), (double)df, err_code);}
   if(err_code!=0)
     {Print("MathQuantileT_TMP() error: ",err_code); return;}

//--- 9. Resize the output stats array and populate metrics for each regressor
   ArrayResize(stats, (int)total_cols);

   for(ulong i = 0; i < total_cols; ++i)
     {
      //--- Store the point estimate and calculate the inflated standard error
      stats[i].estimate  = b_vec[i];
      stats[i].std_error = MathSqrt(variances[i]);

      //--- Calculate the t-statistic and its corresponding one-sided p-value
      if(stats[i].std_error > 0)
        {
         stats[i].t_stat = stats[i].estimate / stats[i].std_error;

         //--- MathCumulativeDistributionT returns the left-tail probability P(T <= t)
         double p_cumulative = MathCumulativeDistributionT(MathAbs(stats[i].t_stat), (double)df, err_code);
         stats[i].p_value = 1.0 - p_cumulative;
        }
      else
        {
         stats[i].t_stat = 0.0;
         stats[i].p_value = 1.0;
        }

      //--- Calculate the confidence interval bounds for the coefficient
      stats[i].conf_low  = stats[i].estimate - (t_crit * stats[i].std_error);
      stats[i].conf_high = stats[i].estimate + (t_crit * stats[i].std_error);
     }

//--- Out-of-sample forecasting stage using the instrumented approach

//--- 1. Construct the complete out-of-sample instrument vector Z0 = [xnew_exo | znew_iv]
   vector z0 = prog.xnew_exo.Concat(prog.znew_iv);

//--- 2. Clean the future: project the future endogenous variables using first-stage coefficients
   prog.xnew_endo = z0.MatMul(B_first);

//--- 3. Form the final cleaned out-of-sample regressor vector for the second-stage equation
   vector x0_hat = prog.xnew_exo.Concat(prog.xnew_endo);

//--- 4. Calculate the symmetric point forecast (mean prediction) using vector dot product
   prog.point_progn = x0_hat.MatMul(b_vec);

//--- 5. Calculate the variance scale factor (g-factor) for the out-of-sample bar
   matrix X0_mat(1, x0_hat.Size());
   X0_mat.Row(x0_hat, 0);
   matrix X0_T = X0_mat.Transpose();
   matrix shift_mat = X0_mat.MatMul(XX_inv).MatMul(X0_T);

//--- Extract the scalar value from the 1x1 matrix using explicit indexing
   double g_factor = shift_mat[0][0];

//--- 6. Calculate the standard errors for both the confidence interval and the prediction interval
   double se_conf = MathSqrt(sigma_sq * g_factor);
   double se_pred = MathSqrt(sigma_sq * (1.0 + g_factor));

//--- 7. Derive half-widths for the interval bounds using the critical t-value
   double conf_width = t_crit * se_conf;
   double pred_width = t_crit * se_pred;

//--- 8. Write upper and lower boundaries into the SPrognose_TSLS structure fields
   prog.conf_low  = prog.point_progn - conf_width;
   prog.conf_high = prog.point_progn + conf_width;
   prog.progn_low  = prog.point_progn - pred_width;
   prog.progn_high = prog.point_progn + pred_width;

  }

//+------------------------------------------------------------------+
//| Computation of the correlation matrix (assumed X1==const)        |
//| Inputs: y - target variable, X - regressors                      |
//| Output: CM - correlation matrix                                  |
//+------------------------------------------------------------------+
void corr_matrix(vector& y, matrix& X, matrix& CM)
  {
//--- Input data validation
   ulong k=X.Cols();
   if(k<1)
     {Print("corr_matrix() error: empty X");  return;}
   ulong n=X.Rows();
   if(y.Size()!=n)
     {Print("corr_matrix() error: wrong size of y"); return;}
   matrix yX=X;
   yX.Col(y,0);
   CM=yX.CorrCoef(false);
  }
//+------------------------------------------------------------------+
//| Residuals plot vs. price bar index                               |
//+------------------------------------------------------------------+
void t_residuals_plot(vector& residuals)
  {
   double e[];
   vector2array(residuals,e);
   t_residuals_plot(e);
  }
//+------------------------------------------------------------------+
//| EPDF of residuals vs. normal density with residual SD            |
//| nofx - number of points on plot                                  |
//+------------------------------------------------------------------+
void epdf_vs_normalpdf(vector& residuals, int nofx = 30)
  {
   double e[];
   vector2array(residuals,e);
   epdf_vs_normalpdf(e,nofx);
  }
//+------------------------------------------------------------------+
//| QQ-plot of residuals vs. normal distribution with residual SD    |
//+------------------------------------------------------------------+
void qq_plot(vector& residuals)
  {
   double e[];
   vector2array(residuals,e);
   qq_plot(e);
  }
//+------------------------------------------------------------------+
//| Correlogram of residuals                                         |
//+------------------------------------------------------------------+
void correlogram(vector& residuals)
  {
   double e[];
   vector2array(residuals,e);
   correlogram(e);
  }
//+------------------------------------------------------------------+
//| Scatter plot of (x, y) points with line y = a * x + b            |
//+------------------------------------------------------------------+
void scatter_plot(vector& x, vector& y, bool add_line = false, double a = 0.0, double b = 0.0)
  {
   double xa[],ya[];
   vector2array(x,xa);
   vector2array(y,ya);
   scatter_plot(xa,ya,add_line,a,b);
  }
//+------------------------------------------------------------------+
//| Scatter plot of (Xi, y) and fitted regression line               |
//| Inputs: X - regressors matrix, y - target variable,              |
//|         i - regressor's index (X column index)                   |
//+------------------------------------------------------------------+
void scatter_plot_Xi_y(matrix& X, ulong i, vector& y)
  {
   if(i>=X.Cols())
     {
      Print("input regressor's index is out of range");
      return;
     }
   double a,b,c;
   vector e;
   regression1(y,X.Col(i),a,b,c,e);
   scatter_plot(X.Col(i),y,true,a,b);
  }
//+------------------------------------------------------------------+
//| Scatter plot of (Xi, Xj) and fitted regression line              |
//| Inputs: X - regressors matrix,                                   |
//|         i, j - regressors' indexes (X columnes indexes)          |
//+------------------------------------------------------------------+
void scatter_plot_Xi_Xj(matrix& X, ulong i, ulong j)
  {
   if(i>=X.Cols()||j>=X.Cols())
     {
      Print("input regressor's index is out of range");
      return;
     }
   double a,b,c;
   vector e;
   regression1(X.Col(j),X.Col(i),a,b,c,e);
   scatter_plot(X.Col(i),X.Col(j),true,a,b);
  }
//+------------------------------------------------------------------+
//| Partial regression plot for Xi                                   |
//| Inputs: X - regressors matrix, y - target variable,              |
//|         i - regressor's index (X column index)                   |
//+------------------------------------------------------------------+
void partial_regression_plot(matrix& X, ulong i, vector& y)
  {
//--- Input data validation
   ulong k=X.Cols();
   if(k<2)
     {
      Print("too few regressors");
      return;
     }
   if(i>=k)
     {
      Print("input regressor's index is out of range");
      return;
     }
   double a,b,c;
   vector x=X.Col(i),ey,ex,e1,bv;
//--- matrix without Xi
   matrix X_1=X;
   if(i!=k-1)
      X_1.Col(X_1.Col(k-1),i);
   X_1.Resize(X_1.Rows(),k-1);
   regression(y,X_1,bv,ey,c);
   regression(x,X_1,bv,ex,c);
   regression1(ey,ex,a,b,c,e1);
   scatter_plot(ex,ey,true,a,b);
  }
//+------------------------------------------------------------------+
//| Scatter plot of (Yfit, Yreal) and line y = x                     |
//| Inputs: X - regressors matrix, b - parameters vector             |
//|         y - real values of dependent variable                    |
//+------------------------------------------------------------------+
void scatter_plot_Yfit_Yreal(matrix& X, vector& b, vector& y)
  {
//--- Input data validation
   ulong k=X.Cols();
   if(k<1)
     {Print("error: empty X");  return;}
   if(b.Size()!=k)
     {Print("error: wrong size of b"); return;}
   ulong n=X.Rows();
   if(y.Size()!=n)
     {Print("error: wrong size of y"); return;}
   scatter_plot(X.MatMul(b),y,true,0.0,1.0);
  }
//+------------------------------------------------------------------+
//| Scatter plot of (Yfit, residuals)                                |
//| Inputs: X - regressors matrix, b - parameters vector             |
//|         e - residuals                                            |
//+------------------------------------------------------------------+
void scatter_plot_Yfit_residuals(matrix& X, vector& b, vector& e)
  {
//--- Input data validation
   ulong k=X.Cols();
   if(k<1)
     {Print("error: empty X");  return;}
   if(b.Size()!=k)
     {Print("error: wrong size of b"); return;}
   ulong n=X.Rows();
   if(e.Size()!=n)
     {Print("error: wrong size of e"); return;}
   scatter_plot(X.MatMul(b),e);
  }
//+------------------------------------------------------------------+
//| R^2, coefficient of determination                                |
//| Inputs: y - real values of dependent variable, e - residuals     |
//+------------------------------------------------------------------+
double R2(vector& y, vector& e)
  {
//--- Input data validation
   ulong n=y.Size();
   if(n<2)
     {Print("R2() error: y too short"); return 0.0;}
   if(e.Size()!=n)
     {Print("R2() error: different sizes of y and e"); return 0.0;}
   double Sy=y.Var(0);
   if(Sy<=DBL_MIN)
     {Print("R2() error: y = const"); return 0.0;}
   return 1.0-e.Var(0)/Sy;
  }
//+------------------------------------------------------------------+
//| R^2_adj, adjusted coefficient of determination                   |
//| Inputs: y - real values of dependent variable, e - residuals     |
//|         k - number of regressors                                 |
//+------------------------------------------------------------------+
double R2_adj(vector& y, vector& e, ulong k)
  {
//--- Input data validation
   ulong n=y.Size();
   if(n<2)
     {Print("R2_adj() error: y too short"); return 0.0;}
   if(n<=k)
     {Print("R2_adj() error: n <= k"); return 0.0;}
   if(e.Size()!=n)
     {Print("R2_adj() error: different sizes of y and e"); return 0.0;}
   double Sy=y.Var(1);
   if(Sy<=DBL_MIN)
     {Print("R2_adj() error: y = const"); return 0.0;}
   return 1.0-e.Var((int)k)/Sy;
  }
//+------------------------------------------------------------------+
//| VIF, Variance Inflation Factor for X2, ..., Xk (assumed X1=const)|
//| Input: X - regressors                                            |
//| Output: vif - (k-1)-size VIF vector                              |
//+------------------------------------------------------------------+
void VIF(matrix& X, vector& vif)
  {
//--- Input data validation
   ulong k=X.Cols(), n=X.Rows();
   if(k<2)
     {Print("VIF() error: not enough regressors"); return;}
   if(n<k)
     {Print("VIF() error: not enough samples"); return;}
   vif.Resize(k-1);
   matrix X_1;
   vector b,e;
   double c;
   for(ulong i=1;i<k;++i)
     {
      X_1=X;
      if(i!=k-1)
         X_1.Col(X_1.Col(k-1),i);
      X_1.Resize(n,k-1);
      regression(X.Col(i),X_1,b,e,c);
      vif[i-1]=1.0/(1.0-R2(X.Col(i),e));
     }
  }
//+------------------------------------------------------------------+
//| Cook's distance                                                  |
//| Input: y - dependent variable, X - regressors                    |
//| Output: cd - n-size distance vector                              |
//+------------------------------------------------------------------+
void Cook_dist(vector& y, matrix& X, vector& D)
  {
//--- Input data validation
   ulong k=X.Cols(), n=X.Rows();
   if(k<2)
     {Print("Cook_dist() error: not enough regressors"); return;}
   if(n<=k)
     {Print("Cook_dist() error: not enough samples"); return;}
   if(X.Rank()<k)
     {Print("Cook_dist() error: low rank of X"); return;}
   if(y.Size()!=n)
     {Print("Cook_dist() error: wrong size of y"); return;}
   D.Resize(n);
   matrix A=X.Transpose().MatMul(X).Inv();
   vector b,e,x;
   double c,h,r;
   regression(y,X,b,e,c);
   for(ulong i=0;i<n;++i)
     {
      x=X.Row(i);
      h=x.MatMul(A.MatMul(x));
      r=e[i]/(c*MathSqrt(1-h));
      D[i]=r*r*h/(k*(1-h));
     }
  }
//+------------------------------------------------------------------+
//| Ljung-Box test for autocorrelation in residuals                  |
//| m - lags number, pq - sum of p and q in ARMA(p,q) model          |
//+------------------------------------------------------------------+
void Ljung_Box_test(vector& residuals, int m = 10, int pq = 0)
  {
   double e[];
   vector2array(residuals,e);
   Ljung_Box_test(e,m,pq);
  }
//+------------------------------------------------------------------+
//| Jarque-Bera test for normality of distribution                   |
//+------------------------------------------------------------------+
void Jarque_Bera_test(vector& residuals)
  {
   double e[];
   vector2array(residuals,e);
   Jarque_Bera_test(e);
  }
//+------------------------------------------------------------------+
//| Pettitt test for structural break (abrupt change in mean)        |
//+------------------------------------------------------------------+
void Pettitt_test(vector& residuals)
  {
   double e[];
   vector2array(residuals,e);
   Pettitt_test(e);
  }
//+------------------------------------------------------------------+
//| Breusch-Pagan test                                               |
//| Input: y - dependent variable, X - regressors                    |
//+------------------------------------------------------------------+
void Breusch_Pagan_test(vector& y, matrix& X)
  {
   Print("Breusch-Pagan test result:");
//--- Input data validation
   if(X.Rank()<X.Cols())
     {Print("error: low rank of X"); return;}
   if(X.Rows()!=y.Size())
     {Print("error: wrong size of y"); return;}
//--- Test statistic calculation
   vector b,e,e2;
   double c;
   regression(y,X,b,e,c);
   e2=e*e;
   regression(e2,X,b,e,c);
   double LM=X.Rows()*R2(e2,e);
//--- p-value computation
   int err=0;
   double p_value=MathCumulativeDistributionChiSquare(LM,X.Cols()-1,false,false,err);
   if(err!=0)
     {Print("error: MathCumulativeDistributionChiSquare() error"); return;}
//--- Print result
   PrintFormat("LM = %.3f, p-value = %.3f", LM, p_value);
  }
//+------------------------------------------------------------------+
//| Ramsey RESET test                                                |
//| Input: y - dependent variable, X - regressors                    |
//+------------------------------------------------------------------+
void RESET_test(vector& y, matrix& X)
  {
   Print("RESET test result:");
//--- Input data validation
   if(X.Rank()<X.Cols())
     {Print("error: low rank of X"); return;}
   if(X.Rows()!=y.Size())
     {Print("error: wrong size of y"); return;}
//--- Test statistic calculation
   matrix X2=X;
   vector b,e,yfit;
   double c;
   regression(y,X,b,e,c);
   yfit=X2.MatMul(b);
   add_column(X2,yfit*yfit);
   add_column(X2,yfit*yfit*yfit);
   Print("based on F-test:");
   ulong irs[2]= {X2.Cols()-2,X2.Cols()-1};
   F_test(y,X2,irs);
  }
//+------------------------------------------------------------------+
//| F-test for all regressors (exclude X1)                           |
//| Input: y - dependent variable, X - regressors                    |
//+------------------------------------------------------------------+
void F_test_all(vector& y, matrix& X)
  {
   Print("F-test for all regressors result:");
//--- Input data validation
   ulong k=X.Cols();
   if(k<2)
     {Print("error: regressors list too short"); return;}
   ulong irs[];
   ArrayResize(irs,(int)k-1);
   for(int i=0;i<(int)k-1;++i)
      irs[i]=i+1;
   Print("based on F-test:");
   F_test(y,X,irs);
  }
//+------------------------------------------------------------------+
//| F-test (assumed X1==const and it not tested)                     |
//| Input: y - dependent variable, X - regressors                    |
//|        irs[] - list of indexes of tested regressors              |
//+------------------------------------------------------------------+
void F_test(vector& y, matrix& X, ulong& irs[])
  {
   Print("F-test result:");
//--- Input data validation
   int n_irs=ArraySize(irs);
   if(n_irs<1)
     {Print("error: empty regressors list"); return;}
   ulong n=X.Rows(),k=X.Cols();
   if(X.Rank()<k)
     {Print("error: low rank of X"); return;}
   if(n!=y.Size())
     {Print("error: wrong size of y"); return;}
   if(!ArraySort(irs))
     {Print("error: ArraySort() error"); return;}
   if(irs[0]==0)
     {Print("error: constant should not be tested"); return;}
   if(irs[n_irs-1]>=k)
     {Print("error: out of regressors range"); return;}
   for(int i=0;i<n_irs-1;++i)
     {
      if(irs[i]==irs[i+1])
        {Print("error: repeating regressors"); return;}
     }
//--- Matrix for short regression calculation
   matrix Xshort=X;
   ulong j=k-1;
   for(int i=n_irs-1;i>=0;--i)
     {
      if(irs[i]<j)
         Xshort.SwapCols(irs[i],j);
      --j;
     }
   ulong kr=k-n_irs;
   Xshort.Resize(n,kr);
//--- Test statistic calculation
   vector b,e;
   double c,cr;
   regression(y,X,b,e,c);
   regression(y,Xshort,b,e,cr);
   double F=cr/c;
   F*=F;
   F=(F*(n-kr)-n+k)/n_irs;
//--- p-value computation
   int err=0;
   double p_value=MathCumulativeDistributionF(F,n_irs,n-k,false,false,err);
   if(err!=0)
     {Print("error: MathCumulativeDistributionF() error ",err); return;}
//--- Print result
   PrintFormat("F = %.3f, p-value = %.3f", F, p_value);
  }
//+------------------------------------------------------------------+
//| helper function for copying a vector to an array                 |
//+------------------------------------------------------------------+
void vector2array(vector& v, double& a[])
  {
   int na=(int)v.Size();
   ArrayResize(a,na);
   for(int i=0;i<na;++i)
      a[i]=v[(ulong)i];
  }
//+------------------------------------------------------------------+
//| helper functions for adding a column to a matrix (as last)       |
//+------------------------------------------------------------------+
void add_column(matrix& m, vector& col)
  {
   ulong k=m.Cols();
   if(k==0)
     {
      m.Assign(col);
      m=m.Transpose();
      return;
     }
   ulong n=m.Rows();
   if(n!=col.Size())
      return;
   m.Resize(n,k+1);
   m.Col(col,k);
  }
//+------------------------------------------------------------------+
void add_column(matrix& m, double& col[])
  {
   ulong k=m.Cols();
   if(k==0)
     {
      m.Assign(col);
      m=m.Transpose();
      return;
     }
   ulong n=m.Rows();
   if(n!=col.Size())
      return;
   m.Resize(n,k+1);
   for(ulong i=0;i<n;++i)
      m[i][k]=col[i];
  }
//+------------------------------------------------------------------+
//| vector scaling function                                          |
//+------------------------------------------------------------------+
void scale(vector& v)
  {
   if(v.Size()<2)
      return;
   double m,s;
   m=v.Mean();
   s=v.Std();
   v=(v-m)/s;
  }
//+------------------------------------------------------------------+
//| matrix columns scaling function                                  |
//+------------------------------------------------------------------+
void scale_cols(matrix& A, ulong from, ulong to)
  {
   if(from>to)
      return;
   if(to>=A.Cols())
      return;
   if(A.Rows()<2)
      return;
   double m, s;
   for(ulong i=from;i<=to; ++i)
     {m=A.Col(i).Mean(); s=A.Col(i).Std(); A.Col((A.Col(i)-m)/s,i);}
  }
//+------------------------------------------------------------------+
//| Helper functions to export vector, matrix and vector+matrix      |
//| to text file (for data analysis in other programs)               |
//+------------------------------------------------------------------+
string vector2string(vector& v)
  {
   ulong k=v.Size();
   if(k<1)
      return "";
   string res=DoubleToString(v[0]);
   for(ulong i=1;i<k;++i)
      res+=" "+DoubleToString(v[i]);
   return res;
  }
//+------------------------------------------------------------------+
void vector2csv(vector& y, string fldr, string fnm)
  {
//--- Initialization and input data validation
   ulong ny = y.Size();
   if(ny < 1)
     {Print("No data for vector2csv()"); return;}
   int ftxt = FileOpen(fldr + "\\" + fnm, FILE_WRITE | FILE_TXT | FILE_ANSI | FILE_COMMON);
   if(ftxt == INVALID_HANDLE)
     {Print("FileOpen() error"); return;}
//--- Writing vector to file
   FileWriteString(ftxt, DoubleToString(y[0]));
   for(ulong i = 1; i < ny; ++i)
      FileWriteString(ftxt, "\n" + DoubleToString(y[i]));
   FileClose(ftxt);
  }
//+------------------------------------------------------------------+
void matrix2csv(matrix& X, string fldr, string fnm)
  {
//--- Initialization and input data validation
   ulong n = X.Rows();
   if(n < 1)
     {Print("No data for matrix2csv()"); return;}
   int ftxt = FileOpen(fldr + "\\" + fnm, FILE_WRITE | FILE_TXT | FILE_ANSI | FILE_COMMON);
   if(ftxt == INVALID_HANDLE)
     {Print("FileOpen() error"); return;}
//--- Writing matrix to file
   FileWriteString(ftxt, vector2string(X.Row(0)));
   for(ulong i = 1; i < n; ++i)
      FileWriteString(ftxt, "\n"+vector2string(X.Row(i)));
   FileClose(ftxt);
  }
//+------------------------------------------------------------------+
void vector_matrix2csv(vector& y, matrix& X, string fldr, string fnm)
  {
//--- Initialization and input data validation
   ulong n = X.Rows();
   if(n < 1)
     {Print("No data for vector_matrix2csv()"); return;}
   if(y.Size()!=n)
     {Print("wrong y size"); return;}
   int ftxt = FileOpen(fldr + "\\" + fnm, FILE_WRITE | FILE_TXT | FILE_ANSI | FILE_COMMON);
   if(ftxt == INVALID_HANDLE)
     {Print("FileOpen() error"); return;}
//--- Writing vector and matrix to file
   FileWriteString(ftxt,  DoubleToString(y[0])+" "+vector2string(X.Row(0)));
   for(ulong i = 1; i < n; ++i)
      FileWriteString(ftxt, "\n"+DoubleToString(y[i])+" "+vector2string(X.Row(i)));
   FileClose(ftxt);
  }
//+------------------------------------------------------------------+