Trading systems

Beyond GARCH (Part IV): Partition Analysis in MQL5 - MQL5 Articles

May 30, 202634 min read

Need a reliable hosting solution for your robots? Contact your broker and find out about available Sponsored MetaTrader VPS offerings Learn more

MetaTrader 5 / Trading systems

Introduction

The first three parts of this series detailed the construction of the entire MMAR (Multifractal Model of Asset Returns) pipeline in Python. We began by loading EURUSD data, subsequently confirming multifractal scaling through partition function analysis, extracting the Hurst exponent, and fitting the singularity spectrum.

Following this, we built a multiplicative cascade and a Fractional Brownian Motion process, performed Monte Carlo simulations, and benchmarked the MMAR against GARCH. The MMAR exhibited superior performance in these comparisons, validating the underlying concept through these tests.

Our current objective presents a significant engineering challenge: Python is excellent for research and rapid prototyping, but it does not execute natively within MetaTrader 5. To enable the MMAR to drive real-time trading decisions, respond to every tick, and integrate seamlessly with the Strategy Tester, a native MQL5 implementation is essential. This requires a self-contained library with no external dependencies, no bridging scripts, and no latency introduced by inter-process communication—just a clean, callable library for any Expert Advisor.

This article marks the beginning of that MQL5 construction. We will develop the first major analytical module of the MMAR library: the Partition Analysis engine.

This component will process raw price data to compute the partition function across various time scales and moment orders. It will then extract the scaling function tau(q), estimate the Hurst exponent H, and perform diagnostic checks to ascertain if the data genuinely exhibits multifractal behavior.

By the end of this article, we will have a functional MQL5 script capable of loading historical bars, processing them through our analysis pipeline, and reporting on the multifractal nature of the market data.

Topics we will cover include:

Library Architecture and Constants
The OLS Foundation
Estimating the Hurst Exponent
Partition Analysis — The Core Engine
Running the Analysis on Live Data
Conclusion

Library Architecture and Constants

Before delving into any analytical code, establishing a shared vocabulary is crucial. The MMAR library is structured as a series of modules, each building upon the other and collectively referencing the same set of enumerations, structures, and type definitions. These foundational elements are centralized in a single header file, MMARConstants.mqh, which is included by every other module.

The architecture employs a layered design. At the foundational layer reside utility modules: an OLS (Ordinary Least Squares) regression class and a Hurst exponent estimator.

Above these, the Partition Analysis engine leverages both utilities to compute the scaling function and assess multifractality. In subsequent articles, we will introduce the Spectrum Fitter (which consumes the tau(q) curve generated by Partition Analysis), the Simulation Engine (responsible for generating synthetic paths based on the fitted distribution), and the Monte Carlo forecaster (which orchestrates multiple simulations).

A facade class, CMMAR, will ultimately unify these components, presenting a simplified Fit-and-Forecast interface. For now, our development proceeds from the bottom up.

Figure 1 illustrates this layered architecture:

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 draw inspiration from the article "Implementing the Generalized Hurst Exponent and the Variance Ratio test in MQL5." That article provided the basis for our OLS implementation (which serves as a self-contained regression container) and the GHE estimation methodology. While we have refined and adapted both modules for our specific application within the MMAR pipeline, credit for establishing these MQL5 implementations rightfully belongs to the original work.

Let us start with the constants file:

mqll

1//+------------------------------------------------------------------+
2//| MMARConstants.mqh |
3//| MMAR Library - Constants, Enumerations, and Data Structures |
4//+------------------------------------------------------------------+
5#property strict
6
7#ifndef MMAR_CONSTANTS_MQH
8#define MMAR_CONSTANTS_MQH
9
10//+------------------------------------------------------------------+
11//| Cascade multiplier distribution types |
12//+------------------------------------------------------------------+
13enum ENUM_MMAR_DISTRIBUTION
14 {
15 MMAR_DIST_NORMAL = 0,
16 MMAR_DIST_BINOMIAL,
17 MMAR_DIST_POISSON,
18 MMAR_DIST_GAMMA
19 };
20
21//+------------------------------------------------------------------+
22//| Fractional Brownian Motion synthesis methods |
23//+------------------------------------------------------------------+
24enum ENUM_MMAR_FBM_METHOD
25 {
26 MMAR_FBM_DAVIES_HARTE = 0,
27 MMAR_FBM_CHOLESKY
28 };
29
30//+------------------------------------------------------------------+
31//| Model fitting status |
32//+------------------------------------------------------------------+
33enum ENUM_MMAR_STATUS
34 {
35 MMAR_STATUS_NOT_FITTED = 0,
36 MMAR_STATUS_FITTED,
37 MMAR_STATUS_PARTITION_ONLY,
38 MMAR_STATUS_ERROR
39 };
40
41//+------------------------------------------------------------------+
42//| Fitted model parameters |
43//+------------------------------------------------------------------+
44struct MMARParams
45 {
46 double H;
47 double tau_q[];
48 double q_values[];
49 ENUM_MMAR_DISTRIBUTION dist_type;
50 double dist_params[4];
51 double alpha_0;
52 double mean_r_squared;
53 double sample_volatility;
54 bool is_multifractal;
55 int multifractal_score;
56 };
57
58//+------------------------------------------------------------------+
59//| Monte Carlo forecast results |
60//+------------------------------------------------------------------+
61struct MMARForecast
62 {
63 double mean_volatility;
64 double std_volatility;
65 double ci_lower;
66 double ci_upper;
67 double median_volatility;
68 int n_simulations;
69 int forecast_horizon;
70 };
71
72#endif
73//+------------------------------------------------------------------+

The ENUM_MMAR_DISTRIBUTION enumeration outlines the four candidate distributions—Normal, Binomial, Poisson, and Gamma—that will later be fitted to the multifractal spectrum. These directly correspond to the theoretical models explored in the Python implementation.

