Trading systems

Beyond GARCH (Part V): Fitting the Multifractal Spectrum in MQL5 - MQL5 Articles

June 9, 202624 min read

Create your own AI for tradingRead our book "Neural Networks in Algo Trading with MQL5"Begin

MetaTrader 5 / Trading systems

Introduction

com/en/articles/21299), we built the partition analysis engine: the module that takes raw price data, computes the scaling function τ(q), estimates the Hurst exponent, and determines whether the data exhibits multifractal behavior. The output of that module is a curve, τ(q), that encodes how different statistical moments scale across time.

But a curve alone is not enough to generate synthetic price paths. To build the multiplicative cascade that creates multifractal trading time, we need a specific distribution and its parameters.

This article bridges that gap. We will implement the Spectrum Fitter, the module that transforms the raw τ(q) curve into a fitted probability distribution ready for simulation.

The process has two stages. First, we apply the Legendre transform to convert τ(q) into the singularity spectrum f(α), a more physically interpretable representation of the multifractal structure.

Second, we fit four candidate distributions to the empirical spectrum using bounded optimization (BLEIC from ALGLIB). We then select the model with the lowest sum of squared errors.

The winner provides the exact parameters that the cascade generator will use in Part 6.

We will cover:

From τ(q) to f(α): The Legendre Transform

The scaling function τ(q) tells us how the q-th moment of returns scales with time. But there is a more intuitive representation of multifractal structure: the singularity spectrum f(α).

Here, α represents the local Hölder exponent (how "rough" or "smooth" the process is at a given point), and f(α) represents the fractal dimension of the set of points that share that exponent. A broad f(α) curve means the process has many different local scaling behaviors, which is the hallmark of multifractality.

The two representations are connected by the Legendre transform:

text

1α(q) = d(τ(q)) / dq
2
3f(α) = q * α(q) - τ(q)

The first equation says that α at a given q is the slope of the τ(q) curve at that point. The second equation computes f(α) from the tangent line relationship.

Figure 1 animates the process. A tangent line sweeps along the τ(q) curve; at each position, its slope gives α, and the tangent-line relationship yields f(α).

Point by point, the singularity spectrum builds up on the right, revealing the characteristic parabolic shape whose width measures the degree of multifractality.

Fig. 1. The Legendre transform in action: each tangent slope on τ(q) maps to a point on the singularity spectrum f(α)

In practice, since we have τ(q) at discrete points rather than as a continuous function, we approximate the derivative using central finite differences:

text

1//+------------------------------------------------------------------+
2//| ComputeSpectrum - Legendre transform: tau(q) → f(alpha)          |
3//+------------------------------------------------------------------+
4bool CSpectrumFitter::ComputeSpectrum(void)
5  {
6   double temp_alpha[];
7   double temp_f[];
8   ArrayResize(temp_alpha, m_n_q);
9   ArrayResize(temp_f, m_n_q);
10   int count = 0;
11
12   for(int i = 1; i < m_n_q - 1; i++)
13     {
14      if(m_q_values[i] < 0.5)
15         continue;
16
17      double q_prev = m_q_values[i - 1];
18      double q_next = m_q_values[i + 1];
19      double tau_prev = m_tau_q[i - 1];
20      double tau_next = m_tau_q[i + 1];
21
22      double dq = q_next - q_prev;
23      if(MathAbs(dq) < 1e-15)
24         continue;
25
26      double alpha = (tau_next - tau_prev) / dq;
27      double f_alpha = alpha * m_q_values[i] - m_tau_q[i];
28
29      if(!MathIsValidNumber(alpha) || !MathIsValidNumber(f_alpha))
30         continue;
31
32      if(f_alpha < -0.5)
33         continue;
34
35      temp_alpha[count] = alpha;
36      temp_f[count] = f_alpha;
37      count++;
38     }
39
40   if(count < 3)
41     {
42      Print("SpectrumFitter::ComputeSpectrum - only ", count, " valid spectrum points");
43      return false;
44     }
45
46   m_n_spectrum = count;
47   ArrayResize(m_alpha, count);
48   ArrayResize(m_f_alpha, count);
49   ArrayCopy(m_alpha, temp_alpha, 0, 0, count);
50   ArrayCopy(m_f_alpha, temp_f, 0, 0, count);
51
52   double a_min = m_alpha[0], a_max = m_alpha[0];
53   double f_max = m_f_alpha[0];
54   int f_max_idx = 0;
55   for(int i = 1; i < m_n_spectrum; i++)
56     {
57      if(m_alpha[i] < a_min)
58         a_min = m_alpha[i];
59      if(m_alpha[i] > a_max)
60         a_max = m_alpha[i];
61      if(m_f_alpha[i] > f_max)
62        {
63         f_max = m_f_alpha[i];
64         f_max_idx = i;
65        }
66     }
67
68   PrintFormat("SpectrumFitter: %d spectrum points, alpha=[%.4f, %.4f], peak at alpha_0=%.4f f=%.4f",
69               m_n_spectrum, a_min, a_max, m_alpha[f_max_idx], f_max);
70
71   return true;
72  }

