Beyond GARCH (Part I): Mandelbrot's MMAR versus Engle's GARCH

Introduction

Volatility forecasting is the backbone of risk management in algorithmic trading. Get it right, and your position sizing adapts to market conditions, your stop losses breathe with price action, and your drawdowns stay controlled. Get it wrong, and a single regime change wipes out months of profit. For decades, the go-to tool for this job has been GARCH, the Generalized Autoregressive Conditional Heteroskedasticity model. It was introduced by Robert Engle in 1982 and generalized by Tim Bollerslev in 1986. GARCH earned Engle a Nobel Prize, and for good reason: it captures volatility clustering, the empirical fact that large price moves tend to follow large moves and small moves follow small ones.

But GARCH has blind spots. It assumes returns follow a normal (or at best, Student-t) distribution, yet real markets produce fat tails far more extreme than any bell curve predicts. It has no mechanism for long memory, the tendency for volatility persistence to stretch across weeks and months, not just a few periods. And it treats time as uniform, even though any trader knows that ten minutes during a news release carry more information than ten minutes at 3 AM. These are not minor quibbles. There are structural limitations baked into the model's assumptions.

Enter Benoit Mandelbrot, the father of fractal geometry. As early as 1963, Mandelbrot observed that financial returns exhibit self-similarity, patterns that repeat across different time scales. A 5-minute candlestick chart of EURUSD looks structurally indistinguishable from a daily chart. Both show the same clustering of volatile and calm periods, the same heavy-tailed return distributions, and the same gradual decay of autocorrelations. This is not a coincidence. It is a signature of an underlying fractal process. In 1997, Mandelbrot, together with Calvet and Fisher, formalized this insight into the Multifractal Model of Asset Returns (MMAR).

The MMAR does not patch GARCH's weaknesses with ad hoc extensions. Instead, it starts from a fundamentally different premise: financial returns are generated by a Fractional Brownian Motion running on multifractal trading time. This is a non-uniform clock that speeds up during volatile periods and slows down during calm ones. This single architectural choice produces long memory, fat tails, volatility clustering, and scale consistency all at once. They emerge as natural consequences of the same generative mechanism, not as separate patches bolted onto a base model.

In this article, we begin the MMAR construction pipeline following the methodology of Zhang (2017). His paper, "Volatility Forecasting with the Multifractal Model of Asset Returns", demonstrated that the MMAR significantly outperforms GARCH across 20 stocks. Using Python and the MetaTrader5 Python API, we will lay the analytical groundwork: loading and preparing the data, and testing whether it exhibits multifractal scaling using the partition function. By the end of this article, we will have a definitive answer to a prerequisite question: is our data fractal at all?

We will cover:

1. GARCH: Engle's Legacy
2. The MMAR: Mandelbrot's Fractal Alternative
3. Setting Up: Data and Configuration
4. Step 1: Testing for Fractality
5. Conclusion

GARCH: Engle's Legacy

Before we can appreciate what the MMAR brings to the table, we need to understand what we are comparing it against. GARCH remains the most widely used approach to volatility forecasting in both academia and industry, and for good reason: it is elegant, well-understood, and fast to compute. We will go through the core idea, mention its most important extensions for context, and then carefully explain what GARCH cannot do. Understanding those limitations motivates the entire rest of this article.

The Core Idea

GARCH models the conditional variance of returns, that is, the variance at time t given everything we know up to time t-1. The standard GARCH(1,1) model consists of two equations. The mean equation states that the return at time t equals a constant plus a noise term:

r_t = mu + epsilon_t

Where epsilon_t = sigma_t * z_t and z_t is drawn from a standard normal distribution. The variance equation is where the action happens:

sigma_t^2 = omega + alpha * epsilon_(t-1)^2 + beta * sigma_(t-1)^2

This says that today's variance is a weighted combination of three things: a long-run baseline (omega), yesterday's squared shock (alpha term), and yesterday's variance (beta term). The alpha term captures reactivity, how quickly variance responds to new information. The beta term captures persistence, how long elevated volatility lingers. Together, they produce volatility clustering. After a large return, sigma_t^2 increases through the alpha term, which raises the conditional variance for the next period and makes another large return more likely. As the large returns subside, the beta term gradually decays the variance back toward omega.

