Hawkes Process Order Flow: Self-Exciting Intensity Models for Market Microstructure

Market orders do not arrive like raindrops. Anyone who has stared at a trade tape knows that a quiet minute can be followed by a burst of fifty prints in two seconds, and that the burst tends to feed on itself: one aggressive buy lifts the offer, algorithms react, and more buys follow. A homogeneous Poisson model cannot represent this, because it assumes the next event is independent of everything that came before. The Hawkes process order flow model fixes exactly that flaw with one idea: every event temporarily raises the probability of further events.

This matters in 2026 because event-driven data is now the default input for execution, market-making and surveillance systems, and the question “how much of this activity is genuinely new information, and how much is the market echoing itself?” has a direct, estimable answer: the branching ratio. This post builds the model from the intensity definition upward, derives the O(n) likelihood recursion, works one numeric example by hand, extends it to buy and sell streams, and ships runnable numpy code that I executed and checked. It closes with the pitfalls that wreck real fits and a production estimation architecture.

What this covers: definition and the exponential kernel, stability and the branching ratio, likelihood and EM estimation, a verified worked example, multivariate buy/sell flow, applications, failure modes, and a deployment blueprint.

Disclaimer: this is quantitative and systems analysis for education. It is not investment advice, and nothing here recommends trading any instrument.

Context and Background

Point-process thinking entered finance through a simple observation: the object of interest in high-frequency data is the sequence of event times, not a regularly sampled price. Early duration models such as the autoregressive conditional duration (ACD) family of Engle and Russell treated inter-arrival times as the variable to forecast. They work well for a single stream, but they are awkward when several event types interact, and they are defined in event time rather than calendar time, which complicates anything that must be priced in seconds.

The Hawkes process, introduced by Alan Hawkes in 1971 as a “self-exciting” and “mutually exciting” point process, models the instantaneous arrival rate directly in continuous time. Ogata’s 1981 paper on Lewis’ simulation method gave the thinning algorithm that most practitioners still use to simulate it, and his earlier likelihood work made maximum likelihood estimation routine. For a modern survey aimed squarely at this audience, see Bacry, Mastromatteo and Muzy, “Hawkes processes in finance”, which reviews definitions, volatility estimation from transaction data, market stability, systemic-risk contagion, optimal execution and limit order book dynamics.

Two empirical papers frame why the model took hold. Filimonov and Sornette, “Quantifying reflexivity in financial markets” (Physical Review E), fitted a Hawkes model to E-mini S&P 500 futures between 1998 and 2010 and reported that the share of price changes resulting from revealed exogenous information fell from about 70% in 1998 to under 30% since 2007, which is to say the endogenous, self-excited share rose sharply (branching ratio around 0.3 early on and roughly 0.7 to 0.8 after 2007, per the paper body). Hardiman, Bercot and Bouchaud, “Critical reflexivity in financial markets” (European Physical Journal B), used a power-law kernel on the same contract over 1998 to 2011 and concluded that the kernel integrates to about one regardless of period, meaning markets sit close to criticality. These two readings disagree on whether reflexivity trended upward, which is a lesson in itself and returns in the pitfalls section: the answer depends heavily on the kernel family you impose.

For readers building the data layer, two sibling posts matter. Hawkes fits are only as good as the event timestamps, so ITCH order book reconstruction systems architecture covers how to get trustworthy nanosecond events with correct aggressor tagging. And order book imbalance features for microstructure models covers the state variables that Hawkes intensities are often combined with.

The state of practice, as I read it, is a split. Academic and research desks use Hawkes models for diagnosis (how reflexive is this market, how fast does a shock decay) and for simulation (stress-testing execution algorithms against realistic clustering). Production trading systems use cheaper proxies, such as exponentially weighted event counts, that are in fact the intensity of a restricted Hawkes process without anyone calling it that. Recognizing that equivalence is the first step toward doing it properly.

How a Hawkes Process Models Order Flow: Intensity, Kernel and Branching Ratio

A Hawkes process is a self-exciting point process whose conditional intensity equals a baseline rate plus a sum of decaying contributions from every past event. Each trade temporarily raises the chance of further trades. The branching ratio, the average number of direct “children” per event, measures how reflexive the flow is, and the process is stable only when it stays below one.

