Beyond GARCH (Part III): Building the MMAR and the Verdict

Introduction

In Part 1 of this series, we covered the theoretical foundations of GARCH and the MMAR, loaded EURUSD 5-minute data, and confirmed via partition function analysis that the data exhibits multifractal scaling. In Part 2, we extracted the scaling function tau(q), estimated the Hurst exponent H, ran a rigorous multifractality decision test, and fitted the multifractal spectrum to four theoretical distributions. The result: our data is genuinely multifractal, and we now have the specific parameters, the Hurst exponent, the winning distribution type, and its fitted coefficients, needed to build the full MMAR process.

In this article, we put those parameters to work. We construct two MMAR components: a multiplicative cascade that creates multifractal trading time (the "deformed clock" that warps uniform time into market time) and Fractional Brownian Motion that provides long-memory structure. We combine them into the compound process X(t) = B_H[theta(t)]. Next, we run 1,000 Monte Carlo simulations to forecast volatility. We then fit a GARCH(1,1) model on the same training data. Finally, we compare both models head-to-head against realized volatility. No tricks, no cherry-picking, just a fair fight.

We will cover:

1. Step 4: Building Multifractal Trading Time
2. Step 5: Generating Fractional Brownian Motion
3. Step 6: The MMAR Equation Comes Alive
4. Step 7: Monte Carlo Forecasting
5. The GARCH Implementation
6. The Showdown: Results
7. Conclusion

Step 4: Building Multifractal Trading Time

With the multifractal spectrum fitted and a distribution type selected, we now construct the first of the two MMAR components: the multifractal trading time theta(t). This is the "deformed clock" that will warp the FBM to produce volatility clustering and fat tails.

The Multiplicative Cascade Algorithm

The multiplicative cascade is an iterative process that creates a multifractal measure, a distribution of "mass" across intervals that exhibits self-similar clustering. The algorithm works as follows:

1. Initialization: Start with a single interval [0, 1] carrying mass = 1.
2. For each cascade level (k = 1 to 10):
   • Every existing interval is divided into b = 2 subintervals.
   • For each pair of child intervals, sample b = 2 random multipliers from the distribution chosen in Step 3.
   • Normalize these multipliers so they sum to b = 2 (conserving total mass at each split).
   • Each child receives its parent's mass times its multiplier.
3. Normalize the final measure to unit mass.
4. Integrate: Cumulative sum gives theta(t).

After k = 10 levels, we have 2^10 = 1,024 intervals. The mass distribution is highly non-uniform: some intervals concentrate large amounts of mass (representing high-volatility periods), while others carry almost none (calm periods). This heterogeneity arises from the multiplicative nature of the cascade: mass that starts large at an early level gets multiplied by potentially large factors at subsequent levels, amplifying the initial asymmetry.

# step4_generate_cascade.py

class CascadeGenerator:
    def __init__(self, fitter, b=2, k=10, verbose=config.VERBOSE):
        self.fitter = fitter
        self.b = b
        self.k = k
        self.H = fitter.H
        self.best_distribution = fitter.best_distribution
        self.best_params = fitter.best_params
        self.measure = None
        self.trading_time = None
        self.n_points = b**k

    def generate_cascade(self):
        measure = np.array([1.0])

        for level in range(self.k):
            n_intervals = len(measure)
            new_measure = np.zeros(n_intervals * self.b)
            multipliers = self.sample_multipliers(n_intervals * self.b)

            for i in range(n_intervals):
                parent_mass = measure[i]
                for j in range(self.b):
                    child_idx = i * self.b + j
                    new_measure[child_idx] = parent_mass * multipliers[child_idx]
            measure = new_measure

        self.measure = measure
        self.measure = self.measure / self.measure.sum()
        return self.measure

Sampling Multipliers

The multiplier sampling is where the distribution choice from Step 3 connects to the cascade construction. The theoretical relationship is that each multiplier M is related to a random variable V by M = b^(-V), where V is drawn from the fitted distribution. The conversion from V to M maps a positive random variable to a multiplier between 0 and infinity, with larger V producing smaller M (less mass assigned to that child interval).