The parameters omega, alpha, and beta are estimated via Maximum Likelihood Estimation (MLE). The likelihood function measures how probable the observed returns are given a particular set of parameters. Starting from an initial guess, the optimizer adjusts omega, alpha, and beta to maximize this probability. The constraints are that all three must be non-negative and that alpha + beta must be less than 1 for the process to be stationary. This stationarity condition ensures the unconditional variance converges to a finite value: omega / (1 - alpha - beta). The MLE approach is well-suited to GARCH because the conditional distribution of each return is fully specified given the past, making the likelihood function straightforward to compute.

Fig. 1 — GARCH variance dynamics: a large shock causes sigma to spike, then gradually decay back to baseline

The Extensions

The standard GARCH(1,1) model treats positive and negative shocks symmetrically: a 2% drop has the same effect on future volatility as a 2% rise. In practice, markets exhibit a leverage effect: negative returns tend to increase volatility more than positive returns of the same magnitude. This asymmetry is well-documented, particularly in equity markets, and has been attributed to the leverage mechanism (falling stock prices increase firm leverage, increasing risk and thus volatility).

Extensions like EGARCH (Nelson, 1991) and GJR-GARCH (Glosten, Jagannathan, and Runkle, 1993) address the leverage effect, where negative returns increase volatility more than positive returns of the same magnitude. These are well-documented improvements, but they do not address GARCH's deeper structural limitations. In our comparison, we use the standard GARCH(1,1) as the benchmark, as it is the most widely used and well-understood variant.

Fig. 2 — The leverage effect: identical-magnitude positive and negative shocks produce asymmetric variance responses

What GARCH Cannot Do

For all its strengths, GARCH makes assumptions that financial data consistently violates. Understanding these limitations is essential because the MMAR addresses each one of them, not through patches, but through its fundamental architecture.

Limitation 1: Thin Tails. Even with a Student-t innovation distribution, GARCH underestimates the frequency of extreme events. The issue is that GARCH's conditional distribution is parametric, meaning it assumes a specific functional form for the tails. Real return distributions have heavier tails than any single parametric family can capture. This matters in practice: a model that underestimates the probability of a 5-sigma move will underestimate Value-at-Risk and Expected Shortfall, leading to inadequate capital reserves.

Limitation 2: Short Memory. GARCH has exponentially decaying memory. The autocorrelation of squared returns in a GARCH(1,1) process decays as (alpha + beta)^k, where k is the lag. For typical parameter values (alpha + beta around 0.95-0.99), this means correlations become negligible within a few dozen periods. But empirical studies consistently find hyperbolic decay in squared returns: the autocorrelation at lag k decays as k^(2d-1), where d is the long-memory parameter, typically between 0.3 and 0.5. This is dramatically slower than exponential decay. A volatility cluster that started three months ago still influences today's volatility, but GARCH has effectively "forgotten" it. The FIGARCH extension partially addresses this, but it does so within the GARCH framework rather than from first principles.

Limitation 3: No Scale Consistency. GARCH operates on a single time scale. A GARCH model fitted to daily returns says nothing about hourly or weekly return dynamics. But empirical evidence proves that the same statistical patterns appear whether you look at 5-minute, hourly, or daily returns: the same fat tails, the same volatility clustering, and the same long memory. This property is known as moment scaling or scale consistency. A good model should produce moment scaling naturally, due to its structure, not require re-fitting at each scale. GARCH does not have this property.

These are precisely the gaps the MMAR was designed to fill.

Fig. 3 — GARCH's three blind spots: thin tails, short memory, and no scale consistency

The MMAR: Mandelbrot's Fractal Alternative

The Multifractal Model of Asset Returns takes a fundamentally different approach to modeling financial returns. Rather than modeling variance as a time series (as GARCH does), the MMAR models returns as a stochastic process running on a deformed clock. The core equation, following Zhang (2017), is deceptively simple:

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

This says that the log-price process X(t) is a Fractional Brownian Motion B_H evaluated at a multifractal trading time theta(t). Two components, each handling different properties of financial returns. Let us examine each one carefully, because understanding these components is essential for understanding the entire pipeline we will build.