The conditional intensity

Let the event times be t_1 < t_2 < … and let N(t) count them. The conditional intensity is the instantaneous arrival rate given the history H_t:

lambda(t) = mu + sum over t_i < t of phi(t – t_i)

Here mu > 0 is the baseline (immigrant) rate, representing flow caused by information and decisions outside the model, and phi >= 0 is the excitation kernel, which says how much one past event lifts today’s rate as a function of its age. The most used choice is the exponential kernel, phi(s) = alpha * exp(-beta * s). Each event adds a jump of size alpha to the intensity, and the excess then decays with rate beta. The half-life of an excitation is ln(2)/beta.

The exponential kernel has a decisive engineering property: it is Markovian. The pair (lambda(t), time) is a sufficient state, so you never need to store the event history to continue computing. That is why the likelihood runs in O(n) rather than O(n^2), and it is why exponentially weighted event counters in production code are secretly Hawkes intensities.

The branching ratio and stability

Integrating the kernel gives the branching ratio, n = integral of phi(s) ds, which equals alpha/beta for the exponential kernel. It has two equivalent readings. In the cluster (Hawkes-Oakes) representation, immigrants arrive as a Poisson process with rate mu, and every event, immigrant or descendant, independently spawns a Poisson number of direct children with mean n. In the intensity representation, n is the fraction of the current event rate that is caused by earlier events.

The process is stationary only if n < 1. The stationary mean rate is mu / (1 – n). The expected total size of a cluster, counting the immigrant itself, is 1 / (1 – n). With n = 0.5 a typical shock yields two events. With n = 0.9 it yields ten. With n = 0.99 it yields one hundred, and at n = 1 the expected cascade diverges, so activity could sustain itself with no outside input. That divergence is what “criticality” means in the empirical literature cited above.

The endogenous fraction of all events is therefore simply n, because in steady state a fraction (1 – n) of events are immigrants. When someone says “the market is 70% reflexive,” they mean an estimate of n of about 0.7, under whichever kernel they fit. Two caveats apply immediately: n is a model parameter, not a measured fact, and its estimate depends on what you count as an event.

Figure 1: How a Hawkes process generates clustered order flow. Baseline flow creates immigrants, each event raises the intensity, and offspring arrive with mean n per parent.

Figure 1 shows the generative story. Exogenous flow produces immigrant events at rate mu. Each event, wherever it comes from, adds a kernel to the intensity, and the intensity then spawns offspring events whose number is Poisson with mean n. The loop from offspring back to the kernel jump is the self-excitation. The two boxes at the bottom contrast the stable regime, where clusters die out, with the near-critical regime, where they run long.

Why order flow, not prices

Applying the model to signed market orders rather than to mid-price changes is a deliberate choice. Order flow has a well-documented long memory: the sign series of market orders is strongly autocorrelated for long lags in many markets, which the metaorder-splitting literature attributes largely to institutions slicing large parents into many child orders. A Hawkes process with a short exponential kernel will not reproduce long memory by itself. A sum of two or three exponentials, or a power-law kernel, can approximate it, and that is the first design decision in practice.

There is also a mechanical reason to use event data. A price change requires a trade that consumes a level, so trade arrivals are the more primitive object, and the arrival times carry information that fixed-interval bars discard. An interval of one second with five trades looks the same as one with five trades in the last 10 milliseconds once you aggregate, but their intensities, and their implied impact on the book, are very different.

Kernels beyond the exponential

Three kernel families cover most use cases. The single exponential is cheap, Markovian and identifiable, but it forces a single memory timescale. A sum of exponentials, phi(s) = sum of alpha_k * exp(-beta_k * s), keeps the Markov property at the cost of one state variable per term, and often captures both a fast reaction (milliseconds, algorithmic) and a slow one (tens of seconds, human or metaorder). The power-law kernel, phi(s) = alpha / (c + s)^(1 + gamma), reproduces the heavy tails seen empirically, as in the roughly -1.15 and -1.45 exponents Hardiman and colleagues report across two time ranges for the E-mini contract. It is not Markovian, so exact computation costs O(n^2) unless you approximate it with a sum of exponentials, which is the standard trick.

