//+------------------------------------------------------------------+
//|                                                        tests.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"bds.mqh"
#include"..\..\cmath.mqh"
#include"..\..\Regression\utils.mqh"
#include"..\..\Regression\ols.mqh"
#include"..\..\np_graphics.mqh"
#include <Math\Alglib\fasttransforms.mqh>
#include <Math\Stat\T.mqh>
#include <Math\Stat\ChiSquare.mqh>
#include <Math\Stat\F.mqh>
#include <Math\Stat\Beta.mqh>
//+------------------------------------------------------------------+
//| statistical significance level calculation                       |
//+------------------------------------------------------------------+
enum ENUM_STAT_CONFIDENCE
{
 CONFIDENCE_99=1,//99%
 CONFIDENCE_95=5,//95%
 CONFIDENCE_90=10//90%
};

//+------------------------------------------------------------------+
//| Volatility proxy  specification                                  |
//+------------------------------------------------------------------+
enum ENUM_VOL_PROXY
  {
   VOL_ABSOLUTE_RETURNS,//Absolute returns
   VOL_SQUARED_RETURNS//Squared returns
  };
//+------------------------------------------------------------------+
//|Calculates the Hurst Exponent using Rescaled Range (R/S) Analysis.|
//+------------------------------------------------------------------+
double hurst_exponent_rs(vector& time_series, int partitions = 20)
  {

   int N = (int)time_series.Size();

   if(N<200)
     {
      Print(__FUNCTION__,": Time series too short for robust R/S analysis.\nMinimum length is 200");
      return EMPTY_VALUE;
     }

   int max_chunks = 0;

   if(partitions<10)
      partitions = 10;

   max_chunks = partitions;

   vector log_space = np::logspace(log10(50),log10(floor(N/2)),ulong(max_chunks));
   vector chunk_sizes = np::unique(log_space);

   vector rs_values,rslist,valid_sizes,chunk,y,z;
   double R,S;
   ulong num_chunks,rsl,rs;

   rs = 0;
   rs_values.Resize(0,chunk_sizes.Size());
   valid_sizes.Resize(0,chunk_sizes.Size());
   ulong size = 0;

   for(ulong j = 0; j<chunk_sizes.Size(); ++j)
     {
      size = (ulong)chunk_sizes[j];
      num_chunks = ulong(N)/ulong(size);
      rslist.Resize(0,num_chunks);
      rsl = 0;
      for(ulong i = 0; i<num_chunks; ++i)
        {
         chunk = np::sliceVector(time_series,long(i*size),long((i+1)*size));
         if(!chunk.Size())
           {
            Print(__FUNCTION__,": vector size ", time_series.Size(), " params: ", long(i*size), ",", long((i+1)*size));
            return EMPTY_VALUE;
           }
         y = chunk - chunk.Mean();
         z = y.CumSum();
         R = z.Max() - z.Min();
         S = chunk.Std(1);
         if(S>0.0)
           {
            if(rslist.Size()<rsl+1)
               rslist.Resize(rsl+1,num_chunks-(rsl+1));
            rslist[rsl++] = R/S;
           }
        }
      if(rslist.Size())
        {
         if(rs_values.Size()<rs+1)
           {
            rs_values.Resize(rs+1,chunk_sizes.Size()-(rs+1));
            valid_sizes.Resize(rs+1,chunk_sizes.Size()-(rs+1));
           }
         rs_values[rs] = rslist.Mean();
         valid_sizes[rs++] = double(size);
        }
     }
   valid_sizes = log(valid_sizes);
   rs_values = log(rs_values);
   vector poly = np::polyfit(valid_sizes,rs_values,1);

   return poly[0];
  }
