//+------------------------------------------------------------------+
//|                                                EconometricsA.mqh |
//|                             Copyright 2000-2026, MetaQuotes Ltd. |
//|                                                     www.mql5.com |
//+------------------------------------------------------------------+
#include <Math\Stat\T.mqh>
#include <Math\Stat\ChiSquare.mqh>
#include <Math\Stat\F.mqh>
#include <Math\Stat\Math.mqh>
#include <Graphics\Graphic.mqh>
//---
//+------------------------------------------------------------------+
//| Computation of parameters and residuals                          |
//| for regression on a constant                                     |
//+------------------------------------------------------------------+
void regression0(double& y[], double& a, double& c, double& e[])
  {
//--- Initialization and input data validation
   a = c = 0.0;
   ArrayResize(e, 0);
   int ny = ArraySize(y);
   if(ny < 1)
     {
      Print("no data for regression0()");
      return;
     }
   ArrayResize(e, ny);
//--- Parameters computation
   a = MathMean(y);
   c = MathStandardDeviation(y);
//--- Residuals computation
   for(int i = 0; i < ny; ++i)
      e[i] = y[i] - a;
  }
//+------------------------------------------------------------------+
//| Calculation of confidence interval for parameter A               |
//| in regression on a constant                                      |
//+------------------------------------------------------------------+
void Aconf0(double a, double c, int n, double conf_level)
  {
//--- Initialization and input data validation
   if(c <= 0.0)
     {
      Print("c must be positive in Aconf0()");
      return;
     }
   if(n < 2)
     {
      Print("n must be greater 1 in Aconf0()");
      return;
     }
   if(conf_level <= 0.0 || conf_level >= 1.0)
     {
      Print("conf_level must be between 0 and 1 in Aconf0()");
      return;
     }
   int err = 0;
//--- t-statistic calculation
//--- double t = MathQuantileT((1.0 - conf_level) / 2.0, n - 1, err);  // faulty library function
   double t = MathQuantileT_TMP((1.0 - conf_level) / 2.0, n - 1, err); // custom implementation
   if(err != 0)
     {
      Print("in MathQuantileT() error ", err);
      return;
     }
   t *= c / MathSqrt(n);
//--- Print result
   PrintFormat("A between %.4f and %.4f with confidence level %.3f", a + t, a - t, conf_level);
  }
//+------------------------------------------------------------------+
//| Test of hypothesis H₀: A = 0 vs H₁: A > 0                        |
//| for regression on a constant                                     |
//+------------------------------------------------------------------+
void Aeq0_test0(double a, double c, int n)
  {
   Print("test H0: A = 0 vs H1: A > 0 (for regession on constant) result:");
//--- Input data validation
   if(c <= 0.0)
     {
      Print("c must be positive in Aeq0_test0()");
      return;
     }
   if(n < 2)
     {
      Print("n must be greater 1 in Aeq0_test0()");
      return;
     }
//--- t-statistic calculation
   double t = a * MathSqrt(n) / c;
   int err = 0;
//--- p-value calculation
   double p_value = MathCumulativeDistributionT(t, n - 1, false, false, err);
   if(err != 0)
     {
      Print("in MathCumulativeDistributionT() error ", err);
      return;
     }
//--- Print result if success
   PrintFormat("t = %.3f, p-value = %.3f", t, p_value);
  }
//+------------------------------------------------------------------+
//| Calculation of prediction interval for Y                         |
//| in regression on a constant                                      |
//+------------------------------------------------------------------+
void prognose_interval0(double a, double c, int n, double conf_level)
  {
//--- Input data validation
   if(c <= 0.0)
     {
      Print("c must be positive in prognose_interval0()");
      return;
     }
   if(n < 2)
     {
      Print("n must be greater 1 in prognose_interval0()");
      return;
     }
   if(conf_level <= 0.0 || conf_level >= 1.0)
     {
      Print("conf_level must be between 0 and 1 in prognose_interval0()");
      return;
     }
   int err = 0;
//--- t-statistic calculation
//--- double t = MathQuantileT((1.0 - conf_level) / 2.0, n - 1, err);  // faulty library function
   double t = MathQuantileT_TMP((1.0 - conf_level) / 2.0, n - 1, err); // custom implementation
   if(err != 0)
     {
      Print("in MathQuantileT() error ", err);
      return;
     }
   t *= c * MathSqrt(1.0 + 1.0 / n);
//--- Print result if success
   PrintFormat("Y between %.4f and %.4f with confidence level %.3f", a + t, a - t, conf_level);
  }
