Beyond GARCH (Part II): Measuring the Fractal Dimension of Markets

Introduction

In Part 1 of this series, we introduced the theoretical foundations of GARCH and the Multifractal Model of Asset Returns (MMAR), loaded EURUSD 5-minute data via the MetaTrader 5 Python API, and ran the partition function analysis to test for multifractal scaling. The result was clear: our data exhibits power-law moment scaling across time intervals, confirming that a fractal structure is present.

But confirming that fractality exists is only the first step. We now need to characterize it. In this article, we extract the scaling function tau(q) from the partition-function slopes in Part 1. We then estimate the Hurst exponent H via R/S analysis, run a three-component multifractality decision test, and transform tau(q) into the multifractal spectrum f(alpha) using the Legendre transform. Finally, we fit the empirical spectrum to four theoretical distributions (Lognormal, Binomial, Poisson, Gamma) and select the multiplicative-cascade type that best matches the observed structure. By the end of this article, we will have the three ingredients needed to build the MMAR process: the Hurst exponent, the winning distribution type, and its fitted parameters.

Before the implementation, we clarify two key concepts: the Hurst exponent and the multifractal spectrum. Both were introduced briefly in Part 1 as components of the MMAR theory, but they are central enough to Part 2's work that understanding them deeply, not just knowing what they are but knowing what they mean, will make every line of code that follows more intelligible.

We will cover:

1. The Hurst Exponent — Memory in Markets
2. Step 2: Extracting the Scaling Function and Hurst Exponent
3. The Multifractal Spectrum — Reading f(alpha)
4. Step 3: Fitting the Multifractal Spectrum
5. Conclusion

The Hurst Exponent — Memory in Markets

In Part 1, we introduced the Hurst exponent H as the parameter that controls the memory structure of Fractional Brownian Motion: H = 0.5 gives a random walk, H > 0.5 gives persistence, H < 0.5 gives anti-persistence. That was enough to understand where H fits in the MMAR equation. Because we now estimate H from real data and use it as a construction parameter, we need a deeper understanding of what H represents, where it came from, and what it costs to get it wrong.

Origins: Hurst and the Nile

The Hurst exponent is named after Harold Edwin Hurst, a British hydrologist who spent decades studying the Nile's hydrology. In the early 1950s, Hurst was working on a practical engineering problem: how large should a reservoir be to guarantee a steady water supply despite the Nile's unpredictable annual flooding? The conventional statistical approach assumed that each year's flood level was independent of previous years, a random walk. Under that assumption, the required reservoir capacity could be computed using standard probability theory.

But Hurst noticed something that the standard theory could not explain. When he examined centuries of Nile flood records (some dating back to the 7th century via nilometer measurements at Roda Island), he found that the range of cumulative deviations from the long-term mean grew faster than the square root of time. For independent observations, the range R divided by the standard deviation S should scale as n^0.5, where n is the number of observations. Hurst found that R/S scaled as n^H where H was approximately 0.73, substantially larger than 0.5. This meant that wet years tended to cluster with wet years and dry years with dry years, far more than random chance would predict. The river had long memory.

Hurst then tested his method on dozens of other geophysical time series, including tree rings, sunspot counts, lake levels, and rainfall records. In nearly every case, he found H values between 0.6 and 0.9. The phenomenon was so pervasive that it became known as the Hurst phenomenon, and the exponent that measured it took his name. Mandelbrot later provided the theoretical foundation by showing that Fractional Brownian Motion with parameter H naturally produces exactly this kind of long-range dependence.

What H Means Geometrically

The Hurst exponent controls the roughness of a path. Consider three Fractional Brownian Motion paths, each 1,000 steps long but with different H values. At H = 0.3 (anti-persistent), the path is extremely jagged. Every upward move is likely followed by a downward move, creating a rapidly oscillating, rough trajectory that constantly reverses direction. At H = 0.5 (the standard random walk), the path has moderate roughness, the familiar irregular wandering of a coin-flip process. At H = 0.7 (persistent), the path is noticeably smoother. Upward moves tend to continue upward, producing longer sweeps and gentler undulations. The path still wanders randomly, but it does so in broader, more coherent strokes.