//+------------------------------------------------------------------+
//| Geweke and Porter-Hudak (GPH) Test for fractional integration    |
//+------------------------------------------------------------------+
vector gph_test(vector& time_series, double alpha=0.5)
  {
   if(time_series.Size()<30)
     {
      Print(__FUNCTION__,": Sample size is too small to compute a meaningful periodogram.");
      return vector::Zeros(0);
     }
   vector x = time_series - time_series.Mean();

   CRowDouble real = x;
   CRowComplex dft;

   CFastFourierTransform::FFTR1D(real,real.Size(),dft);
   vectorc vec = dft.ToVector();
   vector periodogram = pow(c_abs(vec),2.)/(2.*M_PI*double(x.Size()));

   int g = int(floor(pow(x.Size(),alpha)));
   if(g<2)
      g = 2;

   vector frequencies = 2.*M_PI*np::arange(1.0,double(g+1),1.)/double(x.Size());
   vector raw = np::sliceVector(periodogram,1,long(g+1));
   for(ulong i = 0; i<raw.Size(); ++i)
      if(!raw[i])
         raw[i] = 1.e-13;

   vector y = log(raw);
   vector x_reg = -2.*log(2.*sin(frequencies/2.));
   LinRegressResult rr = linregress(x_reg,y);
   double denom = (pow(x_reg-x_reg.Mean(),2.)).Sum();
   if(!denom)
      return vector::Zeros(0);
   double theoretical_std_err = sqrt((pow(M_PI,2.)/6.0)/denom);
   double z_stat = rr.slope/theoretical_std_err;
   double p_val_theoretical = 2.*(1.-CNormalDistr::NormalCDF(fabs(z_stat)));
   vector out(2);
   out[0] = rr.slope;
   out[1] = p_val_theoretical;
   return out;
  }
//+------------------------------------------------------------------+
//|Run bds test                                                      |
//+------------------------------------------------------------------+
void run_bds_test(vector& residuals, ulong max_dim = 4, double epsilon_scale=1.5)
  {
   double epsilon = epsilon_scale*residuals.Std(1);
   BDSResult bdsr = bds(residuals,(int)max_dim,epsilon);
   PrintFormat("%15s %15s %15s","Dimension(m)","BDS Z-Stat","P-value");
   for(ulong i = 2,j = 0; i<max_dim+1; ++i,++j)
      PrintFormat("%15d %15.4f %15.4e",i,bdsr.bds_stats[j],bdsr.bds_pvalues[j]);
   return;
  }
//+------------------------------------------------------------------+
//| Msct result                                                      |
//+------------------------------------------------------------------+
struct MsctResult
  {
   double            avg_coarse_to_fine;
   double            avg_fine_to_coarse;
   double            cascade_delta;
   string            recommended_model;

                     MsctResult(void)
     {
      avg_coarse_to_fine = avg_fine_to_coarse = cascade_delta = 0.0;
      recommended_model = NULL;
     }

                     MsctResult(double avg_coarse,double avg_fine,double delta, string recommendation)
     {
      avg_coarse_to_fine = avg_coarse;
      avg_fine_to_coarse = avg_fine;
      cascade_delta = delta;
      recommended_model = recommendation;
     }

                     MsctResult(MsctResult & other)
     {
      avg_coarse_to_fine = other.avg_coarse_to_fine;
      avg_fine_to_coarse = other.avg_fine_to_coarse;
      cascade_delta = other.cascade_delta;
      recommended_model = other.recommended_model;
     }

   void              operator=(MsctResult &other)
     {
      avg_coarse_to_fine = other.avg_coarse_to_fine;
      avg_fine_to_coarse = other.avg_fine_to_coarse;
      cascade_delta = other.cascade_delta;
      recommended_model = other.recommended_model;
     }
  };
//+----------------------------------------------------------------------------+
//|Performs a cross-correlation analysis across distinct time horizons         |
//|to detect asymmetric volatility cascades (Coarse -> Fine vs Fine -> Coarse).|
//+----------------------------------------------------------------------------+
MsctResult multi_scale_correlation_test(vector& lowertf_close_prices,ulong highertf_window,ulong max_lags, double threshold_diagnostic = 0.05)
  {
   vector highertf_close_prices = np::sliceVector(lowertf_close_prices,long(highertf_window - 1),END,long(highertf_window));

   vector fine_returns = fabs(np::diff(log(lowertf_close_prices),1));
   vector coarse_returns = fabs(np::diff(log(highertf_close_prices),1));

   vector ctf_lags = vector::Zeros(max_lags);
   vector ftc_lags = ctf_lags;
   vector c,f;

   for(long lag = 1; lag<long(max_lags+1); ++lag)
     {
      c = np::sliceVector(coarse_returns,0,long(coarse_returns.Size() - lag));
      f = fine_returns;
      ctf_lags[lag - 1] = f.CorrCoef(c);

      f = np::sliceVector(fine_returns,0,long(fine_returns.Size() - (lag*highertf_window)));
      c = coarse_returns;
      ftc_lags[lag - 1] = c.CorrCoef(f);
     }

   MsctResult result;
   result.avg_coarse_to_fine = ctf_lags.Mean();
   result.avg_fine_to_coarse = ftc_lags.Mean();
   result.cascade_delta = result.avg_coarse_to_fine - result.avg_fine_to_coarse;

   if(result.cascade_delta > threshold_diagnostic)
      result.recommended_model = "HARCH(Asymmetric Mutli-Scale Cascade Detected)";
   else
      result.recommended_model = "FIGARCH(Symmetric Fractionally Integrated Dynamics)";
   return result;
  }
