Beyond GARCH (Part IV): Partition Analysis in MQL5

MetaTrader 5 — Trading systems | 30 May 2026, 11:42

489

Muhammad Minhas Qamar

Introduction

In the first three parts of this series, we built the entire MMAR pipeline in Python. We loaded EURUSD data, confirmed multifractal scaling through partition function analysis, extracted the Hurst exponent, and fitted the singularity spectrum. We then constructed the multiplicative cascade and a Fractional Brownian Motion process, ran Monte Carlo simulations, and put the MMAR head-to-head against GARCH. The MMAR performed better. The concept is supported by these tests.

Now we face a different kind of challenge: engineering. Python is excellent for research and prototyping, but it does not run natively inside MetaTrader 5. If we want the MMAR to drive real-time trading decisions, respond to every tick, and integrate seamlessly with the Strategy Tester, we need a native MQL5 implementation. No external dependencies, no bridge scripts, no latency from inter-process communication. Just a clean library that any Expert Advisor can include and call.

This article begins that construction. We will build the first major analytical module of the MMAR library: the Partition Analysis engine. This component takes raw price data and computes the partition function across multiple time scales and moment orders. It then extracts the scaling function tau(q), estimates the Hurst exponent H, and runs diagnostics to assess whether the data exhibits genuine multifractal behavior. By the end, we will have a working MQL5 script that loads historical bars, feeds them through our analysis pipeline, and reports whether the market data is multifractal.

We will cover:

Library Architecture and Constants

Before writing any analytical code, we need a shared vocabulary. The MMAR library is composed of several modules that build on each other, and they all reference the same set of enumerations, structures, and type definitions. We store these in a single header file, MMARConstants.mqh, which every other module includes.

The architecture follows a layered design. At the bottom sit the utility modules: an OLS regression class and a Hurst exponent estimator. Above them is the Partition Analysis engine, which uses both utilities to compute the scaling function and evaluate multifractality. In later articles, we will add the Spectrum Fitter (which consumes the tau(q) curve produced by Partition Analysis), the Simulation Engine (which uses the fitted distribution to generate synthetic paths), and the Monte Carlo forecaster (which orchestrates many simulations). At the top sits a facade class, CMMAR, that ties everything together with a simple Fit-and-Forecast interface. For now, we build from the bottom up.

Figure 1 shows this layered architecture:

MMAR library architecture with the bottom three layers highlighted

Fig. 1. The MMAR library architecture: this article covers the bottom three modules highlighted in gold

Note: the OLS regression class and the Generalized Hurst Exponent estimator are inspired by the article Implementing the Generalized Hurst Exponent and the Variance Ratio test in MQL5. That article provided the foundation for our OLS implementation (which we use as a self-contained regression container) and the GHE estimation approach. We have trimmed and modified both modules to fit our specific use case within the MMAR pipeline, but credit belongs to the original work for establishing the MQL5 implementations we built upon.

Let us begin with the constants file:

//+------------------------------------------------------------------+
//|                                               MMARConstants.mqh  |
//|     MMAR Library - Constants, Enumerations, and Data Structures  |
//+------------------------------------------------------------------+
#property strict

#ifndef MMAR_CONSTANTS_MQH
#define MMAR_CONSTANTS_MQH

//+------------------------------------------------------------------+
//| Cascade multiplier distribution types                            |
//+------------------------------------------------------------------+
enum ENUM_MMAR_DISTRIBUTION
  {
   MMAR_DIST_NORMAL = 0,
   MMAR_DIST_BINOMIAL,
   MMAR_DIST_POISSON,
   MMAR_DIST_GAMMA
  };

//+------------------------------------------------------------------+
//| Fractional Brownian Motion synthesis methods                     |
//+------------------------------------------------------------------+
enum ENUM_MMAR_FBM_METHOD
  {
   MMAR_FBM_DAVIES_HARTE = 0,
   MMAR_FBM_CHOLESKY
  };

//+------------------------------------------------------------------+
//| Model fitting status                                             |
//+------------------------------------------------------------------+
enum ENUM_MMAR_STATUS
  {
   MMAR_STATUS_NOT_FITTED = 0,
   MMAR_STATUS_FITTED,
   MMAR_STATUS_PARTITION_ONLY,
   MMAR_STATUS_ERROR
  };

//+------------------------------------------------------------------+
//| Fitted model parameters                                          |
//+------------------------------------------------------------------+
struct MMARParams
  {
   double                  H;
   double                  tau_q[];
   double                  q_values[];
   ENUM_MMAR_DISTRIBUTION  dist_type;
   double                  dist_params[4];
   double                  alpha_0;
   double                  mean_r_squared;
   double                  sample_volatility;
   bool                    is_multifractal;
   int                     multifractal_score;
  };

//+------------------------------------------------------------------+
//| Monte Carlo forecast results                                     |
//+------------------------------------------------------------------+
struct MMARForecast
  {
   double            mean_volatility;
   double            std_volatility;
   double            ci_lower;
   double            ci_upper;
   double            median_volatility;
   int               n_simulations;
   int               forecast_horizon;
  };

#endif
//+------------------------------------------------------------------+

The ENUM_MMAR_DISTRIBUTION enumeration defines the four candidate distributions that will later be fitted to the multifractal spectrum: Normal, Binomial, Poisson, and Gamma. These correspond exactly to the four theoretical models we explored in the Python implementation. The ENUM_MMAR_FBM_METHOD gives us two options for generating Fractional Brownian Motion: the fast FFT-based Davies-Harte method and the exact but slower Cholesky decomposition. The ENUM_MMAR_STATUS tracks how far through the pipeline a model instance has progressed, which is important because the analysis can partially succeed (partition analysis works but spectrum fitting fails).