//+------------------------------------------------------------------+
//| Computation of parameters and residuals for                      |
//| simple linear regression                                         |
//+------------------------------------------------------------------+
void regression1(double& y[], double& x[], double& a, double& b, double& c, double& e[])
  {
//--- Initialization and input data validation
   a = b = c = 0.0;
   ArrayResize(e, 0);
   int ny = ArraySize(y);
   if(ny < 3)
     {
      Print("no data for regression1()");
      return;
     }
   if(ny != ArraySize(x))
     {
      Print("arrays y[] and x[] have different lengths in regression1()");
      return;
     }
//--- Check x[] is not constant
   double Sx = MathStandardDeviation(x);
   if(Sx <= DBL_MIN)
     {
      Print("x[]==const, use regression0()");
      return;
     }
//--- Parameters A and B estimation
   if(!MathCorrelationPearson(y, x, b))
     {
      Print("MathCorrelationPearson() error in regression1()");
      return;
     }
   ArrayResize(e, ny);
   b *= MathStandardDeviation(y) / Sx;
   a = MathMean(y) - b * MathMean(x);
//--- Residuals computation
   for(int i = 0; i < ny; ++i)
      e[i] = y[i] - a - b * x[i];
//--- Parameter C estimation
   c = MathStandardDeviation(e) * MathSqrt(1.0 + 1.0 / (ny - 2.0));
  }
//+------------------------------------------------------------------+
//| Test of hypothesis H₀: B = 0 vs H₁: B ≠ 0 for                    |
//| simple linear regression                                         |
//+------------------------------------------------------------------+
void Beq0_test1(double b, double c, double var_x, int n)
  {
   Print("test H0: B = 0 vs H1: B != 0 (for pair regession) result:");
//--- Input data validation
   if(c <= 0.0)
     {
      Print("c must be positive in Beq0_test1()");
      return;
     }
   if(var_x <= 0.0)
     {
      Print("var_x must be positive in Beq0_test1()");
      return;
     }
   if(n < 3)
     {
      Print("n must be 3 or greater in Beq0_test1()");
      return;
     }
//--- t-statistic calculation
   double t = MathAbs(b * MathSqrt(var_x * (n - 1)) / c);
   int err = 0;
//--- p-value computation
   double p_value = 2.0 * MathCumulativeDistributionT(t, n - 2, false, false, err);
   if(err != 0)
     {
      Print("in MathCumulativeDistributionT() error ", err);
      return;
     }
//--- Print result if success
   PrintFormat("t = %.3f, p-value = %.3f", t, p_value);
  }
//+------------------------------------------------------------------+
//| Calculation of prediction interval for Y in                      |
//| simple linear regression                                         |
//+------------------------------------------------------------------+
void prognose_interval1(double x_new, double a, double b, double c, double var_x, double mean_x, int n, double conf_level)
  {
//--- Input data validation
   if(c <= 0.0)
     {
      Print("c must be positive in prognose_interval1()");
      return;
     }
   if(n < 2)
     {
      Print("n must be greater 2 in prognose_interval1()");
      return;
     }
   if(conf_level <= 0.0 || conf_level >= 1.0)
     {
      Print("conf_level must be between 0 and 1 in prognose_interval1()");
      return;
     }
   int err = 0;
//--- t-statistic calculation
//--- double t = MathQuantileT((1.0 - conf_level) / 2.0, n - 1, err);  // faulty library function
   double t = MathQuantileT_TMP((1.0 - conf_level) / 2.0, n - 1, err); // custom implementation
   if(err != 0)
     {
      Print("in MathQuantileT() error ", err);
      return;
     }
   t *= c * MathSqrt(1.0 + 1.0 / n + MathPow((x_new - mean_x), 2) * (n - 1) / var_x);
//--- Print result if success
   PrintFormat("Y between %.4f and %.4f with confidence level %.3f", a + b * x_new + t, a + b * x_new - t, conf_level);
  }
