//+------------------------------------------------------------------+ //| VarianceRatioScript.mq5 | //| Lo-MacKinlay Variance Ratio - Diagnostic Tool | //| | //| Runs a multi-horizon Variance Ratio test on the current chart's | //| symbol/timeframe and prints the results to the Experts log. | //| A script runs once per launch, so this is an on-demand tool. | //+------------------------------------------------------------------+ #property copyright "Quantitative Series" #property version "1.02" #property script_show_inputs #include //--- Inputs input int InpWindowSize = 200; // Rolling sample size (1-period returns) input string InpLags = "2,5,10,20"; // Comma-separated lags to scan input double InpCriticalZ = 1.96; // Two-sided standard-normal critical Z threshold //+-------------------------------------------------------------------+ //| NormSInv | //| Inverse standard normal CDF (Acklam's approximation). See | //| article, Section 4.3, for the Bonferroni correction it supports. | //+-------------------------------------------------------------------+ double NormSInv(const double p) { double a1=-3.969683028665376e+01, a2= 2.209460984245205e+02, a3=-2.759285104469687e+02; double a4= 1.383577518672690e+02, a5=-3.066479806614716e+01, a6= 2.506628277459239e+00; double b1=-5.447609879822406e+01, b2= 1.615858368580409e+02, b3=-1.556989798598866e+02; double b4= 6.680131188771972e+01, b5=-1.328068155288572e+01; double c1=-7.784894002430293e-03, c2=-3.223964580411365e-01, c3=-2.400758277161838e+00; double c4=-2.549732539343734e+00, c5= 4.374664141464968e+00, c6= 2.938163982698783e+00; double d1= 7.784695709041462e-03, d2= 3.224671290700398e-01, d3= 2.445134137142996e+00; double d4= 3.754408661907416e+00; double plow = 0.02425; if(p <= 0.0) return -DBL_MAX; if(p >= 1.0) return DBL_MAX; if(p < plow) { double q = MathSqrt(-2.0 * MathLog(p)); return (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1.0); } if(p > 1.0 - plow) return -NormSInv(1.0 - p); double q = p - 0.5; double r = q * q; return (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1.0); } //+------------------------------------------------------------------+ //| NormCDF | //| Standard normal cumulative distribution function. | //+------------------------------------------------------------------+ double NormCDF(const double z) { //--- Abramowitz-Stegun approximation to the standard normal CDF. double sign = 1.0; if(z < 0.0) sign = -1.0; double x = MathAbs(z); double t = 1.0 / (1.0 + 0.2316419 * x); double d = 0.3989422804014327 * MathExp(-0.5 * x * x); double p = d * t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429)))); if(sign > 0.0) return 1.0 - p; return p; } //+------------------------------------------------------------------+ //| GetInterpretation | //| | //| Classifies the robust Z-statistic against the selected critical | //| threshold. The label describes only whether VR(k) departs | //| significantly from 1 at THIS single lag. It does not identify a | //| market regime, a mean-reversion effect or a trading opportunity. | //+------------------------------------------------------------------+ string GetInterpretation(const double z_hetero, const double threshold) { if(z_hetero < -threshold) return "SIGNIFICANT VR < 1"; if(z_hetero > threshold) return "SIGNIFICANT VR > 1"; return "NO SIGNIFICANT DEPARTURE"; } //+------------------------------------------------------------------+ //| ParseLags | //| Splits the comma-separated InpLags input into an int array. | //| Rejects lags below 2 and drops duplicates so the multiple- | //| testing correction counts each horizon once. | //| Returns the number of valid lags parsed, or -1 on failure. | //+------------------------------------------------------------------+ int ParseLags(const string input_str, int &lags_out[]) { string parts[]; int count = StringSplit(input_str, ',', parts); if(count <= 0) return -1; ArrayResize(lags_out, count); int valid = 0; for(int i = 0; i < count; i++) { string trimmed = parts[i]; StringTrimLeft(trimmed); StringTrimRight(trimmed); if(StringLen(trimmed) == 0) continue; int lag = (int)StringToInteger(trimmed); if(lag < 2) { PrintFormat("VarianceRatioScript: ignoring invalid lag '%s' (must be >= 2).", trimmed); continue; } bool duplicate = false; for(int j = 0; j < valid; j++) if(lags_out[j] == lag) { duplicate = true; break; } if(duplicate) { PrintFormat("VarianceRatioScript: ignoring duplicate lag %d.", lag); continue; } lags_out[valid] = lag; valid++; } ArrayResize(lags_out, valid); return valid; } //+------------------------------------------------------------------+ //| OnStart | //+------------------------------------------------------------------+ void OnStart(void) { //--- 1. Parse the requested lags int lags[]; int lagCount = ParseLags(InpLags, lags); if(lagCount <= 0) { Print("VarianceRatioScript: no valid lags to test. Check the InpLags input (e.g. \"2,5,10,20\")."); return; } //--- 2. Configure the engine, then fetch chronological price history. //--- The engine may clamp the requested window, so its effective value is //--- read back rather than assuming InpWindowSize was accepted. //--- start_pos = 1 skips bar 0, which is still forming, so only fully //--- closed bars enter the return series. Ordering is verified through //--- each bar's .time field instead of being assumed. CVarianceRatio vrEngine; vrEngine.SetWindowSize(InpWindowSize); int effectiveWindow = vrEngine.GetWindowSize(); int needed = effectiveWindow + 1; // window_size + 1 prices -> window_size returns MqlRates rates[]; int copied = CopyRates(_Symbol, _Period, 1, needed + 10, rates); if(copied < needed) // also catches copied == -1 { PrintFormat("VarianceRatioScript: insufficient history on %s (%s) -- copied %d bars, need at least %d.", _Symbol, EnumToString((ENUM_TIMEFRAMES)_Period), copied, needed); return; } bool alreadyChronological = (rates[0].time < rates[copied - 1].time); double closeArr[]; ArrayResize(closeArr, copied); for(int i = 0; i < copied; i++) closeArr[i] = alreadyChronological ? rates[i].close : rates[copied - 1 - i].close; //--- 3. Batch scan: CalculateMultiLag() returns the robust Z-score //--- for each requested lag in a single call, with an explicit //--- per-lag success flag (ok_batch) rather than relying on comparing //--- doubles against EMPTY_VALUE. double vr_batch[], z_batch[]; bool ok_batch[]; bool batchOk = vrEngine.CalculateMultiLag(closeArr, lags, vr_batch, z_batch, ok_batch); int bestIdx = -1; double bestAbsZ = -1.0; int testedCount = 0; if(batchOk) { for(int i = 0; i < lagCount; i++) { if(!ok_batch[i]) continue; testedCount++; double absZ = MathAbs(z_batch[i]); if(absZ > bestAbsZ) { bestAbsZ = absZ; bestIdx = i; } } } //--- 3b. Bonferroni-adjusted threshold for multiple testing. alpha is //--- derived from InpCriticalZ so the naive and adjusted thresholds refer //--- to the same significance level, and only lags that actually //--- computed are counted as tests. double alpha = 2.0 * (1.0 - NormCDF(InpCriticalZ)); double bonferroniAlpha = (testedCount > 0) ? alpha / testedCount : alpha; double bonferroniZ = -NormSInv(bonferroniAlpha / 2.0); // two-sided //--- 4. Detailed per-lag report: call Calculate() individually so we can //--- print both the homoskedastic Z(k) and the robust Z*(k) side by //--- side for comparison, alongside VR(k) and the interpretation. PrintFormat("================ Variance Ratio Diagnostic: %s (%s) ================", _Symbol, EnumToString((ENUM_TIMEFRAMES)_Period)); PrintFormat("Window = %d returns (%d closes) | Naive critical Z = +/-%.2f | \nBonferroni-adjusted (m=%d tests) critical Z ~= +/-%.2f", effectiveWindow, effectiveWindow + 1, InpCriticalZ, testedCount, bonferroniZ); PrintFormat("%-6s | %-10s | %-12s | %-12s | %-26s", "Lag k", "VR(k)", "Z(k) homo", "Z*(k) robust", "Interpretation"); Print("--------------------------------------------------------------------------------"); for(int i = 0; i < lagCount; i++) { double vr, z_homo, z_hetero; bool ok = vrEngine.Calculate(closeArr, lags[i], vr, z_homo, z_hetero); if(!ok) { PrintFormat("%-6d | %-10s | %-12s | %-12s | %-26s", lags[i], "n/a", "n/a", "n/a", "SKIPPED"); continue; } string interpretation = GetInterpretation(z_hetero, InpCriticalZ); PrintFormat("%-6d | %10.4f | %12.4f | %12.4f | %-26s", lags[i], vr, z_homo, z_hetero, interpretation); } Print("--------------------------------------------------------------------------------"); if(bestIdx >= 0) { string bestInterpretation = GetInterpretation(z_batch[bestIdx], InpCriticalZ); bool bonferroniSignificant = (bestAbsZ > bonferroniZ); PrintFormat("Strongest signal: k=%d VR(k)=%.4f Z*(k)=%.4f [%s]", lags[bestIdx], vr_batch[bestIdx], z_batch[bestIdx], bestInterpretation); PrintFormat(" -> descriptive only: %s the Bonferroni-adjusted threshold (+/-%.2f) \nacross %d lags actually tested (of %d requested).", bonferroniSignificant ? "EXCEEDS" : "does NOT exceed", bonferroniZ, testedCount, lagCount); } else Print("Strongest signal: unavailable (all lags failed to compute)."); Print("================================================================================"); } //+------------------------------------------------------------------+