Article-23677-EGARCH-MQL5-Volatility-Modeling/Regression/utils.mqh

//+------------------------------------------------------------------+
//|                                                        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;
  }
//+------------------------------------------------------------------+