alexeynikolaev2/Article-23435-OLS-Generalization
Watch
1
0
Fork
You've already forked Article-23435-OLS-Generalization
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Article-23435-OLS-Generaliz.../EconometricsM.mqh

1345 lines
53 KiB
MQL5
Raw Permalink Normal View History

new files added
2026-08-12 10:37:03 +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);
}
//+-------------------------------------------------------------------+
//| OLS regression with Newey-West HAC standard errors estimation |
//| |
//| ARGUMENTS: |
//| [in] y - Dependent variable vector (size n x 1) |
//| [in] X - Regressors matrix (size n x k, first col is 1s) |
//| [out] stats - Dynamic array to store coefficients statistics |
//| [in/out] prog - Structure containing new data and forecast results|
//| [in] useNWE - Switch: true = Newey-West, false = Classical OLS |
//+-------------------------------------------------------------------+
void regression_Newey_West(const vector& y, const matrix &X, CoefficientStats& stats[], SPrognose& prog, bool useNWE=true)
{
//--- Input data validation
ulong n=y.Size(),k=X.Cols();
if(n<1)
{Print("regression_Newey_West() error: empty y"); return;}
if(n!=X.Rows())
{ Print("regression_Newey_West() error: X rows matrix size must be equal to y size"); return; }
if(n<=k)
{Print("regression_Newey_West() error: not enough samples"); return;}
if(X.Rank()<k)
{Print("regression_Newey_West() error: multicollinearity in the X matrix"); return;}
if(prog.xnew.Size()!=k)
{Print("regression_Newey_West() error: wrong vector xnew size for prognose"); return;}
//--- Coefficients and residuals computation
matrix XT_X_1=X.Transpose().MatMul(X).Inv();
vector b=XT_X_1.MatMul(X.Transpose()).MatMul(y);
vector e=y-X.MatMul(b);
// The dot product of the residual vector with itself yields the Residual Sum of Squares (RSS)
double se2=e.Dot(e)/(n-k);
//--- Vb - coefficients covariation matrix
matrix Vb(k, k);
//--- Newey-West
if(useNWE)
{
matrix Q=matrix::Zeros(k, k), dq;
// 1. First sum: main diagonal of matrix Q (heteroskedasticity adjustment)
for(ulong t=1;t<=n;++t)
Q+=e[t-1]*e[t-1]*X.Row(t-1).Outer(X.Row(t-1));
// 2. Second sum: off-diagonal elements (autocorrelation adjustment)
ulong L=(ulong)MathFloor(4.0*MathPow(n/100.0,2.0/9.0));
double wj,tau;
for(ulong j=1;j<=L;++j)
{
tau=double(j)/(L+1);
// Parzen kernel weights calculation
wj=(tau<=0.5)?1.0-6.0*tau*tau*(1.0-tau):2.0*MathPow(1.0-tau,3.0);
dq=matrix::Zeros(k, k);
for(ulong t=j+1;t<=n;++t)
dq+=e[t-1]*e[t-j-1]*(X.Row(t-1).Outer(X.Row(t-j-1))+X.Row(t-j-1).Outer(X.Row(t-1)));
Q+=wj*dq;
}
Vb=XT_X_1.MatMul(Q).MatMul(XT_X_1);
}
//--- Classical OLS
else
Vb=se2*XT_X_1;
//--- Extract standard errors and compute coefficients statistics
ArrayResize(stats,(int)k);
int ner=0;
double t_crit;
for(ulong i = 0; i < k; ++i)
{
// Standard error is the square root of the diagonal element of Vb
stats[i].estimate = b[i];
stats[i].std_error = MathSqrt(Vb[i][i]);
// T-statistic computation (testing H0: b[i] == 0)
if(stats[i].std_error > 0)
stats[i].t_stat = stats[i].estimate / stats[i].std_error;
else
stats[i].t_stat = 0.0;
// One-sided p-value computation using Student's t-distribution
stats[i].p_value = MathCumulativeDistributionT(MathAbs(stats[i].t_stat),n-k,false,false,ner);
if(ner!=0)
{Print("regression_Newey_West() error: MathCumulativeDistributionT() error: ",ner); return;}
// 95% Confidence interval for the coefficient (critical t-value approx 1.96 for large n)
// For precision, we fetch the exact critical value from Student's inverse CDF
t_crit=MathQuantileT(0.025,n-k,false,false,ner);
if(ner!=0)
{Print("regression_Newey_West() error: MathQuantileT() error: ",ner); return;}
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;
}
//--- Forecast and interval boundaries computation
// 1. Point forecast (mean prediction)
// Vector dot product: prog.xnew (k x 1) and b (k x 1)
prog.point_progn = prog.xnew.Dot(b);
// 2. Variance of the forecast mean (for Confidence Interval)
// Quadratic form: x_new^T * Vb * x_new
double var_ci = prog.xnew.Dot(Vb.MatMul(prog.xnew));
double std_ci = MathSqrt(var_ci);
// 3. Variance of the individual prediction (for Prediction Interval)
double var_pi = var_ci + se2;
double std_pi = MathSqrt(var_pi);
// Get critical t-value for 95% threshold (pre-calculated t_crit from coefficients statistics computing can be reused)
t_crit=MathQuantileT(0.025,n-k,false,false,ner);
if(ner!=0)
{Print("regression_Newey_West() error: MathQuantileT() error: ",ner); return;}
// 4. Compute final boundaries for both intervals
prog.conf_low = prog.point_progn - t_crit * std_ci;
prog.conf_high = prog.point_progn + t_crit * std_ci;
prog.progn_low = prog.point_progn - t_crit * std_pi;
prog.progn_high = prog.point_progn + t_crit * std_pi;
}
//+-------------------------------------------------------------------+
//| Conditional Least Squares (CLS) for ARMA(p,q) model |
//| |
//| ARGUMENTS: |
//| [in] y - Dependent variable vector (size n x 1) |
//| [in] p - Order of the AR component |
//| [in] q - Order of the MA component |
//| [in] y_ - Pre-sample history of dependent variable (p x 1) |
//| [in] e_ - Pre-sample history of residuals/errors (q x 1) |
//| [out] e - residuals/errors (n x 1) |
//| [out] stats - Dynamic array to store coefficients statistics |
//| [out] prog - Structure containing next step forecast results |
//+-------------------------------------------------------------------+
void regression_CLS(const vector& y, const ulong p, const ulong q,
const vector& y_, const vector& e_, vector& e,
CoefficientStats& stats[], SPrognose& prog,
const ulong max_iter=1000, const double min_change=1e-5)
{
//--- Input data validation
ulong n=y.Size(),k=p+q+1;
if(max_iter==0)
{Print("regression_CLS() error: max_iter==0"); return;}
if(min_change<=0.0)
{Print("regression_CLS() error: min_change<=0.0"); return;}
if(n<=k)
{Print("regression_CLS() error: y too short"); return;}
if(y_.Size()!=p)
{Print("regression_CLS() error: wrong y_ size"); return;}
if(e_.Size()!=q)
{Print("regression_CLS() error: wrong e_ size"); return;}
//-- AR-part coefficients initialization with OLS
matrix X0=matrix::Ones(n,p+1);
for(ulong i=0;i<n;++i)
for(ulong j=1;j<=p;++j)
X0[i][j]=(i>=j)?y[i-j]:y_[p+i-j];
if(X0.Rank()<p+1)
{Print("regression_CLS() error: multicollinearity in the X0 matrix"); return;}
vector b=X0.Transpose().MatMul(X0).Inv().MatMul(X0.Transpose()).MatMul(y);
//--- Residuals/errors initialization based on the starting AR regression
e=y-X0.MatMul(b);
X0.Resize(0,0);
//--- Adding MA-part - zeros
b=b.Concat(vector::Zeros(q));
//--- Jacobian matrix
matrix J(n,k),JT_J(k,k);
vector delta_b=vector::Zeros(k);
double change=1.0;
//--- Levenberg-Marquardt parameters
double lambda = 0.01; // Initial damping factor
double v = 2.0; // Lambda multiplication factor
double rss_old = e.Dot(e); // Initial Residual Sum of Squares (RSS)
double rss_new; // New Residual Sum of Squares (RSS)
vector b_old = b; // To store backup of coefficients
vector e_old = e; // To store backup of residuals
//--- Main iteration loop
ulong iter=0;
for(iter=0;iter<max_iter;++iter)
{
for(ulong t=0;t<n;++t)
{
J[t][0]=-1.0;
for(ulong i=1;i<=q;++i)
J[t][0]-=(t>=i)?b[p+i]*J[t-i][0]:0.0;
for(ulong j=1;j<=p;++j)
{
J[t][j]=(t>=j)?-y[t-j]:-y_[p+t-j];
for(ulong i=1;i<=q;++i)
J[t][j]-=(t>=i)?b[p+i]*J[t-i][j]:0.0;
}
for(ulong j=1;j<=q;++j)
{
J[t][p+j]=(t>=j)?-e[t-j]:-e_[q+t-j];
for(ulong i=1;i<=q;++i)
J[t][p+j]-=(t>=i)?b[p+i]*J[t-i][p+j]:0.0;
}
}
//--- delta_b = -(J^T * J)^(-1) * J^T * e
JT_J = J.Transpose().MatMul(J);
// Add Lambda to the main diagonal (Levenberg-Marquardt modification)
for(ulong i = 0; i < k; ++i)
JT_J[i][i] += lambda;
// Compute the damped step
delta_b = -1.0 * JT_J.Inv().MatMul(J.Transpose()).MatMul(e);
//--- 1. Apply the step temporarily
b += delta_b;
//--- 2. Recalculate residuals vector e for the current parameters
for(ulong t = 0; t < n; ++t)
{
e[t] = y[t] - b[0];
for(ulong j = 1; j <= p; ++j)
e[t] -= (t >= j) ? b[j] * y[t - j] : b[j] * y_[p + t - j];
for(ulong j = 1; j <= q; ++j)
e[t] -= (t >= j) ? b[p + j] * e[t - j] : b[p + j] * e_[q + t - j];
}
rss_new = e.Dot(e);
//--- 3. Check if the step improved the model
if(rss_new < rss_old)
{
// Successful iteration: accept updates
lambda /= v; // Decrease damping
rss_old = rss_new;
b_old = b; // Update backups
e_old = e;
// Convergence check based on parameter changes
change = delta_b.Norm(VECTOR_NORM_P, 2.0);
if(change < min_change)
{
PrintFormat("L-M converged successfully after %d iterations.", iter + 1);
break;
}
}
else
{
// Unsuccessful iteration: rollback updates and increase damping
lambda *= v; // Increase damping to behave more like gradient descent
b = b_old; // Rollback coefficients
e = e_old; // Rollback residuals to previous valid state
// Guard against infinite loops with huge lambda
if(lambda > 1e10)
{
Print("L-M stopped: lambda became too large (optimization stuck).");
return;
}
}
}
if(iter>=max_iter)
{Print("L-M not converged"); return;}
//--- Compute unbiased residual variance (sigma^2)
double se2 = e.Dot(e) / double(n - k);
//--- Coefficients covariance matrix for non-linear OLS
//--- Calculated as sigma^2 * (J^T * J)^(-1) using the final iteration's Jacobian matrix
matrix Vb = se2 * J.Transpose().MatMul(J).Inv();
//--- Extract standard errors and compute coefficients statistics
ArrayResize(stats, (int)k);
int ner = 0;
// 95% Confidence interval for the coefficient
double t_crit = MathQuantileT(0.025, double(n - k), false, false, ner);
if(ner != 0)
{ Print("regression_CLS() error: MathQuantileT() failed with code ", ner); return; }
for(ulong i = 0; i < k; ++i)
{
stats[i].estimate = b[i];
// Standard error is the square root of the diagonal element of Vb
stats[i].std_error = MathSqrt(Vb[i][i]);
// T-statistic computation (testing H0: b[i] == 0)
if(stats[i].std_error > 0)
stats[i].t_stat = stats[i].estimate / stats[i].std_error;
else
stats[i].t_stat = 0.0;
// One-sided p-value computation using Student's t-distribution
stats[i].p_value = MathCumulativeDistributionT(MathAbs(stats[i].t_stat), double(n - k), false, false, ner);
if(ner != 0)
{ Print("regression_CLS() error: MathCumulativeDistributionT() failed with code ", ner); return; }
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;
}
//--- Forecast and prediction intervals computation (One-Step-Ahead)
// Create a full regressor vector for the forecast bar (size: p + q + 1)
prog.xnew.Resize(k);
// 1. Set Constant components
prog.xnew[0] = 1.0;
// 2. Fill AR lags from the very end of the sample vector y
for(ulong j = 1; j <= p; ++j)
prog.xnew[j] = (n >= j) ? y[n - j] : y_[p + (n - j)];
// 3. Fill MA lags from the very end of the calculated residuals vector e
for(ulong j = 1; j <= q; ++j)
prog.xnew[p + j] = (n >= j) ? e[n - j] : e_[q + (n - j)];
// 4. Compute Point Forecast (Mean prediction)
prog.point_progn = prog.xnew.Dot(b);
// 5. Variance of the forecast mean (Confidence Interval)
// For non-linear OLS, the gradient vector at the forecast point is equal to prog.xnew
double var_ci = prog.xnew.Dot(Vb.MatMul(prog.xnew));
double std_ci = MathSqrt(var_ci);
// 6. Variance of the individual future value (Prediction Interval)
double var_pi = var_ci + se2;
double std_pi = MathSqrt(var_pi);
// Pre-calculated negative t_crit from step 3 is reused here (with adjusted signs)
prog.conf_low = prog.point_progn - t_crit * std_ci;
prog.conf_high = prog.point_progn + t_crit * std_ci;
prog.progn_low = prog.point_progn - t_crit * std_pi;
prog.progn_high = prog.point_progn + t_crit * std_pi;
}
//+------------------------------------------------------------------+
//| 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