Several decisions are embedded in this implementation. 5 because, at very small q, the partition function is dominated by near-zero returns and the numerical derivative becomes unreliable.

We use the central difference formula (tau[i+1] - tau[i-1]) / (q[i+1] - q[i-1]) for second-order accuracy rather than a forward or backward difference. 5.

Since f(α) represents a fractal dimension, it should theoretically be non-negative, but numerical noise from the discrete derivative can push it slightly negative. 5 threshold is generous enough to keep borderline points while discarding clearly nonsensical values.

The result is an empirical spectrum: a set of (α, f(α)) pairs that form a roughly parabolic curve. The peak of this curve occurs at α₀, where f(α₀) should be close to 1.0 (the fractal dimension of the full time series). The width of the curve measures the degree of multifractality: a wider spectrum means more heterogeneous scaling behavior.

Four Theoretical Spectrum Models

With the empirical spectrum in hand, we need to fit it to a parametric model. The choice of model determines which distribution will generate the cascade multipliers in the simulation stage. Our library implements four candidates, each with a different theoretical basis and number of parameters.

Normal Spectrum

The simplest model assumes the cascade multipliers are drawn from a log-normal distribution. The theoretical spectrum is a downward-opening parabola centered at α₀:

f(α) = 1 - (α - α₀)^2 / [4H(α₀ - H)]

text

1//+------------------------------------------------------------------+
2//| Normal spectrum model: f(alpha) = 1 - (alpha-a0)^2 / [4H(a0-H)]  |
3//+------------------------------------------------------------------+
4double CSpectrumFitter::NormalSpectrum(double alpha, double alpha_0, double H)
5  {
6   double denom = 4.0 * H * (alpha_0 - H);
7   if(denom <= 0.0)
8      return -1e10;
9   return 1.0 - (alpha - alpha_0) * (alpha - alpha_0) / denom;
10  }

This model has a single free parameter: α₀, the location of the spectrum peak. The constraint α₀ > H ensures the denominator is positive and the parabola opens downward. The Normal model produces symmetric spectra and works well when the data has homogeneous scaling at both extremes.

Binomial Spectrum

The binomial model assumes the cascade uses only two possible multiplier values at each split. The spectrum is entropy-based:

text

1f(α) = -[r1 * log(r1) + r2 * log(r2)] / log(2)
2
3where r1 = (α_max - α) / (α_max - α_min)
4and   r2 = (α - α_min) / (α_max - α_min)

text

1//+------------------------------------------------------------------+
2//| Binomial spectrum model: entropy-based                           |
3//+------------------------------------------------------------------+
4double CSpectrumFitter::BinomialSpectrum(double alpha, double alpha_min, double alpha_max)
5  {
6   if(alpha_max <= alpha_min)
7      return -1e10;
8
9   double a = MathMax(alpha_min + 1e-10, MathMin(alpha, alpha_max - 1e-10));
10
11   double r1 = (alpha_max - a) / (alpha_max - alpha_min);
12   double r2 = (a - alpha_min) / (alpha_max - alpha_min);
13
14   r1 = MathMax(r1, 1e-10);
15   r2 = MathMax(r2, 1e-10);
16
17   return -(r1 * MathLog(r1) + r2 * MathLog(r2)) / MathLog(2.0);
18  }

This model has two free parameters: α_min and α_max, defining the endpoints of the spectrum support. 0 when α is at the midpoint and falls to 0 at the extremes.