ENUM_MMAR_FBM_METHOD provides two options for generating Fractional Brownian Motion: the efficient FFT-based Davies-Harte method and the precise but slower Cholesky decomposition. , partition analysis completes, but spectrum fitting fails).

The MMARParams structure serves as a container for all extracted data: the Hurst exponent, the complete tau(q) scaling curve, the selected distribution type and its parameters, and multifractality diagnostics. The MMARForecast structure stores the output from Monte Carlo simulations: the mean, median, standard deviation, and confidence interval bounds for predicted volatility. While these structures will not be fully populated until later articles, their early definition ensures a consistent interface across all modules from the outset.

The OLS Foundation

The partition function method necessitates fitting straight lines to log-log plots. For each moment order q, we compute the partition function Sq(dt) at various time scales dt.

We then take the logarithm of both values and perform a linear fit. The slope of this line yields tau(q) + 1.

To achieve this robustly, complete with R² values and comprehensive diagnostics, a full Ordinary Least Squares regression implementation is required.

Our OLS class utilizes an SVD (Singular Value Decomposition)-based pseudo-inverse for maximal numerical stability. Unlike methods relying on the normal equations (X'X)-1X'y, which can become unstable 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 is particularly important for our log-log regressions, which often span several orders of magnitude.

mqll

1//+------------------------------------------------------------------+
2//| OLS.mqh |
3//| MMAR Library - Ordinary Least Squares Regression |
4//+------------------------------------------------------------------+
5#property strict
6
7#ifndef MMAR_OLS_MQH
8#define MMAR_OLS_MQH
9
10//+------------------------------------------------------------------+
11//| OLS class |
12//+------------------------------------------------------------------+
13class OLS
14 {
15private:
16 bool m_const_prepended;
17 matrix m_exog,
18 m_pinv,
19 m_cov_params,
20 m_m_error,
21 m_norm_cov_params;
22 vector m_endog,
23 m_weights,
24 m_singularvalues,
25 m_params,
26 m_tvalues,
27 m_bse,
28 m_const_cols,
29 m_resid;
30 ulong m_obs,
31 m_model_dof,
32 m_resid_dof,
33 m_kconstant,
34 m_rank;
35 double m_aic,
36 m_bic,
37 m_scale,
38 m_llf,
39 m_sse,
40 m_rsqe,
41 m_centeredtss,
42 m_uncenteredtss;
43 uint m_error;
44
45 //--- Private computation methods
46 ulong countconstants(void);
47 void scale(void);
48 void sse(void);
49 void rsqe(void);
50 void centeredtss(void);
51 void uncenteredtss(void);
52 void aic(void);
53 void bic(void);
54 void bse(void);
55 void llf(void);
56 void tvalues(void);
57 void covariance_matrix(void);
58
59public:
60 OLS(void);
61 ~OLS(void);
62
63 //--- Fitting and prediction
64 bool Fit(vector &y_vars,matrix &x_vars);
65 double Predict(vector &inputs);
66 double Predict(double _input);
67
68 //--- Diagnostics
69 ulong ModelDOF(void);
70 ulong ResidDOF(void);
71 double Scale(void);
72 double Aic(void);
73 double Bic(void);
74 double Sse(void);
75 double Rsqe(void);
76 double Loglikelihood(void);
77 vector Tvalues(void);
78 vector Residuals(void);
79 vector ModelParameters(void);
80 vector Bse(void);
81 matrix CovarianceMatrix(void);
82 };

The class interface mirrors the design found in Python's statsmodels: users 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 to fit an intercept, which is standard practice in econometrics. Let us examine the core fitting method:

mqll

1//+------------------------------------------------------------------+
2//| Fit - Perform OLS regression via SVD pseudo-inverse |
3//+------------------------------------------------------------------+
4bool OLS::Fit(vector &y_vars,matrix &x_vars)
5 {
6 m_endog = y_vars;
7 m_exog = x_vars;
8 m_kconstant = countconstants;
9 m_error = 0;
10 m_model_dof = m_resid_dof = m_rank = m_obs= 0;
11 m_aic = m_bic= m_llf=0;
12 m_obs=m_exog.Rows;
13 m_weights.Resize(m_obs);
14 m_weights.Fill(1.0);
15 m_rank = m_exog.Rank;
16 m_model_dof = m_rank - m_kconstant;
17 m_resid_dof = m_obs - m_rank;
18
19 if(y_vars.Size!=x_vars.Rows)
20 {
21 Print(__FUNCTION__," ",__LINE__," Error Invalid inputs Rows in x_vars not equal to length of y_vars");
22 m_error = 1;
23 return false;
24 }
25
26 matrix U,V;
27 ::ResetLastError;
28
29 if(!m_exog.SVD(U,V,m_singularvalues))
30 {
31 Print(__FUNCTION__," ",__LINE__," SVD operation failed : ",::GetLastError);
32 m_error = 1;
33 return false;
34 }
35
36 m_pinv = m_exog.PInv;
37 matrix pinv_t = m_pinv.Transpose;
38 m_norm_cov_params = m_pinv.MatMul(pinv_t);
39
40 matrix diag;
41 if(!diag.Diag(m_singularvalues))
42 {
43 Print(__FUNCTION__," ",__LINE__," Diag operation failed : ",::GetLastError);
44 m_error = 1;
45 return false;
46 }
47
48 m_rank = diag.Rank;
49 m_params = m_pinv.MatMul(m_endog);
50 m_model_dof = m_rank - m_kconstant;
51 m_resid_dof = m_obs - m_rank;
52 m_resid = m_endog - m_exog.MatMul(m_params);
53 scale;
54 covariance_matrix;
55 bse;
56 sse;
57 llf;
58 centeredtss;
59 uncenteredtss;
60 rsqe;
61 aic;
62 bic;
63 tvalues;
64
65 return true;
66 }

The method performs singular value decomposition on the design matrix and computes the Moore–Penrose pseudo-inverse. The parameter estimates are then derived by multiplying the pseudo-inverse by the dependent variable: params = pinv(X) * y.

Subsequently, it computes various diagnostics including coefficients, scale, covariance, Bounded Standard Error (BSE), Sum of Squared Errors (SSE), log-likelihood, Total Sum of Squares (TSS), R², Akaike Information Criterion (AIC)/Bayesian Information Criterion (BIC), and t-values. For our partition analysis, the R² value is particularly crucial as it indicates how well a straight line represents the log-log relationship at each moment order q, which is a fundamental assumption for scaling behavior.

The R² calculation itself adheres to the standard formula:

mqll

1//+------------------------------------------------------------------+
2//| R-squared |
3//+------------------------------------------------------------------+
4void OLS::rsqe(void)
5 {
6 m_rsqe = (m_kconstant)? 1 - m_sse/m_centeredtss: 1 - m_sse/m_uncenteredtss;
7 }

When a constant column is present in the design matrix (which we always include for our regressions), R² is calculated as 1 minus the ratio of the residual sum of squares to the centered total sum of squares. This metric quantifies the proportion of variance in the dependent variable explained by the regression model. Values approaching 1.0 indicate strong linear scaling, whereas values below our predefined threshold (0.60 by default) suggest that the power-law assumption may not hold for that specific q value.

Estimating the Hurst Exponent

The Hurst exponent (H) is the most critical parameter within the MMAR framework. It governs the long-memory characteristics of the Fractional Brownian Motion component.

5 suggests anti-persistence (trends tend to reverse). Our library offers two estimation methods: a primary method based on the tau(q) zero-crossing and a fallback method utilizing the Generalized Hurst Exponent (GHE).

The GHE method is encapsulated within its own header file and estimates H(q) using q-th order structure functions. It operates by computing lagged differences of the price series across multiple lag scales, fitting a log-log relationship between the structure function and the lag, and then extracting the slope. This implementation requires several helper functions for vector initialization and lagged differencing:

mqll

1#ifndef MMAR_HURST_GHE_MQH
2#define MMAR_HURST_GHE_MQH
3
4//+------------------------------------------------------------------+
5//| Vector arange initialization |
6//+------------------------------------------------------------------+
7template<typename T>
8void mmar_arange(vector<T> &vec, T value = 0.0, T step = 1.0)
9 {
10 for(ulong i = 0; i < vec.Size; i++, value += step)
11 vec[i] = value;
12 }
13
14//+------------------------------------------------------------------+
15//| Lagged difference helper |
16//+------------------------------------------------------------------+
17bool mmar_diff_array(uint st, double &in[], vector &out, vector &other)
18 {
19 ulong nsize = (MathMod(in.Size, st) == 0) ? (ulong)in.Size / st - 1 : (ulong)in.Size / st;
20
21 if(!out.Resize(nsize) || !other.Resize(nsize + 1))
22 {
23 Print("resize error ", GetLastError);
24 return false;
25 }
26
27 uint i, tt = 0;
28 for(i = st; i < in.Size; i += st, tt++)
29 {
30 other[tt] = in[i - st];
31 out[tt] = in[i] - in[i - st];
32 }
33
34 other[tt] = in[i - st];
35
36 return true;
37 }

The mmar_arange function is a straightforward utility that populates a vector with evenly spaced values, akin to NumPy's arange. The mmar_diff_array function calculates lagged differences at a specified step size: for a step j, it generates the series X[j]-X[0], X[2j]-X[j], X[3j]-X[2j], and so forth.

This forms the raw data for computing the structure function at each lag scale. The file also includes a mmar_vecToArray helper, which converts a vector to a double array, along with a vector-input overload of general_hurst that utilizes it.

This allows other modules to call the GHE estimator using either arrays or vectors, depending on convenience.

The central GHE function iterates over lag scales, from a lower to an upper bound. At each scale, it computes the ratio of the q-th moment of detrended differences to the q-th moment of detrended levels. The log-log slope of this quantity, when plotted against the scale, provides an estimate of H(q)/q:

mqll

1//+------------------------------------------------------------------+
2//| Generalized Hurst exponent H(q) - array input |
3//+------------------------------------------------------------------+
4double general_hurst(double &data[], int q, int lower = 2, int upper = 50)
5 {
6 if(data.Size < 100)
7 {
8 Print("data array is of insufficient length");
9 return EMPTY_VALUE;
10 }
11
12 if(lower >= upper || lower < 2 || upper > (int)MathFloor(0.5 * data.Size))
13 {
14 Print("Invalid input for lower and/or upper");
15 return EMPTY_VALUE;
16 }
17
18 uint len = data.Size;
19
20 int k = 0;
21
22 matrix H;
23 if(!H.Resize((ulong)(upper - lower), 1))
24 {
25 Print(__LINE__, " ", __FUNCTION__, " ", GetLastError);
26 return EMPTY_VALUE;
27 }
28
29//--- Main loop over lag scales
30 for(int i = lower; i < upper; i++)
31 {
32 vector x_vector(i);
33 mmar_arange<double>(x_vector, 1.0, 1.0);
34
35 matrix mcord;
36 if(!mcord.Resize((ulong)i, 1))
37 {
38 Print(__LINE__, " ", __FUNCTION__, " ", GetLastError);
39 return EMPTY_VALUE;
40 }
41
42 mcord.Fill(0.0);
43
44 for(int j = 1; j < i + 1; j++)
45 {
46 vector dv, Y;
47 if(!mmar_diff_array(j, data, dv, Y))
48 return EMPTY_VALUE;
49
50 double N = (double)Y.Size;
51
52 vector X(N);
53 mmar_arange<double>(X, 1.0, 1.0);
54
55 double mx = X.Sum / N;
56 vector XP = MathPow(X, 2.0);
57 double SSxx = XP.Sum - N * MathPow(mx, 2.0);
58 double my = Y.Sum / N;
59 vector XY = X * Y;
60 double SSxy = XY.Sum - N * mx * my;
61 double cc1 = SSxy / SSxx;
62 double cc2 = my - cc1 * mx;
63
64 vector ddVd = dv - cc1;
65 vector VVVd = Y - cc1 * X - cc2;
66 vector PddVd = MathAbs(ddVd);
67 PddVd = MathPow(PddVd, q);
68 vector PVVVd = MathAbs(VVVd);
69 PVVVd = MathPow(PVVVd, q);
70
71 mcord[j - 1][0] = PddVd.Mean / PVVVd.Mean;
72 }
73
74 vector Px_vector = MathLog10(x_vector);
75
76 double mx = Px_vector.Mean;
77 vector Sqx = MathPow(Px_vector, 2.0);
78 double SSxx = Sqx.Sum - i * MathPow(mx, 2.0);
79 matrix lmcord = MatMath.Log10(mcord);
80 double my = lmcord.Mean;
81 vector pt = Px_vector * lmcord.Col(0);
82 double SSxy = pt.Sum - i * mx * my;
83
84 H[k][0] = SSxy / SSxx;
85
86 k++;
87 }
88
89 return H.Mean / (double)q;
90 }

The algorithm operates in three nested stages. The outermost loop iterates through lag window sizes, from a specified lower bound to an upper bound.

For each window size i, the inner loop calculates the structure function for lags ranging from 1 to i. At each lag j, the data is differenced at that lag and then detrended using simple linear regression (to remove any underlying drift).

The ratio of the q-th power of the detrended differences to the q-th power of the detrended levels is recorded. Once the inner loop completes, the log-log slope of this ratio against the lag index yields an 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 identify the moment order q* where tau(q*) = 0, and then H = 1/q*.

This approach is theoretically more robust because it directly utilizes the partition function that we are already computing. 9), which can occur if the tau(q) curve does not cross zero cleanly due to noise or insufficient data.