The choice matters for n. Fit a single exponential to data with a heavy-tailed true kernel and the optimizer will pick the dominant short timescale, often understating n. Fit a power law with a poorly chosen cutoff and you can estimate n suspiciously close to one. The disagreement between early reflexivity studies is partly this effect.

Deeper Analysis: Likelihood, the O(n) Recursion and a Worked Example

The log-likelihood

For a point process observed on [0, T] with events t_1..t_n, the log-likelihood is the sum of log-intensities at the events minus the integrated intensity (the compensator):

log L = sum over i of log lambda(t_i-) minus integral from 0 to T of lambda(s) ds

The first term rewards the model for having a high intensity exactly when events occur. The second term penalizes it for having a high intensity when nothing happens. The balance between them is what pins down alpha and beta: a large alpha with a large beta produces spikes that decay before the next event, while a small alpha with a small beta produces a smooth elevated rate.

For the exponential kernel the compensator has a closed form. Each event at t_i contributes (alpha/beta) * (1 – exp(-beta * (T – t_i))), so the integral is mu * T plus the sum of those terms. The intensity just before event i is lambda(t_i-) = mu + alpha * A_i, where A_i = sum over j < i of exp(-beta * (t_i – t_j)). Ogata’s recursion removes the double sum:

A_1 = 0, and A_i = exp(-beta * (t_i – t_{i-1})) * (1 + A_{i-1}).

This is the whole trick. Each event updates one scalar, so a full-sample evaluation costs O(n) time and O(1) memory beyond the data, and gradients with respect to the three parameters follow similarly. A naive double loop for one million events would need on the order of 5 x 10^11 pair evaluations; the recursion needs one million.

Figure 2: The recursion behind O(n) intensity estimation. The state A decays between events, contributes a log-intensity term at each event, then increments by one.

Figure 2 shows the same loop as a message flow. The feed handler supplies an event, the state decays by the exponential of minus beta times the gap, the likelihood accumulates log(mu + alpha * A), the state increments, and the compensator closes the account at the end of the window. The same loop runs online with no change, which is the bridge from research code to a live signal.

A worked example you can check by hand

The numbers below are illustrative parameters, not estimates from any market. Take mu = 0.2 events per second, alpha = 0.9 and beta = 1.5 per second. Then n = alpha/beta = 0.6, the stationary rate is 0.2 / (1 – 0.6) = 0.5 events per second, the expected cluster size is 1 / 0.4 = 2.5 events, and the excitation half-life is ln(2)/1.5 = 0.462 seconds.

Observe three events at t = 0, 1.0 and 1.4 seconds over a window T = 2 seconds.

Event 1 at t = 0: A_1 = 0, so lambda = 0.2 and log lambda = -1.6094. Immediately afterwards the intensity jumps to 0.2 + 0.9 = 1.1.

Event 2 at t = 1.0: the gap is 1.0, so A_2 = exp(-1.5) * (1 + 0) = 0.2231. The intensity just before the event is 0.2 + 0.9 * 0.2231 = 0.4008, and log lambda = -0.9143. A single second after the first trade, the excitation has decayed from 0.9 to 0.2008 above baseline.

Event 3 at t = 1.4: the gap is 0.4, so A_3 = exp(-0.6) * (1 + 0.2231) = 0.5488 * 1.2231 = 0.6713. The intensity is 0.2 + 0.9 * 0.6713 = 0.8041, and log lambda = -0.2180. Two events close together have compounded: the rate is now four times the baseline.

The sum of log-intensities is -2.7418. The compensator is mu * T = 0.4 plus 0.6 times [(1 – e^-3.0) + (1 – e^-1.5) + (1 – e^-0.9)] = 0.6 * (0.9502 + 0.7769 + 0.5934) = 0.6 * 2.3205 = 1.3923, for a total of 1.7923. The log-likelihood is -2.7418 – 1.7923 = -4.5341. I confirmed this figure by running the code below, which returns -4.53398 after rounding differences in the hand arithmetic.