Component 1: Fractional Brownian Motion (FBM)

Fractional Brownian Motion is a generalization of standard Brownian motion parameterized by the Hurst exponent H, a single value between 0 and 1 that controls the memory structure of the process. When H = 0.5, FBM reduces to ordinary Brownian motion, a pure random walk with independent increments and no memory. H > 0.5 indicates the process is persistent: an upward move makes the next upward move more likely, producing trending behavior. When H < 0.5, the process is anti-persistent: moves tend to reverse, producing mean-reverting behavior.

The mathematical signature of FBM is its autocovariance function. For FBM increments separated by k time steps, the autocovariance is:

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

This formula reveals the crucial difference from standard Brownian motion. When H = 0.5, gamma(k) equals zero for all k > 0, meaning increments are uncorrelated (a random walk). But when H is different from 0.5, gamma(k) decays as a power law: gamma(k) is proportional to H(2H-1) * k^(2H-2) for large k. For H > 0.5, this gives positive correlations that decay hyperbolically. This is fundamentally different from GARCH's exponential decay. In GARCH, correlations that are 50 lags old are exponentially tiny. In FBM, they are still meaningfully present, decaying slowly as a power law. This is long memory, and it is precisely what empirical studies find in financial data.

For financial markets, typical Hurst exponent values fall between 0.48 and 0.58. Values near 0.5 suggest behavior close to a random walk; values significantly above 0.5 suggest persistent trending; values below 0.5 suggest anti-persistent mean reversion. FX markets typically cluster around 0.48-0.55, while stock indices sometimes show slightly higher values (0.52-0.58).

Fig. 4 — FBM paths for H=0.3 (anti-persistent), H=0.5 (random walk), and H=0.7 (persistent), with autocorrelation decay below each

Component 2: Multifractal Trading Time

Trading time theta(t) is a cumulative distribution function of a multifractal measure, constructed via a multiplicative cascade. To understand what this means, think about what "time" means in financial markets.

Clock time ticks uniformly. Every minute is the same duration as every other minute. But market activity is not uniform. During high-volatility events, like a central bank announcement or an unexpected geopolitical development, enormous amounts of information arrive in a short clock-time interval. During quiet Asian-session hours, almost no meaningful price discovery happens. The idea of "trading time" captures this non-uniformity: it is a clock that runs fast during active periods and slow during quiet ones.

Formally, trading time theta(t) is a monotonically increasing function from 0 to 1. Unlike a straight line representing uniform clock time, it has steep sections and flat sections. In the steep sections, a small increment of clock time maps to a large increment of trading time. The FBM gets sampled at rapidly changing positions, producing large returns. In the flat sections, clock time advances but trading time barely changes, so the FBM is sampled at nearly the same position, producing tiny returns. This mechanism naturally produces fat tails and volatility clustering without any distributional assumption on the returns themselves. The steep sections generate extreme returns, and the steep and flat sections cluster together due to the cascade's self-similar structure.

The cascade construction that produces this trading time starts with a single interval [0, 1] carrying mass 1. At each level, every interval is divided into b = 2 subintervals (a binary cascade), and each child receives its parent's mass multiplied by a random multiplier. The multipliers are drawn from a specific distribution (Normal, Binomial, Poisson, or Gamma, determined by fitting in Step 3) and normalized so that each pair sums to 2, conserving total mass. After k = 10 levels, we have 2^10 = 1,024 intervals with a rich, heterogeneous mass distribution. The cumulative sum of this distribution gives us theta(t).

Fig. 5 — Multiplicative cascade construction: mass is recursively redistributed across binary splits over 10 levels

Why Both Components Together?

Neither component alone reproduces all the stylized facts of financial returns. FBM alone gives you long memory but Gaussian returns with no fat tails. A multifractal measure alone gives you fat tails and clustering but no long-range dependence in the return series. The compound process X(t) = B_H[theta(t)] produces all four properties simultaneously:

• Long memory: from FBM's hyperbolic autocovariance decay
• Fat tails: from trading time's non-uniform sampling of the FBM
• Volatility clustering: from the cascade's self-similar clumping of mass
• Scale consistency (moment scaling): from the multifractal structure itself, which is, by construction, self-similar across scales