Partition Analysis — The Core Engine

With the necessary utilities in place, we now turn to 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 Sq(dt) over 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 assesses whether the data genuinely exhibits multifractal behavior using three independent diagnostic tests.

Class Structure and Configuration

The class offers a clean, two-step interface: first, initialize with data, then perform the analysis. All configuration parameters have sensible defaults but can be overridden before invoking Analyze:

mqll

1//+------------------------------------------------------------------+
2//| CPartitionAnalysis class |
3//+------------------------------------------------------------------+
4class CPartitionAnalysis
5 {
6private:
7 //--- Input data
8 double m_returns[];
9 int m_n_returns;
10
11 //--- Configuration
12 int m_delta_t_min;
13 int m_delta_t_max;
14 double m_delta_t_spacing;
15 double m_q_min;
16 double m_q_max;
17 double m_q_step;
18 double m_min_r_squared;
19
20 //--- Results
21 double m_q_values[];
22 int m_n_q;
23 double m_tau_q[];
24 double m_slopes[];
25 double m_r_squared[];
26 double m_H;
27 double m_sample_volatility;
28 int m_multifractal_score;
29 bool m_is_multifractal;
30 bool m_analyzed;
31
32 //--- Internal computation
33 double CalculatePartitionFunction(double q, int delta_t);
34 void GenerateDeltaTValues(int &dt_values[], int &n_dt);
35 void GenerateQValues(void);
36 bool FitPartitionPlots(const int &dt_values[], int n_dt);
37 double EstimateHurst(void);
38 double EstimateHurstFromTauQ(void);
39 int CheckMultifractality(void);
40 double ComputeSampleVolatility(void);
41
42public:
43 CPartitionAnalysis(void);
44 ~CPartitionAnalysis(void);
45
46 //--- Configuration
47 void SetPartitionConfig(int dt_min, int dt_max, double spacing);
48 void SetMomentConfig(double q_min, double q_max, double q_step);
49 void SetMinRSquared(double min_r2) { m_min_r_squared = min_r2; }
50
51 //--- Main pipeline
52 bool Init(const double &returns[], int n_returns);
53 bool Analyze(void);
54
55 //--- Results access
56 double GetH(void) { return m_H; }
57 double GetSampleVolatility(void) { return m_sample_volatility; }
58 void GetTauQ(double &tau[], double &q_vals[]);
59 double GetMeanRSquared(void);
60 bool IsMultifractal(void) { return m_is_multifractal; }
61 int GetMultifractalScore(void) { return m_multifractal_score; }
62 int GetNumQ(void) { return m_n_q; }
63 };

