### Install hawkesbook

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Install the package using pip.

```bash
pip install hawkesbook
```

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

### Fit Exponential Hawkes Process with EM Algorithm

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Estimate parameters for an exponential Hawkes process using the Expectation-Maximization algorithm, supporting custom starting values and progress tracking.

```python
import hawkesbook as hawkes
import numpy as np

arrivals = np.array([0.5, 1.1, 1.3, 2.8, 4.2, 5.1, 5.3, 7.9, 9.2])
T = 10.0

# Basic EM fit
theta_em = hawkes.exp_em(arrivals, T, iters=100)
print(f"EM estimate: lambda={theta_em[0]:.4f}, alpha={theta_em[1]:.4f}, beta={theta_em[2]:.4f}")

# EM with custom starting values
theta_start = np.array([0.5, 1.0, 2.0])
theta_em_custom = hawkes.exp_em(arrivals, T, theta_start=theta_start, iters=100)
print(f"EM with custom start: {theta_em_custom}")

# EM with progress tracking and likelihood history
theta_em_verbose, ll_history = hawkes.exp_em(
    arrivals, T,
    theta_start=np.array([1.0, 2.0, 3.0]),
    iters=50,
    verbosity=10,
    calcLikelihoods=True
)
print(f"\nFinal EM estimate: {theta_em_verbose}")
print(f"Likelihood improved from {ll_history[0]:.4f} to {ll_history[-1]:.4f}")

# Compare MLE and EM
theta_mle = hawkes.exp_mle(arrivals, T)
ll_mle = hawkes.exp_log_likelihood(arrivals, T, theta_mle)
ll_em = hawkes.exp_log_likelihood(arrivals, T, theta_em)

print(f"\nMLE log-likelihood: {ll_mle:.4f}")
print(f"EM log-likelihood: {ll_em:.4f}")
print(f"Better fit: {'MLE' if ll_mle > ll_em else 'EM'}")
```

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

### Simulate Hawkes Process Realizations

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Generate Hawkes process arrival times using thinning and composition algorithms for both exponential and power-law kernels.

```python
import hawkesbook as hawkes
import numpy as np

# Set random seed for reproducibility
hawkes.numba_seed(42)

theta = np.array([1.0, 2.0, 3.0])  # lambda, alpha, beta
T = 100.0  # Observation period

# Simulate using thinning algorithm (exponential kernel)
arrivals_thinning = hawkes.exp_simulate_by_thinning(theta, T)
print(f"Thinning simulation: {len(arrivals_thinning)} arrivals")
print(f"First 5 arrivals: {arrivals_thinning[:5]}")

# Simulate using composition method (fixed number of events)
N = 100  # Number of events to simulate
arrivals_composition = hawkes.exp_simulate_by_composition(theta, N)
print(f"\nComposition simulation: {N} arrivals")
print(f"Time span: [0, {arrivals_composition[-1]:.2f}]")

# Alternative composition (stops at time T)
arrivals_comp_alt = hawkes.exp_simulate_by_composition_alt(theta, T)
print(f"\nAlt composition: {len(arrivals_comp_alt)} arrivals in [0, {T}]")

# Simulate power-law Hawkes
theta_power = np.array([1.0, 1.0, 2.0, 3.0])  # lambda, k, c, p
arrivals_power = hawkes.power_simulate_by_thinning(theta_power, T)
print(f"\nPower-law simulation: {len(arrivals_power)} arrivals")
```

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

### Simulate and Fit Mutually-Exciting Hawkes Processes

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Demonstrates multivariate Hawkes process simulation, intensity calculation, log-likelihood estimation, and MLE parameter recovery.