def sample_multipliers(self, n_samples):
    params = self.best_params

    if self.best_distribution == 'Normal':
        alpha_0 = params['alpha_0']
        mu = alpha_0 / self.H
        variance = 2 * (alpha_0 / self.H - 1) * np.log(self.b)
        if variance <= 0:
            variance = 0.01
        sigma = np.sqrt(variance)
        V = np.random.normal(mu, sigma, n_samples)

    elif self.best_distribution == 'Binomial':
        alpha_min = params['alpha_min']
        alpha_max = params['alpha_max']
        V = np.random.choice([alpha_min, alpha_max],
                             size=n_samples, p=[0.5, 0.5])

    elif self.best_distribution == 'Poisson':
        alpha_0 = params['alpha_0']
        lam = alpha_0 / self.H
        V = np.random.poisson(lam, n_samples)

    elif self.best_distribution == 'Gamma':
        alpha_0 = params['alpha_0']
        gamma = params['gamma']
        beta = np.log(self.b) / (self.b**(1/gamma) - 1)
        V = np.random.gamma(gamma, 1/beta, n_samples)

    # Convert V to multipliers: M = b^(-V)
    M = self.b ** (-V)
    M = np.clip(M, 1e-10, 1e10)

    # Normalize: each group of b multipliers sums to b
    M = M.reshape(-1, self.b)
    M = M / M.sum(axis=1, keepdims=True) * self.b
    M = M.flatten()

    return M

The normalization step at the end of the multiplier-sampling code is critical. After converting V to M, we reshape the multipliers into pairs (since b=2) and normalize each pair to sum to b=2. This conserves total mass at each cascade level. Without normalization, the total mass would drift randomly, and the resulting trading time would not be a proper CDF. The np.clip guards against numerical extremes (V could be huge or very negative for the Normal and Gamma distributions, producing essentially zero multipliers or astronomically large).

From Measure to Trading Time

The trading time theta(t) is simply the cumulative sum of the normalized measure, scaled to [0, 1]:

def integrate_measure(self):
    if self.measure is None:
        raise ValueError("Must run generate_cascade() first")

    self.trading_time = np.cumsum(self.measure)
    self.trading_time = self.trading_time / self.trading_time[-1]
    return self.trading_time

The resulting theta(t) is a monotonically increasing function from 0 to 1, but unlike a straight line (which would represent uniform clock time), it has steep sections where mass is concentrated and flat sections where mass is sparse. In Step 6, the steep sections will cause the FBM to be sampled rapidly, producing large returns. The flat sections will cause nearly identical FBM values to be sampled, producing tiny returns.

Fig. 1. Raw cascade measure: mass is concentrated in spiky peaks with large gaps between them

Fig. 2. Trading time theta(t) versus uniform clock time (diagonal): steep sections represent high-activity periods

Fig. 3. Log-scale histogram of cascade mass showing the heavy-tailed distribution of interval weights

Fig. 4. Trading time speed (derivative of theta): clustered bursts mirror real market volatility clustering

Step 5: Generating Fractional Brownian Motion

The second MMAR component is the Fractional Brownian Motion (FBM), which provides the long-memory structure. The FBM is generated using the Davies-Harte method, an exact simulation technique that is both efficient (O(n log n) via FFT) and statistically exact (no approximation errors, unlike truncated methods).

The Autocovariance Function

The starting point is the autocovariance function of FBM increments. For increments separated by k steps:

gamma(k) = (1/2) * [|k+1|^(2H) - 2|k|^(2H) + |k-1|^(2H)]

# step5_generate_fbm.py

class FBMGenerator:
    def __init__(self, H, n_points, verbose=config.VERBOSE):
        self.H = H
        self.n_points = n_points
        self.fbm = None
        self.fbm_increments = None

    def fbm_autocovariance(self, k):
        k = np.abs(k)
        return 0.5 * (np.abs(k + 1)**(2*self.H)
                    - 2*np.abs(k)**(2*self.H)
                    + np.abs(k - 1)**(2*self.H))

When k = 0, this gives gamma(0) = 1 (unit variance for each increment). When H = 0.5, gamma(k) = 0 for all k > 0, confirming that standard Brownian motion has independent increments. For H different from 0.5, the correlations are nonzero: positive for H > 0.5 (persistence) and negative for H < 0.5 (anti-persistence).

The Davies-Harte Algorithm