The configuration is naturally divided into two groups. The "partition config" governs the time scale grid: dt_min and dt_max define the range, and spacing controls the logarithmic density of sampling within that range.

The "moment config" controls the q-value grid: we sweep from q_min to q_max in fixed steps. 5 yields 61 moment orders, each requiring its own log-log regression.

This grid must be sufficiently fine 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 fixed increments. A linear spacing from 1 to 150 would yield 150 points, with most concentrated at smaller scales. Logarithmic spacing, conversely, allocates roughly equal representation to each decade of scale:

mqll

1//+------------------------------------------------------------------+
2//| GenerateDeltaTValues - Build log-spaced time scale array |
3//+------------------------------------------------------------------+
4void CPartitionAnalysis::GenerateDeltaTValues(int &dt_values[], int &n_dt)
5 {
6 int temp[];
7 ArrayResize(temp, 200);
8 int count = 0;
9
10 double current = (double)m_delta_t_min;
11 while(current <= (double)m_delta_t_max)
12 {
13 int dt = (int)MathRound(current);
14 if(dt < 1)
15 dt = 1;
16 if(dt > m_n_returns / 2)
17 break;
18
19 bool duplicate = false;
20 for(int i = 0; i < count; i++)
21 {
22 if(temp[i] == dt)
23 {
24 duplicate = true;
25 break;
26 }
27 }
28
29 if(!duplicate)
30 {
31 if(count >= ArraySize(temp))
32 ArrayResize(temp, count + 100);
33 temp[count++] = dt;
34 }
35
36 current *= m_delta_t_spacing;
37 }
38
39 for(int i = 0; i < count - 1; i++)
40 for(int j = i + 1; j < count; j++)
41 if(temp[i] > temp[j])
42 {
43 int swap = temp[i];
44 temp[i] = temp[j];
45 temp[j] = swap;
46 }
47
48 n_dt = count);
49 ArrayResize(dt_values, count);
50 ArrayCopy(dt_values, temp, 0, 0, count);
51 }

