Hawkes Process

Hawkes Process 霍克斯过程

霍克斯过程（Hawkes Process）是一类自激点过程（Self-exciting Point Process）：历史事件会提高未来事件的发生强度。它在金融高频交易、地震预测、社交网络传播、营销归因、风控欺诈检测等场景中非常常见。

读完本文，你应该能：

1. 写出条件强度函数（Conditional Intensity）和指数核（Exponential Kernel）的公式
2. 手撕 Ogata 稀疏化算法 做模拟
3. 手撕 最大似然估计（MLE） 做参数拟合
4. 理解 分支比（Branching Ratio） 和 Lift 的业务含义

---

从泊松过程说起

齐次泊松过程（Homogeneous Poisson Process）

在区间 $[0, T]$ 上，事件到达时刻记为 ${t_1, t_2, \dots, t_N}$。

齐次泊松过程的核心假设：任意时刻的事件发生率是常数 $\mu$。

• 强度函数：$\lambda(t) = \mu$
• 相邻事件间隔 $W_i = t_i - t_{i-1}$ 服从指数分布：$W_i \sim \text{Exp}(\mu)$
• 事件总数 $N(T) \sim \text{Poisson}(\mu T)$

模拟（手撕版）：不断采样间隔 $W \sim \text{Exp}(\mu)$，累加得到事件时刻。

为什么泊松过程不够用？

真实世界里，事件往往是成簇（clustering）的：

• 一条爆款微博发出后，短时间内转发/评论激增
• 一笔大单成交后，跟风单密集出现
• 一次营销活动后，转化事件在窗口期内扎堆

泊松过程的”无记忆性”无法刻画这种历史依赖。霍克斯过程通过在强度函数里加入”历史激励项”来解决这个问题。

---

霍克斯过程的数学定义

条件强度函数

一维霍克斯过程由基线强度 $\mu > 0$ 和激励核 $\phi(\cdot)$ 定义：

$$
\lambda(t) = \mu + \sum_{t_i < t} \phi(t - t_i)
$$

含义很直观：

符号	含义
$\mu$	外生基线：没有历史事件时，单位时间平均发生率
$\phi(s)$	核函数：一个发生在 $s$ 时间单位之前的事件，对当前强度的贡献
$t_i$	第 $i$ 个历史事件的发生时刻

直觉：每发生一次事件，就在时间轴上”叠加”一个衰减的激励波形；当前强度 = 基线 + 所有历史波形的叠加。

最常用的指数核

$$
\phi(s) = \alpha e^{-\beta s}, \quad s \geq 0
$$

参数	含义
$\alpha > 0$	激励幅度：每个子事件能激发多强的后续事件
$\beta > 0$	衰减速度：激励多快消失（$\beta$ 越大衰减越快）

代入后，条件强度变为：

$$
\lambda(t) = \mu + \sum_{t_i < t} \alpha e^{-\beta(t - t_i)}
$$

分支比（Branching Ratio）—— 系统是否稳定

$$
n = \int_0^{\infty} \phi(s), ds = \int_0^{\infty} \alpha e^{-\beta s}, ds = \frac{\alpha}{\beta}
$$

物理含义：一个事件平均会”繁衍”出 $n$ 个直接子事件。

条件	含义
$n < 1$	稳定：子事件会逐渐湮灭，过程有平稳分布
$n \in [0.5, 1)$	强自激，接近临界爆发边缘
$n \geq 1$	不稳定：事件会无限爆炸（雪崩），现实中通常意味着模型设定有问题

拟合时若 $n \geq 1$，优先检查：数据是否混入了外生冲击、核函数是否选错、时间窗口是否过长。

业务层的 Lift

在营销归因、事件触发分析中，常用 Lift 衡量”历史事件 X 对未来事件 Y 的触发效应”：

$$
\text{Lift} = \frac{P(Y \mid X)}{P(Y)}
$$

即”在发生事件 X 的激励下，短时间内发生事件 Y 的概率”除以”事件 Y 的无条件基线概率”。

经验上，Lift > 3.0 通常被认为是显著的强触发。在霍克斯框架下，Lift 与核函数参数 $\alpha, \beta$ 直接相关：$\alpha$ 越大、衰减越慢（$\beta$ 越小），Lift 越高。

---

核心递推：激励状态的 $O(1)$ 更新

模拟和 MLE 都依赖同一个技巧。定义激励状态：

$$
R(t) = \sum_{t_i < t} e^{-\beta(t - t_i)}
$$

则：

$$
\lambda(t) = \mu + \alpha R(t)
$$

关键递推：当新事件在 $t_k$ 发生时，