//+------------------------------------------------------------------+
//| Residuals plot vs. price bar index                               |
//+------------------------------------------------------------------+
void t_residuals_plot(double& residuals[])
  {
   int n_residuals = ArraySize(residuals);
//--- Input data validation
   if(n_residuals < 1)
     {
      Print("no data for t_residuals_plot()");
      return;
     }
//--- Array T[] as x on plot
   double T[];
   if(!MathSequenceByCount(1, n_residuals, n_residuals, T))
     {
      Print("MathSequenceByCount() error");
      return;
     }
//--- Plot construction
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   graphic.Create(0, "G", 0, 0, 0, 750, 300);
   graphic.CurveAdd(T, residuals, CURVE_LINES, "Residuals");
   graphic.CurvePlotAll();
   graphic.Update();
   MessageBox("Press OK to close graph", "Close graph");
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
  }
//+------------------------------------------------------------------+
//| EPDF of residuals vs. normal density with residual SD            |
//| nofx - number of points on plot                                  |
//+------------------------------------------------------------------+
void epdf_vs_normalpdf(double& residuals[], int nofx = 30)
  {
   int n_residuals = ArraySize(residuals);
//--- Input data validation
   if(n_residuals < 1)
     {
      Print("no data for epdf_vs_normalpdf()");
      return;
     }
//--- Array for x
   double x[];
   if(!MathSequenceByCount(MathMin(residuals), MathMax(residuals), nofx, x))
     {
      Print("sequence error");
      return;
     }
//--- Array for EPDF
   double epdf[];
   if(!MathProbabilityDensityEmpirical(residuals, nofx, x, epdf))
     {
      Print("MathProbabilityDensityEmpirical() error");
      return;
     }
//--- Array for theoretical PDF
   double normal_pdf[];
   if(!MathProbabilityDensityNormal(x, 0.0, MathStandardDeviation(residuals), normal_pdf))
     {
      Print("MathProbabilityDensityNormal() error");
      return;
     }
//--- Plot construction
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   graphic.Create(0, "G", 0, 0, 0, 750, 400);
   graphic.CurveAdd(x, epdf, CURVE_LINES, "Residuals");
   graphic.CurveAdd(x, normal_pdf, CURVE_LINES, "Normal");
   graphic.CurvePlotAll();
   graphic.Update();
   MessageBox("Press OK to close graph", "Close graph");
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
  }
//+------------------------------------------------------------------+
//| QQ-plot of residuals vs. normal distribution with residual SD    |
//+------------------------------------------------------------------+
void qq_plot(double& residuals[])
  {
   int n_residuals = ArraySize(residuals);
//--- Input data validation
   if(n_residuals < 1)
     {
      Print("no data for qq_plot()");
      return;
     }
//--- Copy input array and sort this copy
   double residuals_copy[];
   ArrayResize(residuals_copy, n_residuals);
   ArrayCopy(residuals_copy, residuals);
   ArraySort(residuals_copy);
//--- Probability sequence
   double pseq[];
   if(!MathSequenceByCount(0.5 / n_residuals, 1.0 - 0.5 / n_residuals, n_residuals, pseq))
     {
      Print("MathSequenceByCount() error");
      return;
     }
//--- Theoretical quantiles
   double normal_qs[];
   if(!MathQuantileNormal(pseq, 0.0, MathStandardDeviation(residuals_copy), normal_qs))
     {
      Print("MathQuantileNormal() error");
      return;
     }
//--- Plot construction
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   graphic.Create(0, "G", 0, 0, 0, 750, 550);
   graphic.CurveAdd(normal_qs, residuals_copy, CURVE_LINES, "QQ-plot");
   graphic.CurveAdd(residuals_copy, residuals_copy, CURVE_LINES, "Y=X");
   graphic.CurvePlotAll();
   graphic.Update();
   MessageBox("Press OK to close graph", "Close graph");
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
  }