This is the key insight of Mandelbrot, Calvet, and Fisher (1997), and it is what makes the MMAR architecturally superior to GARCH extensions that bolt on fixes one at a time. GARCH patches the leverage effect with EGARCH/GJR. It patches long memory with FIGARCH. It patches fat tails with Student-t innovations. Each patch adds parameters and complexity without addressing the root cause: returns are generated by a multiscale process. Any model that ignores this structure will always be playing catch-up.

Fig. 6 — The MMAR compound process: FBM evaluated at multifractal trading time produces all four stylized facts simultaneously

Setting Up: Data and Configuration

Our entire implementation follows Zhang (2017) and is built as a 7-step Python pipeline, where each step produces outputs that feed into the next. All configuration lives in a single file, config.py, which centralizes every parameter the pipeline uses. This design makes experiments reproducible and parameter sweeps straightforward.

Let us walk through the configuration file in detail, because every parameter has a reason behind it:

# config.py

# =============================================================================
# DATA PARAMETERS
# =============================================================================

SYMBOL = "EURUSD"
TIMEFRAME_MT5 = "M5"
TIMEFRAME_MINUTES = 5

START_DATE = "2024-08-03"
END_DATE = "2026-03-01"

FORECAST_DAYS = 25
FORECAST_INTERVAL_MINUTES = TIMEFRAME_MINUTES
FORECAST_LENGTH = int(FORECAST_DAYS * (24 * 60) / TIMEFRAME_MINUTES)
TRADING_DAYS_PER_YEAR = 252

We are analyzing EURUSD on the 5-minute timeframe, with a training period from August 2024 to March 2026. This is a substantial period, roughly 18 months of 5-minute bars, giving us tens of thousands of returns to work with. The forecast horizon is 25 calendar days, and the FORECAST_LENGTH is calculated dynamically as the number of 5-minute intervals in 25 days (7,200 returns). The FORECAST_INTERVAL_MINUTES is set equal to TIMEFRAME_MINUTES to ensure that the forecast, the realized volatility, and the annualization factor all use the same return interval. If these diverged, our volatility comparison would be apples to oranges.

# =============================================================================
# STEP 1: PARTITION FUNCTION PARAMETERS
# =============================================================================

DELTA_T_MIN = 1
DELTA_T_MAX = 150
DELTA_T_SPACING_FACTOR = 1.1
NUM_DELTA_T_VALUES = 30

Q_MIN = 0.01
Q_MAX = 30.0
Q_STEP = 0.5

MIN_R_SQUARED = 0.60

The partition function parameters deserve explanation. Delta_t ranges from 1 to 150 observations, spaced logarithmically with a factor of 1.1. Logarithmic spacing is important because we are testing for power-law scaling, which is a linear relationship on a log-log plot. Uniform spacing would over-sample short timescales and under-sample long ones, distorting the regression. The q moment orders range from 0.01 to 30.0 in steps of 0.5 (yielding about 60 q values). Zhang used this range, and it covers both the low-order moments (q < 2, which characterize the center of the distribution) and high-order moments (q > 5, which characterize the tails). The R-squared threshold of 0.60 is our minimum evidence bar for fractality; Zhang achieved an average of 0.66 across his 20 stocks.

The function generate_delta_t_values() produces the log-spaced grid:

def generate_delta_t_values():
    delta_t_values = []
    current = DELTA_T_MIN
    while current <= DELTA_T_MAX:
        delta_t_values.append(int(current))
        current *= DELTA_T_SPACING_FACTOR

    # Remove duplicates and sort
    return np.unique(np.array(delta_t_values))

Starting from 1, each subsequent value is 1.1 times the previous one, giving us a geometrically-spaced sequence. The int() truncation and np.unique() call handle the fact that adjacent values might round to the same integer at small scales. The cascade and Monte Carlo parameters are also defined here:

# =============================================================================
# STEP 4-6: CASCADE AND FBM PARAMETERS
# =============================================================================

CASCADE_B = 2    # Binary cascade (divide into 2 intervals at each step)
CASCADE_K_MAX = 10 # Standalone cascade depth (2^10 = 1024 intervals)

