1345 lines
53 KiB
MQL5
1345 lines
53 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| 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);
|
|
}
|
|
//+------------------------------------------------------------------+
|