Mathematically, this roughness is captured by the Hölder regularity of the path. An FBM path with Hurst exponent H is almost surely Hölder continuous with exponent H, meaning that the increments over a time step delta_t scale as delta_t^H. Larger H means increments grow faster with the time step, which means the path covers more ground in smooth sweeps. Smaller H means increments grow more slowly, which means the path changes direction frequently and the total variation is higher. This geometric interpretation matters for trading: a persistent market (H > 0.5) rewards trend-following strategies, while an anti-persistent market (H < 0.5) rewards mean-reversion strategies. A random walk (H = 0.5) rewards neither, and any apparent pattern is noise.

Fig. 1 — FBM paths for H=0.3 (anti-persistent), H=0.5 (random walk), and H=0.7 (persistent). The roughness of the path decreases as H increases.

The H = 0.5 Assumption and What It Costs You

The entire edifice of classical quantitative finance rests on the assumption that H = 0.5. The Black-Scholes option pricing model, the Efficient Market Hypothesis in its strong form, and standard Value-at-Risk calculations all assume that returns are independent and identically distributed, which is equivalent to assuming standard Brownian motion, which is FBM with H = 0.5. Under this assumption, price changes carry no memory. Yesterday's volatility tells you nothing about tomorrow's. A flash crash is just as likely after a calm week as after a turbulent one.

When H deviates from 0.5, these models break down in specific, quantifiable ways. If H > 0.5 (persistence), the variance of cumulative returns over a horizon T scales as T^(2H) instead of T. For H = 0.55 (a modest deviation), the variance over 100 periods is 100^1.1 = 158, compared to 100^1.0 = 100 under the random walk assumption. That is a 58% underestimation of risk over just 100 periods. Over longer horizons the gap widens further, because the power-law exponent compounds. This means that risk models assuming H = 0.5 systematically underestimate the probability of large drawdowns over medium-to-long horizons, precisely where risk management matters most.

Conversely, if H < 0.5 (anti-persistence), the variance grows slower than T, and models assuming H = 0.5 overestimate risk, leading to excessively conservative position sizing that leaves returns on the table. Either way, getting H wrong means your risk model is miscalibrated, and the error grows with the forecast horizon.

Local versus Global H

There is a subtlety that connects the Hurst exponent directly to the multifractal analysis we are building. When we estimate H from a time series, we get a single number, a global Hurst exponent that represents the average scaling behavior across the entire dataset. But a multifractal process does not have a single scaling exponent. Different regions of the time series have different local regularity: some segments are smooth (locally high H), others are rough (locally low H). The global H is a weighted average of these local exponents, useful but incomplete.

This motivates the multifractal spectrum f(alpha). Where the global H gives you one number, the spectrum gives you the entire distribution of local scaling behaviors. A process with a single H value at every point is monofractal; its spectrum collapses to a single point. A process where H varies from point to point is multifractal; its spectrum is a wide curve. The width of that curve measures how much the local scaling varies, which is the degree of multifractality. This is the bridge between the Hurst exponent (Section 2) and the multifractal spectrum (Section 4): the global H tells you the average memory structure, and the spectrum tells you how much that memory structure varies across the data. Both are needed to build the MMAR.

Step 2: Extracting the Scaling Function and Hurst Exponent

With fractality confirmed and the Hurst exponent understood, we now extract two critical quantities: the scaling function tau(q) and H itself. These are the parameters that define the specific multifractal character of our data. The scaling function tells us whether our data is truly multifractal (as opposed to monofractal), and the Hurst exponent determines the long-memory structure of the FBM component.

Extracting tau(q)

The scaling function comes directly from Step 1's regression slopes. Since the partition function obeys the power law S_q proportional to delta_t^(tau(q) + 1), and the slope of the log-log regression equals tau(q) + 1, the scaling function is simply:

tau(q) = slope - 1

The ScalingExtractor class takes the completed FractalityChecker from Step 1 and immediately computes tau(q) for all q values:

# step2_extract_scaling.py

class ScalingExtractor:
    def __init__(self, checker, verbose=config.VERBOSE):
        self.checker = checker
        self.verbose = verbose

        # Extract from checker
        self.q_values = np.array(sorted(checker.slopes.keys()))
        self.slopes = np.array([checker.slopes[q] for q in self.q_values])
        self.r_squared = np.array([checker.r_squared_values[q]
                                    for q in self.q_values])

        # Calculate tau(q) = slope - 1
        self.tau_q = self.slopes - 1

        self.H = None
        self.tau_q_func = None

The shape of tau(q) tells us what kind of fractal process we are dealing with. For a monofractal process (like a random walk with drift), tau(q) is perfectly linear: tau(q) = qH - 1. This means every moment of the distribution scales at the same rate. For a multifractal process, tau(q) is concave, it curves downward. This means different moments scale at different rates, which is the defining characteristic of multifractal processes. If our data is monofractal, the MMAR offers no advantage over simpler models.

The Multifractality Decision Test

This is one of the most important parts of the pipeline: a rigorous, three-component test that determines whether the MMAR is worth applying. Our implementation runs three independent tests and produces a composite score from 0 to 9:

def check_multifractality(self):
    # Test 1: Concavity (second derivative of tau(q))
    d2_tau = np.diff(self.tau_q, 2)
    concave_pct = 100 * np.sum(d2_tau < 0) / len(d2_tau)
    mean_d2 = np.mean(d2_tau)

    if concave_pct >= 80:
        concavity_score = 3  # STRONGLY CONCAVE
    elif concave_pct >= 60:
        concavity_score = 2  # MODERATELY CONCAVE
    elif concave_pct >= 40:
        concavity_score = 1  # WEAKLY CONCAVE
    else:
        concavity_score = 0  # NOT CONCAVE

Test 1: Concavity. We compute the second differences (discrete approximation of the second derivative) of tau(q). For a multifractal process, these should be predominantly negative (concave). We measure what percentage of second differences are negative. If 80% or more are negative, the process is strongly concave (score 3). If less than 40% are negative, the process is not concave (score 0). This test directly measures the key theoretical prediction: multifractal tau(q) curves downward.

    # Test 2: Compare with theoretical random walk
    q_range = self.q_values[self.q_values > 0]
    tau_range = self.tau_q[self.q_values > 0]

    from scipy.stats import linregress
    slope_linear, intercept_linear, r_value_linear, _, _ = linregress(
        q_range, tau_range
    )
    r2_linear = r_value_linear ** 2

    tau_linear_fit = slope_linear * q_range + intercept_linear
    rmse_from_linear = np.sqrt(np.mean((tau_range - tau_linear_fit)**2))

    if r2_linear > 0.99 and rmse_from_linear < 0.05:
        linearity_score = 0  # HIGHLY LINEAR (like random walk)
    elif r2_linear > 0.95:
        linearity_score = 1  # MOSTLY LINEAR (weakly multifractal)
    elif r2_linear > 0.90:
        linearity_score = 2  # MODERATELY NONLINEAR
    else:
        linearity_score = 3  # STRONGLY NONLINEAR

Test 2: Linearity Deviation. We fit a straight line to tau(q) and measure how well it fits. For a random walk, tau(q) = 0.5q - 1, which is perfectly linear (R-squared = 1.0). The more our tau(q) deviates from linearity, the stronger the evidence for multifractality. An R-squared above 0.99 with RMSE below 0.05 means the data looks essentially like a random walk (score 0). R-squared below 0.90 means strong nonlinearity (score 3).

    # Test 3: Range of local slopes (should vary for multifractals)
    d_tau = np.diff(self.tau_q)
    slope_range = np.ptp(d_tau)  # Peak-to-peak
    slope_std = np.std(d_tau)

    if slope_range > 0.5:
        variation_score = 3  # HIGH VARIATION
    elif slope_range > 0.3:
        variation_score = 2  # MODERATE VARIATION
    elif slope_range > 0.1:
        variation_score = 1  # LOW VARIATION
    else:
        variation_score = 0  # MINIMAL VARIATION (monofractal)

    total_score = concavity_score + linearity_score + variation_score