The algorithm generates correlated Gaussian increments in four stages:

Stage 1: Build the circulant covariance vector. We compute the autocovariance at lags 0, 1, 2, ..., n-1, then mirror it to create a symmetric vector of length 2n. This vector is the first row of a circulant matrix, a special structure that can be diagonalized by the discrete Fourier transform.

Stage 2: Compute eigenvalues via FFT. The eigenvalues of the circulant matrix are the FFT of its first row. For a valid covariance matrix, all eigenvalues must be non-negative. If any are negative (which can happen for certain H values and n), we fall back to Cholesky decomposition.

Stage 3: Generate scaled complex Gaussians. We draw 2n independent complex Gaussian random variables (real and imaginary parts each N(0,1)) and scale each by the square root of the corresponding eigenvalue divided by 2n.

Stage 4: Apply inverse FFT. The inverse FFT produces correlated increments. We take the real part of the first n values, which are the FBM increments. The cumulative sum gives the FBM path.

def generate_fbm_davies_harte(self):
    n = self.n_points

    # Stage 1: Circulant covariance vector
    r = np.zeros(2 * n)
    for k in range(n):
        r[k] = self.fbm_autocovariance(k)
    for k in range(1, n):
        r[2*n - k] = r[k]

    # Stage 2: Eigenvalues via FFT
    lam = np.fft.fft(r).real

    if np.any(lam < -1e-10):
        return self.generate_fbm_cholesky()

    lam = np.maximum(lam, 0)

    # Stage 3: Scaled complex Gaussians
    Z_real = np.random.randn(2*n)
    Z_imag = np.random.randn(2*n)
    Z = Z_real + 1j * Z_imag
    W = np.sqrt(lam / (2*n)) * Z

    # Stage 4: Inverse FFT to get correlated increments
    w = np.fft.ifft(W)
    self.fbm_increments = w[:n].real

    # Cumulative sum gives FBM path
    self.fbm = np.cumsum(self.fbm_increments)
    self.fbm = self.fbm - self.fbm[0]
    return self.fbm

The check for negative eigenvalues (np.any(lam < -1e-10)) is a validity test. Negative eigenvalues mean the circulant embedding does not produce a valid covariance matrix, which can happen when n is small relative to the correlation length or when H is close to 0 or 1. The threshold -1e-10 (rather than exactly 0) handles floating-point noise.

The Cholesky Fallback

When the Davies-Harte method detects negative eigenvalues, it falls back to Cholesky decomposition. This method explicitly constructs the n x n covariance matrix from the autocovariance function, computes its Cholesky factorization (L such that C = L * L^T), then multiplies L by a vector of independent standard Gaussians to produce correlated increments. This is slower (O(n^3) vs. O(n log n)) but always works:

def generate_fbm_cholesky(self):
    n = self.n_points

    cov_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            cov_matrix[i, j] = self.fbm_autocovariance(i - j)

    try:
        L = np.linalg.cholesky(cov_matrix)
    except np.linalg.LinAlgError:
        cov_matrix += 1e-10 * np.eye(n)
        L = np.linalg.cholesky(cov_matrix)

    Z = np.random.randn(n)
    increments = L @ Z

    self.fbm_increments = increments
    self.fbm = np.cumsum(increments)
    self.fbm = self.fbm - self.fbm[0]
    return self.fbm

The tiny diagonal perturbation (1e-10 * np.eye(n)) handles the case where the covariance matrix is positive semi-definite but not positive definite (zero eigenvalues cause Cholesky to fail). Adding a negligible amount to the diagonal makes all eigenvalues strictly positive.

Volatility Scaling

After generation, the FBM is scaled so that its increment standard deviation matches the observed market volatility:

def scale_fbm(self, volatility):
    current_std = np.std(self.fbm_increments)
    scale_factor = volatility / current_std if current_std > 0 else 1.0

    self.fbm = self.fbm * scale_factor
    self.fbm_increments = self.fbm_increments * scale_factor

This is important because without scaling, the FBM produces unit-variance increments that have no relationship to actual market volatility. The scale_factor is the ratio of the target volatility (the standard deviation of historical returns from Step 1) to the FBM's current increment standard deviation.

Fig. 5. Generated FBM path scaled to match historical volatility