The MMARParams structure is the container that will eventually hold everything we extract from the data: the Hurst exponent, the full tau(q) scaling curve, the chosen distribution type and its parameters, and the multifractality diagnostics. The MMARForecast structure holds the output of Monte Carlo simulations: mean, median, standard deviation, and confidence interval bounds for the predicted volatility. These structures will not be fully populated until later articles, but defining them now ensures all modules share a common interface from the start.

The OLS Foundation

The partition function method requires fitting straight lines through log-log plots. For each moment order q, we compute the partition function S_q(dt) at multiple time scales dt, take logarithms of both, and fit a line. The slope of that line gives us tau(q) + 1. To do this rigorously, with R² values and proper diagnostics, we need a full Ordinary Least Squares regression implementation.

Our OLS class uses an SVD-based pseudo-inverse for maximum numerical stability. Unlike the normal equations (X'X)^(-1)X'y, which can blow up when the design matrix is ill-conditioned, the SVD approach decomposes the matrix into orthogonal components and inverts only the non-zero singular values. This matters because our log-log regressions can span several orders of magnitude.

//+------------------------------------------------------------------+
//|                                                         OLS.mqh  |
//|                MMAR Library - Ordinary Least Squares Regression  |
//+------------------------------------------------------------------+
#property strict

#ifndef MMAR_OLS_MQH
#define MMAR_OLS_MQH

//+------------------------------------------------------------------+
//| OLS class                                                        |
//+------------------------------------------------------------------+
class OLS
  {
private:
   bool              m_const_prepended;
   matrix m_exog,
          m_pinv,
          m_cov_params,
          m_m_error,
          m_norm_cov_params;
   vector m_endog,
          m_weights,
          m_singularvalues,
          m_params,
          m_tvalues,
          m_bse,
          m_const_cols,
          m_resid;
   ulong  m_obs,
          m_model_dof,
          m_resid_dof,
          m_kconstant,
          m_rank;
   double m_aic,
          m_bic,
          m_scale,
          m_llf,
          m_sse,
          m_rsqe,
          m_centeredtss,
          m_uncenteredtss;
   uint              m_error;

   //--- Private computation methods
   ulong             countconstants(void);
   void              scale(void);
   void              sse(void);
   void              rsqe(void);
   void              centeredtss(void);
   void              uncenteredtss(void);
   void              aic(void);
   void              bic(void);
   void              bse(void);
   void              llf(void);
   void              tvalues(void);
   void              covariance_matrix(void);

public:
                   OLS(void);
                  ~OLS(void);

   //--- Fitting and prediction
   bool              Fit(vector &y_vars,matrix &x_vars);
   double            Predict(vector &inputs);
   double            Predict(double _input);

   //--- Diagnostics
   ulong             ModelDOF(void);
   ulong             ResidDOF(void);
   double            Scale(void);
   double            Aic(void);
   double            Bic(void);
   double            Sse(void);
   double            Rsqe(void);
   double            Loglikelihood(void);
   vector            Tvalues(void);
   vector            Residuals(void);
   vector            ModelParameters(void);
   vector            Bse(void);
   matrix            CovarianceMatrix(void);
  };

The class interface follows the pattern you would find in Python's statsmodels: call Fit() with a dependent variable vector and an independent variable matrix, then access diagnostics through getter methods. The design matrix must include a constant column if you want an intercept, which is standard practice in econometrics. Let us look at the core fitting method:

//+------------------------------------------------------------------+
//| Fit - Perform OLS regression via SVD pseudo-inverse              |
//+------------------------------------------------------------------+
bool OLS::Fit(vector &y_vars,matrix &x_vars)
  {
   m_endog = y_vars;
   m_exog  = x_vars;
   m_kconstant = countconstants();
   m_error = 0;
   m_model_dof = m_resid_dof = m_rank = m_obs= 0;
   m_aic = m_bic= m_llf=0;
   m_obs=m_exog.Rows();
   m_weights.Resize(m_obs);
   m_weights.Fill(1.0);
   m_rank = m_exog.Rank();
   m_model_dof = m_rank - m_kconstant;
   m_resid_dof = m_obs - m_rank;

   if(y_vars.Size()!=x_vars.Rows())
     {
      Print(__FUNCTION__," ",__LINE__," Error Invalid inputs Rows in x_vars not equal to length of y_vars");
      m_error = 1;
      return false;
     }

   matrix U,V;
   ::ResetLastError();

   if(!m_exog.SVD(U,V,m_singularvalues))
     {
      Print(__FUNCTION__," ",__LINE__," SVD operation failed : ",::GetLastError());
      m_error = 1;
      return false;
     }

   m_pinv = m_exog.PInv();
   matrix pinv_t = m_pinv.Transpose();
   m_norm_cov_params = m_pinv.MatMul(pinv_t);

   matrix diag;
   if(!diag.Diag(m_singularvalues))
     {
      Print(__FUNCTION__," ",__LINE__," Diag operation failed : ",::GetLastError());
      m_error = 1;
      return false;
     }

   m_rank = diag.Rank();
   m_params = m_pinv.MatMul(m_endog);
   m_model_dof = m_rank - m_kconstant;
   m_resid_dof = m_obs - m_rank;
   m_resid = m_endog - m_exog.MatMul(m_params);
   scale();
   covariance_matrix();
   bse();
   sse();
   llf();
   centeredtss();
   uncenteredtss();
   rsqe();
   aic();
   bic();
   tvalues();

   return true;
  }

The method decomposes the design matrix via SVD and computes the Moore–Penrose pseudo-inverse. The parameter estimates are then simply the pseudo-inverse multiplied by the dependent variable: params = pinv(X) * y. It then estimates coefficients and computes diagnostics (scale, covariance, BSE, SSE, log-likelihood, TSS, R², AIC/BIC, and t-values). For our purposes in partition analysis, R² is the critical diagnostic. It tells us how well a straight line describes the log-log relationship at each moment order q, which is the fundamental assumption of scaling behavior.

The R² calculation itself follows the standard formula:

//+------------------------------------------------------------------+
//| R-squared                                                        |
//+------------------------------------------------------------------+
void OLS::rsqe(void)
  {
   m_rsqe = (m_kconstant)? 1 - m_sse/m_centeredtss: 1 - m_sse/m_uncenteredtss;
  }

When a constant column is present in the design matrix (which we always include), R² is computed as 1 minus the ratio of residual sum of squares to the centered total sum of squares. This gives the proportion of variance in the dependent variable that is explained by the regression. Values close to 1.0 indicate strong linear scaling, while values below our threshold (0.60 by default) suggest the power-law assumption does not hold at that particular q value.

Estimating the Hurst Exponent

The Hurst exponent H is the single most important parameter in the MMAR framework. It controls the long-memory structure of the Fractional Brownian Motion component. An H of 0.5 means no memory (standard Brownian motion), H greater than 0.5 means persistence (trends continue), and H less than 0.5 means anti-persistence (trends reverse). Our library provides two estimation methods: a primary method based on the tau(q) zero-crossing, and a fallback based on the Generalized Hurst Exponent (GHE).

The GHE method lives in its own header file and estimates H(q) using qth-order structure functions. It works by computing lagged differences of the price series at multiple lag scales, fitting a log-log relationship between the structure function and the lag, and extracting the slope. The implementation requires several helper functions that handle vector initialization and lagged differencing:

#ifndef MMAR_HURST_GHE_MQH
#define MMAR_HURST_GHE_MQH

//+------------------------------------------------------------------+
//| Vector arange initialization                                     |
//+------------------------------------------------------------------+
template<typename T>
void mmar_arange(vector<T> &vec, T value = 0.0, T step = 1.0)
  {
   for(ulong i = 0; i < vec.Size(); i++, value += step)
      vec[i] = value;
  }

//+------------------------------------------------------------------+
//| Lagged difference helper                                         |
//+------------------------------------------------------------------+
bool mmar_diff_array(uint st, double &in[], vector &out, vector &other)
  {
   ulong nsize = (MathMod(in.Size(), st) == 0) ? int(in.Size() / st) - 1 : int(in.Size() / st);

   if(!out.Resize(nsize) || !other.Resize(nsize + 1))
     {
      Print("resize error ", GetLastError());
      return false;
     }

   uint i, tt = 0;
   for(i = st; i < in.Size(); i += st, tt++)
     {
      other[tt] = in[i - st];
      out[tt] = in[i] - in[i - st];
     }

   other[tt] = in[i - st];

   return true;
  }

The mmar_arange function is a simple utility that fills a vector with evenly spaced values, analogous to numpy's arange. The mmar_diff_array function computes lagged differences at a given step size: for a step of j, it produces the series X[j]-X[0], X[2j]-X[j], X[3j]-X[2j], and so on. This is the raw material for the structure function computation at each lag scale. The file also includes a mmar_vecToArray helper that converts a vector to a double array, along with a vector-input overload of general_hurst that uses it. This lets other modules call the GHE estimator with either arrays or vectors depending on what is convenient.

The main GHE function loops over lag scales from a lower bound to an upper bound, computing at each scale the ratio of the qth moment of detrended differences to the qth moment of detrended levels. The log-log slope of this quantity across scales gives an estimate of H(q)/q:

//+------------------------------------------------------------------+
//| Generalized Hurst exponent H(q) - array input                    |
//+------------------------------------------------------------------+
double general_hurst(double &data[], int q, int lower = 2, int upper = 50)
  {
   if(data.Size() < 100)
     {
      Print("data array is of insufficient length");
      return EMPTY_VALUE;
     }

   if(lower >= upper || lower < 2 || upper > int(floor(0.5 * data.Size())))
     {
      Print("Invalid input for lower and/or upper");
      return EMPTY_VALUE;
     }

   uint len = data.Size();

   int k = 0;

   matrix H;
   if(!H.Resize(ulong(upper - lower), 1))
     {
      Print(__LINE__, " ", __FUNCTION__, " ", GetLastError());
      return EMPTY_VALUE;
     }

//--- Main loop over lag scales
   for(int i = lower; i < upper; i++)
     {
      vector x_vector(ulong(i), mmar_arange<double>, 1.0, 1.0);

      matrix mcord;
      if(!mcord.Resize(ulong(i), 1))
        {
         Print(__LINE__, " ", __FUNCTION__, " ", GetLastError());
         return EMPTY_VALUE;
        }

      mcord.Fill(0.0);

      for(int j = 1; j < i + 1; j++)
        {
         vector dv, Y;
         if(!mmar_diff_array(j, data, dv, Y))
            return EMPTY_VALUE;

         double N = double(Y.Size());

         vector X(ulong(N), mmar_arange<double>, 1.0, 1.0);

         double mx = X.Sum() / N;
         vector XP = MathPow(X, 2.0);
         double SSxx = XP.Sum() - N * pow(mx, 2.0);
         double my = Y.Sum() / N;
         vector XY = X * Y;
         double SSxy = XY.Sum() - N * mx * my;
         double cc1 = SSxy / SSxx;
         double cc2 = my - cc1 * mx;

         vector ddVd = dv - cc1;
         vector VVVd = Y - cc1 * X - cc2;
         vector PddVd = MathAbs(ddVd);
         PddVd = pow(PddVd, q);
         vector PVVVd = MathAbs(VVVd);
         PVVVd = pow(PVVVd, q);

         mcord[j - 1][0] = PddVd.Mean() / PVVVd.Mean();
        }

      vector Px_vector = MathLog10(x_vector);

      double mx = Px_vector.Mean();
      vector Sqx = MathPow(Px_vector, 2.0);
      double SSxx = Sqx.Sum() - i * pow(mx, 2.0);
      matrix lmcord = log10(mcord);
      double my = lmcord.Mean();
      vector pt = Px_vector * lmcord.Col(0);
      double SSxy = pt.Sum() - i * mx * my;

      H[k][0] = SSxy / SSxx;

      k++;
     }

   return H.Mean() / double(q);
  }

The algorithm works in three nested stages. The outer loop iterates over lag window sizes from lower to upper. For each window size i, the inner loop computes the structure function at lags 1 through i. At each lag j, the data is differenced at that lag and detrended using a simple linear regression (removing any drift). The ratio of the qth power of the detrended differences to the qth power of the detrended levels is stored. After the inner loop completes, the log-log slope of this ratio against the lag index gives one estimate of H(q). The final result is the average across all window sizes, divided by q.

This estimator serves as a fallback. The primary method for estimating H in our library is the tau(q) zero-crossing: we find the moment order q* where tau(q*) = 0, and then H = 1/q*. This is theoretically cleaner because it uses the partition function we are already computing. The GHE fallback activates only when the zero-crossing produces an implausible result (outside the range 0.1 to 0.9), which can happen if the tau(q) curve does not cross zero cleanly due to noise or insufficient data.

Partition Analysis — The Core Engine

With the utilities in place, we arrive at the central module of this article: CPartitionAnalysis. This class encapsulates the complete multifractal partition function analysis pipeline. It takes an array of log-returns, computes the partition function S_q(dt) across a grid of moment orders q and time scales dt, fits log-log regression lines to extract the scaling exponents tau(q), estimates the Hurst exponent, and evaluates whether the data is genuinely multifractal through three independent diagnostic tests.

Class Structure and Configuration

The class exposes a clean two-step interface: initialize with data, then analyze. All configuration parameters have sensible defaults but can be overridden before calling Analyze():

//+------------------------------------------------------------------+
//| CPartitionAnalysis class                                         |
//+------------------------------------------------------------------+
class CPartitionAnalysis
  {
private:
   //--- Input data
   double            m_returns[];
   int               m_n_returns;

   //--- Configuration
   int               m_delta_t_min;
   int               m_delta_t_max;
   double            m_delta_t_spacing;
   double            m_q_min;
   double            m_q_max;
   double            m_q_step;
   double            m_min_r_squared;

   //--- Results
   double            m_q_values[];
   int               m_n_q;
   double            m_tau_q[];
   double            m_slopes[];
   double            m_r_squared[];
   double            m_H;
   double            m_sample_volatility;
   int               m_multifractal_score;
   bool              m_is_multifractal;
   bool              m_analyzed;

   //--- Internal computation
   double            CalculatePartitionFunction(double q, int delta_t);
   void              GenerateDeltaTValues(int &dt_values[], int &n_dt);
   void              GenerateQValues(void);
   bool              FitPartitionPlots(const int &dt_values[], int n_dt);
   double            EstimateHurst(void);
   double            EstimateHurstFromTauQ(void);
   int               CheckMultifractality(void);
   double            ComputeSampleVolatility(void);

public:
                   CPartitionAnalysis(void);
                  ~CPartitionAnalysis(void);

   //--- Configuration
   void              SetPartitionConfig(int dt_min, int dt_max, double spacing);
   void              SetMomentConfig(double q_min, double q_max, double q_step);
   void              SetMinRSquared(double min_r2) { m_min_r_squared = min_r2; }

   //--- Main pipeline
   bool              Init(const double &returns[], int n_returns);
   bool              Analyze(void);

   //--- Results access
   double            GetH(void)                 { return m_H; }
   double            GetSampleVolatility(void)  { return m_sample_volatility; }
   void              GetTauQ(double &tau[], double &q_vals[]);
   double            GetMeanRSquared(void);
   bool              IsMultifractal(void)       { return m_is_multifractal; }
   int               GetMultifractalScore(void) { return m_multifractal_score; }
   int               GetNumQ(void)              { return m_n_q; }
  };

The configuration splits naturally into two groups. The partition config controls the time scale grid: dt_min and dt_max define the range, and spacing controls how densely we sample between them (logarithmically). The moment config controls the q-value grid: we sweep from q_min to q_max in fixed steps. The default range of q from 0.01 to 30.0 in steps of 0.5 gives us 61 moment orders, each requiring its own log-log regression. This grid must be fine enough to capture the curvature of tau(q) without being so dense that computation becomes prohibitive.

Generating Time Scales

The time scales are spaced logarithmically rather than linearly. This is critical because scaling behavior manifests across orders of magnitude, not across fixed increments. A linear spacing from 1 to 150 would give 150 points, most of them clustered at small scales. Logarithmic spacing gives roughly equal representation to each decade of scale:

//+------------------------------------------------------------------+
//| GenerateDeltaTValues - Build log-spaced time scale array         |
//+------------------------------------------------------------------+
void CPartitionAnalysis::GenerateDeltaTValues(int &dt_values[], int &n_dt)
  {
   int temp[];
   ArrayResize(temp, 200);
   int count = 0;

   double current = (double)m_delta_t_min;
   while(current <= (double)m_delta_t_max)
     {
      int dt = (int)MathRound(current);
      if(dt < 1)
         dt = 1;
      if(dt > m_n_returns / 2)
         break;

      bool duplicate = false;
      for(int i = 0; i < count; i++)
        {
         if(temp[i] == dt)
           {
            duplicate = true;
            break;
           }
        }

      if(!duplicate)
        {
         if(count >= ArraySize(temp))
            ArrayResize(temp, count + 100);
         temp[count++] = dt;
        }

      current *= m_delta_t_spacing;
     }

   for(int i = 0; i < count - 1; i++)
      for(int j = i + 1; j < count; j++)
         if(temp[i] > temp[j])
           {
            int swap = temp[i];
            temp[i] = temp[j];
            temp[j] = swap;
           }

   n_dt = count;
   ArrayResize(dt_values, count);
   ArrayCopy(dt_values, temp, 0, 0, count);
  }

Starting from dt_min (1), we multiply by the spacing factor (1.1) at each step. This produces a geometric sequence: 1, 1.1, 1.2, 1.3, 1.5, 1.6, 1.8, 2, 2.2, ... up to dt_max. Because these are time scales in integer bar counts, we round each value and eliminate duplicates. The guard condition dt > m_n_returns / 2 prevents us from creating time windows larger than half the data, which would give too few non-overlapping blocks for a meaningful partition function estimate. The resulting array is sorted ascending for the log-log regression to make physical sense.

Computing the Partition Function

The heart of multifractal analysis is the partition function S_q(dt). For a given moment order q and time scale dt, it answers the question: what is the average qth power of absolute returns measured over non-overlapping windows of size dt?

//+------------------------------------------------------------------+
//| CalculatePartitionFunction - S_q(dt) for given q and time scale  |
//+------------------------------------------------------------------+
double CPartitionAnalysis::CalculatePartitionFunction(double q, int delta_t)
  {
   int N = m_n_returns / delta_t;
   if(N <= 0)
      return -1.0;

   double sum_sq = 0.0;
   for(int i = 0; i < N; i++)
     {
      double interval_return = 0.0;
      int start = i * delta_t;
      for(int j = 0; j < delta_t; j++)
         interval_return += m_returns[start + j];

      double abs_r = MathAbs(interval_return);
      if(abs_r < 1e-300)
         abs_r = 1e-300;
      sum_sq += MathPow(abs_r, q);
     }

   return sum_sq / N;
  }

The function divides the return series into N non-overlapping blocks of size delta_t. For each block, it sums the individual log-returns within that block to get the block return, takes the absolute value, raises it to the power q, and accumulates the result. The final output is the average across all blocks. The floor value of 1e-300 prevents taking the logarithm of zero in the subsequent log-log fitting step. This is a defensive measure for quiet market periods where multiple consecutive returns might be exactly zero (especially on higher timeframes or illiquid instruments).

The key insight is that if the data exhibits scaling behavior, then S_q(dt) follows a power law: S_q(dt) is proportional to dt^(tau(q)+1). Taking logarithms linearizes this relationship: log(S_q) = (tau(q)+1) * log(dt) + constant. The slope of this log-log plot gives us tau(q) + 1.

Figure 2 animates this process: returns are grouped into blocks, the partition function is computed at multiple scales, and the log-log slopes reveal different scaling rates for different moment orders q:

$Log-returns grouped by scale to show partition function and multifractal scaling$

Fig. 2. Partition function scaling: different slopes for different q values reveal multifractal behavior

Fitting the Log-Log Plots

For each moment order q in our grid, we compute S_q(dt) at every time scale, take base-10 logarithms of both quantities, and run an OLS regression:

//+------------------------------------------------------------------+
//| FitPartitionPlots - Log-log regression for each q value          |
//+------------------------------------------------------------------+
bool CPartitionAnalysis::FitPartitionPlots(const int &dt_values[], int n_dt)
  {
   ArrayResize(m_tau_q, m_n_q);
   ArrayResize(m_slopes, m_n_q);
   ArrayResize(m_r_squared, m_n_q);

   int valid_q = 0;

   for(int qi = 0; qi < m_n_q; qi++)
     {
      double q = m_q_values[qi];

      double log_dt[];
      double log_sq[];
      ArrayResize(log_dt, n_dt);
      ArrayResize(log_sq, n_dt);
      int n_valid = 0;

      for(int di = 0; di < n_dt; di++)
        {
         double s_q = CalculatePartitionFunction(q, dt_values[di]);
         if(s_q > 0.0)
           {
            log_dt[n_valid] = MathLog10((double)dt_values[di]);
            log_sq[n_valid] = MathLog10(s_q);
            n_valid++;
           }
        }

      if(n_valid < 3)
        {
         m_slopes[qi] = 0.0;
         m_tau_q[qi] = -1.0;
         m_r_squared[qi] = 0.0;
         continue;
        }

      vector y_vec;
      y_vec.Resize(n_valid);
      matrix x_mat;
      x_mat.Resize(n_valid, 2);

      for(int i = 0; i < n_valid; i++)
        {
         y_vec[i] = log_sq[i];
         x_mat[i][0] = log_dt[i];
         x_mat[i][1] = 1.0;
        }

      OLS ols;
      if(!ols.Fit(y_vec, x_mat))
        {
         m_slopes[qi] = 0.0;
         m_tau_q[qi] = -1.0;
         m_r_squared[qi] = 0.0;
         continue;
        }

      vector params = ols.ModelParameters();
      m_slopes[qi] = params[0];
      m_tau_q[qi] = params[0] - 1.0;
      m_r_squared[qi] = ols.Rsqe();
      valid_q++;
     }

   if(valid_q < 3)
     {
      Print("PartitionAnalysis::FitPartitionPlots - only ", valid_q, " valid q values");
      return false;
     }

   return true;
  }

For each q value, we build the design matrix with two columns: log10(dt) as the predictor and a column of ones as the intercept. The dependent variable is log10(S_q). After fitting, params[0] gives us the slope, which equals tau(q) + 1. We subtract 1 to get tau(q) directly. The R² value tells us how well the power-law model fits at this particular q. If fewer than 3 valid data points exist (which can happen at extreme q values where numerical overflow or underflow occurs), we skip that q value entirely rather than fitting a meaningless regression.

Estimating the Hurst Exponent from tau(q)

The relationship between the scaling function and the Hurst exponent is elegant: tau(q) crosses zero at q* = 1/H. This means we can estimate H by finding where the tau(q) curve crosses zero and inverting that value. The implementation uses linear interpolation between the two q values that straddle the zero crossing:

//+------------------------------------------------------------------+
//| EstimateHurstFromTauQ - Find q* where tau(q*)=0, then H=1/q*     |
//+------------------------------------------------------------------+
double CPartitionAnalysis::EstimateHurstFromTauQ(void)
  {
   int idx_below = -1;
   int idx_above = -1;

   for(int i = 0; i < m_n_q - 1; i++)
     {
      if(m_tau_q[i] * m_tau_q[i + 1] <= 0.0)
        {
         idx_below = i;
         idx_above = i + 1;
         break;
        }
     }

   if(idx_below < 0)
     {
      double min_abs = DBL_MAX;
      int min_idx = 0;
      for(int i = 0; i < m_n_q; i++)
        {
         if(MathAbs(m_tau_q[i]) < min_abs)
           {
            min_abs = MathAbs(m_tau_q[i]);
            min_idx = i;
           }
        }
      double q_zero = m_q_values[min_idx];
      if(q_zero <= 0.0)
         return 0.5;
      return MathMax(0.1, MathMin(0.9, 1.0 / q_zero));
     }

   double q1 = m_q_values[idx_below];
   double q2 = m_q_values[idx_above];
   double t1 = m_tau_q[idx_below];
   double t2 = m_tau_q[idx_above];

   double q_zero;
   if(MathAbs(t2 - t1) > 1e-15)
      q_zero = q1 + (0.0 - t1) * (q2 - q1) / (t2 - t1);
   else
      q_zero = (q1 + q2) / 2.0;

   if(q_zero <= 0.0)
      return 0.5;
   return MathMax(0.1, MathMin(0.9, 1.0 / q_zero));
  }

The sign-change detection (tau[i] * tau[i+1] <= 0) finds adjacent q values where tau transitions from negative to positive (or vice versa). Linear interpolation between these two points gives us the precise q* where tau(q*) = 0. Then H = 1/q*. The result is clamped to [0.1, 0.9] because values outside this range are physically implausible for financial data and likely indicate a numerical artifact. If no zero crossing exists at all (which happens when the data has very weak scaling), the fallback finds the q value closest to zero and uses that approximation.

The parent EstimateHurst() method tries the tau(q) method first, and if it returns an implausible value, falls back to the GHE estimator we discussed in the previous section.

The Multifractality Decision

Extracting tau(q) is only half the job. We also need to determine whether the scaling is genuinely multifractal (the curve is nonlinear) or merely monofractal (the curve is a straight line, meaning all moments scale identically). Our library uses three independent diagnostic tests, each scored 0-3, for a maximum total of 9:

//+------------------------------------------------------------------+
//| CheckMultifractality - Score 0-9 based on three diagnostics      |
//+------------------------------------------------------------------+
int CPartitionAnalysis::CheckMultifractality(void)
  {
   if(m_n_q < 5)
      return 0;

//--- Test 1: Concavity of tau(q)
   int n_concave = 0;
   int n_d2 = m_n_q - 2;
   for(int i = 0; i < n_d2; i++)
     {
      double d2 = m_tau_q[i + 2] - 2.0 * m_tau_q[i + 1] + m_tau_q[i];
      if(d2 < 0.0)
         n_concave++;
     }
   double concave_pct = 100.0 * n_concave / n_d2;

   int concavity_score;
   if(concave_pct >= 80.0)
      concavity_score = 3;
   else
      if(concave_pct >= 60.0)
         concavity_score = 2;
      else
         if(concave_pct >= 40.0)
            concavity_score = 1;
         else
            concavity_score = 0;

//--- Test 2: Linearity (high R^2 → monofractal)
   double sum_q = 0.0, sum_t = 0.0, sum_qq = 0.0, sum_qt = 0.0;
   int n_pos = 0;
   for(int i = 0; i < m_n_q; i++)
     {
      if(m_q_values[i] > 0.0)
        {
         sum_q += m_q_values[i];
         sum_t += m_tau_q[i];
         sum_qq += m_q_values[i] * m_q_values[i];
         sum_qt += m_q_values[i] * m_tau_q[i];
         n_pos++;
        }
     }

   double mean_q = sum_q / n_pos;
   double mean_t = sum_t / n_pos;
   double ss_qq = sum_qq - n_pos * mean_q * mean_q;
   double ss_qt = sum_qt - n_pos * mean_q * mean_t;
   double slope_lin = ss_qt / ss_qq;
   double intercept_lin = mean_t - slope_lin * mean_q;

   double ss_res = 0.0, ss_tot = 0.0;
   for(int i = 0; i < m_n_q; i++)
     {
      if(m_q_values[i] > 0.0)
        {
         double fitted = slope_lin * m_q_values[i] + intercept_lin;
         ss_res += (m_tau_q[i] - fitted) * (m_tau_q[i] - fitted);
         ss_tot += (m_tau_q[i] - mean_t) * (m_tau_q[i] - mean_t);
        }
     }
   double r2_linear = (ss_tot > 0.0) ? 1.0 - ss_res / ss_tot : 1.0;
   double rmse_lin = MathSqrt(ss_res / n_pos);

   int linearity_score;
   if(r2_linear > 0.99 && rmse_lin < 0.05)
      linearity_score = 0;
   else
      if(r2_linear > 0.95)
         linearity_score = 1;
      else
         if(r2_linear > 0.90)
            linearity_score = 2;
         else
            linearity_score = 3;

//--- Test 3: Variation in dtau/dq slopes
   double d_tau_min = DBL_MAX, d_tau_max = -DBL_MAX;
   for(int i = 0; i < m_n_q - 1; i++)
     {
      double d = m_tau_q[i + 1] - m_tau_q[i];
      if(d < d_tau_min)
         d_tau_min = d;
      if(d > d_tau_max)
         d_tau_max = d;
     }
   double slope_range = d_tau_max - d_tau_min;

   int variation_score;
   if(slope_range > 0.5)
      variation_score = 3;
   else
      if(slope_range > 0.3)
         variation_score = 2;
      else
         if(slope_range > 0.1)
            variation_score = 1;
         else
            variation_score = 0;

   int total = concavity_score + linearity_score + variation_score;

   PrintFormat("Multifractality check: concavity=%d linearity=%d variation=%d total=%d/9",
               concavity_score, linearity_score, variation_score, total);

   return total;
  }

Test 1 — Concavity: A multifractal tau(q) curve must be concave (curving downward). We compute the second finite difference at every interior point. If d2 = tau[i+2] - 2*tau[i+1] + tau[i] is negative, that point is concave. A score of 3 requires 80% or more of points to be concave. A perfectly linear tau(q) would have d2 approximately zero everywhere, scoring 0.

Test 2 — Linearity: This is the inverse test. We fit a straight line directly to the tau(q) curve and measure how well it fits. If R² exceeds 0.99 and the RMSE is below 0.05, the curve is essentially linear and the data is monofractal (score 0). The dual condition prevents false negatives: a high R² alone can occur when the tau(q) curve has a slight systematic curvature that happens to be well-correlated with a line, but the RMSE catches cases where the residuals are still meaningfully large. Lower R² means the curve deviates significantly from a line, suggesting genuine multifractal structure (score up to 3).

Test 3 — Slope Variation: For a multifractal process, the local slope dtau/dq varies across the q range. The local slope at small q reflects the behavior of large fluctuations, while the slope at large q reflects the behavior of small fluctuations. If these differ substantially (range > 0.5), it scores 3. A monofractal process has nearly constant slope throughout.

A total score of 5 or higher (out of 9) classifies the data as multifractal. This threshold was calibrated against known multifractal and monofractal processes during development. It is deliberately conservative: we would rather miss some weakly multifractal data than falsely classify random noise as multifractal.

Figure 3 visualizes the distinction between monofractal and multifractal tau(q) curves, and shows how the three diagnostic tests detect the difference:

$Monofractal vs. multifractal τ(q) curves with key diagnostic tests and Hurst exponent$

Fig. 3. Monofractal vs multifractal tau(q): a straight line means monofractal, a concave curve means multifractal

Running the Analysis on Live Data

With all modules in place, we can now write a test script that loads real market data from the MetaTrader terminal and runs the full partition analysis pipeline. This script will be our first end-to-end validation that the MQL5 implementation produces sensible results:

//+------------------------------------------------------------------+
//|                                              MMAR Partition Test |
//|                      Copyright 2026, MMQ - Muhammad Minhas Qamar |
//|                           https://www.mql5.com/en/articles/21299 |
//+------------------------------------------------------------------+
#property copyright "Copyright 2026, MMQ - Muhammad Minhas Qamar"
#property link      "https://www.mql5.com/en/articles/21299"
#property version   "1.00"
#property strict
#property script_show_inputs

#include <MMAR\PartitionAnalysis.mqh>

input int    InpBars       = 50000;
input double InpQMin       = 0.01;
input double InpQMax       = 30.0;
input double InpQStep      = 0.5;
input int    InpDtMin      = 1;
input int    InpDtMax      = 150;
input double InpDtSpacing  = 1.1;
input double InpMinR2      = 0.60;

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   PrintFormat("=== MMAR Partition Analysis Test ===");
   PrintFormat("Symbol: %s | Period: %s", _Symbol, EnumToString(_Period));

//--- Load close prices and compute log returns
   double close[];
   int copied = CopyClose(_Symbol, _Period, 0, InpBars, close);
   if(copied < 1000)
     {
      PrintFormat("Not enough bars: %d (need 1000+)", copied);
      return;
     }
   PrintFormat("Loaded %d bars", copied);

   double returns[];
   int n_returns = copied - 1;
   ArrayResize(returns, n_returns);
   for(int i = 0; i < n_returns; i++)
      returns[i] = MathLog(close[i + 1] / close[i]);

//--- Run partition analysis
   CPartitionAnalysis partition;
   partition.SetPartitionConfig(InpDtMin, InpDtMax, InpDtSpacing);
   partition.SetMomentConfig(InpQMin, InpQMax, InpQStep);
   partition.SetMinRSquared(InpMinR2);

   if(!partition.Init(returns, n_returns))
     {
      Print("Init failed");
      return;
     }

   uint start = GetTickCount();
   if(!partition.Analyze())
     {
      Print("Analyze failed");
      return;
     }
   uint elapsed = GetTickCount() - start;
   PrintFormat("Analysis time: %u ms", elapsed);

//--- Print results
   PrintFormat("\n--- Results ---");
   PrintFormat("H (Hurst exponent): %.4f", partition.GetH());
   PrintFormat("Sample volatility:  %.6f", partition.GetSampleVolatility());
   PrintFormat("Mean R-squared:     %.4f", partition.GetMeanRSquared());
   PrintFormat("Multifractal score: %d/9", partition.GetMultifractalScore());
   PrintFormat("Is multifractal:    %s", partition.IsMultifractal() ? "YES" : "NO");

//--- Print tau(q) curve
   double tau[], q_vals[];
   partition.GetTauQ(tau, q_vals);
   int n_q = partition.GetNumQ();

   PrintFormat("\n--- tau(q) curve (every 5th value) ---");
   PrintFormat("%8s %10s", "q", "tau(q)");
   for(int i = 0; i < n_q; i += 5)
      PrintFormat("%8.2f %10.4f", q_vals[i], tau[i]);

   Print("\n=== Partition Analysis Test Complete ===");
  }