//+----------------------------------------------------------------------------------------+
//| Plot graph of GPH Fractional integration estimates over a range of bandwiths           |
//+----------------------------------------------------------------------------------------+
void plot_gph_test_stats(vector& time_series, double min_alpha = 0.4, double max_alpha = 0.6, double step = 0.01)
  {
   if(time_series.Size()<30)
     {
      Print(__FUNCTION__,": Sample size is too small to compute a meaningful stats.");
      return;
     }

   vector xvals = np::arange(min_alpha,max_alpha+step,step);
   vector yvals = vector::Zeros(xvals.Size());
   vector stats;
   for(ulong i = 0; i<yvals.Size(); ++i)
     {
      stats = gph_test(time_series,xvals[i]);
      yvals[i] = stats[0];
     }

   np::plotxy(xvals,yvals,"GPH test stats","Bandwidth(alpha)","Fractional Integration estimate(D)",false,0,0,10,10,750,500,true,1,CURVE_LINES,30);
  }
//+------------------------------------------------------------------+
//|Plot Hurst exponent values over a range of partition sizes        |
//+------------------------------------------------------------------+
void plot_hurst_exponet_estimates(vector& time_series, int min_partion_size = 20, int max_partion_size = 50, int step = 5)
  {
   if(time_series.Size()<200)
     {
      Print(__FUNCTION__,": Time series too short for robust R/S analysis.\nMinimum length is 200");
      return;
     }
   vector xvals = np::range(max_partion_size+step,min_partion_size,step);
   vector yvals = vector::Zeros(xvals.Size());
   for(ulong i = 0; i<yvals.Size(); ++i)
      yvals[i] = hurst_exponent_rs(time_series,(int)xvals[i]);

   np::plotxy(xvals,yvals,"Hurst Exponent stats","Min Partition Size","Hurst Exponent (H)",false,0,0,10,10,750,500,true,1,CURVE_LINES,30);
  }
//+------------------------------------------------------------------+
//|Computes and plots the log-periodogram vs log-frequency           |
//+------------------------------------------------------------------+
void plot_frequency_domain_fit(vector& time_series, double alpha=0.5, string plotname = NULL)
  {
   if(time_series.Size()<30)
     {
      Print(__FUNCTION__,": Time series too short for meaningful spectal analysis");
      return;
     }

   vector x = time_series - time_series.Mean();
   CRowDouble real = x;
   CRowComplex dft;

   CFastFourierTransform::FFTR1D(real,real.Size(),dft);
   vectorc vec = dft.ToVector();
   vector periodogram = pow(c_abs(vec),2.)/(2.*M_PI*double(x.Size()));

   ulong half_n = x.Size()/2;
   vector frequencies = 2.*M_PI * np::arange(1.0,double(half_n+1),1.0)/double(x.Size());
   vector power_spectrum = np::sliceVector(periodogram,1,long(half_n+1));

   if(!frequencies.Clip(1.e-13,frequencies.Max()) ||
      !power_spectrum.Clip(1.e-13,power_spectrum.Max()))
     {
      Print(__FUNCTION__,": Clip() vector method error ", GetLastError());
      return;
     }

   vector log_freq = log(frequencies);
   vector log_power = log(power_spectrum);

   int g = int(pow(x.Size(),alpha));
   g = fmax(2,fmin(g,int(log_freq.Size())));
   vector log_freq_low = np::sliceVector(log_freq,0,long(g));
   vector log_power_low = np::sliceVector(log_power,0,long(g));

   vector pf = np::polyfit(log_freq_low,log_power_low,1);
   vector trend_line_low = pf[0] * log_freq_low + pf[1];
   Print("Fraction Decay Fit, slope/d :", pf[0]);
   string plots[] = {"Log-Periodogram (Full Spectrum)","Long Memory Zone (Low Freq)","Fractional Decay Fit (Slope/d:"+DoubleToString(pf[0],2)+")"};
   vector yplots[3];
   yplots[0] = log_power;
   yplots[1] = log_power_low;
   yplots[2] = trend_line_low;
   np::plotxys(log_freq,yplots,plots,"Frequency Domain Periodogram Analysis: "+plotname,"Log frequency","Log Power(Spectral Density)",false,0,0,20,20,750,500,true,2,1,30);
   return;
  }