Two facts fall out of this example that matter later. First, the intensity at an event depends on all the earlier gaps only through one scalar, which is why streaming implementations are trivial. Second, a burst such as events two and three raises the intensity at the moment events arrive, which is exactly the signal the optimizer uses to push alpha upward, while the compensator term stops it from pushing alpha without limit.

Runnable code: simulation and maximum likelihood

The following numpy and scipy code simulates a Hawkes process by Ogata thinning and estimates the three parameters. It uses a log parametrization so mu, alpha and beta stay positive. I ran it with seed 7 on a 5,000-second window of the parameters above.

import numpy as np
from scipy.optimize import minimize

def simulate(mu, alpha, beta, T, rng):
    """Ogata thinning for a univariate exponential Hawkes process."""
    t, events, lam_ex = 0.0, [], 0.0   # lam_ex = excitation part at time t
    while True:
        lam_bar = mu + lam_ex            # upper bound until next event
        w = rng.exponential(1.0 / lam_bar)
        t_new = t + w
        if t_new > T:
            return np.array(events)
        lam_ex *= np.exp(-beta * w)
        if rng.uniform() * lam_bar <= mu + lam_ex:
            events.append(t_new)
            lam_ex += alpha
        t = t_new

def neg_loglik(theta, t, T):
    mu, alpha, beta = np.exp(theta)
    A, ll, prev = 0.0, np.log(mu), t[0]
    for i in range(1, len(t)):
        A = np.exp(-beta * (t[i] - prev)) * (1.0 + A)
        ll += np.log(mu + alpha * A)
        prev = t[i]
    comp = mu * T + (alpha / beta) * np.sum(1.0 - np.exp(-beta * (T - t)))
    return -(ll - comp)

rng = np.random.default_rng(7)
t = simulate(0.2, 0.9, 1.5, 5000.0, rng)
res = minimize(neg_loglik, np.log([0.1, 0.5, 1.0]), args=(t, 5000.0),
               method="Nelder-Mead")
mu_h, a_h, b_h = np.exp(res.x)
print(len(t), mu_h, a_h, b_h, a_h / b_h)

On my run this produced 2,509 events (the stationary expectation is 2,500) and estimates mu = 0.193, alpha = 0.891, beta = 1.448, giving n = 0.615 against a true 0.6. That is one seed on a synthetic sample, not a benchmark, but it shows the estimator behaves when the model is correctly specified. The pure Python loop is slow for millions of events; in production compile the loop with numba or write it in C++ or Rust.

EM estimation and the latent-parent view

Maximum likelihood on the Hawkes surface can be flat in the direction of high alpha and high beta, so many practitioners prefer expectation-maximization, the nonparametric flavor popularized for point processes by Lewis and Mohler. The latent variable is the parent of each event. In the E-step, event i is an immigrant with probability mu / lambda(t_i), and has event j as its parent with probability alpha * exp(-beta * (t_i – t_j)) / lambda(t_i). In the M-step, mu is the expected immigrant count over T, n is the expected number of parent links divided by the number of events, and beta is the parent-link count divided by the sum of probability-weighted lags.

def em_exp(t, T, iters=200, window=30.0, init=(0.1, 0.5, 1.0)):
    mu, alpha, beta = init
    n = len(t)
    lo = np.searchsorted(t, t - window)          # truncate old parents
    I = np.repeat(np.arange(n), np.arange(n) - lo)
    J = np.concatenate([np.arange(lo[i], i) for i in range(n)])
    dt = t[I] - t[J]
    for _ in range(iters):
        g = alpha * np.exp(-beta * dt)
        denom = mu + np.bincount(I, g, minlength=n)
        p_bg = mu / denom                          # immigrant probability
        p_ij = g / denom[I]                        # parent probabilities
        mu = p_bg.sum() / T
        beta = p_ij.sum() / np.sum(p_ij * dt)
        alpha = (p_ij.sum() / n) * beta
    return mu, alpha, beta