The script follows a straightforward pattern that will recur throughout this series. First, we load historical close prices using CopyClose. We request up to InpBars bars (default 50,000) and require at least 1,000 for meaningful analysis. Next, we compute log-returns: the natural logarithm of consecutive price ratios. These log-returns are the input to all subsequent analysis. Finally, we configure the partition analysis parameters, initialize with the return data, call Analyze(), and print the results.

Let us look at the output when we run this on EURUSD M10 with InpBars set to 500,000 (overriding the default of 50,000 to stress-test the pipeline with a larger sample):

=== MMAR Partition Analysis Test ===
Symbol: EURUSD | Period: PERIOD_M10
Loaded 500000 bars
Multifractality check: concavity=3 linearity=1 variation=1 total=5/9
Analysis time: 3859 ms

--- Results ---
H (Hurst exponent): 0.4892
Sample volatility:  0.000422
Mean R-squared:     0.8867
Multifractal score: 5/9
Is multifractal:    YES

--- tau(q) curve (every 5th value) ---
       q     tau(q)
    0.01    -0.9930
    2.51     0.1724
    5.01     0.7441
    7.51     1.2193
   10.01     1.7181
   12.51     2.2156
   15.01     2.7076
   17.51     3.1952
   20.01     3.6796
   22.51     4.1617
   25.01     4.6422
   27.51     5.1215