Despite its simplicity (only two multiplier values), the binomial model can produce excellent fits because the entropy function closely mimics the shape of many empirical spectra. The clamping of α to [α_min, α_max] and the floor on r1, r2 prevent log(0) errors at the boundaries.

Poisson Spectrum

The Poisson model uses multipliers derived from a Poisson-distributed random variable:

text

1//+------------------------------------------------------------------+
2//| Poisson spectrum model                                           |
3//+------------------------------------------------------------------+
4double CSpectrumFitter::PoissonSpectrum(double alpha, double alpha_0, double H, double b)
5  {
6   if(alpha_0 <= 0.0)
7      return -1e10;
8
9   double a = MathMax(alpha, 1e-10);
10   double term1 = 1.0 - alpha_0 / (H * MathLog(b));
11   double term2 = (a / H) * (MathLog(alpha_0 * M_E / a) / MathLog(b));
12   return term1 + term2;
13  }

The Poisson model has one free parameter (α₀) but produces an asymmetric spectrum, unlike the Normal model's symmetric parabola. The asymmetry comes from the Poisson distribution's inherent skewness. The cascade base b is fixed at 2 (binary cascade).

Gamma Spectrum

The most flexible model uses multipliers from a Gamma distribution:

text

1//+------------------------------------------------------------------+
2//| Gamma spectrum model                                             |
3//+------------------------------------------------------------------+
4double CSpectrumFitter::GammaSpectrum(double alpha, double alpha_0, double gamma, double b)
5  {
6   if(alpha_0 <= 0.0 || gamma <= 0.0)
7      return -1e10;
8
9   double a = MathMax(alpha, 1e-10);
10   double term1 = gamma * MathLog(a / alpha_0) / MathLog(b);
11   double term2 = gamma * (alpha_0 - a) / (alpha_0 * MathLog(b));
12   return 1.0 + term1 + term2;
13  }

With two free parameters (α₀ and γ), the Gamma model can produce both symmetric and asymmetric spectra. The γ shape parameter controls how peaked or flat the distribution is.

When γ is large, the spectrum narrows; when γ is small, it broadens. This flexibility often makes it competitive with the Binomial model in terms of fit quality.

With all four models defined, we can identify which one best matches the empirical spectrum. Figure 2 puts them head to head.

Each theoretical curve is drawn over the same set of empirical data points, and the sum of squared errors (SSE) reveals a clear hierarchy. The Normal parabola, constrained to symmetry, misses badly.

The Poisson improves but still drifts. The Gamma gets close.

And the Binomial, with its entropy-based shape, hugs the empirical dots almost exactly.

Fig. 2. Four distributions compete to fit the empirical spectrum: the Binomial wins decisively with SSE = 0.148

Bounded Optimization with BLEIC

Fitting each distribution model to the empirical spectrum is a nonlinear optimization problem: find the parameter values that minimize the sum of squared errors (SSE) between the theoretical spectrum model and the empirical points. We use the BLEIC (Boundary L-BFGS with Equality/Inequality Constraints) algorithm from the ALGLIB library, which is included with MetaTrader 5. BLEIC supports box constraints (upper and lower bounds on each parameter), which is essential because our parameters have physical constraints: α₀ must exceed H, α_max must exceed α_min, γ must be positive, and so on.

The optimizer works by evaluating the objective function (SSE) at different parameter values and iteratively moving toward the minimum. Since we use the derivative-free variant (MinBLEICCreateF), it approximates gradients via finite differences. This is slightly slower than providing analytical gradients but much simpler to implement and debug.

Figure 3 shows BLEIC's behavior for the Binomial fit. The SSE surface over the two-dimensional parameter space (α_min, α_max) forms a smooth bowl.

The optimizer starts at the midpoint of the feasible region and moves downhill while respecting the box constraints, with the step size shrinking as the minimum approaches. The purple dashed boundaries show the box constraints that keep the parameters physically meaningful.

Without them, the optimizer could wander into regions where α_min exceeds α_max, producing nonsensical results.

Fig. 3. BLEIC navigates the SSE surface: bounded optimization finds the minimum while respecting physical parameter constraints

Let us look at the Normal fit as an example of the code pattern. It optimizes the single parameter α₀:

text