```python
import hawkesbook as hawkes
import numpy as np

# Set random seed
np.random.seed(42)

# Define 2-dimensional mutually-exciting Hawkes parameters
m = 2  # Number of processes

# Base intensities for each process
lambda_vec = np.array([0.5, 0.3])

# Alpha matrix: alpha[i,j] = jump in process j when event occurs in process i
alpha = np.array([
    [0.3, 0.1],  # Effect of process 0 on [process 0, process 1]
    [0.2, 0.4]   # Effect of process 1 on [process 0, process 1]
])

# Decay rates for each process
beta = np.array([2.0, 1.5])

theta = (lambda_vec, alpha, beta)

# Simulate mutually-exciting process
T = 100.0
events = hawkes.mutual_exp_simulate_by_thinning(theta, T)
times = np.array([e[0] for e in events])
ids = np.array([e[1] for e in events])

print(f"Simulated {len(events)} total events")
print(f"Process 0: {np.sum(ids == 0)} events")
print(f"Process 1: {np.sum(ids == 1)} events")

# Compute intensity at a specific time
t_query = 50.0
mask = times < t_query
intensity = hawkes.mutual_exp_hawkes_intensity(
    t_query, times[mask], ids[mask], theta
)
print(f"\nIntensity at t={t_query}: process 0 = {intensity[0]:.4f}, process 1 = {intensity[1]:.4f}")

# Compute log-likelihood
ll = hawkes.mutual_exp_log_likelihood(times, ids, T, theta)
print(f"Log-likelihood: {ll:.4f}")

# Compute compensators for model validation
compensators = hawkes.mutual_exp_hawkes_compensators(times, ids, theta)
print(f"\nCompensators shape: {compensators.shape}")

# MLE for mutual Hawkes (starting from true parameters for demo)
theta_start = (lambda_vec * 1.1, alpha * 0.9, beta * 1.1)
theta_mle, ll_mle = hawkes.mutual_exp_mle(times, ids, T, theta_start)
print(f"\nMLE log-likelihood: {ll_mle:.4f}")
print(f"Estimated lambda: {theta_mle[0]}")
print(f"Estimated alpha:\n{theta_mle[1]}")
print(f"Estimated beta: {theta_mle[2]}")
```

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

### Fit and analyze Hawkes processes for earthquake data

Source: https://github.com/pat-laub/hawkesbook/blob/main/README.md

Demonstrates loading earthquake data, performing MLE and EM parameter estimation, comparing models using BIC, and generating a Q-Q plot for model validation.

```python
import hawkesbook as hawkes

import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
from statsmodels.graphics.gofplots import qqplot

# Load data to fit
quakes = pd.read_csv("japanese-earthquakes.csv")
quakes.index = pd.to_datetime(quakes.Day.astype(str) + "/" + quakes.Month.astype(str) + "/" + quakes.Year.astype(str) + " " + quakes.Time, dayfirst=True)
quakes.sort_index(inplace=True)

# Calculate each arrival as a (fractional) number of days since the
# beginning of the observation period
timeToQuake = quakes.index - pd.Timestamp("1/1/1973")
ts = np.array(timeToQuake.total_seconds() / 60 / 60 / 24)

# Calculate the length of the observation period
obsPeriod = pd.Timestamp("31/12/2020") - pd.Timestamp("1/1/1973")
T = obsPeriod.days

# Calculate the maximum likelihood estimate for the Hawkes process
# with an exponentially decaying intensity
𝛉_exp_mle = hawkes.exp_mle(ts, T)
print("Exp Hawkes MLE fit: ", 𝛉_exp_mle)

# Calculate the EM estimate or the same type of Hawkes process
𝛉_exp_em = hawkes.exp_em(ts, T, iters=100)
print("Exp Hawkes EM fit: ", 𝛉_exp_mle)

# Get the likelihoods of each fit to find the better one
ll_mle = hawkes.exp_log_likelihood(ts, T, 𝛉_exp_mle)
ll_em = hawkes.exp_log_likelihood(ts, T, 𝛉_exp_em)

if ll_mle > ll_em:
	print("MLE was a better fit than EM in this case")
	𝛉_exp = 𝛉_exp_mle
	ll_exp = ll_mle
else:
	print("EM was a better fit than MLE in this case")
	𝛉_exp = 𝛉_exp_em
	ll_exp = ll_em

# Fit instead the Hawkes with a power-law decay
𝛉_pl = hawkes.power_mle(ts, T)
ll_pl = hawkes.power_log_likelihood(ts, T, 𝛉_pl)

# Compare the BICs
BIC_exp = 3 * np.log(len(ts)) - 2 * ll_exp
BIC_pl = 4 * np.log(len(ts)) - 2 * ll_pl
if BIC_exp < BIC_pl:
	print(f"The exponentially-decaying Hawkes was the better fit with BIC={BIC_exp:.2f}.")
	print(f"The power-law Hawkes had BIC={BIC_pl:.2f}.")
else:
	print(f"The power-law Hawkes was the better fit with BIC={BIC_pl:.2f}.")
	print(f"The exponentially-decaying Hawkes had BIC={BIC_exp:.2f}.")

# Create a Q-Q plot for the exponential-decay fit by
# first transforming the points to a unit-rate Poisson
# process as outlined by the random time change theorem
tsShifted = hawkes.exp_hawkes_compensators(ts, 𝛉_exp)
iat = np.diff(np.insert(tsShifted, 0, 0))
qqplot(iat, dist=stats.expon, fit=False, line="45")
plt.show()
```

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