at each step. , extending up to dt_max.

Since these are time scales representing integer bar counts, each value is rounded, and duplicates are removed. The condition dt > m_n_returns / 2 serves as a safeguard, preventing the creation of time windows larger than half the available data.

This ensures a sufficient number of non-overlapping blocks for a meaningful partition function estimate. The resulting array is sorted in ascending order for correct physical interpretation during log-log regression.

Computing the Partition Function

At the heart of multifractal analysis lies the partition function Sq(dt). For given moment order q and time scale dt, this function answers: what is the average q-th power of absolute returns measured across non-overlapping windows of size dt?

mqll

1//+------------------------------------------------------------------+
2//| CalculatePartitionFunction - S_q(dt) for given q and time scale |
3//+------------------------------------------------------------------+
4double CPartitionAnalysis::CalculatePartitionFunction(double q, int delta_t)
5 {
6 int N = m_n_returns / delta_t;
7 if(N <= 0)
8 return -1.0;
9
10 double sum_sq = 0.0;
11 for(int i = 0; i < N; i++)
12 {
13 double interval_return = 0.0;
14 int start = i * delta_t;
15 for(int j = 0; j < delta_t; j++)
16 interval_return += m_returns[start + j];
17
18 double abs_r = MathAbs(interval_return);
19 if(abs_r < 1e-300)
20 abs_r = 1e-300;
21 sum_sq += MathPow(abs_r, q);
22 }
23
24 return sum_sq / N;
25 }

This function divides the return series into N non-overlapping blocks, each of length delta_t. For each block, it sums the individual log-returns to compute the block's aggregate return, takes its absolute value, raises it to the power q, and accumulates the result.

The final output is the average of these q-th power sums across all blocks. A floor value of 1e-300 is applied to prevent taking the logarithm of zero during the subsequent log-log fitting step.

This acts as a defensive measure, particularly important for quiet market periods where multiple consecutive returns might be exactly zero (especially on higher timeframes or with illiquid instruments).

The crucial insight is that if the data exhibits scaling behavior, then Sq(dt) follows a power law: Sq(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 provides tau(q) + 1.

Figure 2 illustrates this process dynamically: returns are aggregated into blocks, the partition function is calculated at various scales, and the log-log slopes reveal distinct scaling rates for different moment orders q:

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

Fitting the Log-Log Plots

For each moment order q within our defined grid, we compute Sq(dt) at every time scale, then take the base-10 logarithm of both quantities, and perform an OLS regression:

mqll

1//+------------------------------------------------------------------+
2//| FitPartitionPlots - Log-log regression for each q value |
3//+------------------------------------------------------------------+
4bool CPartitionAnalysis::FitPartitionPlots(const int &dt_values[], int n_dt)
5 {
6 ArrayResize(m_tau_q, m_n_q);
7 ArrayResize(m_slopes, m_n_q);
8 ArrayResize(m_r_squared, m_n_q);
9
10 int valid_q = 0;
11
12 for(int qi = 0; qi < m_n_q; qi++)
13 {
14 double q = m_q_values[qi];
15
16 double log_dt[];
17 double log_sq[];
18 ArrayResize(log_dt, n_dt);
19 ArrayResize(log_sq, n_dt);
20 int n_valid = 0;
21
22 for(int di = 0; di < n_dt; di++)
23 {
24 double s_q = CalculatePartitionFunction(q, dt_values[di]);
25 if(s_q > 0.0)
26 {
27 log_dt[n_valid] = MathLog10((double)dt_values[di]);
28 log_sq[n_valid] = MathLog10(s_q);
29 n_valid++;
30 }
31 }
32
33 if(n_valid < 3)
34 {
35 m_slopes[qi] = 0.0;
36 m_tau_q[qi] = -1.0;
37 m_r_squared[qi] = 0.0;
38 continue;
39 }
40
41 vector y_vec;
42 y_vec.Resize(n_valid);
43 matrix x_mat;
44 x_mat.Resize(n_valid, 2);
45
46 for(int i = 0; i < n_valid; i++)
47 {
48 y_vec[i] = log_sq[i];
49 x_mat[i][0] = log_dt[i];
50 x_mat[i][1] = 1.0;
51 }
52
53 OLS ols;
54 if(!ols.Fit(y_vec, x_mat))
55 {
56 m_slopes[qi] = 0.0;
57 m_tau_q[qi] = -1.0;
58 m_r_squared[qi] = 0.0;
59 continue;
60 }
61
62 vector params = ols.ModelParameters;
63 m_slopes[qi] = params[0];
64 m_tau_q[qi] = params[0] - 1.0;
65 m_r_squared[qi] = ols.Rsqe;
66 valid_q++;
67 }
68
69 if(valid_q < 3)
70 {
71 Print("CPartitionAnalysis::FitPartitionPlots - detected only ", valid_q, " valid q values for analysis");
72 return false;
73 }
74
75 return true;
76 }

For each q value, we construct the design matrix with two columns: log10(dt) as the predictor and a column of ones for the intercept. The dependent variable is log10(S_q).

After fitting, params[0] provides the slope, which equals tau(q) + 1. We subtract 1 to directly obtain tau(q).

The R² value indicates how well the power-law model fits for this particular q. If fewer than 3 valid data points are available (which can occur at extreme q values due to numerical overflow or underflow), we skip that q value entirely to avoid 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 implies that we can estimate H by locating where the tau(q) curve intersects the zero axis and then inverting that q* value. The implementation uses linear interpolation between the two q values that bracket the zero crossing:

mqll

1//+------------------------------------------------------------------+
2//| EstimateHurstFromTauQ - Find q* where tau(q*)=0, then H=1/q* |
3//+------------------------------------------------------------------+
4double CPartitionAnalysis::EstimateHurstFromTauQ(void)
5 {
6 int idx_below = -1;
7 int idx_above = -1;
8
9 for(int i = 0; i < m_n_q - 1; i++)
10 {
11 if(m_tau_q[i] * m_tau_q[i + 1] <= 0.0)
12 {
13 idx_below = i;
14 idx_above = i + 1;
15 break;
16 }
17 }
18
19 if(idx_below < 0)
20 {
21 double min_abs = DBL_MAX;
22 int min_idx = 0;
23 for(int i = 0; i < m_n_q; i++)
24 {
25 if(MathAbs(m_tau_q[i]) < min_abs)
26 {
27 min_abs = MathAbs(m_tau_q[i]);
28 min_idx = i;
29 }
30 }
31 double q_zero = m_q_values[min_idx];
32 if(q_zero <= 0.0)
33 return 0.5;
34 return MathMax(0.1, MathMin(0.9, 1.0 / q_zero));
35 }
36
37 double q1 = m_q_values[idx_below];
38 double q2 = m_q_values[idx_above];
39 double t1 = m_tau_q[idx_below];
40 double t2 = m_tau_q[idx_above];
41
42 double q_zero;
43 if(MathAbs(t2 - t1) > 1e-15)
44 q_zero = q1 + (0.0 - t1) * (q2 - q1) / (t2 - t1);
45 else
46 q_zero = (q1 + q2) / 2.0;
47
48 if(q_zero <= 0.0)
49 return 0.5;
50 return MathMax(0.1, MathMin(0.9, 1.0 / q_zero));
51 }

The sign-change detection (tau[i] * tau[i+1] <= 0) identifies adjacent q values where tau transitions from negative to positive (or vice versa). Linear interpolation between these two points provides the precise q* where tau(q*) = 0.

Consequently, *H = 1/q*. 9] because values outside this range are physically implausible for financial data and often indicate a numerical anomaly.