Test 3: Slope Variation. We compute the first differences of tau(q) (the local slope) and measure their range. For a monofractal process, all local slopes are identical (the function is linear), so the range is zero. For a multifractal, the local slope varies as q changes, producing a nonzero range. A range above 0.5 indicates high variation (score 3), while a range below 0.1 indicates near-monofractal behavior (score 0).

The composite score is interpreted as follows:

• Score 7-9: Strongly multifractal. MMAR highly recommended, proceed with confidence.
• Score 5-6: Moderately multifractal. MMAR worth trying, may outperform GARCH.
• Score 3-4: Weakly multifractal. MMAR uncertain, GARCH may win.
• Score 0-2: Not multifractal. MMAR not recommended, data is too close to a random walk.

This acts as a strict applicability gate. The MMAR is not a universal tool, and applying it to data that is not multifractal would produce worse forecasts than GARCH while being more expensive to compute.

Estimating the Hurst Exponent H

The Hurst exponent H is estimated using R/S (rescaled range) analysis via the nolds library. This is the same technique Hurst himself developed for the Nile data, modernized with corrections for finite-sample bias. We use R/S rather than the alternative approach of numerically finding where tau(1/H) = 0, which can be unstable with noisy tau(q) estimates. The R/S method directly measures the scaling of the range of the cumulative deviation from the mean, relative to the standard deviation:

def estimate_H_nolds(self):
    if not NOLDS_AVAILABLE:
        return None

    returns = self.checker.returns

    # corrected=True applies Anis-Lloyd-Peters correction
    # fit='RANSAC' is robust to outliers in R/S statistics
    H_nolds = nolds.hurst_rs(returns, fit='RANSAC', corrected=True)
    return H_nolds

The corrected=True flag applies the Anis-Lloyd-Peters correction, which adjusts for finite-sample bias in the R/S statistic. Without this correction, short time series tend to overestimate H. The fit='RANSAC' flag uses Random Sample Consensus instead of ordinary least squares for the regression, making the estimation robust to outliers in the R/S statistics at individual subsample sizes.

If the nolds library is not available, a fallback method finds where tau(q) crosses zero. Theoretically, tau(1/H) = 0, so H = 1/q_zero where q_zero is the value of q at the crossing. The fallback uses Brent's root-finding method on a cubic spline interpolation of tau(q). If no sign change exists in the computed range (tau(q) does not cross zero), it uses the closest point as an approximation:

# Fallback: find where tau(q) = 0
self.tau_q_func = interp1d(self.q_values, self.tau_q,
                            kind='cubic',
                            fill_value='extrapolate',
                            bounds_error=False)

tau_min = self.tau_q_func(self.q_values.min())
tau_max = self.tau_q_func(self.q_values.max())

if tau_min * tau_max > 0:
    # No sign change - use closest point
    idx_closest = np.argmin(np.abs(self.tau_q))
    q_zero = self.q_values[idx_closest]
    self.H = 1.0 / q_zero
else:
    # Brent's method for root finding
    q_zero = brentq(self.tau_q_func, self.q_values.min(), self.q_values.max())
    self.H = 1.0 / q_zero

Fig. 2 — The scaling function tau(q) with the Hurst exponent marked

The Multifractal Spectrum — Reading f(alpha)

We are about to transform the scaling function tau(q) into the multifractal spectrum f(alpha) and fit it to theoretical distributions. But before we write that code, we need to understand what the spectrum actually tells us. The Legendre transform and the f(alpha) curve were mentioned in Part 1, but only as a mathematical step in the MMAR pipeline. Here we explain what the spectrum means: what its shape reveals about the data, why its width matters, and how to read it like a diagnostic tool.

Alpha: Local Regularity at a Point