Fig. 6. FBM increment distribution: Gaussian as expected, with normal overlay for comparison

Fig. 7. Autocorrelation function of FBM increments showing the slow power-law decay characteristic of long memory

Step 6: The MMAR Equation Comes Alive

This is the step where everything connects. We take the Fractional Brownian Motion from Step 5 and evaluate it at the trading time positions from Step 4. The result is the MMAR process: a synthetic price series that exhibits long memory, fat tails, volatility clustering, and scale consistency all at once.

The core equation is:

X(j) = B_H[theta(j)]

The implementation uses linear interpolation because trading time theta(j) maps to fractional indices in the FBM array. As Zhang describes: "if the value was 23.45, the relevant points of the FBM are the 24th and 25th entries. The price for this entry j of the series is a linear interpolation between the values of the FBM for the 24th and 25th entries."

# step6_combine_model.py

class MMARCombiner:
    def __init__(self, fbm_generator, cascade_generator, verbose=config.VERBOSE):
        self.fbm = fbm_generator.fbm
        self.trading_time = cascade_generator.trading_time
        self.mmar_process = None
        self.mmar_returns = None

    def combine_fbm_and_trading_time(self):
        n_points = len(self.trading_time)

        # Create interpolation function mapping [0,1] to FBM values
        fbm_grid = np.linspace(0, 1, len(self.fbm))
        fbm_interp = interp1d(fbm_grid, self.fbm,
                                 kind='linear',
                                 bounds_error=False,
                                 fill_value=(self.fbm[0], self.fbm[-1]))

        self.mmar_process = np.zeros(n_points)
        for j in range(n_points):
            theta_j = self.trading_time[j]
            self.mmar_process[j] = fbm_interp(theta_j)

        self.mmar_returns = np.diff(self.mmar_process)
        return self.mmar_process, self.mmar_returns

Let us trace through the intuition of what this warping does.

Consider two adjacent positions j and j+1 in the cascade. If the measure has a large mass concentration between them, theta(j+1) - theta(j) is large, meaning we jump a big distance along the FBM. The interpolated FBM values at theta(j) and theta(j+1) will be far apart, producing a large return. If the measure has almost no mass between them, theta(j+1) - theta(j) is tiny, and we sample nearly the same FBM point twice, producing a near-zero return. Because the cascade's mass distribution is self-similar and clustered, large-mass regions tend to be adjacent to other large-mass regions, producing volatility clustering. And because some mass concentrations are extreme, the resulting returns have fat tails.

Meanwhile, the FBM provides long memory through its correlated increments. The combination X(t) = B_H[theta(t)] preserves both the FBM's memory structure and the cascade's heterogeneity, producing a process that captures all four stylized facts simultaneously.

The choice of linear interpolation (kind='linear') over higher-order methods (like cubic spline) is deliberate. Linear interpolation is exact for the mapping Zhang describes and avoids oscillation artifacts that cubic splines can introduce at sharp transitions in trading time. The bounds_error=False and fill_value parameters handle the edge case where theta(j) falls exactly at 0 or 1 (rounding errors could push it slightly outside the [0,1] range).

Fig. 8. The MMAR compound process: a single synthetic price path generated by evaluating FBM at multifractal trading time

Fig. 9. MMAR returns: volatility clustering is clearly visible, with bursts of large moves followed by calm periods

Step 7: Monte Carlo Forecasting

A single MMAR simulation produces one realization of the process, one possible future. Because the cascade is random (different multipliers each time) and the FBM is random (different Gaussian draws), each simulation produces a different volatility estimate. To generate a robust forecast, we run 1,000 independent simulations, each with a freshly generated cascade and FBM. The volatility forecast is the mean of these 1,000 volatility estimates.

The Simulation Loop

Each simulation repeats Steps 4 through 6 with fresh random draws:

# step7_monte_carlo.py

