Volatility Clustering with Merton-Hawkes Jump-Diffusion Simulations in Python

Standard Geometric Brownian Motion assumes constant volatility and normally distributed returns — a clean mathematical convenience that real markets routinely violate. In practice, large price moves cluster together: a spike in volatility today makes another spike tomorrow significantly more likely. This phenomenon, known as volatility clustering, is a defining feature of intraday microstructure across equities, ETFs, indices, and crypto assets. Capturing it accurately is not an academic exercise — it directly impacts risk management, option pricing, and Monte Carlo scenario generation.

This article implements a Merton jump-diffusion model enhanced with a Hawkes self-exciting point process to simulate realistic intraday price dynamics across a multi-asset universe including tech stocks, ETFs, major indices, and BTC-USD. We combine return bootstrapping from live market data fetched via yfinance with calibrated jump intensities that respond to their own history, producing simulation paths that replicate the fat tails and volatility bursts observed in one-minute return data. All code is runnable end-to-end in Python.

Most algo trading content gives you theory.

This gives you the code.

3 Python strategies. Fully backtested. Colab notebook included.

Plus a free ebook with 5 more strategies the moment you subscribe.

5,000 quant traders already run these:

This article covers:

• Section 1 — Core Concepts:** Explains GBM's limitations, the Merton jump-diffusion framework, and how the Hawkes process introduces self-exciting jump clustering without heavy mathematical prerequisite


• Section 2 — Python Implementation:** Full setup with configurable parameters (2.1), bootstrapped return calibration and Hawkes intensity simulation (2.2), multi-asset Monte Carlo path generation (2.3), and visualization of simulated price paths and jump intensities (2.4)


• Section 3 — Results and Analysis:** What the simulation reveals about cross-asset volatility behavior, jump frequency calibration, and realistic scenario coverage


• Section 4 — Use Cases:** Practical applications in risk management, options desk scenario generation, intraday strategy stress testing, and crypto volatility modeling


• Section 5 — Limitations and Edge Cases:** Honest discussion of calibration sensitivity, Hawkes stationarity requirements, data frequency constraints, and overfitting risk

1. From GBM to Jump-Diffusion with Self-Exciting Intensity

Geometric Brownian Motion is the workhorse of classical quantitative finance. It models an asset price as a continuous process driven by a constant drift and a Wiener process scaled by constant volatility. The famous Black-Scholes formula is built entirely on this foundation. The problem is structural: GBM produces log returns that are normally distributed and serially independent. Real intraday returns are neither. They exhibit fat tails — extreme moves far more frequently than a Gaussian would predict — and they cluster, meaning periods of turbulence beget more turbulence.

Robert Merton's jump-diffusion model addresses the fat-tail problem directly by adding a compound Poisson jump component to the standard GBM. Between jumps, the asset diffuses smoothly. Jumps arrive at a rate governed by a Poisson intensity parameter (lambda), and each jump size is drawn from a log-normal distribution with its own mean and variance. This immediately produces heavier tails in the return distribution. However, the classical Merton model still treats jumps as memoryless: the probability of a jump in the next instant is the same regardless of whether a jump just occurred. That is not what we observe intraday.

The Hawkes process solves the memory problem. It is a self-exciting point process in which each arriving event temporarily increases the rate at which future events arrive. Think of it like an earthquake model: a main shock raises the probability of aftershocks, which themselves raise the probability of further aftershocks, with the excitation decaying exponentially over time. Applied to financial jumps, a large price shock increases the intensity of subsequent jumps — precisely the clustering behavior we observe around earnings releases, macro announcements, and liquidity crises.

Combining Merton's jump sizes with Hawkes-driven jump timing produces a simulation framework that is both analytically tractable and empirically realistic. The conditional jump intensity at time t follows: λ(t) = λ₀ + α · Σ exp(−β · (t − tᵢ)) where λ₀ is the baseline intensity, α is the excitation magnitude, β is the decay rate, and the sum runs over all past jump times tᵢ. When α/β < 1, the process is stationary — critical for stable long-run simulations. This is the model we calibrate and simulate below.

2. Python Implementation

2.1 Setup and Parameters

The parameters below control every aspect of the simulation. ASSETS defines the multi-asset universe pulled from Yahoo Finance. INTERVAL and PERIOD set the intraday data granularity — one-minute bars over five days, the maximum window Yahoo permits at this resolution. The Hawkes parameters ALPHA and BETA govern excitation strength and decay; keeping ALPHA/BETA < 1 ensures stationarity. LAMBDA_0 is the baseline jump arrival rate per minute. N_STEPS and N_SIMS control the Monte Carlo workload.

CODE Kopieren

import numpy as np

import pandas as pd

import yfinance as yf

import matplotlib.pyplot as plt

import matplotlib.gridspec as gridspec

from scipy.stats import norm



# --- Universe and data config ---

ASSETS   = ["AAPL", "QQQ", "SPY", "^VIX", "BTC-USD"]

INTERVAL = "1m"

PERIOD   = "5d"



# --- Hawkes process parameters ---

LAMBDA_0 = 0.05    # baseline jump intensity (jumps per minute)

ALPHA    = 0.6     # excitation magnitude — how much each jump raises intensity

BETA     = 0.8     # decay rate — how quickly excitation fades

# Stationarity check: ALPHA / BETA must be < 1

assert ALPHA / BETA < 1, "Hawkes process is non-stationary — reduce ALPHA or increase BETA"



# --- Merton jump-diffusion parameters ---

JUMP_MEAN = 0.0    # mean log jump size

JUMP_STD  = 0.015  # std of log jump size (1.5% per jump)



# --- Monte Carlo config ---

N_STEPS = 390      # minutes in a standard US trading session

N_SIMS  = 200      # number of simulation paths per asset



np.random.seed(42)

: longer, more detailed articles + full Google Colab implementations for every approach.

Or get everything in one powerful package with or