1//+------------------------------------------------------------------+
2//| FitNormal - Optimize alpha_0 for Normal spectrum model           |
3//+------------------------------------------------------------------+
4double CSpectrumFitter::FitNormal(void)
5  {
6   double alpha_min_data = m_alpha[0];
7   double alpha_max_data = m_alpha[0];
8   for(int i = 1; i < m_n_spectrum; i++)
9     {
10      if(m_alpha[i] < alpha_min_data)
11         alpha_min_data = m_alpha[i];
12      if(m_alpha[i] > alpha_max_data)
13         alpha_max_data = m_alpha[i];
14     }
15
16   double lb = MathMax(m_H + 1e-6, alpha_min_data);
17   double ub = alpha_max_data;
18
19   if(lb >= ub)
20     {
21      PrintFormat("SpectrumFitter::FitNormal - no valid range (need alpha_0 > H=%.4f, alpha_max=%.4f)", m_H, alpha_max_data);
22      return DBL_MAX;
23     }
24
25   double x0[];
26   ArrayResize(x0, 1);
27   x0[0] = (lb + ub) / 2.0;
28
29   double bndl[];
30   ArrayResize(bndl, 1);
31   bndl[0] = lb;
32
33   double bndu[];
34   ArrayResize(bndu, 1);
35   bndu[0] = ub;
36
37   CMinBLEICState state;
38   CMinBLEIC::MinBLEICCreateF(1, x0, 1e-6, state);
39   CMinBLEIC::MinBLEICSetBC(state, bndl, bndu);
40   CMinBLEIC::MinBLEICSetCond(state, 1e-8, 1e-8, 1e-8, 200);
41
42   while(CMinBLEIC::MinBLEICIteration(state))
43     {
44      if(state.m_needf)
45        {
46         double alpha_0 = state.m_x[0];
47         state.m_f = ComputeSSE_Normal(alpha_0);
48        }
49     }
50
51   double result[];
52   CMinBLEICReport rep;
53   CMinBLEIC::MinBLEICResults(state, result, rep);
54
55   if(rep.m_terminationtype <= 0)
56     {
57      PrintFormat("SpectrumFitter::FitNormal - optimizer failed (type=%d)", rep.m_terminationtype);
58      return DBL_MAX;
59     }
60
61   m_params[MMAR_DIST_NORMAL][0] = result[0];
62   m_n_params[MMAR_DIST_NORMAL] = 1;
63   double sse = ComputeSSE_Normal(result[0]);
64   m_sse[MMAR_DIST_NORMAL] = sse;
65
66   PrintFormat("  Normal:   SSE=%.6f  alpha_0=%.6f", sse, result[0]);
67   return sse;
68  }

The pattern is consistent across all four distributions. We (1) derive parameter bounds from the data, (2) set the initial guess to the feasible midpoint, (3) configure BLEIC with convergence tolerances and maximum iterations, (4) run the loop that evaluates the SSE at each candidate point, and (5) extract the best parameters from the result. The termination type check ensures we do not accept results from a failed optimization (which can happen if the feasible region is too narrow or the objective is flat).

The SSE evaluation function is straightforward. It loops over all empirical spectrum points, computes the theoretical f(α) at each, and accumulates the squared differences:

text

1//+------------------------------------------------------------------+
2//| ComputeSSE_Normal                                                |
3//+------------------------------------------------------------------+
4double CSpectrumFitter::ComputeSSE_Normal(double alpha_0)
5  {
6   double sse = 0.0;
7   for(int i = 0; i < m_n_spectrum; i++)
8     {
9      double f_fit = NormalSpectrum(m_alpha[i], alpha_0, m_H);
10      double diff = m_f_alpha[i] - f_fit;
11      sse += diff * diff;
12     }
13   return sse;
14  }

The Binomial and Gamma fits follow the same pattern but with two-dimensional parameter spaces. The Binomial optimizes α_min and α_max simultaneously, with the constraint that α_max > α_min enforced both through the bounds and an explicit check inside the objective function (returning a penalty of 1e10 if violated). The Gamma optimizes α₀ and the shape parameter γ, bounded between 0.1 and 10.0.

The CSpectrumFitter Class

The full class brings together the Legendre transform and the four fitting routines into a clean pipeline. Let us first look at the class interface to understand the overall structure:

text