$$
R(t_k^+) = R(t_k^-) + 1
$$

在两次事件之间（$t \in [t_k, t_{k+1})$），只有指数衰减：

$$
R(t) = R(t_k), e^{-\beta(t - t_k)}
$$

这个 $O(1)$ 递推是后面所有代码的核心，避免每次计算 $\lambda(t)$ 都遍历全部历史事件。

---

模拟：Ogata 稀疏化算法（Thinning）

算法思想

霍克斯过程的强度 $\lambda(t)$ 随时间变化。Ogata 稀疏化的思路：

1. 找一个上界 $\bar{\lambda}$，使得 $\lambda(t) \leq \bar{\lambda}$ 在剩余时间内成立
2. 按速率 $\bar{\lambda}$ 生成一个候选事件时刻 $t^*$
3. 以概率 $\lambda(t^*) / \bar{\lambda}$ 接受，否则拒绝并继续

因为真实接受率被”稀释”了，所以叫 thinning（稀疏化）。

指数核下的上界

在事件 $t_k$ 之后、下一个事件到来之前，$R(t)$ 单调递减，因此 $\lambda(t)$ 也单调递减。上界可以取：

$$
\bar{\lambda} = \lambda(t_k^+) = \mu + \alpha R(t_k^+)
$$

手撕伪代码

输入: μ, α, β, T_max
t = 0, R = 0, events = []

while t < T_max:
    λ_bar = μ + α * R          # 当前上界
    if λ_bar <= 0: break

    # 按上界生成候选间隔
    dt = Exp(λ_bar)
    t = t + dt
    if t >= T_max: break

    # 衰减到候选时刻
    R = R * exp(-β * dt)
    λ = μ + α * R              # 真实强度

    # 接受-拒绝
  if Uniform(0,1) < λ / λ_bar:
        events.append(t)
        R = R + 1              # 新事件叠加激励

分支表示法（理解用，非必须实现）

霍克斯过程等价于一个分支过程：

• 以速率 $\mu$ 产生”移民”（immigrant）事件，它们不再繁衍
• 每个事件以核 $\phi$ 繁衍子事件，子事件继续繁衍……

当 $n < 1$ 时，总事件数有限。这个视角对理解”簇”的结构很有帮助，但工程上 Ogata 更直接。

---

参数估计：最大似然（MLE）

给定观测到的事件时刻 ${t_1, \dots, t_N}$，$T$ 为观测窗口右端点。

对数似然

点过程的对数似然：

$$
\ell(\mu, \alpha, \beta) = \sum_{i=1}^{N} \log \lambda(t_i) - \int_0^T \lambda(s), ds
$$

两项含义：

• $\sum \log \lambda(t_i)$：事件确实在”高概率时刻”发生 → 越大越好
• $-\int \lambda$：惩罚过高的总体强度 → 防止过拟合

指数核下的积分项（手推结论）

利用 $R(t)$ 的递推，积分项可以解析计算：

$$
\int_0^T \lambda(s), ds = \mu T + \frac{\alpha}{\beta} \sum_{i=1}^{N} \left(1 - e^{-\beta(T - t_i)}\right)
$$

MLE 实现要点

遍历事件，维护 $R(t)$：

# 似然的第一项
log_lik = sum(log(μ + α * R_at_event_i))

# 第二项（补偿项）
compensator = μ * T + (α / β) * sum(1 - exp(-β * (T - t_i)))
log_lik -= compensator

对 $\mu, \alpha, \beta$ 用 scipy.optimize.minimize 最小化 $-\ell$ 即可。注意约束：$\mu > 0, \alpha > 0, \beta > 0$，且建议检查 $n = \alpha/\beta < 1$。

---

多元霍克斯过程（简述）

当有 $D$ 种事件类型时，强度向量：

$$
\lambda_k(t) = \mu_k + \sum_{j=1}^{D} \sum_{t_i^{(j)} < t} \alpha_{kj}, e^{-\beta_{kj}(t - t_i^{(j)})}
$$

参数	含义
$\mu_k$	类型 $k$ 的基线强度
$\alpha_{kj}$	类型 $j$ 的事件对类型 $k$ 的激励强度
$\beta_{kj}$	对应的衰减速度

解读 $\alpha_{kj}$：$\alpha_{kj}$ 大 → “做了 j 之后容易发生 k”。营销上就是归因矩阵：渠道 j 对转化 k 的触发效应。

模拟和 MLE 的思路与一元情况相同，只是状态从标量 $R$ 变成矩阵，计算量 $O(D^2)$ per event。

---

应用场景与案例