At the heart of the multifractal spectrum is the Hölder exponent alpha, which measures how regular (smooth) a function is at a specific point. Consider a price series. At some moments, the price changes gradually and predictably, and the chart looks smooth. At other moments, the price jumps or spikes, and the chart looks rough, almost discontinuous. The Hölder exponent quantifies this: a large alpha means the price process is locally smooth at that point (small, gradual changes), while a small alpha means it is locally rough (sharp, abrupt changes).

More precisely, alpha at a point t_0 is defined by the condition that the price increment over a small interval delta_t around t_0 scales as delta_t^alpha. If alpha = 1, the price is locally Lipschitz continuous (as smooth as a differentiable function). If alpha = 0.5, the price increments scale like a Brownian motion. If alpha < 0.5, the price is rougher than Brownian motion at that point, meaning abrupt changes are happening at an unusually fast rate relative to the time scale. In a real financial time series, different points have different alpha values: quiet overnight sessions have high alpha (smooth), while news-driven minutes have low alpha (rough).

f(alpha): The Fractal Dimension of Regularity Sets

If alpha tells us the local regularity at a single point, f(alpha) tells us how prevalent that regularity level is across the entire time series. Specifically, f(alpha) is the Hausdorff dimension of the set of all points t where the local Hölder exponent equals alpha. Think of it this way: collect every point in the time series that has regularity alpha, and measure the "size" of that collection using fractal dimension. A Hausdorff dimension of 1 means the set fills a one-dimensional interval (these alpha values are everywhere). A dimension of 0 means the set is just a scattering of isolated points (these alpha values are extremely rare). Values between 0 and 1 indicate increasingly sparse sets.

This is why the peak of the spectrum always touches or approaches f = 1. The alpha value at the peak is the most common local regularity in the data, the one that fills a full one-dimensional interval of time. All other alpha values are less common (f < 1), and alpha values far from the peak are progressively rarer, forming thinner and thinner fractal dust on the time axis.

The Shape of the Spectrum: Monofractal versus Multifractal

The shape and width of the f(alpha) curve carry direct information about the complexity of the scaling structure. A monofractal process has a single scaling exponent everywhere. Every point in the time series has the same local Hölder exponent alpha = H. The spectrum degenerates to a single point at (H, 1): one alpha value, and the set of points with that alpha fills the entire time axis (f = 1). Standard Brownian motion is monofractal with alpha = H = 0.5. Simple FBM without multifractal trading time is also monofractal.

A multifractal process has a range of local scaling exponents. Different points in the time series have different alpha values, and the spectrum opens up into an inverted parabola (or more precisely, a concave curve). The width of this curve, measured as alpha_max - alpha_min, is a direct measure of multifractality strength. A narrow spectrum (width close to zero) means the data is nearly monofractal: all regions scale similarly, volatility is relatively homogeneous. A wide spectrum means the data contains both very smooth regions (high alpha) and very rough regions (low alpha), indicating strong heterogeneity in the scaling behavior, which is the hallmark of complex volatility clustering.

For financial data, typical spectrum widths range from 0.2 to 0.8. A width of 0.2 indicates mild multifractality (slightly heterogeneous volatility), while a width of 0.8 indicates strong multifractality (the data contains dramatically different local behaviors at different times). FX markets like EURUSD tend to have moderate widths (0.3 to 0.5), reflecting the mixture of quiet Asian sessions, volatile London opens, and news-driven spikes that characterize currency trading.

Asymmetry: Left-Skew versus Right-Skew

The spectrum is not always symmetric. A left-skewed spectrum (the left tail, toward small alpha, is longer than the right tail) means that the extreme rough events are more diverse than the smooth ones. In market terms, the data has a rich variety of different crash intensities but relatively uniform calm behavior. A right-skewed spectrum (longer right tail) means the opposite: smooth periods come in many flavors but extreme events are more uniform. Most financial time series exhibit slight left-skewness, reflecting the empirical observation that volatility bursts come in a wider range of intensities than quiet periods.

The asymmetry can also indicate which regime dominates the data's multifractal structure. If the left branch of the spectrum (governed by large positive q values in the Legendre transform) is more developed, it means the scaling behavior is more varied during high-moment events (extreme returns). If the right branch (governed by small q values near zero) is more developed, the scaling diversity lives in the small fluctuations.