1//+------------------------------------------------------------------+
2//| CSpectrumFitter class                                            |
3//+------------------------------------------------------------------+
4class CSpectrumFitter
5  {
6private:
7   //--- Input data
8   double            m_q_values[];
9   double            m_tau_q[];
10   double            m_H;
11   int               m_n_q;
12
13   //--- Computed spectrum
14   double            m_alpha[];
15   double            m_f_alpha[];
16   int               m_n_spectrum;
17
18   //--- Fitting results
19   double            m_sse[4];
20   double            m_params[4][4];
21   int               m_n_params[4];
22   ENUM_MMAR_DISTRIBUTION m_best_dist;
23   bool              m_fitted;
24
25   //--- Spectrum model functions
26   double            NormalSpectrum(double alpha, double alpha_0, double H);
27   double            BinomialSpectrum(double alpha, double alpha_min, double alpha_max);
28   double            PoissonSpectrum(double alpha, double alpha_0, double H, double b);
29   double            GammaSpectrum(double alpha, double alpha_0, double gamma, double b);
30
31   //--- Distribution fitting
32   double            FitNormal(void);
33   double            FitBinomial(void);
34   double            FitPoisson(void);
35   double            FitGamma(void);
36
37   //--- SSE evaluation helpers
38   double            ComputeSSE_Normal(double alpha_0);
39   double            ComputeSSE_Binomial(double alpha_min, double alpha_max);
40   double            ComputeSSE_Poisson(double alpha_0);
41   double            ComputeSSE_Gamma(double alpha_0, double gamma);
42
43public:
44                     CSpectrumFitter(void);
45                    ~CSpectrumFitter(void);
46
47   //--- Initialization
48   bool              Init(const double &q_values[], const double &tau_q[], int n_q, double H);
49
50   //--- Computation
51   bool              ComputeSpectrum(void);
52   bool              FitAllDistributions(void);
53
54   //--- Results access
55   ENUM_MMAR_DISTRIBUTION GetBestDistribution(void) { return m_best_dist; }
56   void              GetBestParams(double &params[]);
57   void              GetSpectrum(double &alpha[], double &f_alpha[]);
58   int               GetSpectrumSize(void) { return m_n_spectrum; }
59   double            GetBestSSE(void)     { return m_sse[(int)m_best_dist]; }
60   bool              IsFitted(void)       { return m_fitted; }
61  };

The private section splits into three logical groups. The input data stores the τ(q) curve and Hurst exponent received from the partition analysis module.

The computed spectrum holds the (α, f(α)) pairs produced by the Legendre transform. mqh.

The two-dimensional m_params[4][4] array allows up to four parameters per distribution (though our models use at most two). The public interface follows a two-step pattern: call Init() and ComputeSpectrum() to produce the empirical spectrum, then FitAllDistributions() to find the best model.

Result accessors let the caller retrieve the winning distribution, its parameters, and the raw spectrum data.

After initialization with the τ(q) curve from Part 4, calling ComputeSpectrum() followed by FitAllDistributions() runs the entire process and selects the best model. Here is the orchestration method:

sql

1//+------------------------------------------------------------------+
2//| FitAllDistributions - Fit all 4 models and select best           |
3//+------------------------------------------------------------------+
4bool CSpectrumFitter::FitAllDistributions(void)
5  {
6   if(m_n_spectrum < 3)
7     {
8      Print("SpectrumFitter::FitAllDistributions - spectrum not computed");
9      return false;
10     }
11
12   Print("Fitting spectrum to 4 distributions...");
13
14   ArrayInitialize(m_sse, DBL_MAX);
15   int n_success = 0;
16
17   double sse_n = FitNormal();
18   if(sse_n < DBL_MAX)
19      n_success++;
20
21   double sse_b = FitBinomial();
22   if(sse_b < DBL_MAX)
23      n_success++;
24
25   double sse_p = FitPoisson();
26   if(sse_p < DBL_MAX)
27      n_success++;
28
29   double sse_g = FitGamma();
30   if(sse_g < DBL_MAX)
31      n_success++;
32
33   if(n_success == 0)
34     {
35      Print("SpectrumFitter::FitAllDistributions - all fits failed");
36      return false;
37     }
38
39//--- Select best fit (lowest SSE)
40   m_best_dist = MMAR_DIST_NORMAL;
41   double best_sse = m_sse[MMAR_DIST_NORMAL];
42
43   for(int d = 1; d < 4; d++)
44     {
45      if(m_sse[d] < best_sse)
46        {
47         best_sse = m_sse[d];
48         m_best_dist = (ENUM_MMAR_DISTRIBUTION)d;
49        }
50     }
51
52   string dist_names[] = {"Normal", "Binomial", "Poisson", "Gamma"};
53   double rmse = MathSqrt(best_sse / m_n_spectrum);
54   PrintFormat("BEST FIT: %s (SSE=%.6f, RMSE=%.6f)", dist_names[(int)m_best_dist], best_sse, rmse);
55
56   m_fitted = true;
57   return true;
58  }