金融高频交易

场景：限价单成交、撤单、价格跳动等事件在毫秒级成簇出现。

建模：

• 事件类型：买单成交 / 卖单成交
• $\alpha$ 大 → 存在明显的”订单流自相关”（order flow persistence）
• $n$ 接近 1 → 市场处于高波动、临界状态

用途：预测短期成交强度、检测异常交易簇、做市商报价策略。

营销归因与转化分析

场景：用户先”点击广告”（事件 A），再”加购”（事件 B），再”下单”（事件 C）。

建模：

• 多元霍克斯，$\alpha_{CB}$ 衡量”加购对下单的触发”
• Lift $= P(\text{下单} \mid \text{加购}) / P(\text{下单})$，与 $\alpha_{CB}$ 正相关

用途：判断哪个触点真正”激发”了转化，而非仅仅”先于”转化发生。

社交网络与舆情

场景：一条推文引发转发链，转发又引发更多转发。

建模：分支比 $n$ 衡量传播力；$n \to 1$ 时话题接近”病毒式”爆发。

用途：预测话题热度衰减、识别关键种子节点、评估干预效果。

地震余震

场景：主震之后余震在数天/数周内成簇衰减——这是霍克斯过程最经典的应用之一。

建模：$\beta$ 小 → 余震持续时间长；$\alpha/\beta$ 衡量主震”触发”余震的能力。

风控欺诈检测

场景：盗刷交易往往短时间内密集出现。

建模：对用户交易时间序列拟合霍克斯，$\lambda(t)$ 突然飙升 → 异常簇。

用途：实时风险评分，而非固定时间窗口计数。

---

完整手撕代码

以下代码只依赖 numpy 和 scipy，约 150 行，可直接运行。

import numpy as np
from scipy.optimize import minimize


# ============ 1. 模拟：Ogata Thinning ============

def simulate_hawkes(mu, alpha, beta, T_max, seed=42):
    """
    一维霍克斯过程（指数核）模拟。
    返回: 事件时刻数组
    """
    rng = np.random.default_rng(seed)
    t, R = 0.0, 0.0
    events = []

    while t < T_max:
        lam_bar = mu + alpha * R
        if lam_bar <= 0:
            break

        dt = rng.exponential(1.0 / lam_bar)
        t += dt
        if t >= T_max:
            break

        R *= np.exp(-beta * dt)
        lam = mu + alpha * R

        if rng.uniform() < lam / lam_bar:
            events.append(t)
            R += 1.0

    return np.array(events)


# ============ 2. 拟合：MLE ============

def log_likelihood(params, events, T):
    """负对数似然（用于最小化）"""
    mu, alpha, beta = params
    if mu <= 0 or alpha <= 0 or beta <= 0:
        return 1e10

    n = len(events)
    if n == 0:
        return mu * T

    R = 0.0
    log_sum = 0.0
    t_prev = 0.0

    for t_i in events:
        R *= np.exp(-beta * (t_i - t_prev))
        lam = mu + alpha * R
        if lam <= 0:
            return 1e10
        log_sum += np.log(lam)
        R += 1.0
        t_prev = t_i

    compensator = mu * T + (alpha / beta) * np.sum(
        1.0 - np.exp(-beta * (T - events))
    )
    return -(log_sum - compensator)


def fit_hawkes(events, T, init=(0.5, 0.3, 1.0)):
    """MLE 拟合，返回 μ, α, β"""
    result = minimize(
        log_likelihood,
        x0=np.array(init),
        args=(events, T),
        method="L-BFGS-B",
        bounds=[(1e-8, None), (1e-8, None), (1e-8, None)],
    )
    mu, alpha, beta = result.x
    return mu, alpha, beta, result.fun


# ============ 3. 预测：未来强度 ============

def intensity_at(t, events, mu, alpha, beta):
    """计算时刻 t 的条件强度 λ(t)"""
    R = 0.0
    for t_i in events:
        if t_i < t:
            R += np.exp(-beta * (t - t_i))
    return mu + alpha * R


# ============ 4. 业务指标 ============

def branching_ratio(alpha, beta):
    return alpha / beta


def empirical_lift(events, window, baseline_rate):
    """
    简化版 Lift：统计每个事件后 window 内是否还有事件。
    baseline_rate: 无条件的基线发生率（如 1/平均间隔）
    """
    if len(events) < 2:
        return 0.0
    triggered = 0
    for i, t_i in enumerate(events[:-1]):
        has_followup = np.any((events[i+1:] - t_i) <= window)
        if has_followup:
            triggered += 1
    p_cond = triggered / (len(events) - 1)
    return p_cond / baseline_rate if baseline_rate > 0 else 0.0