On a 2,000-second synthetic sample (1,010 events) with the same true parameters, this returned mu = 0.191, alpha = 0.893 and beta = 1.434, so n = 0.623. The window truncation caps cost at the number of events within the window times n; with a 30-second window and a half-life under half a second the truncation error is negligible, but with a heavy-tailed kernel it is not. EM is monotone, needs no learning rate, and gives the cluster interpretation for free, which is why I reach for it first when exploring unfamiliar data. Its drawbacks are slow convergence near the boundary n close to 1 and a quadratic worst case when the window is large.

Checking the fit: time rescaling

Any estimate needs a goodness-of-fit test, and the standard one is the random time change theorem. If the model is correct, the transformed gaps tau_i = integral of lambda between t_{i-1} and t_i are independent unit-rate exponentials. Compute them with the same recursion, then draw a QQ plot against Exp(1) and run a Kolmogorov-Smirnov test plus a Ljung-Box test on the tau_i for residual autocorrelation. Misspecified kernels typically show up as excess short residual gaps (the model under-predicts bursts) or as serial correlation at lags one through ten. A fit that passes the marginal KS test but fails the Ljung-Box test is a model that has the right average clustering but the wrong timing.

Multivariate Hawkes for Buy and Sell Order Flow

From one stream to a matrix of influences

A single stream of unsigned trades throws away the most interesting question in order flow: does a buy beget buys, or does it beget sells? The multivariate extension answers it by giving each event type its own intensity. For buy and sell market orders indexed 1 and 2:

lambda_i(t) = mu_i + sum over j of the sum over events of type j before t of phi_ij(t – t_k)

With exponential kernels phi_ij(s) = alpha_ij * exp(-beta_ij * s), the model has a 2 x 2 matrix of excitation strengths. Dividing by the decay rates gives the branching matrix N, whose entry n_ij = alpha_ij / beta_ij is the expected number of type-i children produced by one type-j parent. The diagonal captures self-excitation (buys following buys, the footprint of order splitting and momentum chasing). The off-diagonal captures cross-excitation (sells reacting to buys, typically market makers and arbitrageurs leaning against the flow, or stop-loss cascades after a move).

Stability generalizes cleanly: the process is stationary if and only if the spectral radius of N is below one. The stationary mean-rate vector is Lambda = (I – N)^-1 mu.

Figure 3: A two-dimensional Hawkes process for buy and sell market orders. Self-excitation loops keep each side clustered, cross terms let one side provoke the other, and the spectral radius of N sets stability.

Figure 3 draws the two loops and the two cross links. Both streams feed an imbalance and price-impact stage, which is where the model connects to the signals described in the order book imbalance post.

A symmetric numeric example

Again these are illustrative values. Let N = [[0.45, 0.15], [0.15, 0.45]] and mu = (0.1, 0.1) events per second for buys and sells. The eigenvalues of N are 0.60 (eigenvector (1, 1)) and 0.30 (eigenvector (1, -1)), so the spectral radius is 0.60 and the system is stable. On the symmetric vector, (I – N) acts as 1 – 0.60 = 0.40, so each side has a stationary rate of 0.1 / 0.4 = 0.25 events per second, a 2.5-fold amplification of the baseline.

The eigen-decomposition says more than the stability test. The (1, 1) mode, with eigenvalue 0.6, is total activity: buys and sells rising together. The (1, -1) mode, with eigenvalue 0.3, is imbalance: a buy surplus. Because the imbalance mode has the smaller eigenvalue, a one-sided burst dies out faster than a two-sided one. In this example, imbalance shocks decay with an effective multiplier of 1 / (1 – 0.3) = 1.43 events per parent, while activity shocks multiply by 1 / (1 – 0.6) = 2.5. In real data it is plausible, though I have not verified it for any specific venue, that the imbalance mode is the more persistent one, because persistent one-sided execution is what produces long sign autocorrelation. If a fit shows that ordering, it is the quantitative fingerprint of metaorder splitting.

A caution on interpretation: sums of cross terms are not causal. A positive n_21 says that sells tend to follow buys within the kernel timescale. It does not tell you why, and in a market with a single large metaorder being worked, the pattern can come from the counterparty side of the same sequence.

Fitting the multivariate model