Fig. 3 — Monofractal spectrum (single point) versus multifractal spectrum (wide inverted parabola). The width measures multifractality strength.

This interpretation is used in Section 4 to fit the empirical spectrum to theoretical cascades. We are not just curve-fitting for its own sake. We are determining which theoretical cascade model produces the specific shape of f(alpha) that our data exhibits, its width, its peak location, its asymmetry, because that shape is the fingerprint of how volatility organizes itself across scales.

Step 3: Fitting the Multifractal Spectrum

Now we need to characterize what type of multifractal our data exhibits. Not all multifractals are the same. They differ in how their local singularities are distributed, and this distribution determines how we sample multipliers in the cascade construction (Step 4). We transform the scaling function tau(q) into the multifractal spectrum f(alpha) using the Legendre transform, then fit the empirical spectrum to four theoretical distributions.

The Legendre Transform: From tau(q) to f(alpha)

The Legendre transform converts from q-space (moment orders) to alpha-space (singularity exponents). The transformation is:

alpha(q) = d(tau)/dq

f(alpha) = alpha * q - tau(q)

Here, alpha represents the local Hölder exponent we described in Section 3: a measure of local regularity at a point, where small alpha means rough behavior and large alpha means smooth behavior. The function f(alpha) is the Hausdorff dimension of the set of points sharing that regularity level. Differentiating tau(q) with respect to q yields the local scaling exponent alpha, indicating which singularity strength dominates at each moment order. The relation f(alpha) = alpha * q - tau(q) then gives the dimension of the corresponding singularity set.

The implementation computes alpha numerically using central differences of the actual neighboring tau(q) values:

# step3_fit_spectrum.py

def compute_multifractal_spectrum(self):
    alpha_values = []
    f_alpha_values = []
    q_for_spectrum = []

    q_values = self.extractor.q_values
    tau_q = self.extractor.tau_q

    for i, q in enumerate(q_values):
        if q < 0.5:  # Skip very small q (numerical issues)
            continue
        if i == 0 or i == len(q_values) - 1:  # Skip edge points
            continue

        # Central difference: d(tau)/dq
        q_prev = q_values[i-1]
        q_next = q_values[i+1]
        tau_prev = tau_q[i-1]
        tau_next = tau_q[i+1]

        alpha = (tau_next - tau_prev) / (q_next - q_prev)
        f_alpha = alpha * q - tau_q[i]

        if np.isfinite(alpha) and np.isfinite(f_alpha):
            alpha_values.append(alpha)
            f_alpha_values.append(f_alpha)

    self.alpha_values = np.array(alpha_values)
    self.f_alpha_values = np.array(f_alpha_values)

Note two important details. First, we skip q values below 0.5 because the derivative of tau(q) near q = 0 is numerically unstable (small denominators amplify noise). Second, we skip edge points where we cannot compute a central difference. These precautions prevent NaN and Inf values from contaminating the spectrum.

Fitting to Four Theoretical Distributions

Following Zhang's Appendix B, we fit the empirical spectrum to four theoretical distributions, each corresponding to a different type of multiplicative cascade. Each distribution has a closed-form expression for f(alpha), and we select the best fit by minimizing the sum of squared errors. Let us examine each distribution and understand what it means for the market's multifractal structure.

1. Lognormal Distribution:

def normal_spectrum(self, alpha, alpha_0, H):
    denominator = 4 * H * (alpha_0 - H)
    if denominator <= 0:
        return np.full_like(alpha, -np.inf)
    return 1 - (alpha - alpha_0)**2 / denominator

This produces a parabolic spectrum centered at alpha_0, the most common singularity exponent. The single parameter alpha_0 controls both the center and width of the parabola. The constraint alpha_0 > H is required to keep the denominator positive; if alpha_0 were equal to or less than H, the spectrum would degenerate. A lognormal cascade means the multiplier magnitudes are drawn from a lognormal distribution, implying that trading activity is the product of many independent random factors at each cascade level, a natural consequence of the central limit theorem applied to log-multipliers.