### Analyze Earthquake Data with Hawkes Models

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Full workflow for fitting exponential and power-law Hawkes models to earthquake data, comparing models via BIC, and validating with Q-Q plots.

```python
import hawkesbook as hawkes
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
from statsmodels.graphics.gofplots import qqplot

# Load earthquake data
quakes = pd.read_csv("japanese-earthquakes.csv")
quakes.index = pd.to_datetime(
    quakes.Day.astype(str) + "/" + quakes.Month.astype(str) + "/" +
    quakes.Year.astype(str) + " " + quakes.Time,
    dayfirst=True
)
quakes.sort_index(inplace=True)

# Convert to arrival times (days since start)
timeToQuake = quakes.index - pd.Timestamp("1/1/1973")
ts = np.array(timeToQuake.total_seconds() / 60 / 60 / 24)

# Observation period
obsPeriod = pd.Timestamp("31/12/2020") - pd.Timestamp("1/1/1973")
T = obsPeriod.days

print(f"Analyzing {len(ts)} earthquakes over {T} days")

# Fit exponential Hawkes via MLE
theta_exp_mle = hawkes.exp_mle(ts, T)
print(f"\nExp Hawkes MLE: lambda={theta_exp_mle[0]:.6f}, alpha={theta_exp_mle[1]:.6f}, beta={theta_exp_mle[2]:.6f}")

# Fit exponential Hawkes via EM
theta_exp_em = hawkes.exp_em(ts, T, iters=100)
print(f"Exp Hawkes EM:  lambda={theta_exp_em[0]:.6f}, alpha={theta_exp_em[1]:.6f}, beta={theta_exp_em[2]:.6f}")

# Compare fits
ll_mle = hawkes.exp_log_likelihood(ts, T, theta_exp_mle)
ll_em = hawkes.exp_log_likelihood(ts, T, theta_exp_em)

if ll_mle > ll_em:
    print("MLE was a better fit")
    theta_exp = theta_exp_mle
    ll_exp = ll_mle
else:
    print("EM was a better fit")
    theta_exp = theta_exp_em
    ll_exp = ll_em

# Fit power-law Hawkes
theta_pl = hawkes.power_mle(ts, T)
ll_pl = hawkes.power_log_likelihood(ts, T, theta_pl)
print(f"\nPower-law Hawkes: lambda={theta_pl[0]:.6f}, k={theta_pl[1]:.6f}, c={theta_pl[2]:.6f}, p={theta_pl[3]:.6f}")

# Model comparison via BIC
BIC_exp = 3 * np.log(len(ts)) - 2 * ll_exp
BIC_pl = 4 * np.log(len(ts)) - 2 * ll_pl

print(f"\nModel Comparison:")
print(f"  Exponential BIC: {BIC_exp:.2f}")
print(f"  Power-law BIC: {BIC_pl:.2f}")
print(f"  Better model: {'Exponential' if BIC_exp < BIC_pl else 'Power-law'}")

# Goodness-of-fit via Q-Q plot
tsShifted = hawkes.exp_hawkes_compensators(ts, theta_exp)
iat = np.diff(np.insert(tsShifted, 0, 0))

fig, ax = plt.subplots(figsize=(8, 6))
qqplot(iat, dist=stats.expon, fit=False, line="45", ax=ax)
ax.set_title("Q-Q Plot: Transformed Inter-arrival Times vs Exp(1)")
plt.savefig("earthquake_qqplot.png", dpi=150)
plt.show()

print("\nIf model fits well, points should lie on the 45-degree line")
```

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

### Simulate Hawkes Process via Inverse Compensator

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Uses the inverse compensator method to generate arrival times for a generic Hawkes process.

```python
mu = lambda x: 2.0 * np.exp(-3.0 * x)
M = lambda t: (2.0/3.0) * (1 - np.exp(-3.0 * t))
theta_generic = (1.0, mu, M)

import numpy.random as rnd
rnd.seed(42)
arrivals_generic = hawkes.simulate_inverse_compensator(theta_generic, hawkes.hawkes_compensator, N=50)
print(f"\nInverse compensator simulation: {len(arrivals_generic)} arrivals")
print(f"First 5 arrivals: {arrivals_generic[:5]}")
```

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