The O(n) recursion extends by keeping one state variable per (target, source) pair. With K types and one exponential per pair, the state is K x K numbers, updated by decay at every event and incremented in the column of the event type that occurred. The cost is O(n * K^2) per pass, which for K = 2 is trivial and for K = 8 (buy and sell at four book levels, say) is still tractable. The parameter count is 3K^2 + K for free (alpha, beta) pairs, which grows quickly, so most production fits share a single beta per target row or per source column.

import numpy as np

def branching_matrix(alpha, beta):
    """alpha, beta: KxK arrays. Returns N, spectral radius, mean rates helper."""
    N = alpha / beta
    rho = np.max(np.abs(np.linalg.eigvals(N)))
    return N, rho

def stationary_rates(N, mu):
    return np.linalg.solve(np.eye(len(mu)) - N, mu)

N = np.array([[0.45, 0.15], [0.15, 0.45]])
mu = np.array([0.1, 0.1])
print(np.linalg.eigvals(N))          # [0.6 0.3]
print(stationary_rates(N, mu))       # [0.25 0.25]

Model selection across K is where overfitting hides. Likelihood always improves when you add types or exponentials, so use the Bayesian information criterion or out-of-sample likelihood on later windows, never in-sample likelihood. Sparse penalties (an L1 penalty on the off-diagonal terms) and low-rank structure on N are the two regularizers that work in practice.

Applications: What the Intensity Is Used For

Order flow clustering and regime flags

The intensity itself is the cleanest definition of order flow clustering. A threshold on the ratio lambda(t) / Lambda, where Lambda is the stationary rate, flags bursts adaptively across quiet and busy sessions. Because lambda is a function of the event history alone, it is computable the instant a trade prints and does not wait for a bar to close. A market-making system can use a rising intensity to widen quotes or reduce size before realized volatility has registered, which is the same logic as the volatility estimators described in the finance survey cited earlier.

Volatility and price-impact modeling

The finance survey by Bacry and coauthors reviews how mutually exciting models for upward and downward price moves connect the branching structure to realized volatility and to microstructure effects such as the signature plot. A standard result for the event count helps with intuition: the long-horizon variance of the number of events per unit time is mu / (1 – n)^3. That cubic dependence is why a small change in n near one has an outsized effect on how erratic the activity is. It is a statement about counts; mapping it to price variance needs the signs and sizes of the moves, so treat it as intuition, not a calibrated volatility model.

Liquidity, execution and stress testing

On the execution side, a Hawkes baseline lets you simulate realistic impact of your own orders: your child order is an event that raises the intensity of other participants’ orders, so a schedule that bunches child orders is self-penalizing. This gives a principled reason to randomize and spread slices, and the survey above covers optimal-execution formulations built on exactly this mechanism. For stress tests, a simulator calibrated at n = 0.7 can be pushed to n = 0.9 to see how a strategy behaves in a more reflexive world, which is a controlled way to test the tail rather than waiting for one.

Surveillance and quote stuffing

For surveillance, a multivariate fit on messages (adds, cancels, executes) gives a baseline of normal cross-excitation. Bursts that violate the fitted cross-excitation structure, such as high add intensity that does not provoke any matching execution response, are candidates for review. This is an anomaly flag, not a legal conclusion, and any such flag needs a human review path.

Production Estimation Architecture

A research notebook and a production fit differ mainly in three things: timestamps, windows and monitoring. The architecture in Figure 4 separates a slow path that fits parameters from a fast path that evaluates intensity with those parameters.

Figure 4: A production layout for intensity estimation. A windowed fit job writes parameters to a registry, an online recursion consumes them, and time-rescaling diagnostics trigger refits.

The feed handler and book reconstruction stage, covered in the ITCH architecture post, produce the event store. Two details here decide model quality. Aggressor tagging must be correct, because mislabelled signs corrupt the cross-excitation matrix first. And timestamps must come from one clock domain, since a few hundred microseconds of skew between venues can create spurious cross-excitation at the shortest lags.