2. Binomial Distribution:

def binomial_spectrum(self, alpha, alpha_min, alpha_max):
    if alpha_max <= alpha_min:
        return np.full_like(alpha, -np.inf)

    alpha = np.clip(alpha, alpha_min + 1e-10, alpha_max - 1e-10)

    term1_ratio = (alpha_max - alpha) / (alpha_max - alpha_min)
    term2_ratio = (alpha - alpha_min) / (alpha_max - alpha_min)

    term1_ratio = np.maximum(term1_ratio, 1e-10)
    term2_ratio = np.maximum(term2_ratio, 1e-10)

    f = -(term1_ratio * np.log2(term1_ratio)
          + term2_ratio * np.log2(term2_ratio))
    return f

The binomial spectrum has two parameters: alpha_min and alpha_max. It represents the simplest possible cascade, where at each level the multiplier takes one of exactly two values with equal probability. This produces a symmetric spectrum bounded between alpha_min and alpha_max. The np.clip and epsilon guards prevent log(0) errors at the boundaries. A binomial cascade implies that market activity switches between exactly two regimes (high and low), a crude but sometimes effective model.

3. Poisson Distribution:

def poisson_spectrum(self, alpha, alpha_0, H, b=2):
    if alpha_0 <= 0:
        return np.full_like(alpha, -np.inf)

    alpha = np.maximum(alpha, 1e-10)

    term1 = 1 - alpha_0 / (H * np.log(b))
    term2 = (alpha / H) * (np.log(alpha_0 * np.e / alpha) / np.log(b))
    return term1 + term2

The Poisson spectrum has one parameter (alpha_0) and produces an asymmetric spectrum with a sharp cutoff on one side. This corresponds to a cascade where extreme multipliers are rare but possible, following a Poisson process. This can model markets where extreme volatility events are discrete, identifiable shocks (earnings announcements, policy changes) rather than continuous fluctuations.

4. Gamma Distribution:

def gamma_spectrum(self, alpha, alpha_0, gamma, b=2):
    if alpha_0 <= 0 or gamma <= 0:
        return np.full_like(alpha, -np.inf)

    alpha = np.maximum(alpha, 1e-10)

    term1 = gamma * np.log(alpha / alpha_0) / np.log(b)
    term2 = gamma * (alpha_0 - alpha) / (alpha_0 * np.log(b))
    return 1 + term1 + term2

The Gamma spectrum has two parameters (alpha_0 and gamma) and is the most flexible of the four. The gamma parameter controls the shape of the spectrum independently of the center, allowing both symmetric and asymmetric forms. This flexibility means the Gamma distribution often provides the best fit, but with two parameters it requires more data to fit reliably.

Distribution Selection

All four distributions are fitted using L-BFGS-B optimization with parameter bounds, minimizing the sum of squared errors between the theoretical and empirical spectra. The distribution with the lowest SSE wins:

def fit_all_distributions(self):
    distributions = {
        'Normal': self.fit_normal_distribution,
        'Binomial': self.fit_binomial_distribution,
        'Poisson': self.fit_poisson_distribution,
        'Gamma': self.fit_gamma_distribution
    }

    for name, fit_func in distributions.items():
        try:
            params, sse = fit_func()
            self.fitted_distributions[name] = {
                'params': params,
                'sse': sse,
                'rmse': np.sqrt(sse / len(self.alpha_values))
            }
        except Exception as e:
            self.fitted_distributions[name] = None

    valid_dists = {k: v for k, v in self.fitted_distributions.items()
                   if v is not None}
    self.best_distribution = min(
        valid_dists.keys(), key=lambda k: valid_dists[k]['sse']
    )
    self.best_params = valid_dists[self.best_distribution]['params']

Each fit function is wrapped in a try/except block. If a particular distribution fails to fit (for example, the Lognormal cannot fit when alpha_0 falls below H, making the denominator negative), it is simply excluded from the comparison. This robustness ensures the pipeline does not crash due to a single poorly-conditioned distribution.