# ============ 5. 跑通示例 ============

if __name__ == "__main__":
    # --- 真实参数 ---
    mu_true, alpha_true, beta_true = 0.2, 0.6, 1.5
    T = 500.0

    # 模拟
    events = simulate_hawkes(mu_true, alpha_true, beta_true, T)
    print(f"模拟事件数: {len(events)}")
    print(f"真实分支比 n = {branching_ratio(alpha_true, beta_true):.3f}")

    # 拟合
    mu_hat, alpha_hat, beta_hat, nll = fit_hawkes(events, T)
    n_hat = branching_ratio(alpha_hat, beta_hat)
    print(f"\n拟合结果:")
    print(f"  μ: {mu_hat:.4f} (真实 {mu_true})")
    print(f"  α: {alpha_hat:.4f} (真实 {alpha_true})")
    print(f"  β: {beta_hat:.4f} (真实 {beta_true})")
    print(f"  分支比 n: {n_hat:.4f}")
    print(f"  负对数似然: {nll:.2f}")

    # 预测强度
    t_query = events[-1] + 0.5
    lam = intensity_at(t_query, events, mu_hat, alpha_hat, beta_hat)
    print(f"\n在 t={t_query:.2f} 的预测强度 λ(t) = {lam:.4f}")

    # Lift
    avg_gap = T / len(events) if len(events) > 0 else T
    baseline = 1.0 / avg_gap
    lift = empirical_lift(events, window=2.0, baseline_rate=baseline)
    print(f"经验 Lift (window=2.0): {lift:.2f}")

预期输出（随机种子 42，$T=500$）大致为：

模拟事件数: ~120
μ ≈ 0.2, α ≈ 0.6, β ≈ 1.5
分支比 n ≈ 0.4（稳定）

---

实战踩坑与进阶

常见坑

问题	原因	建议
$n \geq 1$	数据含外生冲击 / 核函数不匹配	换幂律核、分段建模、剔除异常簇
拟合不稳定	事件太少	至少几百个事件；用正则化 MLE
$\beta$ 极大	激励衰减太快，退化成泊松	检查时间单位是否统一（秒 vs 毫秒）
非平稳	基线 $\mu$ 随时间变化	用时变基线 $\mu(t)$ 或分时段拟合

核函数扩展

指数核简单但衰减太快。实务中常用：

• 幂律核：$\phi(s) = \frac{\alpha c}{(s+c)^{1+p}}$，适合长尾簇（地震、舆情）
• 和核 / 多指数核：$\phi(s) = \sum_k \alpha_k e^{-\beta_k s}$，多时间尺度叠加

模拟和 MLE 的递推公式会变复杂，但 Ogata 稀疏化的框架不变。

与相关方法的关系

方法	与霍克斯的关系
Cox 过程	强度本身是随机场，不是历史事件叠加
ARIMA / 计数回归	离散时间；霍克斯是连续时间
生存分析	霍克斯可看作”被历史激励的 hazard rate”
Transformer 点过程	用神经网络替代手工核 $\phi$，适合复杂场景

推荐库（不必重复造轮子）

• Python: tick（Hawkes 专用，含多种核和求解器）
• Python: hawkeslib
• R: hawkes

手撕一遍再使用库，理解会深很多。

---

小结

flowchart LR
    A[历史事件 t_i] --> B[核函数 φ 叠加]
    B --> C[条件强度 λ(t)]
    C --> D{Ogata 稀疏化}
    D --> E[模拟未来事件]
    C --> F{MLE 对数似然}
    F --> G[拟合 μ, α, β]
    G --> H[分支比 n = α/β]
    H --> I{稳定性判断}

核心就四句话：

1. 强度 = 基线 + 历史激励的叠加
2. $R(t)$ 递推让模拟和 MLE 都是 $O(N)$
3. 分支比 $n = \alpha/\beta$ 判断系统是否稳定
4. Ogata 稀疏化是模拟的标准算法，MLE 是拟合的标准算法

把第 8 节的代码跑通，你就已经手撕了霍克斯过程最核心的部分。接下来按业务选核函数、上多元、加时变基线，都是在这个骨架上长肉。

---

参考资料

• Hawkes, A. G. (1971). Spectra of some self-exciting and mutually exciting point processes. Biometrika.
• Ogata, Y. (1981). On Lewis’ simulation method for point processes. IEEE TIT.
• Laub, P. J. et al. (2015). Hawkes Processes. arXiv:1507.02822 — 非常好的综述与代码索引
• Python 库: tick, hawkeslib