=== Partition Analysis Test Complete ===

The results tell a clear story. The Hurst exponent of 0.4892 is very close to 0.5, which means the overall drift structure is near a random walk. This is consistent with the efficient market hypothesis at low frequencies. However, the multifractal score of 5/9 crosses our threshold, meaning the data exhibits genuine multifractal scaling. The concavity test scored 3/3 (the tau(q) curve is strongly concave), confirming that different moment orders scale at different rates. The mean R² of 0.8867 tells us the power-law model fits well across the q range, giving us confidence that the scaling exponents are meaningful.

The tau(q) curve itself reveals the multifractal signature: it starts negative (around -0.99 at q=0.01), crosses zero near q=2 (consistent with H approximately 0.5, since 1/H approximately 2), and grows sublinearly. If the data were monofractal, tau(q) would be a straight line with slope H. The visible curvature, especially the fact that the increments between successive values decrease slightly as q increases, is what the concavity test detects.

Now compare this to a different scenario. Running the same analysis on EURUSD M30 with 50,000 bars produces a score of only 2/9, which does not pass our multifractality threshold. This is not a failure of the code. It is the code correctly identifying that at this particular combination of timeframe and sample size, the multifractal structure is too weak to be statistically meaningful. The partition analysis serves as an honest gatekeeper: if the data does not support multifractal modeling, the library tells you so before you waste computation on cascade generation and Monte Carlo simulation.

