Chapter 13: Advanced Simulation Techniques

The simulation algorithms presented in earlier chapters are correct but not always efficient. This chapter covers advanced techniques: the branching representation for Hawkes processes, vectorized O(n) computation, strategies for large-scale simulation, and a validation framework that ensures your code is correct.

13.1 Branching Simulation for Hawkes Processes

An alternative to Ogata’s thinning is the cluster (branching) simulation, which directly exploits the Galton-Watson branching structure of the Hawkes process:

Algorithm:

1. Simulate background (immigrant) events: t_bg ~ HPP(μ) on [0, T]
2. For each immigrant, generate offspring:
   • Number of offspring: K ~ Poisson(n*)
   • Delays from parent: d ~ φ(·)/n* (normalized kernel as a density)
3. For each offspring, generate their offspring recursively
4. Collect all events, sort by time

For the exponential kernel φ(t) = α·exp(−βt):

• n* = α/β
• Offspring delays: d ~ Exponential(β) (since φ(t)/n* = β·exp(−βt))

def simulate_hawkes_branching(mu, alpha, beta, T):
    n_star = alpha / beta
    all_events = []

    def generate_offspring(t_parent):
        k = np.random.poisson(n_star)
        delays = np.random.exponential(1.0 / beta, size=k)
        for d in delays:
            t_child = t_parent + d
            if t_child < T:
                all_events.append(t_child)
                generate_offspring(t_child)

    # Background events
    n_bg = np.random.poisson(mu * T)
    bg_times = np.random.uniform(0, T, n_bg)
    for t in bg_times:
        all_events.append(t)
        generate_offspring(t)

    return np.sort(all_events)

Comparison with Ogata thinning:

• Branching simulation: exact, easy to implement, reveals cluster structure
• Ogata thinning: exact, works for any kernel (not just exponential)
• Both have the same asymptotic cost; branching avoids the rejection sampling overhead

13.2 Vectorized O(n) Implementation

The recursive computation of Rᵢ = exp(−β·Δtᵢ) · (1 + Rᵢ₋₁) is inherently sequential — each step depends on the previous. However, the compensator computation and MLE gradient can be parallelized.

For large n (millions of events), avoid Python loops by using numba JIT compilation:

from numba import njit

@njit
def hawkes_loglik_fast(events, T, mu, alpha, beta):
    n = len(events)
    R = 0.0
    loglik = 0.0
    comp = 0.0

    for i in range(n):
        if i > 0:
            dt = events[i] - events[i-1]
            R = np.exp(-beta * dt) * (1.0 + R)
        intensity = mu + alpha * R
        loglik += np.log(intensity)
        comp += 1.0 - np.exp(-beta * (T - events[i]))

    loglik -= mu * T + (alpha / beta) * comp
    return loglik

With numba, this runs ~100x faster than a pure Python loop for n = 100,000 events.

13.3 Computational Complexity

The O(n²) vs O(n) difference is critical: for n = 10,000, naive computation takes ~100× longer.

Truncation for long histories: If the event sequence is very long, the exponential kernel decays rapidly. For events more than τ_cutoff = 5/β time units in the past, exp(−β·τ_cutoff) < 0.007 — negligible. We can safely truncate the sum:

# Only include events within a window of 5/beta
def hawkes_intensity_truncated(t, past_events, mu, alpha, beta):
    window = 5.0 / beta
    recent = past_events[past_events > t - window]
    return mu + alpha * np.sum(np.exp(-beta * (t - recent)))

13.4 Simulation of Very Long Sequences

For very long T (e.g., simulating years of seismicity):

1. Stationary initialization: Don’t start from empty history (there is a transient). Run for a burn-in period and discard the first T_burn events.

2. Memory management: Store only the last 5/β seconds of history at any time — older events don’t contribute.

3. Parallel independent runs: For uncertainty estimation, simulate K independent realizations in parallel using multiprocessing or joblib.

from joblib import Parallel, delayed

results = Parallel(n_jobs=4)(
    delayed(simulate_hawkes_branching)(mu, alpha, beta, T)
    for _ in range(K)
)

13.5 Validation Framework

Every simulation function should be validated by checking that estimates from simulated data recover the known parameters:

def validate_simulator(true_params, n_sim=100, T=200):
    mu, alpha, beta = true_params
    estimates = []
    for _ in range(n_sim):
        events = simulate_hawkes_branching(mu, alpha, beta, T)
        est = fit_hawkes(events, T).x
        estimates.append(est)
    estimates = np.array(estimates)

    print("True:    ", true_params)
    print("Mean est:", estimates.mean(axis=0))
    print("Bias:    ", estimates.mean(axis=0) - true_params)
    print("Std err: ", estimates.std(axis=0))
    # Check ~95% coverage
    z = (estimates - true_params) / estimates.std(axis=0)
    print("Coverage:", (np.abs(z) < 1.96).mean(axis=0))

This detects bugs in the simulation or estimation code and measures finite-sample bias.

See code/13_advanced_simulation.py for benchmarks comparing Ogata vs. branching simulation, vectorized implementations, and a full validation harness.

13.6 Key Takeaways

• The branching simulation is an exact alternative to Ogata thinning, revealing cluster structure directly.
• The recursive Rᵢ reduces log-likelihood computation from O(n²) to O(n).
• Use numba for large-scale simulation; truncate distant history for efficiency.
• Always validate simulations with a recover-parameters test: simulate → estimate → check bias and coverage.