alexeynikolaev2/Article-22235-Multiple-Regression-And-Diagnostics
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Article-22235-Multiple-Regr.../EconometricsM.mqh

1025 lines
40 KiB
MQL5
Raw Permalink Normal View History

new files added
2026-07-09 07:53:14 +03:00
//+------------------------------------------------------------------+
//| 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);
}
//+------------------------------------------------------------------+
Reference in a new issue Copy permalink