### Maximum Likelihood Estimation

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Fit Hawkes process parameters using MLE. Includes variants with gradient and Hessian support for faster convergence.

```python
import hawkesbook as hawkes
import numpy as np

# Observed arrival times and observation period
arrivals = np.array([0.5, 1.1, 1.3, 2.8, 4.2, 5.1, 5.3, 7.9, 9.2])
T = 10.0  # Observation period length

# MLE for exponential decay Hawkes process
# Returns [lambda, alpha, beta]
theta_exp_mle = hawkes.exp_mle(arrivals, T)
print(f"Exponential Hawkes MLE: lambda={theta_exp_mle[0]:.4f}, alpha={theta_exp_mle[1]:.4f}, beta={theta_exp_mle[2]:.4f}")

# MLE with gradient (faster convergence)
theta_grad = hawkes.exp_mle_with_grad(arrivals, T)
print(f"MLE with gradient: {theta_grad}")

# MLE with Hessian (even faster for large datasets)
theta_hess = hawkes.exp_mle_with_hess(arrivals, T)
print(f"MLE with Hessian: {theta_hess}")

# MLE for power-law decay Hawkes process
# Returns [lambda, k, c, p]
theta_power_mle = hawkes.power_mle(arrivals, T)
print(f"Power-law Hawkes MLE: lambda={theta_power_mle[0]:.4f}, k={theta_power_mle[1]:.4f}, c={theta_power_mle[2]:.4f}, p={theta_power_mle[3]:.4f}")

# Custom starting values
theta_start = np.array([0.5, 1.0, 2.0])
theta_custom = hawkes.exp_mle(arrivals, T, theta_start=theta_start)
print(f"MLE with custom start: {theta_custom}")
```

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

### Compute Log-Likelihoods and Model Comparison

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Calculate log-likelihoods for exponential and power-law Hawkes processes and compare models using the Bayesian Information Criterion (BIC).

```python
ll_exp = hawkes.exp_log_likelihood(arrivals, T, theta_exp)
ll_power = hawkes.power_log_likelihood(arrivals, T, theta_power)

print(f"Exponential Hawkes log-likelihood: {ll_exp:.4f}")
print(f"Power-law Hawkes log-likelihood: {ll_power:.4f}")

# Model comparison using BIC
n = len(arrivals)
BIC_exp = 3 * np.log(n) - 2 * ll_exp    # 3 parameters
BIC_power = 4 * np.log(n) - 2 * ll_power  # 4 parameters

print(f"\nBIC (exponential): {BIC_exp:.4f}")
print(f"BIC (power-law): {BIC_power:.4f}")
print(f"Better model: {'Exponential' if BIC_exp < BIC_power else 'Power-law'}")

# Generic log-likelihood function (for custom kernels)
mu = lambda x: theta_exp[1] * np.exp(-theta_exp[2] * x)
M = lambda t: (theta_exp[1]/theta_exp[2]) * (1 - np.exp(-theta_exp[2]*t))
theta_generic = (theta_exp[0], mu, M)

ll_generic = hawkes.log_likelihood(
    arrivals, T, theta_generic,
    hawkes.hawkes_intensity, hawkes.hawkes_compensator
)
print(f"\nGeneric log-likelihood: {ll_generic:.4f}")
```

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

### Fit Exponential Hawkes Process with GMM

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Perform parameter estimation using the Generalized Method of Moments by matching empirical and theoretical moments.

```python
import hawkesbook as hawkes
import numpy as np

# Generate arrivals for testing
arrivals = np.array([0.5, 1.1, 1.3, 2.8, 4.2, 5.1, 5.3, 7.9, 9.2, 10.5, 12.1, 15.3, 18.2, 20.0])
T = 25.0

# Compute empirical moments
tau = 5.0  # Bin width
lag = 2    # Lag for autocovariance
emp_moments = hawkes.empirical_moments(arrivals, T, tau, lag)
print(f"Empirical moments (mean, var, autocov): {emp_moments.flatten()}")

# Theoretical moments for given parameters
theta_test = np.array([1.0, 2.0, 3.0])
theo_moments = hawkes.exp_moments(theta_test, tau, lag)
print(f"Theoretical moments for theta={theta_test}: {theo_moments.flatten()}")

# GMM estimation
theta_gmm = hawkes.exp_gmm(arrivals, T, tau=5, lag=2, iters=2)
print(f"\nGMM estimate: lambda={theta_gmm[0]:.4f}, alpha={theta_gmm[1]:.4f}, beta={theta_gmm[2]:.4f}")

# Compare with MLE
theta_mle = hawkes.exp_mle(arrivals, T)
print(f"MLE estimate: lambda={theta_mle[0]:.4f}, alpha={theta_mle[1]:.4f}, beta={theta_mle[2]:.4f}")
```

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