class MonteCarloForecaster:
    def __init__(self, fitter, n_simulations=None,
                 forecast_length=None,
                 sample_volatility=None,
                 random_seed=config.RANDOM_SEED,
                 verbose=config.VERBOSE):
        self.fitter = fitter
        self.n_simulations = n_simulations or config.NUM_SIMULATIONS
        self.forecast_length = forecast_length or config.FORECAST_LENGTH
        self.sample_volatility = sample_volatility or 1.0
        self.cascade_b = config.CASCADE_B
        self.cascade_levels = self._calculate_cascade_levels(
            self.forecast_length
        )
        self.n_points = self.cascade_b ** self.cascade_levels

    def _calculate_cascade_levels(self, forecast_length):
        required_points = forecast_length + 1
        levels = int(np.ceil(np.log(required_points) / np.log(self.cascade_b)))
        return max(1, levels)

The _calculate_cascade_levels() method dynamically determines how many cascade levels are needed. Since the cascade produces 2^k points, we need k = ceil(log2(forecast_length + 1)). For forecast_length = 7,200 (25 days of 5-minute bars), this gives k = 13 (since 2^13 = 8,192 > 7,201). The extra points beyond the forecast length are simply truncated.

Each individual simulation runs a complete Step 4–6 sequence:

def run_single_simulation(self):
    # Step 4: Fresh trading time
    cascade_gen = CascadeGenerator(
        self.fitter, b=self.cascade_b,
        k=self.cascade_levels, verbose=False
    )
    cascade_gen.generate_cascade()
    cascade_gen.integrate_measure()

    # Step 5: Fresh FBM
    fbm_gen = FBMGenerator(
        H=self.fitter.H, n_points=self.n_points, verbose=False
    )
    fbm_gen.generate_fbm_davies_harte()
    fbm_gen.scale_fbm(self.sample_volatility)

    # Step 6: Combine
    combiner = MMARCombiner(fbm_gen, cascade_gen, verbose=False)
    _, returns = combiner.combine_fbm_and_trading_time()
    returns = returns[:self.forecast_length]

    return np.std(returns)

Note that verbose=False is set on all sub-components to prevent 1,000 copies of progress messages. The returns are truncated to forecast_length because the cascade produces more points than needed (2^k > forecast_length for the chosen k).

The main Monte Carlo loop runs all simulations and collects the volatility estimates:

def run_monte_carlo(self):
    if self.random_seed is not None:
        np.random.seed(self.random_seed)

    self.volatility_forecasts = []

    for i in range(self.n_simulations):
        volatility = self.run_single_simulation()
        self.volatility_forecasts.append(volatility)

        if self.verbose and (i + 1) % max(1, self.n_simulations // 10) == 0:
            progress = 100 * (i + 1) / self.n_simulations
            mean_so_far = np.mean(self.volatility_forecasts)
            print(f"  Progress: {progress:.1f}% - Mean vol: {mean_so_far:.6f}")

    self.mean_forecast = np.mean(self.volatility_forecasts)
    self.std_forecast = np.std(self.volatility_forecasts)

The output includes a point estimate (mean volatility across all simulations), a measure of uncertainty (standard deviation of the volatility estimates), and a 95% confidence interval (mean +/- 1.96 * std). The runner script loads the actual sample volatility from the training data to ensure the FBM scaling matches observed market conditions:

# run_step7.py

# Load historical returns from Step 1 to get actual sample volatility
with open(step1_path, 'rb') as f:
    checker = pickle.load(f)

sample_volatility = np.std(checker.returns)

forecaster = run_monte_carlo_forecast(fitter,
    sample_volatility=sample_volatility,
    n_simulations=config.NUM_SIMULATIONS,
    forecast_length=config.FORECAST_LENGTH)

Convergence

The convergence plot, showing the cumulative mean stabilizing as simulations accumulate, serves as a sanity check that 1,000 simulations are sufficient. If the cumulative mean is still drifting significantly at 1,000, more simulations would be needed. In practice, the mean typically stabilizes within the first 500-800 simulations, with subsequent simulations only refining the estimate. The Q-Q plot checks whether the forecast distribution is approximately normal, as suggested by the central limit theorem. Each forecast is an aggregate statistic (the standard deviation of thousands of returns).

Fig. 10.  Distribution of 1,000 Monte Carlo volatility forecasts: approximately Gaussian as expected

Fig. 11. Convergence of the cumulative mean: the estimate stabilizes well before 1,000 simulations

Fig. 12. Individual forecast realizations with 95% confidence interval band

Fig. 13. Q-Q plot confirming normality of the forecast distribution

The GARCH Implementation

For a fair comparison, we fit GARCH models on the same training data used for the MMAR. This is not just good practice, it is essential: any difference in training data would confound the comparison, making it impossible to attribute performance differences to the models themselves.

Our implementation uses the arch Python library, a well-established package for ARCH/GARCH modeling. We fit a standard GARCH(1,1) with normal innovations:

# garch_model.py

class GARCHForecaster:
    def __init__(self, returns, p=1, q=1,
                 dist='normal', random_seed=config.RANDOM_SEED,
                 verbose=True):
        self.returns = returns * 100  # arch expects percentage returns
        self.p = p
        self.q = q
        self.dist = dist

The multiplication by 100 on the first line is a common requirement of the arch library: it expects returns expressed as percentages rather than decimals. This scaling is reversed after forecasting to keep units consistent with the MMAR output.

Model fitting uses rescale=True, which is critical for numerical stability. EURUSD returns on a 5-minute timeframe are tiny numbers (on the order of 0.0001), and even after the 100x scaling, the conditional variance values can be tiny. The arch library's rescale feature automatically scales the data internally for numerical stability during optimization, then converts the results back to the original scale. Without it, the optimizer can fail to converge or produce garbage parameters.

def fit(self):
    self.model = arch_model(
        self.returns,
        vol='GARCH',
        p=self.p,
        q=self.q,
        dist=self.dist,
        rescale=True,
    )
    self.fit_result = self.model.fit(
        disp='off',
        options={'ftol': 1e-6, 'maxiter': 1000}
    )

The Forecasting Methodology

This is where a critical and potentially confusing design decision must be explained. Zhang does not use GARCH's standard analytic multi-step forecast (which would give a variance path). Instead, he simulates one path of FORECAST_LENGTH returns from the fitted model and takes the standard deviation. This mirrors the R function ugarchsim. This approach ensures the GARCH forecast is measured in the same units and over the same horizon as the MMAR forecast:

def forecast(self, n_sim=None):
    if n_sim is None:
        n_sim = config.FORECAST_LENGTH

    if self.random_seed is not None:
        np.random.seed(self.random_seed)

    # Simulate ONE path of n_sim returns
    sim = self.fit_result.forecast(
        horizon=n_sim, method='simulation', simulations=1
    )
    sim_returns = sim.simulations.values[0, 0, :]

    # Unscale if model was rescaled
    if hasattr(self.model, 'scale') and self.model.scale is not None:
        sim_returns = sim_returns / self.model.scale

    # Convert from percentage back to decimal
    self.forecast_volatility = np.std(sim_returns) / 100
    return self.forecast_volatility

There are two scaling corrections happening here. First, if the model used internal rescaling (rescale=True), the simulated returns are in rescaled units and must be divided by self.model.scale to get back to the percentage-return space. Second, we divide by 100 to convert from percentage returns back to decimal returns (reversing the * 100 from the constructor).

Why the Different Methodology?

This is not a methodological mismatch between MMAR and GARCH; it is precisely how Zhang (2017) compared the models. The MMAR needs many simulations (1,000 in our case; Zhang used 10,000) because each cascade realization is random: different multiplier draws produce different mass distributions, leading to different volatility estimates. Averaging over many simulations smooths out this stochastic variability. GARCH needs only one simulation because the model's dynamics are fully determined by the fitted parameters (omega, alpha, beta). Given the same fitted model, the conditional variance path is deterministic; only the innovation draws are random. A single long simulation already averages over enough innovations to produce a stable volatility estimate.

The comparison runner loads the same training data from Step 1, ensuring perfect consistency:

# garch_model.py

def run_garch_comparison(p=1, q=1, dist='normal',
                         output_dir=None, verbose=True):
    # Load training data from Step 1 (same as MMAR uses)
    step1_path = Path(config.OUTPUT_DIR) / "step1_checker.pkl"
    with open(step1_path, 'rb') as f:
        checker = pickle.load(f)

    returns = checker.returns

    forecaster = GARCHForecaster(
        returns=returns, p=p, q=q, dist=dist, verbose=verbose
    )
    forecaster.fit()
    forecast_vol = forecaster.forecast(n_sim=config.FORECAST_LENGTH)

    return {'GARCH': {
        'forecaster': forecaster,
        'forecast_volatility': forecast_vol,
        'aic': forecaster.fit_result.aic,
        'bic': forecaster.fit_result.bic,
        'loglik': forecaster.fit_result.loglikelihood
    }}

The Showdown: Results

With both models fitted and forecasts generated, we load the out-of-sample data (the 25 days following our training period) and compute realized volatility as the standard deviation of actual returns. This is our ground truth. The run_garch_comparison.py script orchestrates the comparison. It loads the Step 7 MMAR forecast, fits the GARCH model, computes realized volatility via ForecastComparison, and ranks both models by percentage error:

# run_garch_comparison.py

# Step 1: Fit GARCH model
garch_results = run_garch_comparison(p=1, q=1, dist='normal', verbose=True)

# Step 2: Load MMAR forecast from Step 7
mmar_path = Path(config.OUTPUT_DIR) / "step7_forecaster.pkl"
with open(mmar_path, 'rb') as f:
    mmar_forecaster = pickle.load(f)
mmar_forecast = mmar_forecaster.mean_forecast

# Step 3: Get realized volatility via ForecastComparison
comparison = run_forecast_comparison(save_plots=False)
realized_vol = comparison.realized_volatility

The ForecastComparison class handles loading the forecast period data (from the day after the training period ends) and computing realized volatility as the sample standard deviation with Bessel's correction (ddof=1), the unbiased estimator:

# compare_forecast.py

class ForecastComparison:
    def load_forecast_period_data(self, start_date=None, end_date=None,
                                 days_ahead=None):
        if start_date is None:
            training_end = pd.to_datetime(config.END_DATE)
            start_date = (training_end + timedelta(days=1)).strftime("%Y-%m-%d")

        if days_ahead is None:
            days_ahead = config.FORECAST_DAYS

        loader = DataLoader(
            symbol=config.SYMBOL,
            start_date=start_date,
            end_date=end_date
        )
        loader.load_from_mt5()
        loader.calculate_returns(price_column='close', method='log')
        self.realized_returns = loader.get_returns_array()

    def calculate_realized_volatility(self):
        self.realized_volatility = np.std(self.realized_returns, ddof=1)

The error metrics include absolute error, percentage error, MAPE, RMSE, confidence interval bounds, and whether the realized volatility falls within the MMAR's 95% confidence interval:

def compute_error_metrics(self):
    self.error = self.forecast_volatility - self.realized_volatility
    self.percent_error = 100 * self.error / self.realized_volatility
    self.rmse = np.sqrt((self.error)**2)
    mape = np.abs(self.percent_error)

    lower_bound = self.forecast_volatility - 1.96 * self.forecast_std
    upper_bound = self.forecast_volatility + 1.96 * self.forecast_std
    within_ci = lower_bound <= self.realized_volatility <= upper_bound

    return {
        'forecast': self.forecast_volatility,
        'realized': self.realized_volatility,
        'error': self.error,
        'percent_error': self.percent_error,
        'abs_percent_error': mape,
        'rmse': self.rmse,
        'forecast_std': self.forecast_std,
        'forecast_95ci_lower': lower_bound,
        'forecast_95ci_upper': upper_bound,
        'within_95ci': within_ci
    }

The comparison also computes annualized volatility at multiple timeframes (5-minute, hourly, daily, annual) using the square root of time rule. This provides context for interpreting the results in terms familiar to practitioners:

# run_comparison.py

periods_per_year = config.periods_per_year()
ann_forecast = comparison.forecast_volatility * (periods_per_year ** 0.5) * 100
ann_realized = comparison.realized_volatility * (periods_per_year ** 0.5) * 100

The final comparison table ranks MMAR and GARCH by percentage error relative to realized volatility:

# terminal output

======================================================================

FULL COMPARISON: MMAR vs GARCH vs REALIZED

======================================================================

Model           Forecast        Error           Error %    AIC          BIC

-------------------------------------------------------------------------------------

MMAR            0.0005768527  0.0002260410    64.43%       N/A       N/A

GARCH           0.0010309505  0.0006801388   193.88%    36880.79    36912.30

-------------------------------------------------------------------------------------

REALIZED        0.0003508117

=====================================================================================

BEST MODEL: MMAR (Error: 64.43%)

======================================================================

MMAR vs GARCH

======================================================================

  MMAR reduces percentage error by 129.44 pp

These results are consistent with Zhang (2017), who demonstrated across 20 US equities that the MMAR produced lower root mean squared forecast error than all GARCH variants tested. The key reason is structural: GARCH models conditional variance as a single time-series process, so it can only capture one scaling behavior. The MMAR, by decomposing returns into a long-memory component (FBM) and a multifractal time-deformation component (trading time), captures the full spectrum of scaling behaviors present in real market data. Where GARCH sees a single volatility process, the MMAR sees a rich multiscale structure.

Conclusion

Across Part 1, Part 2, and this article, we have built and tested the complete MMAR pipeline following Zhang (2017). In Parts 1 and 2, we confirmed that EURUSD 5-minute data is genuinely multifractal and extracted the parameters needed for construction. In this article, we put those parameters to work and ran a fair head-to-head comparison. Here is what the full picture tells us:

• The MMAR captures more. By decomposing returns into Fractional Brownian Motion (long memory) and multifractal trading time (fat tails, clustering), the MMAR reproduces all four stylized facts of financial returns from first principles, not through parameter patching. Long memory comes from FBM's power-law autocovariance. Fat tails come from non-uniform time sampling. Volatility clustering comes from the cascade's self-similar mass distribution. Scale consistency comes from the multifractal structure itself.
• GARCH is not wrong, but limited. GARCH models are fast, well understood, and effective within their assumptions. But those assumptions, thin tails, short memory, uniform time—are systematically violated by real market data. GARCH extensions (EGARCH for asymmetry, FIGARCH for long memory, Student-t for fat tails) each patch one limitation while leaving the others unaddressed. The MMAR handles all of them through a single architectural choice: the compound process X(t) = B_H[theta(t)].
• The comparison is fair. Both models used the same training data, the same forecast horizon, and Zhang's exact methodology for generating forecasts.

That said, the MMAR is not without limitations. It is computationally expensive, as 1,000 Monte Carlo simulations, each running a full cascade + FBM + combination pipeline, take significantly longer than fitting a GARCH model (which is a single optimization). The cascade construction introduces stochastic variability that requires many simulations to average out. And the model's performance depends on the data being genuinely multifractal; for assets or timeframes that fail the multifractality test in Part 2's Step 2, GARCH may well be the better choice.

The computational cost is real but manageable. On a modern machine, 1,000 simulations take under a minute. And the cost is entirely in the Monte Carlo step; the analysis and fitting from Parts 1 and 2 are fast. For production use, one would run the analysis once, cache the parameters (H, distribution type, distribution parameters), and run Monte Carlo forecasts on demand.

In Part 4 of this series, we begin translating the MMAR into native MQL5 code, starting with the partition analysis—the scaling function, Hurst exponent, and multifractality testing that form the estimation foundation. Part 5 tackles fitting the multifractal spectrum to theoretical distributions using bounded optimization. Part 6 builds the simulation engine: Fractional Brownian Motion via the Davies-Harte algorithm and the multiplicative cascade that generates multifractal trading time. Part 7 wraps the simulation engine in a Monte Carlo loop to produce volatility forecasts with confidence intervals. Finally, Part 8 assembles the complete MMAR library and puts it to work in an Expert Advisor.

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 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
step4_generate_cascade.py	Multiplicative cascade with distribution-specific multiplier sampling, mass normalization, trading time integration
step5_generate_fbm.py	FBM generation via Davies-Harte (FFT) with Cholesky fallback, volatility scaling
step6_combine_model.py	MMAR compound process: X(t) = B_H[theta(t)] via linear interpolation
step7_monte_carlo.py	Monte Carlo volatility forecasting with dynamic cascade levels, convergence monitoring, 1,000 simulations
garch_model.py	GARCH(1,1) implementation using arch library with Zhang's simulation methodology
compare_forecast.py	MMAR forecast validation: error metrics, rolling realized volatility, annualized comparison plots
run_step[1-7].py	Runner scripts for each step
run_comparison.py	Runner: MMAR vs realized validation with annualized metrics and directional analysis
run_garch_comparison.py	Runner: fits GARCH, loads MMAR forecast, produces head-to-head comparison table and ranking
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