# =============================================================================
# STEP 7: MONTE CARLO SIMULATION PARAMETERS
# =============================================================================

NUM_SIMULATIONS = 1000  # Paper used 10,000
RANDOM_SEED = 42

The cascade uses a binary split (b=2) with k=10 levels, producing 2^10 = 1,024 grid points per cascade. The Monte Carlo uses 1,000 simulations with a fixed random seed for reproducibility. Zhang used 10,000 in his paper; we use 1,000 as a practical compromise that converges sufficiently while keeping computation time manageable.

Configuration also includes a validation function that catches parameter errors at load time:

def validate_config():
    assert DELTA_T_MIN > 0, "DELTA_T_MIN must be positive"
    assert DELTA_T_MAX > DELTA_T_MIN, "DELTA_T_MAX must be greater than DELTA_T_MIN"
    assert Q_MIN > 0, "Q_MIN must be positive"
    assert Q_MAX > Q_MIN, "Q_MAX must be greater than Q_MIN"
    assert 0 <= MIN_R_SQUARED <= 1, "MIN_R_SQUARED must be between 0 and 1"
    assert FORECAST_INTERVAL_MINUTES == TIMEFRAME_MINUTES, (
        "FORECAST_INTERVAL_MINUTES should match TIMEFRAME_MINUTES so forecasts, "
        "realized volatility, and annualization use the same return interval"
    )
    assert CASCADE_B >= 2, "CASCADE_B must be at least 2"
    assert CASCADE_K_MAX > 0, "CASCADE_K_MAX must be positive"

The assertion that FORECAST_INTERVAL_MINUTES matches TIMEFRAME_MINUTES is important. If you are training on 5-minute bars but forecasting at a different interval, the volatility units will not match, and the comparison with realized volatility will be meaningless.

Loading Data from MetaTrader 5

Data is loaded via the DataLoader class, which wraps the MetaTrader5 Python package. The class supports two data sources: a primary MetaTrader 5 connection and a CSV fallback for environments where MetaTrader 5 is not available. The MetaTrader 5 initialization handles connection gracefully, reporting the terminal, account, and server on success:

# data_loader.py

try:
    import MetaTrader5 as mt5
    MT5_AVAILABLE = True
except ImportError:
    MT5_AVAILABLE = False
    print("Warning: MetaTrader5 module not installed. Only CSV loading available.")

class DataLoader:
    def _initialize_mt5(self):
        if not MT5_AVAILABLE:
            return False

        if not mt5.initialize():
            if self.verbose:
                print(f"MT5 initialization failed: {mt5.last_error()}")
            return False

        self.mt5_initialized = True
        return True

Once connected, load_from_mt5() fetches OHLCV bars for the configured symbol and date range using mt5.copy_rates_range(), then converts the result into a pandas DataFrame with standardized column names:

def load_from_mt5(self, symbol=None, timeframe=None,
                 start_date=None, end_date=None):
    symbol = symbol or self.symbol
    timeframe = timeframe or config.TIMEFRAME_MT5

    if not self._initialize_mt5():
        raise RuntimeError("Failed to initialize MetaTrader 5.")

    mt5_tf = self._get_mt5_timeframe(timeframe)
    rates = mt5.copy_rates_range(symbol, mt5_tf, start_date, end_date)

    if rates is None or len(rates) == 0:
        raise ValueError(f"Failed to fetch data: {mt5.last_error()}")

    df = pd.DataFrame(rates)
    df['time'] = pd.to_datetime(df['time'], unit='s')
    df.set_index('time', inplace=True)
    df.columns = ['open', 'high', 'low', 'close',
                 'tick_volume', 'spread', 'real_volume']
    self.data = df
    return df

The CSV loader provides flexibility for multiple datetime column formats, handling both separate date/time columns and combined datetime columns:

def load_csv(self, filepath):
    df = pd.read_csv(filepath)
    df.columns = df.columns.str.lower()

    if 'date' in df.columns and 'time' in df.columns:
        df['datetime'] = pd.to_datetime(
            df['date'].astype(str) + ' ' + df['time'].astype(str)
        )
    elif 'datetime' in df.columns:
        df['datetime'] = pd.to_datetime(df['datetime'])
    elif 'time' in df.columns:
        df['datetime'] = pd.to_datetime(df['time'])

    df.set_index('datetime', inplace=True)
    df.sort_index(inplace=True)
    df = df.loc[self.start_date:self.end_date]
    self.data = df
    return df