The fit job runs on a rolling or expanding window, for example the previous five trading days restricted to continuous trading hours and with the opening and closing auctions removed. It emits mu, alpha, beta, n and the spectral radius into a versioned parameter registry. The online process, which holds only the K x K state described earlier, reads the latest parameters and updates intensity per event in constant time. A single core handles millions of updates per second in compiled code, so intensity evaluation is almost never the bottleneck; the feed is.

The diagnostics stage reruns time rescaling on the most recent hour of data using the live parameters. When the KS statistic or Ljung-Box statistic crosses a threshold, or when the instantaneous ratio of observed to predicted event counts drifts, a drift alarm triggers an early refit. This closed loop is the practical answer to non-stationarity, discussed next, because a daily fit alone will be wrong at 09:31 and at 15:55 every day.

Intraday seasonality

Market activity follows a U-shaped daily pattern, with heavy volumes at the open and close and a lull around midday. If you fit a constant mu to a full day, the model attributes the seasonal rise to self-excitation and inflates n. Two fixes are standard. One is to replace mu with a deterministic piecewise-constant or spline function mu(t) estimated from multi-day averages. The other is to fit within narrow time-of-day buckets. Either way, report n alongside the baseline treatment you used, because two papers with different seasonal handling are not comparable.

Sizing and cost

The cost arithmetic is straightforward. A liquid futures contract might print on the order of tens of thousands to a few hundred thousand trades per day, depending on the product and the day, which is well within a single-core O(n) fit that completes in seconds when compiled. An equity universe of several thousand names fits embarrassingly in parallel. The expensive cases are multivariate fits across many types with power-law kernels, and message-level fits on full depth-of-book feeds with hundreds of millions of events a day. For those, the sum-of-exponentials approximation and a GPU or vectorized batch likelihood become necessary. These sizes are order-of-magnitude judgments, not measurements from a specific venue.

Trade-offs, Gotchas, and What Goes Wrong

Kernel misspecification

The most common failure is imposing a single exponential on data with several timescales. The fit will land on a compromise: a beta that is too slow for the algorithmic reaction and too fast for the metaorder tail. The residual signature is serial correlation in the rescaled gaps and a branching ratio that moves a lot when you change the window length. The test for it is cheap: refit with two exponentials and see whether the likelihood improves enough to survive the information criterion, and whether n shifts materially. If n moves from 0.55 to 0.8 when you allow a second timescale, the original number was not a measurement.

The near-critical identifiability problem

When the true n is close to one, the likelihood is nearly flat along directions that trade off mu against n. Estimates of n then have a wide confidence band even with millions of events, and small amounts of non-stationarity or a mis-specified baseline push the estimate toward one. Hardiman and colleagues’ conclusion that the kernel integrates to about one is simultaneously the most cited and the most fragile finding, because it sits exactly where the estimator is weakest. Report confidence intervals by block bootstrap or by the observed Fisher information, and avoid headlines like “the market is at criticality” from a point estimate.

Endogeneity is a modeling label

Branching structure is a statistical decomposition, not causal evidence. Two independent traders reacting to the same public news at slightly different latencies look like a parent and a child. An exchange matching-engine rule that splits one order into several prints creates fake clusters. A common exogenous driver, such as an index future moving, will appear as excitation in a single-instrument model. If you need causal claims, you need instruments: exogenous shocks such as scheduled announcements, or cross-venue designs that separate the common factor.

Non-stationarity and event definition

Parameters drift across days, volatility regimes and tick-size changes. A model fitted in a calm week and applied unchanged to a news day will under-predict. Event definition matters just as much: counting every trade print versus every aggressive order (merging fills from one sweep) changes n substantially, because a sweep of five levels is one decision but five prints. Pick the unit that corresponds to a decision, and state it.

Discretization and timestamp resolution

Exchange timestamps are rounded, and many events share an identical timestamp. A continuous-time likelihood assumes distinct times. Break ties by sequence number, add small randomized jitter for sensitivity checks, and be wary of any fit whose parameters depend on that choice. Latency between your capture point and the matching engine also adds a shift that can blur the fastest self-excitation timescale.

Computation traps