This behavior has practical implications. When deploying the MMAR in an Expert Advisor, you will want to select your timeframe and lookback window based on where the data actually exhibits multifractal scaling. The M10 timeframe with a large sample (500,000 bars representing several years) passes the test. Shorter timeframes with fewer bars may not. The library gives you the diagnostic tools to make this decision empirically rather than by guessing.

Conclusion

In this article, we built the first major module of the MMAR library for MetaTrader 5. Starting from shared type definitions, we constructed an OLS regression utility, a Generalized Hurst Exponent estimator, and the full Partition Analysis engine. The result is a self-contained MQL5 class that can take any return series, compute its multifractal scaling properties, and make a rigorous determination about whether multifractal modeling is appropriate for that data.

MMARConstants.mqh provides the shared vocabulary: distribution types, model states, and parameter containers used by all modules.
OLS.mqh delivers SVD-based linear regression with full diagnostics, used internally for log-log slope extraction.
HurstGHE.mqh provides a robust fallback estimator for the Hurst exponent when the primary tau(q) method produces implausible results.
PartitionAnalysis.mqh orchestrates the full pipeline: time scale generation, partition function computation, log-log fitting, Hurst estimation, and the three-test multifractality diagnostic.

The analysis runs in under 4 seconds on 500,000 bars, which is fast enough for periodic recalibration in a live trading environment. More importantly, it produces an honest assessment: not all market data is multifractal, and the library correctly identifies when the assumption does not hold.