The selection criterion is simple: lowest SSE wins. We do not penalize for the number of parameters (as AIC or BIC would) because the parameter counts are small (1 or 2) and the differences in fit quality are typically large enough that parsimony corrections would not change the outcome. The RMSE (root mean squared error) is printed alongside the SSE for interpretability, since it gives the average error in f(α) units.

The results are accessible through getter methods. GetBestDistribution() returns the winning enum value, GetBestParams() fills an array with the distribution parameters, and GetSpectrum() provides the empirical (α, f(α)) pairs for external visualization or validation.

Running the Full Pipeline

We can now write a test script that chains the partition analysis from Part 4 with the new spectrum fitting module. This script loads market data, runs partition analysis to extract τ(q), then passes the result to the spectrum fitter:

cpp

1//+------------------------------------------------------------------+
2//|                                               MMAR Spectrum Test |
3//|                      Copyright 2026, MMQ - Muhammad Minhas Qamar |
4//|                           https://www.mql5.com/en/articles/22519 |
5//+------------------------------------------------------------------+
6#property copyright "Copyright 2026, MMQ - Muhammad Minhas Qamar"
7#property link      "https://www.mql5.com/en/articles/22519"
8#property version   "1.00"
9#property strict
10#property script_show_inputs
11
12#include <MMAR\PartitionAnalysis.mqh>
13#include <MMAR\SpectrumFitter.mqh>
14
15input int    InpBars       = 50000;
16input double InpQMin       = 0.01;
17input double InpQMax       = 30.0;
18input double InpQStep      = 0.5;
19input int    InpDtMin      = 1;
20input int    InpDtMax      = 150;
21input double InpDtSpacing  = 1.1;
22input double InpMinR2      = 0.60;
23
24//+------------------------------------------------------------------+
25//| Script program start function                                    |
26//+------------------------------------------------------------------+
27void OnStart()
28  {
29   PrintFormat("=== MMAR Spectrum Fitting Test ===");
30   PrintFormat("Symbol: %s | Period: %s", _Symbol, EnumToString(_Period));
31
32//--- Load close prices and compute log returns
33   double close[];
34   int copied = CopyClose(_Symbol, _Period, 0, InpBars, close);
35   if(copied < 1000)
36     {
37      PrintFormat("Not enough bars: %d (need 1000+)", copied);
38      return;
39     }
40   PrintFormat("Loaded %d bars", copied);
41
42   double returns[];
43   int n_returns = copied - 1;
44   ArrayResize(returns, n_returns);
45   for(int i = 0; i < n_returns; i++)
46      returns[i] = MathLog(close[i + 1] / close[i]);
47
48//--- Phase 2: Partition Analysis
49   CPartitionAnalysis partition;
50   partition.SetPartitionConfig(InpDtMin, InpDtMax, InpDtSpacing);
51   partition.SetMomentConfig(InpQMin, InpQMax, InpQStep);
52   partition.SetMinRSquared(InpMinR2);
53
54   if(!partition.Init(returns, n_returns))
55     {
56      Print("Partition Init failed");
57      return;
58     }
59
60   uint t0 = GetTickCount();
61   if(!partition.Analyze())
62     {
63      Print("Partition Analyze failed");
64      return;
65     }
66   uint t_partition = GetTickCount() - t0;
67
68   double H = partition.GetH();
69   PrintFormat("\n--- Partition Results (%.0f ms) ---", (double)t_partition);
70   PrintFormat("H = %.4f | Vol = %.6f | R2 = %.4f | Score = %d/9",
71               H, partition.GetSampleVolatility(),
72               partition.GetMeanRSquared(), partition.GetMultifractalScore());
73
74//--- Phase 3: Spectrum Fitting
75   double tau[], q_vals[];
76   partition.GetTauQ(tau, q_vals);
77   int n_q = partition.GetNumQ();
78
79   CSpectrumFitter fitter;
80   if(!fitter.Init(q_vals, tau, n_q, H))
81     {
82      Print("SpectrumFitter Init failed");
83      return;
84     }
85
86   t0 = GetTickCount();
87   if(!fitter.ComputeSpectrum())
88     {
89      Print("ComputeSpectrum failed");
90      return;
91     }
92
93   if(!fitter.FitAllDistributions())
94     {
95      Print("FitAllDistributions failed");
96      return;
97     }
98   uint t_spectrum = GetTickCount() - t0;
99
100//--- Print spectrum results
101   PrintFormat("\n--- Spectrum Results (%.0f ms) ---", (double)t_spectrum);
102
103   string dist_names[] = {"Normal", "Binomial", "Poisson", "Gamma"};
104   ENUM_MMAR_DISTRIBUTION best = fitter.GetBestDistribution();
105   PrintFormat("Best distribution: %s", dist_names[(int)best]);
106
107   double params[];
108   fitter.GetBestParams(params);
109   PrintFormat("Parameters:");
110   switch(best)
111     {
112      case MMAR_DIST_NORMAL:
113         PrintFormat("  alpha_0 = %.6f", params[0]);
114         break;
115      case MMAR_DIST_BINOMIAL:
116         PrintFormat("  alpha_min = %.6f", params[0]);
117         PrintFormat("  alpha_max = %.6f", params[1]);
118         break;
119      case MMAR_DIST_POISSON:
120         PrintFormat("  alpha_0 = %.6f", params[0]);
121         break;
122      case MMAR_DIST_GAMMA:
123         PrintFormat("  alpha_0 = %.6f", params[0]);
124         PrintFormat("  gamma   = %.6f", params[1]);
125         break;
126     }
127
128   double best_sse = fitter.GetBestSSE();
129   double rmse = MathSqrt(best_sse / fitter.GetSpectrumSize());
130   PrintFormat("SSE = %.6f | RMSE = %.6f", best_sse, rmse);
131
132//--- Print empirical spectrum (every 5th point)
133   double alpha[], f_alpha[];
134   fitter.GetSpectrum(alpha, f_alpha);
135   int n_spec = fitter.GetSpectrumSize();
136
137   PrintFormat("\n--- Empirical Spectrum (every 5th point) ---");
138   PrintFormat("%10s %10s", "alpha", "f(alpha)");
139   for(int i = 0; i < n_spec; i += 5)
140      PrintFormat("%10.4f %10.4f", alpha[i], f_alpha[i]);
141
142   PrintFormat("\n--- Summary ---");
143   PrintFormat("Partition time: %u ms", t_partition);
144   PrintFormat("Spectrum time:  %u ms", t_spectrum);
145   PrintFormat("Total time:     %u ms", t_partition + t_spectrum);
146
147   Print("\n=== Spectrum Fitting Test Complete ===");
148  }