The winning distribution and its parameters are stored in the SpectrumFitter object and passed to Step 4, where they determine how cascade multipliers are sampled.

Fig. 4 — Empirical multifractal spectrum f(alpha) with all four theoretical distribution fits overlaid

Fig. 5 — Best-fit distribution selected by minimum sum of squared errors

Conclusion

In this article, we took the partition function results from Part 1 and extracted the quantities that define our data's specific multifractal character. Here is what we established:

• The Hurst exponent is more than a parameter. It encodes the memory structure of the market, from Hurst's original Nile river studies to its role as the exponent that controls path roughness in FBM. The H = 0.5 assumption embedded in classical finance (Black-Scholes, EMH) systematically miscalibrates risk when the true H differs, and the error grows with the forecast horizon.
• Multifractality has degrees. The three-component scoring system in Step 2 (concavity, linearity deviation, slope variation) does not just give a binary yes/no answer. It quantifies how strongly multifractal the data is, allowing practitioners to make informed decisions about model selection.
• The spectrum is a fingerprint, not just a curve. The multifractal spectrum f(alpha) reveals the full distribution of local scaling behaviors: its width measures multifractality strength, its peak locates the dominant regularity, and its asymmetry tells you whether the diversity lives in the extreme events or the quiet periods. Fitting it to four theoretical distributions (Lognormal, Binomial, Poisson, Gamma) determines which cascade type will best reproduce this fingerprint.
• We now have everything we need. The Hurst exponent H, the best-fit distribution type, and its parameters are the three ingredients required to build the MMAR process. These are not arbitrary choices; they are derived directly from the data through the analytical pipeline we built across Parts 1 and 2.

The analytical phase spanning Parts 1 and 2 is not optional overhead. It serves as a strict applicability gate: it tells you whether the MMAR is worth applying, and if so, provides the exact parameters for construction, before you invest the computational effort of building it.

In Part 3 of this series, we will put these fitted parameters to work. We will construct the multiplicative cascade to build multifractal trading time, generate Fractional Brownian Motion via the Davies-Harte algorithm, combine them into the MMAR compound process, and run 1,000 Monte Carlo simulations to produce a volatility forecast. We will then implement a GARCH(1,1) on the same training data and compare both models head-to-head against realized volatility. The question we set up here, whether fractals beat econometrics at forecasting volatility, will finally get its answer.

Disclaimer: This article is for educational and research purposes only. The models and code presented here are based on academic research and should not be interpreted as financial advice. Volatility forecasting carries inherent uncertainty. Always validate models on out-of-sample data and test on demo accounts before risking live funds.

Getting the Source Code via MQL5 Algo Forge

All source files for this article 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
config.py	Configuration parameters, helper functions, and validation for the entire pipeline
data_loader.py	MetaTrader 5 data loading, CSV fallback, return calculation, and data quality checks
utils.py	Shared utility functions: rolling windows, volatility, autocorrelation, Timer, JSON export
step1_check_fractality.py	Partition function analysis with non-overlapping intervals, R-squared assessment, crossover detection
step2_extract_scaling.py	Scaling function extraction, Hurst estimation (nolds R/S + Brent fallback), three-component multifractality test
step3_fit_spectrum.py	Legendre transform, four distribution fits (Lognormal, Binomial, Poisson, Gamma), L-BFGS-B optimization
run_step1.py	Runner: loads data, runs fractality check, saves FractalityChecker for Step 2
run_step2.py	Runner: loads Step 1 output, extracts scaling function and Hurst exponent, runs multifractality test
run_step3.py	Runner: loads Step 2 output, computes multifractal spectrum, fits four distributions
requirements.txt	Python dependencies

Further Reading:

• A Multifractal Model of Asset Returns — Benoit Mandelbrot, Adlai Fisher, Laurent Calvet, 1997
• Generalized autoregressive conditional heteroskedasticity — Tim Bollerslev, 1986
• Volatility Forecasting with the Multifractal Model of Asset Returns — Terrence Y. Zhang, 2017
• Dynamic Conditional Correlation - a Simple Class of Multivariate GARCH Models — Robert F. Engle, 2000