//+------------------------------------------------------------------+
//| Correlogram of residuals                                         |
//+------------------------------------------------------------------+
void correlogram(double& residuals[])
  {
   int n_residuals = ArraySize(residuals);
//--- Input data validation
   if(n_residuals < 2)
     {
      Print("no data for correlogram()");
      return;
     }
//--- Number of lags
   int ncg = (int)MathSqrt(n_residuals);
//--- Significance levels
   double bound = 1.96 / MathSqrt(n_residuals);
   double cg0 = 0.0, xb[2] = {0, ncg + 1}, upb[2] = {bound, bound}, lowb[2] = {-bound, -bound};
//--- Lags array creation
   double lags[];
   if(!MathSequenceByCount(1, ncg, ncg, lags))
     {
      Print("MathSequenceByCount() error");
      return;
     }
//--- ACF array calculation
   double cg[];
   ArrayResize(cg, ncg);
   ArrayFill(cg, 0, ncg, 0.0);
   for(int i = 0; i < n_residuals; ++i)
      cg0 += residuals[i] * residuals[i];
   for(int k = 1; k <= ncg; ++k)
     {
      for(int i = 0; i < n_residuals - k; ++i)
         cg[k - 1] += residuals[i] * residuals[i + k];
      cg[k - 1] /= cg0;
     }
//--- Plot construction
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   graphic.Create(0, "G", 0, 0, 0, 750, 300);
   graphic.CurveAdd(lags, cg, CURVE_HISTOGRAM, "Correlogram");
   graphic.CurveAdd(xb, upb, CURVE_LINES, "Upper bound");
   graphic.CurveAdd(xb, lowb, CURVE_LINES, "Lower bound");
   graphic.CurvePlotAll();
   graphic.Update();
   MessageBox("Press OK to close graph", "Close graph");
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
  }
//+------------------------------------------------------------------+
//| Scatter plot of (x, y) points with line y = a * x + b            |
//+------------------------------------------------------------------+
void scatter_plot(double& x[], double& y[], bool add_line = false, double a = 0.0, double b = 0.0)
  {
   int ny = ArraySize(y);
//--- Input data validation
   if(ny < 1)
     {
      Print("no data for regression1()");
      return;
     }
   if(ny != ArraySize(x))
     {
      Print("arrays y[] and x[] have different lengths in scatter_plot()");
      return;
     }
//--- Regression line points calculation
   double xreg[2] = {MathMin(x), MathMax(x)};
   double yreg[2] = {a + b*xreg[0], a + b*xreg[1]};
//--- Plot construction
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   graphic.Create(0, "G", 0, 0, 0, 750, 550);
   graphic.CurveAdd(x, y, CURVE_POINTS, "scatter plot");
   if(add_line)
      graphic.CurveAdd(xreg, yreg, CURVE_LINES, "y = a * x + b");
   graphic.CurvePlotAll();
   graphic.Update();
   MessageBox("Press OK to close graph", "Close graph");
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
  }
//+------------------------------------------------------------------+
//| 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(double& residuals[], int m = 10, int pq = 0)
  {
   Print("Ljung-Box test result:");
//--- Input data validation
   if(pq < 0)
     {
      Print("pq must be nonnegative");
      return;
     }
   if(m <= pq)
     {
      Print("m must be greater pq");
      return;
     }
   int n_residuals = ArraySize(residuals);
   if(n_residuals < m + 1)
     {
      Print("not enough data");
      return;
     }
//--- Initialization
   double Q = 0.0, cov0 = 0.0, cov, ro;
//--- Test statistic calculation
   for(int i = 0; i < n_residuals; ++i)
      cov0 += residuals[i] * residuals[i];
   for(int k = 1; k <= m; ++k)
     {
      cov = 0.0;
      for(int i = 0; i < n_residuals - k; ++i)
         cov += residuals[i] * residuals[i + k];
      ro = cov / cov0;
      Q += ro * ro / (n_residuals - k);
     }
   Q *= n_residuals * (n_residuals + 2);
//--- p-value computation
   int err = 0;
   double p_value = MathCumulativeDistributionChiSquare(Q, m - pq, false, false, err);
   if(err != 0)
     {
      Print("MathCumulativeDistributionChiSquare() error: ", err);
      return;
     }
//--- Print result
   PrintFormat("Q = %.3f, p-value = %.3f", Q, p_value);
  }
//+------------------------------------------------------------------+
//| Jarque-Bera test for normality of distribution                   |
//+------------------------------------------------------------------+
void Jarque_Bera_test(double& residuals[])
  {
   Print("Jarque-Bera test result:");
//--- Input data validation
   int n_residuals = ArraySize(residuals);
   if(n_residuals < 1)
     {
      Print("not enough data");
      return;
     }
//--- Skewness calculation
   double S = MathSkewness(residuals);
//--- Kurtosis calculation
   double K = MathKurtosis(residuals);
//--- Test statistic calculation
   double JB = n_residuals * (S * S + K * K / 4.0) / 6.0;
//--- p-value computation
   double p_value = MathExp(-JB / 2.0);
//--- Print result
   PrintFormat("JB = %.3f, p-value = %.3f", JB, p_value);
  }