Running this on EURUSD M10 with 50,000 bars produces:

text

1=== MMAR Spectrum Fitting Test ===
2Symbol: EURUSD | Period: PERIOD_M10
3Loaded 50000 bars
4Multifractality check: concavity=0 linearity=1 variation=1 total=2/9
5
6--- Partition Results (438 ms) ---
7H = 0.4819 | Vol = 0.000419 | R2 = 0.8712 | Score = 2/9
8SpectrumFitter: 58 spectrum points, alpha=[0.1308, 0.5148], peak at alpha_0=0.5148 f=0.9924
9Fitting spectrum to 4 distributions...
10  Normal:   SSE=76.783833  alpha_0=0.514816
11  Binomial: SSE=0.148225  alpha_min=0.183507  alpha_max=0.772224
12  Poisson:  SSE=8.972094  alpha_0=0.514816
13  Gamma:    SSE=0.192775  alpha_0=0.514816  gamma=1.195666
14BEST FIT: Binomial (SSE=0.148225, RMSE=0.050553)
15
16--- Spectrum Results (0 ms) ---
17Best distribution: Binomial
18Parameters:
19  alpha_min = 0.183507
20  alpha_max = 0.772224
21SSE = 0.148225 | RMSE = 0.050553
22
23--- Empirical Spectrum (every 5th point) ---
24     alpha   f(alpha)
25    0.5148     0.9924
26    0.2886     0.5433
27    0.1385    -0.0312
28    0.1317    -0.0696
29    0.1360    -0.0303
30    0.1386     0.0005
31    0.1399     0.0185
32    0.1405     0.0282
33    0.1407     0.0330
34    0.1408     0.0349
35    0.1408     0.0353
36    0.1408     0.0347
37
38--- Summary ---
39Partition time: 438 ms
40Spectrum time:  0 ms
41Total time:     438 ms
42
43=== Spectrum Fitting Test Complete ===

