//+------------------------------------------------------------------+
//|                            TestScalingLambdaMaxBrownMovement.mq5 |
//|                             Copyright 2000-2026, MetaQuotes Ltd. |
//|                                              http://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2000-2026, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <Graphics\Graphic.mqh>
//+------------------------------------------------------------------+
//| Generate Brownian motion                                         |
//+------------------------------------------------------------------+
void GenerateBrownian(int N,vector<double> &data)
  {
   data.Resize(N);
   data[0] = 0.0;
   for(int i = 1; i < N; i++)
      data[i] = data[i - 1] + (MathRand()/32767.0 - 0.5);
  }
//+------------------------------------------------------------------+
//| LinearRegression                                                 |
//+------------------------------------------------------------------+
void LinearRegression(const double &x[], const double &y[], int n, double &a, double &b)
  {
   double sx = 0.0, sy = 0.0, sxx = 0.0, sxy = 0.0;
   for(int i = 0; i < n; i++)
     {
      sx  += x[i];
      sy  += y[i];
      sxx += x[i] * x[i];
      sxy += x[i] * y[i];
     }
   double denom = n * sxx - sx * sx;
   a = (n * sxy - sx * sy) / denom;
   b = (sy - a * sx) / n;
  }
//+------------------------------------------------------------------+
//| TestScaling with statistics                                      |
//+------------------------------------------------------------------+
void TestScalingStatistics()
  {
   MathSrand(42);
   int RUNS = 10;     //
   int MC   = 10;   // Monte Carlo
   double alpha_values[];
   ArrayResize(alpha_values, RUNS);
// --- geometric grid of T
   int nT = 8;
   int Tvals[];
   ArrayResize(Tvals, nT);
   int T0 = 64;
   for(int i = 0; i < nT; i++)
      Tvals[i] = T0 << i;
   Print("Scaling test with statistics");
//---
   double logT[];
   double logLambda[];
   vector<double> bm;
   vector<double> l1_trend;
   for(int run = 0; run < RUNS; run++)
     {
      ArrayResize(logT, nT);
      ArrayResize(logLambda, nT);
      //---
      for(int i = 0; i < nT; i++)
        {
         int T = Tvals[i];
         double lambda_sum = 0.0;
         l1_trend.Resize(T);
         for(int k = 0; k < MC; k++)
           {
            GenerateBrownian(T, bm);
            lambda_sum += bm.L1TrendFilterLambdaMax();
            bm.L1TrendFilter(l1_trend,0.2,true);
           }
         double lambda_avg = lambda_sum / MC;
         logT[i]      = MathLog((double)T);
         logLambda[i] = MathLog(lambda_avg);
        }
      // --- regression
      double alpha, c;
      LinearRegression(logT, logLambda, nT, alpha, c);
      alpha_values[run] = alpha;
      PrintFormat("run %d -> alpha = %.6f", run+1, alpha);
     }
//--- statistics
   double mean = 0.0;
   for(int i=0;i<RUNS;i++)
      mean += alpha_values[i];
   mean /= RUNS;
// --- standard deviation
   double var = 0.0;
   for(int i=0;i<RUNS;i++)
      var += (alpha_values[i] - mean)*(alpha_values[i] - mean);
   var /= (RUNS - 1);
   double stddev = MathSqrt(var);
// --- standard error of mean
   double sem = stddev / MathSqrt((double)RUNS);
// --- theoretical comparison
   double alpha_theory = 2.5;
   double percent_error = MathAbs(mean - alpha_theory) / alpha_theory * 100.0;
//--- results
   PrintFormat("mean alpha = %.6f", mean);
   PrintFormat("std deviation = %.6f", stddev);
   PrintFormat("standard error = %.6f", sem);
   PrintFormat("theory = %.4f", alpha_theory);
   PrintFormat("percent error from theory = %.4f %%", percent_error);
  }
//+------------------------------------------------------------------+
//| TestScaling                                                      |
//+------------------------------------------------------------------+
void TestScaling()
  {
   MathSrand(1);
// --- geometric grid of T
   int nT = 8;
   int Tvals[];
   ArrayResize(Tvals, nT);
//---
   int T0 = 64;
   for(int i = 0; i < nT; i++)
      Tvals[i] = T0 << i;   // 64 * 2^i
//---
   double logT[], logLambda[];
   ArrayResize(logT, nT);
   ArrayResize(logLambda, nT);
//---
   Print("scaling test for lambda_max");
   for(int i = 0; i < nT; i++)
     {
      int T = Tvals[i];
      //--- Monte-Carlo simulations
      int MC=1000;
      double lambda_sum = 0.0;
      for(int k = 0; k < MC; k++)
        {
         vector<double> bm;
         GenerateBrownian(T, bm);
         lambda_sum += bm.L1TrendFilterLambdaMax();
        }
      double lambda_avg = lambda_sum / MC;
      logT[i]      = MathLog((double)T);
      logLambda[i] = MathLog(lambda_avg);
      PrintFormat("T=%5d   <lambda_max>=%.6f", T, lambda_avg);
     }
// --- linear regression in log-log
   double alpha, c;
   LinearRegression(logT, logLambda, nT, alpha, c);
//---
   PrintFormat("estimated scaling exponent alpha = %.4f", alpha);
   double alpha_theory = 2.5;
   PrintFormat("theoretical value = %.4f", alpha_theory);
//--- plot scaling law
   CGraphic g;
   g.Create(0, "ScalingLaw", 0, 0, 0, 1000, 600);
   g.BackgroundMain("Scaling law of lambda_max (Brownian motion)");
   g.BackgroundMainSize(16);
   g.CurveAdd(logT, logLambda, CURVE_POINTS, "Simulation");
//---
   double xfit[2], yfit[2];
   xfit[0] = logT[0];
   xfit[1] = logT[nT - 1];
//---
   yfit[0] = alpha * xfit[0] + c;
   yfit[1] = alpha * xfit[1] + c;
//---least squares fit
   g.CurveAdd(xfit, yfit, CURVE_LINES, "LS fit");
   g.CurvePlotAll();
   g.Update();
   DebugBreak();
  }
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//--- calculate scaling with statistics
   TestScalingStatistics();
//--- show sample results
   TestScaling();
  }
//+------------------------------------------------------------------+