If no zero crossing is found (which can occur if the data exhibits very weak scaling), the fallback mechanism identifies the q value closest to zero and uses that as an approximation.

The parent EstimateHurst method first attempts the tau(q) method. If this method yields an implausible value, it falls back to the GHE estimator discussed in the previous section.

The Multifractality Decision

Extracting tau(q) is only one part of the task. We also need to ascertain whether the scaling is genuinely multifractal (implying a nonlinear curve) or merely monofractal (a straight line, indicating identical scaling for all moments). Our library employs three independent diagnostic tests, each scored from 0 to 3, for a maximum total score of 9:

mqll

1//+------------------------------------------------------------------+
2//| CheckMultifractality - Score 0-9 based on three diagnostics |
3//+------------------------------------------------------------------+
4int CPartitionAnalysis::CheckMultifractality(void)
5 {
6 if(m_n_q < 5)
7 return 0;
8
9//--- Test 1: Concavity of tau(q)
10 int n_concave = 0;
11 int n_d2 = m_n_q - 2;
12 for(int i = 0; i < n_d2; i++)
13 {
14 double d2 = m_tau_q[i + 2] - 2.0 * m_tau_q[i + 1] + m_tau_q[i];
15 if(d2 < 0.0)
16 n_concave++;
17 }
18 double concave_pct = 100.0 * n_concave / n_d2;
19
20 int concavity_score;
21 if(concave_pct >= 80.0)
22 concavity_score = 3;
23 else
24 if(concave_pct >= 60.0)
25 concavity_score = 2;
26 else
27 if(concave_pct >= 40.0)
28 concavity_score = 1;
29 else
30 concavity_score = 0;
31
32//--- Test 2: Linearity (high R^2 -> monofractal)
33 double sum_q = 0.0, sum_t = 0.0, sum_qq = 0.0, sum_qt = 0.0;
34 int n_pos = 0;
35 for(int i = 0; i < m_n_q; i++)
36 {
37 if(m_q_values[i] > 0.0)
38 {
39 sum_q += m_q_values[i];
40 sum_t += m_tau_q[i];
41 sum_qq += m_q_values[i] * m_q_values[i];
42 sum_qt += m_q_values[i] * m_tau_q[i];
43 n_pos++;
44 }
45 }
46
47 double mean_q = sum_q / n_pos;
48 double mean_t = sum_t / n_pos;
49 double ss_qq = sum_qq - n_pos * mean_q * mean_q;
50 double ss_qt = sum_qt - n_pos * mean_q * mean_t;
51 double slope_lin = ss_qt / ss_qq;
52 double intercept_lin = mean_t - slope_lin * mean_q;
53
54 double ss_res = 0.0, ss_tot = 0.0;
55 for(int i = 0; i < m_n_q; i++)
56 {
57 if(m_q_values[i] > 0.0)
58 {
59 double fitted = slope_lin * m_q_values[i] + intercept_lin;
60 ss_res += (m_tau_q[i] - fitted) * (m_tau_q[i] - fitted);
61 ss_tot += (m_tau_q[i] - mean_t) * (m_tau_q[i] - mean_t);
62 }
63 }
64 double r2_linear = (ss_tot > 0.0) ? 1.0 - ss_res / ss_tot : 1.0;
65 double rmse_lin = MathSqrt(ss_res / n_pos);
66
67 int linearity_score;
68 if(r2_linear > 0.99 && rmse_lin < 0.05)
69 linearity_score = 0;
70 else
71 if(r2_linear > 0.95)
72 linearity_score = 1;
73 else
74 if(r2_linear > 0.90)
75 linearity_score = 2;
76 else
77 linearity_score = 3;
78
79//--- Test 3: Variation in dtau/dq slopes
80 double d_tau_min = DBL_MAX, d_tau_max = -DBL_MAX;
81 for(int i = 0; i < m_n_q - 1; i++)
82 {
83 double d = m_tau_q[i + 1] - m_tau_q[i];
84 if(d < d_tau_min)
85 d_tau_min = d;
86 if(d > d_tau_max)
87 d_tau_max = d;
88 }
89 double slope_range = d_tau_max - d_tau_min;
90
91 int variation_score;
92 if(slope_range > 0.5)
93 variation_score = 3;
94 else
95 if(slope_range > 0.3)
96 variation_score = 2;
97 else
98 if(slope_range > 0.1)
99 variation_score = 1;
100 else
101 variation_score = 0;
102
103 int total = concavity_score + linearity_score + variation_score;
104
105 PrintFormat("Multifractality check: concavity=%d linearity=%d variation=%d total=%d/9",
106 concavity_score, linearity_score, variation_score, total);
107
108 return total;
109 }