The results are instructive on multiple levels. 78.

and a flat right tail. The Binomial's entropy-based formula handles this shape naturally.

Second, the spectrum peak at α₀ = 0.5148 with f(α₀) = 0.9924 is very close to the theoretical ideal of f(α₀) = 1.0. This confirms that the Legendre transform is working correctly, because the dominant scaling exponent has a fractal dimension near 1.0, meaning it characterizes nearly the entire time series.

define the cascade's behavior. In Part 6, when we build the multiplicative cascade, each split will assign multipliers derived from one of these two exponents.

59 units) indicates substantial heterogeneity in the scaling, meaning some intervals will be very smooth (high α) while others are very rough (low α). This heterogeneity is exactly what produces realistic volatility clustering in the final MMAR process.

Finally, notice that the spectrum fitting itself takes essentially zero milliseconds. All the computational cost is in the partition analysis (438 ms).

The BLEIC optimization converges quickly because the SSE surfaces are smooth and the parameter spaces are low-dimensional (1 or 2 parameters). This means adding spectrum fitting to a live system incurs negligible additional cost beyond the partition analysis we already computed in Part 4.

Conclusion

In this article, we built the second analytical module of the MMAR library: the CSpectrumFitter class. Starting from the τ(q) scaling function produced by the partition analysis engine in Part 4, we applied the Legendre transform to obtain the empirical singularity spectrum f(α), then fitted four theoretical distribution models using bounded nonlinear optimization.

The Legendre transform. Converts the scaling function into a physically interpretable representation, where α gives the local smoothness and f(α) gives the fractal dimension of points sharing that smoothness.
Four candidate models. (Normal, Binomial, Poisson, Gamma) compete based on SSE. Each encodes different assumptions about the cascade multiplier distribution.
BLEIC optimization. From ALGLIB handles the parameter constraints and converges in negligible time for these low-dimensional problems.
The Binomial model. Won on our EURUSD M10 data, providing α_min and α_max parameters that define the cascade's heterogeneity.

The two analytical stages are now complete. We can take any return series, determine whether it is multifractal, and if so, extract the exact distribution and parameters needed for simulation. In Part 6, we will put these parameters to work: building the multiplicative cascade that warps uniform time into multifractal trading time, generating Fractional Brownian Motion via Davies–Harte and Cholesky methods, and composing them into the full MMAR process X(t) = B H [θ(t)].

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 spectrum 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 below. The full repository is also available on MQL5 Algo Forge, a 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 (Part 4)
MQL5\Include\MMAR\OLS.mqh	Ordinary Least Squares regression via SVD pseudo-inverse with full diagnostics (Part 4)
MQL5\Include\MMAR\HurstGHE.mqh	Generalized Hurst Exponent estimator, fallback method (Part 4)
MQL5\Include\MMAR\PartitionAnalysis.mqh	Multifractal partition function analysis with scaling diagnostics (Part 4)
MQL5\Include\MMAR\SpectrumFitter.mqh	Multifractal spectrum computation and four-distribution fitting via BLEIC optimization (Part 5)
MQL5\Scripts\MMAR\MMAR_PartitionTest.mq5	Test script for partition analysis on live market data (Part 4)
MQL5\Scripts\MMAR\MMAR_SpectrumTest.mq5	Test script for the full partition + spectrum pipeline on live market data (Part 5)

Attached files |

Download ZIP

MQL5.zip(84.71 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

Developer by Profession, Trader by Hobby

Gmail: ayanminhasshayar@gmail.com

Beyond GARCH (Part V): Fitting the Multifractal Spectrum in MQL5 - MQL5 Articles

Introduction

From τ(q) to f(α): The Legendre Transform

Four Theoretical Spectrum Models

Bounded Optimization with BLEIC

The CSpectrumFitter Class

Running the Full Pipeline

Conclusion

Other articles by this author