//+------------------------------------------------------------------+
//|Custom structure for model free EngleNg test                      |
//+------------------------------------------------------------------+
struct EngleNgResult
  {
   // Array order: [0]=const, [1]=sign_bias (a1), [2]=neg_size_bias (a2), [3]=pos_size_bias (a3)
   double            coefficients[4];
   double            t_stats[4];
   double            p_values[4];

   double            joint_F_stat;
   double            joint_F_pvalue;
   bool              significant_asymmetry;
   ENUM_STAT_CONFIDENCE confidence_level;

   double            a2_minus_a3;
   double            a2_vs_a3_t_stat;
   double            a2_vs_a3_p_value;
   bool              significant_asymmetry_purged_of_ARCH;

                     EngleNgResult(void)
     {
      joint_F_pvalue = joint_F_stat = a2_minus_a3 = a2_vs_a3_p_value = a2_vs_a3_t_stat = EMPTY_VALUE;
      ArrayFill(coefficients,0,0,EMPTY_VALUE);
      ArrayFill(t_stats,0,0,EMPTY_VALUE);
      ArrayFill(p_values,0,0,EMPTY_VALUE);
      significant_asymmetry_purged_of_ARCH = significant_asymmetry = false;
      confidence_level = CONFIDENCE_95;
     }
                     EngleNgResult(EngleNgResult& other)
     {
      joint_F_pvalue = other.joint_F_pvalue;
      joint_F_stat = other.joint_F_stat;
      a2_minus_a3 = other.a2_minus_a3;
      a2_vs_a3_p_value = other.a2_vs_a3_p_value;
      a2_vs_a3_t_stat = other.a2_vs_a3_t_stat;
      significant_asymmetry = other.significant_asymmetry;
      significant_asymmetry_purged_of_ARCH = other.significant_asymmetry_purged_of_ARCH;
      confidence_level = other.confidence_level;
      ArrayCopy(coefficients,other.coefficients);
      ArrayCopy(t_stats,other.t_stats);
      ArrayCopy(p_values,other.p_values);
     }
   void              operator=(EngleNgResult& other)
     {
      joint_F_pvalue = other.joint_F_pvalue;
      joint_F_stat = other.joint_F_stat;
      a2_minus_a3 = other.a2_minus_a3;
      a2_vs_a3_p_value = other.a2_vs_a3_p_value;
      a2_vs_a3_t_stat = other.a2_vs_a3_t_stat;
      significant_asymmetry = other.significant_asymmetry;
      significant_asymmetry_purged_of_ARCH = other.significant_asymmetry_purged_of_ARCH;
      confidence_level = other.confidence_level;
      ArrayCopy(coefficients,other.coefficients);
      ArrayCopy(t_stats,other.t_stats);
      ArrayCopy(p_values,other.p_values);
     }
  };
//+------------------------------------------------------------------+
//| Leverage correlation test results struct                         |
//+------------------------------------------------------------------+

struct LeverageCorrelationResult
  {
   int               lag;
   double            correlation;
   double            p_value;
   ENUM_STAT_CONFIDENCE confidence_level;
   bool              significant_leverage_effect;

                     LeverageCorrelationResult(void)
     {
      lag = -1;
      correlation = p_value = EMPTY_VALUE;
      significant_leverage_effect = false;
      confidence_level = CONFIDENCE_95;
     }
                     LeverageCorrelationResult(LeverageCorrelationResult& other)
     {
      lag = other.lag;
      correlation = other.correlation;
      p_value = other.p_value;
      significant_leverage_effect = other.significant_leverage_effect;
      confidence_level = other.confidence_level;
     }
   void              operator=(LeverageCorrelationResult& other)
     {
      lag = other.lag;
      correlation = other.correlation;
      p_value = other.p_value;
      significant_leverage_effect = other.significant_leverage_effect;
      confidence_level = other.confidence_level;
     }
  };