Test 1 — Concavity: A multifractal tau(q) curve must exhibit concavity (curve downwards). We compute the second finite difference at each interior point.

If d2 = tau[i+2] - 2*tau[i+1] + tau[i] is negative, that point is concave. Achieving a score of 3 requires 80% or more of the points to be concave.

A perfectly linear tau(q) would have d2 approximately zero everywhere, resulting in a score of 0.

Test 2 — Linearity: This test acts as an inverse check. We fit a straight line directly to the tau(q) curve and evaluate its goodness of fit.

05, the curve is considered essentially linear, indicating monofractal data (scoring 0). This dual condition helps prevent false negatives: a high R² alone might occur if the tau(q) curve has a subtle, but systematic, curvature that happens to correlate well with a line; the RMSE then correctly identifies such cases where residuals remain meaningfully large.

Lower R² values signify a significant deviation from a linear relationship, suggesting genuine multifractal structure (scoring 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, whereas the slope at large q reflects the behavior of small fluctuations.

5), it scores 3. A monofractal process, in contrast, maintains a nearly constant slope throughout.

A total score of 5 or higher (out of a possible 9) classifies the data as multifractal. This threshold was empirically calibrated during development against known multifractal and monofractal processes. It is intentionally conservative: we prefer to err on the side of caution by potentially missing some weakly multifractal data rather than falsely classifying random noise as multifractal.

Figure 3 visually differentiates between monofractal and multifractal tau(q) curves and demonstrates how the three diagnostic tests identify this distinction:

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 now in place, we can develop a test script to load real market data from the MetaTrader terminal and execute the full partition analysis pipeline. This script will serve as our initial end-to-end validation, confirming that the MQL5 implementation yields sensible results:

mqll

1//+------------------------------------------------------------------+
2//| MMAR Partition Test |
3//| Copyright 2026, MMQ - Muhammad Minhas Qamar |
4//| |
5//+------------------------------------------------------------------+
6#property copyright "Copyright 2026, MMQ - Muhammad Minhas Qamar"
7#property link "
8#property version "1.00"
9#property strict
10#property script_show_inputs
11
12#include <MMAR\PartitionAnalysis.mqh>
13
14input int InpBars = 50000;
15input double InpQMin = 0.01;
16input double InpQMax = 30.0;
17input double InpQStep = 0.5;
18input int InpDtMin = 1;
19input int InpDtMax = 150;
20input double InpDtSpacing = 1.1;
21input double InpMinR2 = 0.60;
22
23//+------------------------------------------------------------------+
24//| Script program start function |
25//+------------------------------------------------------------------+
26void OnStart
27 {
28 PrintFormat("=== MMAR Partition Analysis Test ===");
29 PrintFormat("Symbol: %s | Period: %s", _Symbol, EnumToString(_Period));
30
31//--- Load close prices and compute log returns
32 double close[];
33 int copied = CopyClose(_Symbol, _Period, 0, InpBars, close);
34 if(copied < 1000)
35 {
36 PrintFormat("Not enough bars: %d (need 1000+)", copied);
37 return;
38 }
39 PrintFormat("Loaded %d bars", copied);
40
41 double returns[];
42 int n_returns = copied - 1;
43 ArrayResize(returns, n_returns);
44 for(int i = 0; i < n_returns; i++)
45 returns[i] = MathLog(close[i + 1] / close[i]);
46
47//--- Run partition analysis
48 CPartitionAnalysis partition;
49 partition.SetPartitionConfig(InpDtMin, InpDtMax, InpDtSpacing);
50 partition.SetMomentConfig(InpQMin, InpQMax, InpQStep);
51 partition.SetMinRSquared(InpMinR2);
52
53 if(!partition.Init(returns, n_returns))
54 {
55 Print("Init failed");
56 return;
57 }
58
59 uint start = GetTickCount;
60 if(!partition.Analyze)
61 {
62 Print("Analyze failed");
63 return;
64 }
65 uint elapsed = GetTickCount - start;
66 PrintFormat("Analysis time: %u ms", elapsed);
67
68//--- Print results
69 PrintFormat("\n--- Results ---");
70 PrintFormat("H (Hurst exponent): %.4f", partition.GetH);
71 PrintFormat("Sample volatility: %.6f", partition.GetSampleVolatility);
72 PrintFormat("Mean R-squared: %.4f", partition.GetMeanRSquared);
73 PrintFormat("Multifractal score: %d/9", partition.GetMultifractalScore);
74 PrintFormat("Is multifractal: %s", partition.IsMultifractal ? "YES" : "NO");
75
76//--- Print tau(q) curve
77 double tau[], q_vals[];
78 partition.GetTauQ(tau, q_vals);
79 int n_q = partition.GetNumQ;
80
81 PrintFormat("\n--- tau(q) curve (every 5th value) ---");
82 PrintFormat("%8s %10s", "q", "tau(q)");
83 for(int i = 0; i < n_q; i += 5)
84 PrintFormat("%8.2f %10.4f", q_vals[i], tau[i]);
85
86 Print("\n=== Partition Analysis Test Complete ===");
87 }