### Log-Likelihood Computation

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Evaluate the log-likelihood of observed data under different Hawkes models.

```python
import hawkesbook as hawkes
import numpy as np

arrivals = np.array([0.5, 1.1, 1.3, 2.8, 4.2, 5.1, 5.3, 7.9, 9.2])
T = 10.0

# Fit both models
theta_exp = hawkes.exp_mle(arrivals, T)
theta_power = hawkes.power_mle(arrivals, T)
```

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

### Compute Compensators

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Calculate integrated intensity for model validation using the random time change theorem.

```python
import hawkesbook as hawkes
import numpy as np

# Arrival times
arrivals = np.array([0.5, 1.2, 2.1, 3.5, 4.8])
theta_exp = np.array([1.0, 2.0, 3.0])  # lambda, alpha, beta

# Compute compensators at each arrival time (for exponential Hawkes)
compensators = hawkes.exp_hawkes_compensators(arrivals, theta_exp)
print("Compensators at each arrival:")
for i, (t, comp) in enumerate(zip(arrivals, compensators)):
    print(f"  t={t:.1f}: Lambda={comp:.4f}")

# For power-law Hawkes
theta_power = np.array([1.0, 1.0, 2.0, 3.0])  # lambda, k, c, p
power_compensators = hawkes.power_compensators(arrivals, theta_power)
print(f"\nPower-law compensators: {power_compensators}")

# Inter-arrival times in transformed space should be Exp(1) if model is correct
iat = np.diff(np.insert(compensators, 0, 0))
print(f"\nTransformed inter-arrival times (should be ~Exp(1)): {iat}")
```

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

### Compute Conditional Intensity

Source: https://context7.com/pat-laub/hawkesbook/llms.txt

Calculate the conditional intensity of a Hawkes process at a specific time using event history. Supports generic, exponential, and power-law kernels.

```python
import hawkesbook as hawkes
import numpy as np

# Event history (arrival times before time t)
event_history = [0.5, 1.2, 2.1]

# Generic Hawkes intensity with custom excitation function
# Parameters: (lambda, mu_function, M_function)
lambda_base = 1.0
mu = lambda x: 2.0 * np.exp(-3.0 * x)  # Exponential decay kernel
M = lambda t: (2.0/3.0) * (1 - np.exp(-3.0*t))  # Integrated kernel

theta = (lambda_base, mu, M)
intensity = hawkes.hawkes_intensity(t=2.5, H_t=event_history, theta=theta)
print(f"Generic Hawkes intensity at t=2.5: {intensity:.4f}")

# Exponential Hawkes intensity (optimized version)
# Parameters: [lambda, alpha, beta] where mu(t) = alpha * exp(-beta * t)
theta_exp = np.array([1.0, 2.0, 3.0])  # lambda=1, alpha=2, beta=3
exp_intensity = hawkes.exp_hawkes_intensity(t=2.5, H_t=event_history, theta=theta_exp)
print(f"Exponential Hawkes intensity at t=2.5: {exp_intensity:.4f}")

# Power-law Hawkes intensity
# Parameters: [lambda, k, c, p] where mu(t) = k / (c + t)^p
theta_power = np.array([1.0, 1.0, 2.0, 3.0])
power_intensity = hawkes.power_hawkes_intensity(t=2.5, H_t=np.array(event_history), theta=theta_power)
print(f"Power-law Hawkes intensity at t=2.5: {power_intensity:.4f}")
```

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

### Define Hawkes process conditional intensity in Python

Source: https://github.com/pat-laub/hawkesbook/blob/main/README.md

Uses Unicode characters to represent mathematical symbols directly in the function definition.

```python
def hawkes_intensity(t, ℋ_t, 𝛉):
    λ, μ, _ = 𝛉
    λˣ = λ
    for t_i in ℋ_t:
        λˣ += μ(t - t_i)
    return λˣ
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.