In Part 5, we will take the tau(q) curve produced here and fit it to four theoretical distributions using the Legendre transform and bounded optimization. That will give us the distribution parameters needed to generate the multiplicative cascade, the next step toward the complete MMAR simulation engine.

The programs presented in this article are intended for educational purposes only. They are not designed for use in live trading and are not optimized for performance in real-market conditions. The code serves as a foundation for understanding multifractal analysis; always validate thoroughly on demo accounts before considering any adaptation for live trading.

Getting the Source Code via MQL5 Algo Forge

All source files are attached to this article below, but the full repository is also available on MQL5 Algo Forge, the community's Git-based platform for sharing and collaborating on trading projects. Algo Forge stores version history both locally and in the cloud, so you will always have the latest version of the code, including any updates or fixes made after this article was published.

To clone the repository, open a terminal (Command Prompt, PowerShell, or the integrated terminal in your editor) and run:

git clone https://forge.mql5.io/ayantrader/MMAR.git

This requires Git to be installed on your machine. If you are using Visual Studio Code, you can open the integrated terminal with Ctrl+` and run the command directly. The repository will be cloned into an MMAR folder in your current directory.

Alternatively, you can browse the project at forge.mql5.io/ayantrader/MMAR, where you can view the source files, commit history, and any updates directly in your browser. The Algo Forge advantage over the attached files is that you get the full Git history and can pull future updates with a single command:

git pull

File name	Description
MQL5\Include\MMAR\MMARConstants.mqh	Shared enumerations, structures, and type definitions for the MMAR library
MQL5\Include\MMAR\OLS.mqh	Ordinary Least Squares regression via SVD pseudo-inverse with full diagnostics
MQL5\Include\MMAR\HurstGHE.mqh	Generalized Hurst Exponent estimator (fallback method)
MQL5\Include\MMAR\PartitionAnalysis.mqh	Multifractal partition function analysis with scaling diagnostics
MQL5\Scripts\MMAR\MMAR_PartitionTest.mq5	Test script demonstrating partition analysis on live market data

Attached files |

Download ZIP

MQL5.zip (79.09 KB)

Warning: All rights to these materials are reserved by MetaQuotes Ltd. Copying or reprinting of these materials in whole or in part is prohibited.

This article was written by a user of the site and reflects their personal views. MetaQuotes Ltd is not responsible for the accuracy of the information presented, nor for any consequences resulting from the use of the solutions, strategies or recommendations described.