Computing Log Returns

With the raw OHLCV data loaded, we compute log returns, the standard representation for multifractal analysis:

def calculate_returns(self, price_column='close', method='log'):
    if self.data is None:
        raise ValueError("No data loaded. Call load_from_mt5() or load_csv() first.")

    prices = self.data[price_column].values
    self.prices = prices

    if method == 'log':
        # Log returns: ln(P_t / P_{t-1}) = ln(P_t) - ln(P_{t-1})
        returns = np.diff(np.log(prices))
    elif method == 'simple':
        # Simple returns: (P_t - P_{t-1}) / P_{t-1}
        returns = np.diff(prices) / prices[:-1]

    returns_series = pd.Series(returns, index=self.data.index[1:])
    self.returns = returns_series
    return returns_series

Log returns have a key property that makes them ideal for multifractal analysis: they are additive across time scales. A 30-minute log return is precisely the sum of the six 5-minute log returns within it. Simple returns do not have this property (they are multiplicative). This additivity is what allows the partition function in Step 1 to test moment scaling across different aggregation levels. The np.diff(np.log(prices)) computation is numerically equivalent to np.log(prices[1:] / prices[:-1]) but avoids creating a temporary ratio array.

The DataLoader also performs data quality checks via summary_statistics(). It flags NaN values and high percentages of zero returns (which can indicate data gaps or illiquidity) and computes skewness and kurtosis. Excess kurtosis greater than zero signals fat tails, which is what we expect to find:

def summary_statistics(self):
    nan_count = self.returns.isna().sum()
    zero_count = (self.returns == 0).sum()
    zero_pct = 100 * zero_count / len(self.returns)

    print(f"  NaN values: {nan_count}")
    print(f"  Zero returns: {zero_count} ({zero_pct:.2f}%)")

    if nan_count > 0:
        print("  WARNING: Data contains NaN values!")
    if zero_pct > 5:
        print(f"  WARNING: High percentage of zero returns ({zero_pct:.1f}%)")

A zero-return percentage above 5% often indicates data quality issues: weekend gaps not properly filtered, or a symbol with very low liquidity where the price does not move between ticks.

The Pipeline Orchestration

Each step is run via a dedicated runner script (run_step1.py, run_step2.py, etc.) that loads the previous step's output from a pickle file, runs the current step, and saves its output for the next step. This design allows you to re-run any individual step without repeating the entire pipeline. For example, run_step1.py loads data, runs the fractality check, and saves the FractalityChecker object:

# run_step1.py

def main():
    # Load data
    loader = DataLoader(
        symbol=config.SYMBOL,
        start_date=config.START_DATE,
        end_date=config.END_DATE,
        verbose=True
    )

    if config.MT5_ENABLED:
        loader.load_from_mt5()
    else:
        csv_path = Path(config.DATA_DIR) / f"{config.SYMBOL}_{config.START_DATE}_{config.END_DATE}.csv"
        loader.load_csv(csv_path)

    loader.calculate_returns(price_column='close', method='log')
    loader.summary_statistics()

    returns = loader.get_returns_array()
    checker = run_fractality_check(returns, save_plots=True)

    # Save for Step 2
    checker_path = Path(config.OUTPUT_DIR) / "step1_checker.pkl"
    with open(checker_path, 'wb') as f:
        pickle.dump(checker, f)

Step 1: Testing for Fractality

Before we can build an MMAR model, we must answer a prerequisite question: does our data exhibit multifractal behavior at all? If the answer is no, the entire construction falls apart. Applying the MMAR to non-fractal data would be like fitting a polynomial to a straight line: you introduce unnecessary complexity with no benefit. The test uses the partition function, a tool from statistical physics that measures how the moments of returns scale across different time intervals.

The Partition Function

The partition function, following Zhang's equation (4.1), is defined as:

S_q(T, delta_t) = (1/N) * sum |r(i * delta_t, delta_t)|^q

