794 lines
27 KiB
MQL5
794 lines
27 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| utils.mqh |
|
|
//| Copyright 2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2025, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include<Math\Alglib\solvers.mqh>
|
|
#include<Math\Stat\T.mqh>
|
|
#include<Math\Stat\Beta.mqh>
|
|
//+------------------------------------------------------------------+
|
|
//| 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<typename T>
|
|
matrix<T> toeplitz(vector<T>& x)
|
|
{
|
|
ulong n = x.Size();
|
|
matrix<T> 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<typename T>
|
|
matrix<T> toeplitz(vector<T>& x, vector<T>& y)
|
|
{
|
|
ulong nx = x.Size();
|
|
ulong ny = y.Size();
|
|
|
|
// Create a unified vector z representing the complete sequence of diagonals
|
|
vector<T> z = vector<T>::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<T> 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;
|
|
}
|
|
//+------------------------------------------------------------------+
|