//+------------------------------------------------------------------+
//| Model free Engle-Ng Sign Bias Test                               |
//+------------------------------------------------------------------+
EngleNgResult EngleNgSignBiasTest(const vector &returns, ENUM_STAT_CONFIDENCE confidence_level = CONFIDENCE_95)
  {
   EngleNgResult res;
   res.confidence_level = confidence_level;
   ulong total_n = returns.Size();
   if(total_n < 5)
     {
      Print("Error: Insufficient data for Engle-Ng Test.");
      return res;
     }

   ulong n = total_n - 1;
   ulong k = 4; // number of regressors (const, S_neg, S_neg*r_lag, S_pos*r_lag)
   double dof = (double)(n - k);

// Form dependent variable y = r_t^2 (for t = 1 to N-1)
   vector y;
   y.Init(n);
   for(ulong i = 0; i < n; i++)
     {
      y[i] = MathPow(returns[i + 1], 2.0);
     }

// Form design matrix X
   matrix X;
   X.Init(n, k);
   for(ulong i = 0; i < n; i++)
     {
      double r_lag = returns[i];
      double s_neg = (r_lag < 0.0) ? 1.0 : 0.0;
      double s_pos = 1.0 - s_neg;

      X[i, 0] = 1.0;            // const
      X[i, 1] = s_neg;          // S_neg
      X[i, 2] = s_neg * r_lag;  // S_neg * r_lag
      X[i, 3] = s_pos * r_lag;  // S_pos * r_lag
     }

// OLS: beta = (X^T * X)^-1 * X^T * y
   matrix XT = X.Transpose();
   matrix XTX = XT.MatMul(X);
   matrix XTX_inv = XTX.Inv();
   vector XTy = XT.MatMul(y);
   vector beta = XTX_inv.MatMul(XTy);

// Residuals: resid = y - X * beta
   vector y_hat = X.MatMul(beta);
   vector resid = y - y_hat;

// White (HC1) Heteroskedasticity-Robust Covariance Matrix
// meat = X.T * diag(resid^2) * X
   matrix diag_resid_sq;
   diag_resid_sq.Init(n, n);
   diag_resid_sq.Fill(0);
   for(ulong i = 0; i < n; i++)
     {
      diag_resid_sq[i, i] = resid[i] * resid[i];
     }

   matrix meat = XT.MatMul(diag_resid_sq).MatMul(X);
   double scale = (double)n / dof;
   matrix cov_beta = XTX_inv.MatMul(meat).MatMul(XTX_inv) * scale;

// Standard Errors, t-stats and p-values
   int err_code = 0;
   for(ulong j = 0; j < k; j++)
     {
      res.coefficients[j] = beta[j];
      double se = MathSqrt(cov_beta[j, j]);
      res.t_stats[j] = (se == 0.0) ? 0.0 : (beta[j] / se);

      double cdf = MathCumulativeDistributionT(MathAbs(res.t_stats[j]), dof, err_code);
      if(err_code)
        {
         Print(__FUNCTION__,": CDF calculation error.");
         return EngleNgResult();
        }
      res.p_values[j] = 2.0 * (1.0 - cdf);
     }

// Joint F-test (Wald test) that a1 = a2 = a3 = 0
// R matrix (3 x 4)
   matrix R;
   R.Init(3, 4);
   R.Fill(0.0);
   R[0, 1] = 1.0;
   R[1, 2] = 1.0;
   R[2, 3] = 1.0;

   vector r_beta = R.MatMul(beta);
   matrix RT = R.Transpose();
   matrix R_cov_RT = R.MatMul(cov_beta).MatMul(RT);
   matrix R_cov_RT_inv = R_cov_RT.Inv();

   double wald = r_beta.Dot(R_cov_RT_inv.MatMul(r_beta));
   res.joint_F_stat = wald / 3.0;

   double f_cdf = MathCumulativeDistributionF(res.joint_F_stat, 3.0, dof, err_code);
   if(err_code)
     {
      Print(__FUNCTION__,": CDF calculation error.");
      return EngleNgResult();
     }
   res.joint_F_pvalue = 1.0 - f_cdf;
   res.significant_asymmetry = (res.joint_F_pvalue < (double(res.confidence_level)/100.0));

// a2 vs a3 comparison: c = [0, 0, 1, -1]
   vector c_vec;
   c_vec.Init(4);
   c_vec[0] = 0.0;
   c_vec[1] = 0.0;
   c_vec[2] = 1.0;
   c_vec[3] = -1.0;

   res.a2_minus_a3 = beta[2] - beta[3];
   double se_diff = MathSqrt(c_vec.Dot(cov_beta.MatMul(c_vec)));
   res.a2_vs_a3_t_stat = (se_diff == 0.0) ? 0.0 : (res.a2_minus_a3 / se_diff);

   double t_diff_cdf = MathCumulativeDistributionT(MathAbs(res.a2_vs_a3_t_stat), dof, err_code);
   if(err_code)
     {
      Print(__FUNCTION__,": CDF calculation error.");
      return EngleNgResult();
     }
   res.a2_vs_a3_p_value = 2.0 * (1.0 - t_diff_cdf);
   res.significant_asymmetry_purged_of_ARCH = (res.a2_vs_a3_p_value < (double(res.confidence_level)/100.0));

   return res;
  }