where r(i * delta_t, delta_t) is the return over a non-overlapping interval of length delta_t, q is the moment order, and N is the number of complete non-overlapping intervals. Let us unpack what this actually measures.

For a given time scale delta_t, we divide the entire return series into non-overlapping blocks of delta_t consecutive returns and sum the returns within each block (using log-return additivity) to get one "aggregated return" per block. We then compute the mean of the absolute aggregated returns raised to the power q. This gives us a single number S_q for each combination of q and delta_t. We then sweep across many delta_t values (from 1 to 150 observations) and many q values (from 0.01 to 30.0).

The fractal hypothesis is this: if the data is fractal, then log(S_q) should be a linear function of log(delta_t) for each q. That is, S_q should follow a power law in delta_t. The slope of this linear relationship (on the log-log plot) gives us tau(q) + 1, where tau(q) is the scaling function we extract in Step 2.

Why Non-Overlapping Intervals?

The use of non-overlapping intervals is critical, and Zhang explicitly emphasizes this. Overlapping windows would introduce artificial autocorrelation: adjacent windows share most of their data points, so their aggregated returns are mechanically correlated. This inflates the R-squared of the log-log regression, making any data look more fractal than it actually is. Non-overlapping windows ensure statistical independence between observations, giving us an honest test.

Why Mean Instead of Sum?

Another subtle but important choice is using the mean of |r|^q rather than the sum. As delta_t increases, N (the number of intervals that fit) decreases. If we used the sum, S_q would decline simply because there are fewer terms, confounding the true scaling relationship with a sample-size effect. Using the mean isolates the genuine scaling behavior of the moments.

Here is the complete implementation:

# step1_check_fractality.py

class FractalityChecker:
    def __init__(self, returns, delta_t_values=None, q_values=None,
                 verbose=config.VERBOSE):
        self.returns = returns
        self.verbose = verbose

        self._validate_returns()

        if delta_t_values is None:
            self.delta_t_values = config.generate_delta_t_values()
        else:
            self.delta_t_values = delta_t_values

        if q_values is None:
            self.q_values = config.generate_q_values()
        else:
            self.q_values = q_values

        # Storage for results
        self.partition_values = {}  # {q: {delta_t: S_q value}}
        self.r_squared_values = {}  # {q: R-squared}
        self.slopes = {}  # {q: slope (which equals tau(q) + 1)}
        self.intercepts = {}  # {q: intercept}

The constructor accepts the return series and optionally custom delta_t and q grids. The _validate_returns() method checks for NaN/Inf values, flags short series (less than 1,000 observations), and warns about high zero-return percentages. It raises a ValueError if the data contains NaN or infinite values, because these would silently corrupt every downstream computation.

The partition function calculation itself is vectorized for performance:

def calculate_partition_function(self, q, delta_t):
    # Number of complete NON-overlapping intervals (floor division)
    N = len(self.returns) // delta_t

    if N <= 0:
        return np.nan

    # Trim returns to fit complete intervals
    trimmed_returns = self.returns[:N * delta_t]

    # Reshape into (N, delta_t) and sum along axis 1
    # This gives us N non-overlapping interval returns
    interval_returns = trimmed_returns.reshape(N, delta_t).sum(axis=1)

    # MEAN of |r|^q (removes N-dependency as delta_t changes)
    S_q = np.mean(np.abs(interval_returns) ** q)

    return S_q

The key vectorization trick is the reshape(N, delta_t).sum(axis=1) pattern. Instead of looping over each interval and manually summing returns, we reshape the trimmed return array into a 2D matrix where each row is one non-overlapping interval. We then sum across columns in a single NumPy operation. This is both cleaner and orders of magnitude faster than explicit loops for large datasets.

Let us trace through a concrete example. Suppose we have 1,000 returns and delta_t = 50. Then N = 1000 // 50 = 20. We take the first 1,000 returns, reshape them into a 20 x 50 matrix, sum each row to get 20 aggregated returns, take their absolute values, raise to the power q, and compute the mean. That mean is S_q(T, 50) for that particular q.

After computing S_q for all combinations of q and delta_t (about 60 x 30 = 1,800 combinations), we fit a linear regression on the log-log plot for each q:

def fit_partition_plots(self):
    for q in self.q_values:
        delta_t_list = []
        S_q_list = []

        for delta_t in sorted(self.partition_values[q].keys()):
            S_q = self.partition_values[q][delta_t]
            if not np.isnan(S_q) and S_q > 0:
                delta_t_list.append(delta_t)
                S_q_list.append(S_q)

        log_delta_t = np.log10(delta_t_list)
        log_S_q = np.log10(S_q_list)

        # OLS regression
        slope, intercept, r_value, p_value, std_err = stats.linregress(
            log_delta_t, log_S_q
        )

        # slope = tau(q) + 1
        self.slopes[q] = slope
        self.intercepts[q] = intercept
        self.r_squared_values[q] = r_value ** 2

The slope of each regression line gives us tau(q) + 1, where tau(q) is the scaling function we extract in Step 2. The R-squared tells us how well the power-law model fits the data for that particular moment order q. The key output at this stage is the R-squared values across all q. Zhang achieved a mean R-squared of 0.66 for his stock data; our threshold is 0.60. A mean R-squared above 0.60 confirms that moment scaling exists, which is the hallmark of a fractal process.

The checker also includes a crossover detection method, which finds the time scale where high-frequency scaling behavior breaks down. At very short timescales, microstructure effects (bid-ask bounce, discrete tick sizes) can distort the scaling relationship. The crossover detector systematically tests truncated scaling regions, starting from each delta_t and computing the mean R-squared for all q values in that truncated region:

def detect_crossover_point(self, r_squared_threshold=0.85):
    self.crossover_scores = {}
    sorted_dt = sorted(self.delta_t_values)

    for start_idx, candidate_dt in enumerate(sorted_dt):
        candidate_dts = sorted_dt[start_idx:]
        if len(candidate_dts) < 3:
            break

        r2_values = []
        for q in self.q_values:
            # Re-fit regression using only delta_t >= candidate_dt
            # ...
            _, _, r_value, _, _ = stats.linregress(log_delta_t, log_S_q)
            r2_values.append(r_value ** 2)

        if r2_values:
            self.crossover_scores[candidate_dt] = np.mean(r2_values)

    # Find first delta_t where truncated R-squared exceeds threshold
    for dt, mean_r2 in self.crossover_scores.items():
        if mean_r2 >= r_squared_threshold:
            return dt

    return None

If a crossover is detected at, say, delta_t = 5, it means the scaling relationship is strongest for delta_t >= 5 and breaks down below that. This is informative for understanding the data's fractal structure at different scales.

Conclusion

In this article, we laid the groundwork for the Multifractal Model of Asset Returns. We introduced the theoretical foundations of both GARCH and the MMAR, loaded EURUSD 5-minute data via the MetaTrader5 Python API, and ran the first critical test: does our data exhibit multifractal scaling at all? Here is what we established:

• GARCH has structural blind spots. Despite its elegance and widespread use, GARCH assumes thin tails, has only short (exponentially decaying) memory, and operates on a single time scale. These are not minor quibbles but fundamental limitations that real market data consistently violates.
• The MMAR addresses all of them at once. The compound process X(t) = B_H[theta(t)] decomposes returns into Fractional Brownian Motion (long memory) and multifractal trading time (fat tails, clustering). This produces all four stylized facts of financial returns from a single architectural choice, not through parameter patching.
• Fractality is testable, not assumed. The partition function analysis in Step 1 provides a quantitative test for multifractal scaling. If the data fails this test, the MMAR should not be used. This is an honest prerequisite that GARCH-based approaches do not have (and arguably should).

The key takeaway is that this testing phase is not optional overhead. Applying the MMAR to data that is not fractal would produce worse forecasts than GARCH while being more expensive to compute. The partition function serves as an honest gate: it tells you whether there is any fractal structure worth modeling before you invest the effort.

In Part 2 of this series, we will extract the scaling function tau(q) from these partition function results, estimate the Hurst exponent H, run a rigorous multifractality decision test, and fit the multifractal spectrum to four theoretical distributions. Those fitted parameters will be the bridge between the analysis we started here and the MMAR construction that follows in Part 3.

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
run_step1.py	Runner: loads data, runs fractality check, saves FractalityChecker for Step 2
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