The recursion is numerically stable, but exp(-beta * gap) underflows harmlessly to zero while log(mu + alpha * A) is fine. The real traps are optimizer related: unconstrained parameters wandering into n >= 1 (use a log parametrization or a penalty), local optima in multi-exponential fits (initialize from a grid or run EM first), and degeneracy where beta grows without bound and alpha with it, collapsing the kernel to a delta spike. When that last case happens, the data resolution is coarser than the excitation timescale, and the fit is telling you the dynamics are faster than your timestamps.

Practical Recommendations

Start simple and escalate only when diagnostics demand it. Fit a univariate exponential model on a clean, seasonally adjusted event stream, run the time-rescaling checks, and look at the stability of n across windows before adding anything. Only if the residuals show structure should you move to two or three exponentials, then to a signed bivariate model. Treat each added component as a hypothesis with a cost.

Report uncertainty every time. A branching ratio without an interval is a trivia number. Compare in-sample and out-of-sample likelihood, and keep the parameter registry versioned so any signal derived from intensity can be reproduced. Keep the simulator next to the estimator: simulating from fitted parameters and re-estimating is the fastest bug detector for the full pipeline, and the check I used for every figure quoted in this post.

Finally, be clear about what the model is for. As a diagnostic and a simulator, it is robust and informative. As a stand-alone forecasting signal, it competes with simpler decayed counters and with richer state-dependent models, and the extra machinery has to earn its place on out-of-sample data.

Clean timestamps, one clock domain, correct aggressor signs, auctions removed.
Model the intraday baseline mu(t) before reading n.
Fit with a log parametrization; add EM as an initializer or cross-check.
Run time-rescaling QQ, KS and Ljung-Box on held-out windows.
Test a second timescale; keep it only if the information criterion agrees.
Report n with a bootstrap interval, and state the event definition.
Monitor drift and refit on alarm, not on the calendar alone.

Frequently Asked Questions

What is a Hawkes process in simple terms?

A Hawkes process is a model for events that cluster in time. It has a baseline arrival rate, and every event adds a temporary boost to the rate of future events that fades with a decay time. Because past events raise the future rate, the process is called self-exciting. In finance it describes bursts of trades, quote updates, or price jumps far better than a Poisson process, which treats every moment as independent of the last.

What is the branching ratio and why must it be below one?

The branching ratio n is the average number of direct follow-on events triggered by one event, equal to the integral of the kernel. If n is below one, each cluster eventually dies out, the average cluster has 1 / (1 – n) events, and the process has a finite stationary rate of mu / (1 – n). At n of one or more, the expected cascade does not terminate and the rate would grow without outside input, so the model has no stable long-run behavior.

How do you estimate a Hawkes process from order flow data?

Two methods dominate. Maximum likelihood maximizes the sum of log intensities minus the integrated intensity, and with an exponential kernel it runs in O(n) time through a recursion on a single state variable. Expectation-maximization treats each event’s parent as hidden, alternates between parent probabilities and parameter updates, and is stable and monotone. In both cases you should adjust for intraday seasonality and validate with time-rescaling residuals.

Is a Hawkes process better than a Poisson model for trades?

For clustered data, almost always, because a Poisson process cannot generate bursts, and real trade arrivals are heavily clustered with long memory. The comparison to make is on out-of-sample log-likelihood and on rescaled-residual diagnostics, not on in-sample fit, since Hawkes has more parameters. Be aware that part of what Hawkes calls self-excitation can be seasonality or common news that a time-varying Poisson baseline would also capture.

Can Hawkes processes predict price direction?

Not by themselves in any dependable way I can verify. The intensity predicts how active the market will be in the next moments, and a bivariate buy and sell fit predicts the likely direction of the next order flow imbalance over short horizons, but turning that into a profitable trade must overcome spreads, fees and adverse selection. Treat Hawkes outputs as inputs to risk and execution decisions, not as an alpha claim. This is analysis, not advice.

Which Python libraries implement Hawkes processes?

The tick library (from the Paris-based quantitative finance community behind much of the Hawkes literature) provides simulation and several estimators, including nonparametric ones, though its maintenance status should be checked before adopting it. For small projects the hand-rolled numpy and scipy code above is enough, and for large production loads a compiled implementation in numba, C++ or Rust is preferable. I did not benchmark any library for this post.