The script follows a straightforward pattern that will be reiterated throughout this series. First, it loads historical close prices using CopyClose, requesting up to InpBars (defaulting to 50,000) and requiring a minimum of 1,000 bars for meaningful analysis.

Next, it computes log-returns: the natural logarithm of consecutive price ratios. These log-returns serve as the input for all subsequent analysis.

Finally, the script configures the partition analysis parameters, initializes the CPartitionAnalysis object with the return data, calls the Analyze method, and then outputs the results.

Let's examine the output when this script is run on EURUSD M10 with InpBars set to 500,000 (overriding the default of 50,000 to stress-test the pipeline with a larger sample):

text

1=== MMAR Partition Analysis Test ===
2Symbol: EURUSD | Period: PERIOD_M10
3Loaded 500000 bars
4Multifractality check: concavity=3 linearity=1 variation=1 total=5/9
5Analysis time: 3859 ms
6
7--- Results ---
8H (Hurst exponent): 0.4892
9Sample volatility: 0.000422
10Mean R-squared: 0.8867
11Multifractal score: 5/9
12Is multifractal: YES
13
14--- tau(q) curve (every 5th value) ---
15 q tau(q)
16 0.01 -0.9930
17 2.51 0.1724
18 5.01 0.7441
19 7.51 1.2193
20 10.01 1.7181
21 12.51 2.2156
22 15.01 2.7076
23 17.51 3.1952
24 20.01 3.6796
25 22.51 4.1617
26 25.01 4.6422
27 27.51 5.1215
28
29=== Partition Analysis Test Complete ===

The results reveal a clear narrative. 5, indicates that the overall drift structure is akin to a random walk.

This aligns with the efficient market hypothesis at lower frequencies. However, with a multifractal score of 5/9, the data surpasses our threshold, signaling genuine multifractal scaling.

The concavity test awarded a perfect 3/3 score, confirming the strong concavity of the tau(q) curve, which implies that different moment orders scale at distinct rates. 8867 reinforces confidence in the power-law model's fit across the q range, validating the meaningfulness of the scaling exponents.

The tau(q) curve itself unveils the multifractal signature: it begins negatively (approximately -0.99 at q=0.01), crosses zero near q=2 (consistent with H approximately 0.5, as 1/H is roughly 2), and then exhibits sublinear growth. If the data were monofractal, tau(q) would trace a straight line with a slope equal to H. The visible curvature, particularly the slight decrease in increments between successive values as q increases, is precisely what the concavity test detects.

Consider a contrasting scenario. Running the identical analysis on EURUSD M30 with 50,000 bars yields a score of only 2/9, failing to meet our multifractality threshold.

This is not a failure of the code; rather, it correctly identifies that, at this specific combination of timeframe and sample size, the multifractal structure is too weak to be statistically significant. The partition analysis acts as an honest gatekeeper: if the data does not support multifractal modeling, the library informs you as such, preventing fruitless computation on cascade generation and Monte Carlo simulations.

This behavior has practical implications. When deploying the MMAR in an Expert Advisor, it is advisable to select your timeframe and lookback window based on where the data genuinely exhibits multifractal scaling.

The M10 timeframe with a substantial sample (500,000 bars, representing several years) passes the test. Shorter timeframes with fewer bars may not.

The library provides the necessary diagnostic tools to make this decision empirically, rather than relying on speculation.

Conclusion

In this article, we successfully developed the first major module of the MMAR library for MetaTrader 5. Beginning with shared type definitions, we progressively built an OLS regression utility, a Generalized Hurst Exponent estimator, and the complete Partition Analysis engine. The outcome is a self-contained MQL5 class capable of processing any return series, calculating its multifractal scaling properties, and making a rigorous determination regarding the appropriateness of multifractal modeling for that particular dataset.

MMARConstants.mqh establishes the shared vocabulary, encompassing distribution types, model states, and parameter containers utilized across all modules.
OLS.mqh provides SVD-based linear regression with comprehensive diagnostics, which is internally employed for extracting log-log slopes.
HurstGHE.mqh offers a robust fallback estimator for the Hurst exponent, to be used when the primary tau(q) method yields improbable results.
PartitionAnalysis.mqh orchestrates the entire pipeline: generating time scales, computing the partition function, performing log-log fitting, estimating the Hurst exponent, and executing the three-test multifractality diagnostic.

The analysis completes in under 4 seconds for 500,000 bars, a speed suitable for periodic recalibration within a live trading environment. Crucially, it provides an honest assessment: not all market data is multifractal, and the library accurately identifies instances where this assumption does not hold.

In Part 5, we will take the tau(q) curve generated here and fit it to four theoretical distributions using the Legendre transform and bounded optimization. This will provide the distribution parameters necessary for generating the multiplicative cascade, representing the next step toward a complete MMAR simulation engine.

The programs presented in this article are solely for educational purposes. They are not designed for live trading and are not optimized for real-market performance. The code serves as a foundation for understanding multifractal analysis; always conduct thorough validation on demo accounts before considering any adaptation for live trading.

Getting the Source Code via MQL5 Algo Forge

All source files are provided below; however, the complete repository is also accessible on MQL5 Algo Forge, the community’s Git-based platform for sharing and collaborating on trading projects. Algo Forge maintains version history both locally and in the cloud, ensuring you always have access to the latest code, including any updates or fixes implemented after this article’s publication.

To clone the repository, open a terminal (e.g., Command Prompt, PowerShell, or your editor's integrated terminal) and execute:

bash

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

This command requires Git to be installed on your system. 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 within your current directory.

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

bash

1git 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.

Muhammad Minhas Qamar

Pakistan
3304

Developer by Profession, Trader by Hobby

Gmail: ayanminhasshayar@gmail.com

Beyond GARCH (Part IV): Partition Analysis in MQL5 - MQL5 Articles

Introduction

Library Architecture and Constants

The OLS Foundation

Estimating the Hurst Exponent

Partition Analysis — The Core Engine

Running the Analysis on Live Data

Conclusion

Other articles by this author