//+------------------------------------------------------------------+
//| Leverage Correlation Test                                        |
//+------------------------------------------------------------------+
LeverageCorrelationResult LeverageCorrelationTest(const vector &returns, int lag = 1, ENUM_STAT_CONFIDENCE confidence_level = CONFIDENCE_95)
  {
   LeverageCorrelationResult res;
   res.lag = lag;
   res.confidence_level = confidence_level;
   ulong total_n = returns.Size();
   if(total_n <= (ulong)lag)
     {
      Print("Error: Data length is shorter than the requested lag.");
      ZeroMemory(res);
      return res;
     }

   ulong n = total_n - (ulong)lag;
   vector x, y;
   x.Init(n);
   y.Init(n);

   for(ulong i = 0; i < n; i++)
     {
      x[i] = returns[i];
      y[i] = MathPow(returns[i + (ulong)lag], 2.0);
     }

// Pearson correlation calculation
   double mean_x = x.Mean();
   double mean_y = y.Mean();

   double num = 0.0;
   double den_x = 0.0;
   double den_y = 0.0;

   for(ulong i = 0; i < n; i++)
     {
      double diff_x = x[i] - mean_x;
      double diff_y = y[i] - mean_y;
      num += diff_x * diff_y;
      den_x += diff_x * diff_x;
      den_y += diff_y * diff_y;
     }

   if(den_x == 0.0 || den_y == 0.0)
     {
      res.correlation = 0.0;
      res.p_value = 1.0;
     }
   else
     {
      res.correlation = num / MathSqrt(den_x * den_y);

      // Exact p-value for Pearson correlation using Beta Cumulative Distribution
      // t = r * sqrt((n-2) / (1-r^2)). Equivalent beta relationship:
      // p_value = MathCumulativeDistributionBeta(1 - r^2, (n-2)/2, 1/2)
      double r_sq = MathPow(res.correlation, 2.0);
      if(r_sq >= 1.0)
        {
         res.p_value = 0.0;
        }
      else
        {
         int err_code = 0;
         double a = ((double)n - 2.0) / 2.0;
         double b = 0.5;
         res.p_value = MathCumulativeDistributionBeta(1.0 - r_sq, a, b, err_code);
         if(err_code)
           {
            Print(__FUNCTION__,": Beta Distribution CDF calculation error.");
            return LeverageCorrelationResult();
           }
        }
     }

   res.significant_leverage_effect = (res.p_value < (double(res.confidence_level)/100.0)) && (res.correlation < 0.0);
   return res;
  }