//+------------------------------------------------------------------+
//| Pettitt test for structural break (abrupt change in mean)        |
//+------------------------------------------------------------------+
void Pettitt_test(double& residuals[])
  {
   Print("Pettitt test result:");
//--- Input data validation
   int n_residuals = ArraySize(residuals);
   if(n_residuals < 2)
     {
      Print("not enough data");
      return;
     }
//--- Initialization
   double K = 0.0, U, W = 0.0;
   int t_change = 0;
//--- Sample ranking computation
   double rank[];
   MathRank(residuals, rank);
//--- Test statistic calculation
   for(int t = 1; t < n_residuals; ++t)
     {
      W += rank[t - 1];
      U = MathAbs(2.0 * W - t * (n_residuals + 1.0));
      if(U > K)
        {
         K = U;
         t_change = t;
        }
     }
//--- p-value computation
   double p_value = MathMin(2.0 * MathExp(-6.0 * K * K / (n_residuals * n_residuals * (n_residuals + 1.0))), 1.0);
//--- Print result
   PrintFormat("K = %.3f, p-value = %.3f, change-point = %d", K, p_value, t_change);
  }
//+------------------------------------------------------------------+
//| Student's t-distribution quantile                                |
//| (custom implementation replacing faulty library function)        |
//| JUST A TEMPORARY REPLACEMENT! (Based on library variant)         |
//+------------------------------------------------------------------+
double MathQuantileT_TMP(const double prob, const double nu, int &error_code)
  {
   if(!MathIsValidNumber(prob) || !MathIsValidNumber(nu))
     {
      error_code = ERR_ARGUMENTS_NAN;
      return QNaN;
     }
   if(nu <= 0.0 || nu != MathRound(nu))
     {
      error_code = ERR_ARGUMENTS_INVALID;
      return QNaN;
     }
   if(prob < 0.0 || prob > 1.0)
     {
      error_code = ERR_ARGUMENTS_INVALID;
      return QNaN;
     }
   error_code = ERR_RESULT_INFINITE;
   if(prob == 0.0)
      return QNEGINF;
   if(prob == 1.0)
      return QPOSINF;
   error_code = ERR_OK;
   if(prob == 0.5)
      return 0.0;
   if(nu == 1.0)
      return MathTan(M_PI * (prob - 0.5));
   int err_code = 0;
   double x_low = -1.0;
   double x_high = 1.0;
   while(MathCumulativeDistributionT(x_low, nu, err_code) > prob && x_low > -1e6)
      x_low *= 2.0;
   while(MathCumulativeDistributionT(x_high, nu, err_code) < prob && x_high < 1e6)
      x_high *= 2.0;
   double x = 0.0;
   int max_iterations = 100;
   double precision = 1e-12; // Высокая точность
   for(int i = 0; i < max_iterations; i++)
     {
      x = (x_low + x_high) / 2.0;
      double cdf = MathCumulativeDistributionT(x, nu, err_code);
      if(MathAbs(x_high - x_low) < precision)
         break;
      if(cdf < prob)
         x_low = x;
      else
         x_high = x;
      if(i == max_iterations - 1)
         error_code = ERR_NON_CONVERGENCE;
     }
   return x;
  }
//+------------------------------------------------------------------+
//| Helper function to export array to text file                     |
//| (for data analysis in other programs)                            |
//+------------------------------------------------------------------+
void array2csv(double& a[], string fldr, string fnm)
  {
//--- Initialization and input data validation
   int ftxt = FileOpen(fldr + "\\" + fnm, FILE_WRITE | FILE_TXT | FILE_ANSI | FILE_COMMON);
   if(ftxt == INVALID_HANDLE)
     {
      Print("FileOpen() error");
      return;
     }
   int na = ArraySize(a);
   if(na < 1)
     {
      Print("No data for array2csv()");
      return;
     }
//--- Writing array to file
   FileWriteString(ftxt, DoubleToString(a[0]));
   for(int i = 1; i < na; ++i)
      FileWriteString(ftxt, "\n" + DoubleToString(a[i]));
   FileClose(ftxt);
  }
//+------------------------------------------------------------------+