//+------------------------------------------------------------------+ //| utils.mqh | //| Copyright 2025, MetaQuotes Ltd. | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Copyright 2025, MetaQuotes Ltd." #property link "https://www.mql5.com" #include #include #include //+------------------------------------------------------------------+ //| Alternate hypothesis specification for statistical significance | //| Determines the directionality of the test (e.g., in unit root tests) | //+------------------------------------------------------------------+ enum ENUM_ALTERNATE_HYPOTHESIS { TWO_SIDED=0, // Reject if slope is non-zero (Two-tailed test) LESS, // Reject if slope is strictly less than zero GREATER // Reject if slope is strictly greater than zero }; //+------------------------------------------------------------------+ //| Deterministic components (trend/intercept) in regression model | //+------------------------------------------------------------------+ enum ENUM_TREND { TREND_NONE=0, // No deterministic terms TREND_CONST_ONLY, // Intercept only TREND_LINEAR_ONLY, // Linear trend only TREND_LINEAR_CONST, // Linear trend + Intercept TREND_QUAD_LINEAR_CONST // Quadratic + Linear + Intercept }; //+------------------------------------------------------------------+ //| Options for handling existing constant columns in design matrices| //+------------------------------------------------------------------+ enum ENUM_HAS_CONST { HAS_CONST_RAISE=0, // Error if constant exists HAS_CONST_SKIP, // Ignore existing constant HAS_CONST_ADD // Append as an additional term }; //+------------------------------------------------------------------+ //| Options for trimming invalid or NA observations from datasets | //+------------------------------------------------------------------+ enum ENUM_TRIM { TRIM_NONE=0, TRIM_FORWARD, // Trim leading invalid observations TRIM_BACKWARD, // Trim trailing invalid observations TRIM_BOTH // Trim from both ends }; //+------------------------------------------------------------------+ //| Options for handling/integrating the original input data set | //+------------------------------------------------------------------+ enum ENUM_ORIGINAL { ORIGINAL_EX=0, // Exclude original data from result ORIGINAL_IN, // Include original data in output ORIGINAL_SEP // Separate original data for distinct processing }; //+------------------------------------------------------------------+ //| Yule-Walker estimation method variants | //+------------------------------------------------------------------+ enum ENUM_YW_METHOD { YW_ADJUSTED=0, // Use unbiased Yule-Walker estimates YW_MLE // Use Maximum Likelihood estimation }; //+------------------------------------------------------------------+ //| Structure to hold results of Yule-Walker autoregressive analysis | //+------------------------------------------------------------------+ struct YWResult { double sigma; // Estimated innovation variance vector rho; // Autoregressive coefficients matrix inv; // Covariance or inverse information matrix // Initialize results to empty states YWResult(void) { sigma = EMPTY_VALUE; rho = vector::Zeros(0); inv = matrix::Zeros(0, 0); } ~YWResult(void) { } // Copy constructor for result persistence YWResult(YWResult& other) { sigma = other.sigma; rho = other.rho; inv = other.inv; } // Assignment operator for updating results void operator=(YWResult& other) { sigma = other.sigma; rho = other.rho; inv = other.inv; } }; //+------------------------------------------------------------------+ //| Calculates the p-value for a given t-statistic | //| Determines the probability of observing the test result | //| under the null hypothesis given degrees of freedom. | //+------------------------------------------------------------------+ double _get_pvalue(double statistic, double df, ENUM_ALTERNATE_HYPOTHESIS alternative, bool symmetric = true) { int error = 0; double pvalue = EMPTY_VALUE; switch(alternative) { // Left-tailed test: Probability mass in the lower tail case LESS: pvalue = MathCumulativeDistributionT(statistic, df, error); break; // Right-tailed test: Probability mass in the upper tail (1 - CDF) case GREATER: pvalue = 1.0 - MathCumulativeDistributionT(statistic, df, error); break; // Two-sided test: Considers extremes in both directions case TWO_SIDED: // If symmetric, double the tail probability of the absolute statistic // Otherwise, take the minimum tail and multiply for a non-symmetric comparison pvalue = symmetric ? 2. * (1.0 - MathCumulativeDistributionT(fabs(statistic), df, error)) : fmin(MathCumulativeDistributionT(statistic, df, error), 1.0 - MathCumulativeDistributionT(statistic, df, error)); break; } // Handle potential math library errors if(error) pvalue = EMPTY_VALUE; return pvalue; } //+------------------------------------------------------------------+ //| Sequence | //| Generates an arithmetic progression of double-precision numbers. | //| | //| Arguments: | //| from : Starting value of the sequence | //| to : Upper bound of the sequence | //| step : The increment between consecutive values | //| result[] : Output array to store the generated values | //| | //| Return value: true if successful, otherwise false. | //+------------------------------------------------------------------+ bool Sequence(const double from, const double to, const double step, double &result[]) { //--- Validate inputs: Ensure all parameters are valid numeric values if(!MathIsValidNumber(from) || !MathIsValidNumber(to) || !MathIsValidNumber(step)) return(false); //--- Validate logic: Start must be less than or equal to the end value if(to < from) return(false); //--- Calculate total elements: Number of steps that fit in the range [from, to] int count = 1 + int((to - from) / step); //--- Memory management: Ensure the output array is large enough to hold the sequence if(ArraySize(result) < count) if(ArrayResize(result, count) != count) return(false); //--- Generate sequence: Populate array using the arithmetic progression formula: a_i = a_0 + i * d for(int i = 0; i < count; i++) result[i] = from + i * step; return(true); } //+------------------------------------------------------------------+ //| Helper function: Augments the input matrix with a trend and/or constant column //+------------------------------------------------------------------+ bool addtrend(matrix &in, matrix &out, ENUM_TREND trend = TREND_CONST_ONLY, bool prepend = false, ENUM_HAS_CONST has_const = HAS_CONST_SKIP) { // Handle base case: No trend requested, simply copy input to output if(trend == TREND_NONE) return out.Copy(in); // Determine the degree of the polynomial trend based on the requested ENUM_TREND ulong trendorder = 0; if(trend == TREND_CONST_ONLY) trendorder = 0; else if(trend == TREND_LINEAR_CONST || trend == TREND_LINEAR_ONLY) trendorder = 1; else if(trend == TREND_QUAD_LINEAR_CONST) trendorder = 2; ulong nobs = in.Rows(); vector tvector, temp; // Prepare time vector for trend calculation if(!tvector.Resize(nobs) || !temp.Resize(nobs)) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", ::GetLastError()); return false; } double sequence[]; if(!Sequence(1, nobs, 1, sequence) || !tvector.Assign(sequence)) { Print(__FUNCTION__, " ", __LINE__, " Error ", ::GetLastError()); return false; } // Construct the trend matrix: columns represent t^0, t^1, t^2... matrix trendmat; if(!trendmat.Resize(tvector.Size(), (trend == TREND_LINEAR_ONLY) ? 1 : trendorder + 1)) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", ::GetLastError()); return false; } for(ulong i = 0; i < trendmat.Cols(); i++) { // If only linear trend, skip the intercept (t^0) temp = MathPow(tvector, (trend == TREND_LINEAR_ONLY) ? double(i + 1) : double(i)); if(!trendmat.Col(temp, i)) { Print(__FUNCTION__, " ", __LINE__, " Matrix Assign Error ", ::GetLastError()); return false; } } // Detect if input matrix already contains a constant to prevent redundant column creation vector ptp = in.Ptp(0); // Peak-to-peak (max-min) if(!ptp.Min()) // If min range is 0, column is constant { if(has_const == HAS_CONST_RAISE) { Print("Input matrix contains one or more constant columns"); return false; } if(has_const == HAS_CONST_SKIP) { matrix vsplitted[]; ulong parts[] = {1, trendmat.Cols() - 1}; trendmat.Vsplit(parts, vsplitted); // Remove constant column from trendmat if(!trendmat.Copy(vsplitted[1])) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", ::GetLastError()); return false; } } } // Merge original data and trend matrix into the output matrix if(!out.Resize(trendmat.Rows(), trendmat.Cols() + in.Cols())) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", ::GetLastError()); return false; } ulong j = 0; matrix mtemp; // Arrange columns based on whether the trend should be prepended or appended if(prepend) mtemp.Copy(trendmat); else mtemp.Copy(in); for(j = 0; j < mtemp.Cols(); j++) { vector col = mtemp.Col(j); if(!out.Col(col, j)) { Print(__FUNCTION__, " ", __LINE__, " Matrix Assign Error ", ::GetLastError()); return false; } } ulong minus = j; if(prepend) mtemp.Copy(in); else mtemp.Copy(trendmat); for(; j < out.Cols(); j++) { vector col = mtemp.Col(j - minus); if(!out.Col(col, j)) { Print(__FUNCTION__, " ", __LINE__, " Matrix Assign Error ", ::GetLastError()); return false; } } return true; } //+------------------------------------------------------------------+ //| Helper function: Transforms RHS matrix into a lagged representation | //| Used for VAR/ARDL models to create a matrix of lagged predictors | //+------------------------------------------------------------------+ bool lagmat(matrix &in, matrix &out[], ulong mlag, ENUM_TRIM trim = TRIM_BOTH, ENUM_ORIGINAL original = ORIGINAL_IN) { ulong nobs = in.Rows(); ulong nvars = in.Cols(); ulong dropidx = 0; // Validation: Lag depth must be less than the available observations if(mlag >= nobs) { Print("mlag should be < number of rows of the input matrix"); return false; } // Determine if we should drop or separate the original (non-lagged) variables if(original == ORIGINAL_EX || original == ORIGINAL_SEP) dropidx = nvars; // Create a padded work matrix to hold the shifted lagged versions matrix nn; if(!nn.Resize(nobs + mlag, nvars * (mlag + 1))) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", ::GetLastError()); return false; } nn.Fill(0.0); ulong maxlag = mlag; ulong row_end, row_start, col_end, col_start, j, z; row_end = row_start = col_end = col_start = j = z = 0; // Populate the matrix with shifted values (lag structure) for(ulong k = 0; k < (maxlag + 1); k++) { row_start = maxlag - k; row_end = nobs + maxlag - k; col_start = nvars * (maxlag - k); col_end = nvars * (maxlag - k + 1); j = 0; for(ulong irow = row_start; irow < row_end; irow++, j++) { z = 0; for(ulong icol = col_start; icol < col_end; icol++, z++) nn[irow][icol] = in[j][z]; } } // Define observation window based on requested trimming ulong startobs, stopobs; if(trim == TRIM_NONE || trim == TRIM_FORWARD) startobs = 0; else startobs = maxlag; if(trim == TRIM_NONE || trim == TRIM_BACKWARD) stopobs = nn.Rows(); else stopobs = nobs; // Prepare output array (1 for combined, 2 if separating original data) if(dropidx && original == ORIGINAL_SEP) ArrayResize(out, 2); else ArrayResize(out, 1); // Extract lagged variables into the first output index if(!out[0].Resize(stopobs - startobs, nn.Cols() - dropidx)) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", ::GetLastError()); return false; } for(ulong irow = startobs; irow < stopobs; irow++) for(ulong icol = dropidx; icol < nn.Cols(); icol++) out[0][irow - startobs][icol - dropidx] = nn[irow][icol]; // Extract original variables into the second output index if requested if(out.Size() > 1) { if(!out[1].Resize(stopobs - startobs, dropidx)) { Print(__FUNCTION__, " ", __LINE__, " Resize Error ", GetLastError()); return false; } for(ulong irow = startobs; irow < stopobs; irow++) for(ulong icol = 0; icol < out[1].Cols(); icol++) out[1][irow - startobs][icol] = nn[irow][icol]; } return true; } //+------------------------------------------------------------------+ //| Toeplitz matrix construction (symmetric) | //| Generates a square matrix where each descending diagonal from | //| left to right is constant, based on a single input vector. | //+------------------------------------------------------------------+ template matrix toeplitz(vector& x) { ulong n = x.Size(); matrix out(n, n); // Fill matrix: out[i,j] depends on the distance between indices for(ulong i = 0; i < n; ++i) for(ulong j = 0; j < n; ++j) // Symmetric property: diagonals are defined by the absolute difference if(i <= j) out[i, j] = x[j - i]; else out[i, j] = x[i - j]; return out; } //+------------------------------------------------------------------+ //| Toeplitz matrix construction (asymmetric) | //| Generates an (nx by ny) matrix using two vectors to define | //| the first column and the first row respectively. | //+------------------------------------------------------------------+ template matrix toeplitz(vector& x, vector& y) { ulong nx = x.Size(); ulong ny = y.Size(); // Create a unified vector z representing the complete sequence of diagonals vector z = vector::Zeros(nx + ny - 1); ulong start = 1; // Construct the sequence z: [reversed x, followed by y elements excluding the first] for(ulong i = 0; i < z.Size(); ++i) { if(i < nx) z[i] = x[(nx - i) - 1]; else z[i] = y[start++]; } // Map the sequence z into the matrix based on index offsets matrix out(nx, ny); for(ulong i = 0; i < nx; ++i) { for(ulong j = 0; j < ny; ++j) { if(i < j) // Upper triangle entries out[i, j] = z[nx + (j - 1) - i]; else // Lower triangle entries out[i, j] = z[nx - (i + 1) + j]; } } return out; } //+---------------------------------------------------------------------------+ //| Estimate AR(p) parameters from a sequence using the Yule-Walker equations.| //| Solves the system R * rho = r, where R is the Toeplitz autocorrelation matrix | //+---------------------------------------------------------------------------+ YWResult yule_walker(vector& x, long order = 1, ENUM_YW_METHOD method = YW_ADJUSTED, long df = -1, bool inv = false, bool demean = true) { YWResult out; vector in = x; if(!in.Size()) { Print(__FUNCTION__, " input sequence is empty "); return out; } // Optional: Center data around the mean (standard for AR estimation) if(demean) in = in - in.Mean(); long n = df > 0 ? df : long(in.Size()); vector r = vector::Zeros(order + 1); // Estimate the variance (lag 0 autocovariance) r[0] = pow(in, 2.0).Sum() / double(n); double adj_needed = (double)(method == YW_ADJUSTED); // Estimate autocovariances for lags 1 to order for(long k = 1; k < (order + 1); ++k) { double sum = 0; for(long i = 0; i < (long(in.Size()) - k); ++i) sum += (in[i] * in[k + i]); // Apply adjustment factor if YW_ADJUSTED is selected r[k] = sum / double(n - k * adj_needed); } // Prepare the Toeplitz matrix R (autocorrelation matrix) vector rslice = vector::Zeros(r.Size() - 1); for(ulong i = 0; i < rslice.Size(); rslice[i] = r[i], ++i); matrix R = toeplitz(rslice); // Prepare the target vector (r_1 to r_p) for(ulong i = 0; i < rslice.Size(); rslice[i] = r[i + 1], ++i); double b[], sol[]; ArrayResize(b, (int)rslice.Size()); for(uint i = 0; i < b.Size(); b[i] = rslice[i], ++i); // Solve the linear system R * rho = r_slice CMatrixDouble Rr = R; CDenseSolverLSReport ls; int opinfo; CDenseSolver::RMatrixSolveLS(Rr, Rr.Rows(), Rr.Cols(), b, 0, opinfo, ls, sol); // Fallback to pseudo-inverse if the matrix is singular or ill-conditioned if(opinfo < 0) { Print(__FUNCTION__, " possible singular matrix. Using pinv "); out.rho = (R.PInv()).MatMul(rslice); } else out.rho.Assign(sol); // Estimate innovation variance (sigma^2 = r_0 - rho' * r_1:p) double sigmasq = r[0] - (rslice * out.rho).Sum(); if(MathClassify(sigmasq) == FP_NORMAL && sigmasq > 0) out.sigma = sqrt(sigmasq); else out.sigma = double("nan"); // Optionally return the inverse covariance matrix (information matrix) if(inv) out.inv = R.Inv(); return out; } //+------------------------------------------------------------------+ //| LinRegressResult: Container for linear regression statistics | //| Holds coefficients, correlation, significance, and error metrics | //+------------------------------------------------------------------+ struct LinRegressResult { double slope; // The estimated coefficient of the independent variable double intercept; // The estimated Y-axis intercept double rvalue; // Pearson correlation coefficient (R-value) double pvalue; // Statistical significance of the slope double stderr; // Standard error of the slope estimate double intercept_stderror; // Standard error of the intercept estimate // Constructor: Initializes all statistics to zero LinRegressResult(void) { slope = intercept = rvalue = pvalue = stderr = intercept_stderror = 0.0; } // Constructor: Parametric initialization for immediate result storage LinRegressResult(double slope_, double intercept_, double rvalue_, double pvalue_, double stderr_, double inter_stderr_) { slope = slope_; intercept = intercept_; rvalue = rvalue_; pvalue = pvalue_; stderr = stderr_; intercept_stderror = inter_stderr_; } // Copy constructor: Enables object duplication LinRegressResult(LinRegressResult &other) { slope = other.slope; intercept = other.intercept; rvalue = other.rvalue; pvalue = other.pvalue; stderr = other.stderr; intercept_stderror = other.intercept_stderror; } // Assignment operator: Allows deep copying of results between objects void operator=(LinRegressResult &other) { slope = other.slope; intercept = other.intercept; rvalue = other.rvalue; pvalue = other.pvalue; stderr = other.stderr; intercept_stderror = other.intercept_stderror; } }; //+------------------------------------------------------------------------+ //| Calculate a linear least-squares regression for two sets of measurements | //| Fits y = intercept + slope * x by minimizing the sum of squared residuals | //+------------------------------------------------------------------------+ LinRegressResult linregress(vector& x, vector& y, ENUM_ALTERNATE_HYPOTHESIS alternative = TWO_SIDED) { double tiny = 1.e-20; // Small epsilon to prevent division by zero // Validate inputs: Ensure vectors are non-empty if(!x.Size() || !y.Size()) { Print(__FUNCTION__, ": One or more of the inputs is empty"); return LinRegressResult(); } // Validate inputs: Regression requires variation in the independent variable if(x.Max() == x.Min() && x.Size() > 1) { Print(__FUNCTION__, ": Cannot calculate regression if all x values are identical"); return LinRegressResult(); } // Compute means for centralized calculations double xmean = x.Mean(); double ymean = y.Mean(); // Sum of squares: variance/covariance components double ssxm = (pow(x - xmean, 2.)).Mean(); // Variance of X double ssxym = ((x - xmean) * (y - ymean)).Mean(); // Covariance of X and Y double ssym = (pow(y - ymean, 2.)).Mean(); // Variance of Y // Calculate Pearson correlation coefficient (r) double r; if(!ssxm || !ssym) r = 0.0; else { r = ssxym / sqrt(ssxm * ssym); // Clamp r within the valid range [-1, 1] due to potential floating point errors if(r > 1.) r = 1.; else if(r < -1.) r = -1.; } // Calculate regression coefficients: slope = Cov(X,Y)/Var(X), intercept = Mean(Y) - slope * Mean(X) double slope = ssxym / ssxm; double intercept = ymean - slope * xmean; double prob, slope_stderr, intercept_stderr; int df; // Case: Two points determine a perfect line if(x.Size() == 2) { prob = (y[0] == y[1]) ? 1. : 0.0; slope_stderr = intercept_stderr = 0.0; } else { // Statistical inference: degrees of freedom and t-statistic for significance testing df = int(x.Size() - 2); double t = r * sqrt(df / ((1.0 - r + tiny) * (1. + r + tiny))); prob = _get_pvalue(t, double(df), alternative); // Calculate standard errors for slope and intercept coefficients slope_stderr = sqrt((1. - pow(r, 2.)) * ssym / ssxm / double(df)); intercept_stderr = slope_stderr * sqrt(ssxm + pow(xmean, 2.)); } return LinRegressResult(slope, intercept, r, prob, slope_stderr, intercept_stderr); } //+------------------------------------------------------------------+ //|Pearson correlation result | //+------------------------------------------------------------------+ struct PearsonResult { double correlation; double pvalue; ulong n; vector x,y; ENUM_ALTERNATE_HYPOTHESIS alternate_hypothesis; PearsonResult(void) { correlation = pvalue = EMPTY_VALUE; n = 0; x = y = vector::Zeros(0); alternate_hypothesis = TWO_SIDED; } PearsonResult(vector& x_, vector& y_, double correlation_,double pvalue_,ENUM_ALTERNATE_HYPOTHESIS alt_hypo) { correlation = correlation_; x = x_; y = y_; pvalue = pvalue_; alternate_hypothesis = alt_hypo; n = x.Size(); } PearsonResult(PearsonResult& other) { correlation = other.correlation; pvalue = other.pvalue; n = other.n; x = other.x; y = other.y; alternate_hypothesis = other.alternate_hypothesis; } void operator=(PearsonResult& other) { correlation = other.correlation; pvalue = other.pvalue; n = other.n; x = other.x; y = other.y; alternate_hypothesis = other.alternate_hypothesis; } }; //+------------------------------------------------------------------+ //|Pearson correlation coefficient and p-value for testing non-correlation. | //+------------------------------------------------------------------+ PearsonResult pearsonr(vector& x, vector& y, ENUM_ALTERNATE_HYPOTHESIS alt_hypothesis=TWO_SIDED) { if(x.Size()!=y.Size() || x.Size()<2) { Print(__FUNCTION__," : X and Y must have the same length and must have length at least 2."); return PearsonResult(); } double threshold = pow(DBL_EPSILON,0.75); bool const_xy = (x.Max()-x.Min() == 0.0 || y.Max()-y.Min() == 0.0); if(const_xy) { Print(__FUNCTION__," : An input array is constant; the correlation coefficient is not defined."); return PearsonResult(); } vector xm = x - x.Mean(); vector ym = y - y.Mean(); double xmax = (fabs(xm)).Max(); double ymax = (fabs(ym)).Max(); vector nxm = xm/xmax; vector nym = ym/ymax; double normxm = xmax*nxm.Norm(VECTOR_NORM_P); double normym = ymax*nym.Norm(VECTOR_NORM_P); bool nconst_x = (normxm < threshold*fabs(x.Mean())); bool nconst_y = (normym < threshold*fabs(y.Mean())); bool nconst_xy = (nconst_x|nconst_y); if(nconst_xy && (~const_xy)) Print(__FUNCTION__," : WARNING : An input array is nearly constant; the computed correlation coefficient may be inaccurate."); vector xy = xm/normxm * ym/normym; double r = xy.Sum(); PearsonResult result; result.correlation = r; if(r<-1.0) r = -1; else if(r>1.) r = 1.; if(x.Size() == 2) { r = MathRound(r,8); result.pvalue = double("nan"); } else { int errcode = 0; double ab = (double(x.Size())/2.0) - 1.; switch(alt_hypothesis) { case LESS: result.pvalue = MathCumulativeDistributionBeta((r+1.)/2.0,ab,ab,errcode); break; case GREATER: result.pvalue = 1. - MathCumulativeDistributionBeta((r+1.)/2.0,ab,ab,errcode); break; case TWO_SIDED: result.pvalue = 2.* (1. - MathCumulativeDistributionBeta((fabs(r) + 1.)/2.0,ab,ab,errcode)); break; } if(errcode) { Print(__FUNCTION__," : Error evaluating CDF of Beta function. ", GetLastError()); result.pvalue = double("nan"); } } result.n = x.Size(); result.x = x; result.y = y; result.alternate_hypothesis = alt_hypothesis; return result; } //+------------------------------------------------------------------+