//+------------------------------------------------------------------+
//|Realized semivariance test results struct                         |
//+------------------------------------------------------------------+
struct RealizedSemivarianceResult
  {
   int               window;
   int               n_blocks;
   double            mean_RS_minus;
   double            mean_RS_plus;
   double            mean_diff;
   double            t_stat;
   double            p_value;
   bool              significant_asymmetry;
   ENUM_STAT_CONFIDENCE confidence_level;

                     RealizedSemivarianceResult(void)
     {
      window = n_blocks = -1;
      mean_RS_minus = mean_RS_plus = mean_diff = t_stat = p_value = EMPTY_VALUE;
      significant_asymmetry = false;
      confidence_level = CONFIDENCE_95;
     }
                     RealizedSemivarianceResult(RealizedSemivarianceResult& other)
     {
      window = other.window;
      n_blocks = other.n_blocks;
      mean_RS_minus = other.mean_RS_minus;
      mean_RS_plus = other.mean_RS_plus;
      mean_diff = other.mean_diff;
      t_stat = other.t_stat;
      p_value = other.p_value;
      significant_asymmetry = other.significant_asymmetry;
      confidence_level = other.confidence_level;
     }
   void              operator=(RealizedSemivarianceResult& other)
     {
      window = other.window;
      n_blocks = other.n_blocks;
      mean_RS_minus = other.mean_RS_minus;
      mean_RS_plus = other.mean_RS_plus;
      mean_diff = other.mean_diff;
      t_stat = other.t_stat;
      p_value = other.p_value;
      significant_asymmetry = other.significant_asymmetry;
      confidence_level = other.confidence_level;
     }
  };

//+------------------------------------------------------------------+
//| Volatility runs test results container                           |
//+------------------------------------------------------------------+

struct VolatilityRunsResult
  {
   matrix            contingency_table;
   double            p_high_vol_prev_down;
   double            p_high_vol_prev_up;
   double            chi2_stat;
   double            p_value;
   ENUM_STAT_CONFIDENCE confidence_level;
   bool              significant_asymmetry;

                     VolatilityRunsResult(void)
     {
      contingency_table = matrix::Zeros(2,2);
      p_high_vol_prev_down = p_high_vol_prev_up = chi2_stat = p_value = EMPTY_VALUE;
      significant_asymmetry = false;
      confidence_level = CONFIDENCE_95;
     }
                     VolatilityRunsResult(VolatilityRunsResult& other)
     {
      contingency_table = other.contingency_table;
      p_high_vol_prev_down = other.p_high_vol_prev_down;
      p_high_vol_prev_up = other.p_high_vol_prev_up;
      chi2_stat = other.chi2_stat;
      p_value = other.p_value;
      significant_asymmetry = other.significant_asymmetry;
      confidence_level = other.confidence_level;
     }
   void              operator=(VolatilityRunsResult& other)
     {
      contingency_table = other.contingency_table;
      p_high_vol_prev_down = other.p_high_vol_prev_down;
      p_high_vol_prev_up = other.p_high_vol_prev_up;
      chi2_stat = other.chi2_stat;
      p_value = other.p_value;
      significant_asymmetry = other.significant_asymmetry;
      confidence_level = other.confidence_level;
     }
  };

//+------------------------------------------------------------------+
//| Realized Semivariance Asymmetry Test                              |
//+------------------------------------------------------------------+
RealizedSemivarianceResult RealizedSemivarianceAsymmetry(const vector &returns, int window = 22, ENUM_STAT_CONFIDENCE confidence_level = CONFIDENCE_95)
  {
   RealizedSemivarianceResult res;
   res.window = window;
   res.confidence_level = confidence_level;
   ulong total_elements = returns.Size();
   res.n_blocks = (int)(total_elements / (ulong)window);

   if(res.n_blocks <= 1)
     {
      Print("Error: Not enough data blocks for t-test.");
      ZeroMemory(res);
      return res;
     }

   vector rs_neg, rs_pos, diff;
   rs_neg.Init(res.n_blocks);
   rs_pos.Init(res.n_blocks);
   diff.Init(res.n_blocks);

// Loop through non-overlapping blocks
   for(int i = 0; i < res.n_blocks; i++)
     {
      double sum_neg = 0.0;
      double sum_pos = 0.0;

      for(int j = 0; j < window; j++)
        {
         double r = returns[i * window + j];
         if(r < 0.0)
            sum_neg += r * r;
         if(r > 0.0)
            sum_pos += r * r;
        }

      rs_neg[i] = sum_neg;
      rs_pos[i] = sum_pos;
      diff[i]  = sum_neg - sum_pos;
     }

// Statistics calculations
   res.mean_RS_minus = rs_neg.Mean();
   res.mean_RS_plus  = rs_pos.Mean();
   res.mean_diff     = diff.Mean();

   double variance = diff.Var();
   double std_err  = MathSqrt(variance / res.n_blocks);

   if(std_err == 0.0)
     {
      res.t_stat = 0.0;
      res.p_value = 1.0;
     }
   else
     {
      res.t_stat = res.mean_diff / std_err;

      // Compute 2-tailed p-value via MQL5 standard library
      int err_code = 0;
      double cdf = MathCumulativeDistributionT(MathAbs(res.t_stat), res.n_blocks - 1, err_code);
      if(err_code)
        {
         Print(__FUNCTION__,": T-distribution CDF calculation error.");
         return RealizedSemivarianceResult();
        }
      res.p_value = 2.0 * (1.0 - cdf);
     }

   res.significant_asymmetry = (res.p_value < (double(res.confidence_level)/100.0)) && (res.mean_diff > 0.0);
   return res;
  }

//+------------------------------------------------------------------+
//| Volatility Runs / Proportion Test Asymmetry Test                 |
//+------------------------------------------------------------------+
VolatilityRunsResult VolatilityRunsAsymmetryTest(const vector &returns_or_resids, ENUM_STAT_CONFIDENCE confidence_level = CONFIDENCE_95 )
  {
   VolatilityRunsResult res;
   res.confidence_level = confidence_level;
   ulong n = returns_or_resids.Size();
   if(n <= 1)
      return res;

// 1. Calculate absolute returns_or_resids and extract median
   vector abs_r = fabs(returns_or_resids);
// MQL5 native vector sorting to fetch the median cleanly
   double median_abs = abs_r.Median();
// 2. Build Contingency Table elements
   double a = 0, b = 0, c = 0, d = 0;
   for(ulong i = 0; i < n - 1; i++)
     {
      bool prev_sign_neg = (returns_or_resids[i] < 0.0);
      bool next_high_vol  = (abs_r[i + 1] > median_abs);

      if(prev_sign_neg && next_high_vol)
         a++;
      else
         if(prev_sign_neg && !next_high_vol)
            b++;
         else
            if(!prev_sign_neg && next_high_vol)
               c++;
            else
               if(!prev_sign_neg && !next_high_vol)
                  d++;
     }

   res.contingency_table[0][0] = a;
   res.contingency_table[0][1] = b;
   res.contingency_table[1][0] = c;
   res.contingency_table[1][1] = d;

// Proportions
   res.p_high_vol_prev_down = (a + b > 0) ? (a / (a + b)) : double("nan");
   res.p_high_vol_prev_up   = (c + d > 0) ? (c / (c + d)) : double("nan");

// 3. Chi-Square Test with Yates' continuity correction
   double total = a + b + c + d;
   if(total == 0 || (a+b) == 0 || (c+d) == 0 || (a+c) == 0 || (b+d) == 0)
     {
      res.chi2_stat = 0.0;
      res.p_value = 1.0;
      return res;
     }

// Shortcut formula for 2x2 Chi-Square with Yates' Continuity Correction
   double numerator = MathAbs(a * d - b * c) - (total / 2.0);
   if(numerator < 0)
      numerator = 0; // cannot be negative before squaring

   double denominator = (a + b) * (c + d) * (a + c) * (b + d);
   res.chi2_stat = (total * MathPow(numerator, 2.0)) / denominator;

// Compute p-value for 1 degree of freedom (2x2 table)
   int err_code = 0;
// MathCumulativeDistributionChiSquare returns_or_resids P(X <= x). We want upper tail P(X >= x)
   double cdf = MathCumulativeDistributionChiSquare(res.chi2_stat, 1.0, err_code);
// Check for errors
   if(err_code)
     {
      Print(__FUNCTION__, ": Chisquare CDF calculation error.");
      return VolatilityRunsResult();
     }
//---
   res.p_value = 1.0 - cdf;
//---
   res.significant_asymmetry = (res.p_value < (double(res.confidence_level)/100.0)) && (res.p_high_vol_prev_down > res.p_high_vol_prev_up);
//---
   return res;
  }

//+------------------------------------------------------------------+