💡 Hawkes過程をざっくり理解する

Hawkes過程は、イベントの発生が後続のイベントの発生頻度を増加させるという自己励起性を持つ点過程モデルです。金融市場の取引活動、地震の発生、社会的活動の拡散など、自己励起的な特性を持つ現象をモデル化するのに広く用いられています。

Hawkes過程のこの自己励起性は、その強度関数によって特徴付けられています。強度関数への理解を深めるために、まず定常ポアソン過程、次に非定常ポアソン過程、そして点過程の一般理論に沿った発展の順を追って説明していきます。

強度関数とは

強度関数 $\lambda(t)$ は、ある時刻 $t$ におけるイベントの発生率を表します。

定常ポアソン過程の強度関数

定常ポアソン過程では、イベントは時間が経つにつれて独立してランダムに発生します。ここでは、全ての時間間隔でイベントの発生率（強度） $\lambda$ は一定です。時間が経っても発生率は変わらないため、「定常」と呼ばれます。

非定常ポアソン過程の強度関数

対照的に、非定常ポアソン過程ではイベントの発生率 $\lambda(t)$ が時間と共に変化します。これは、時間に応じて変動するイベントの発生率をモデル化するのに適しています。例えば、日中の異なる時間帯で電話のコール数が変わる場合などが該当します。

一般化された点過程の強度関数

さらに一般化すると、イベント発生の現象は時間だけでなく、過去のイベントに依存して発生率が変化する可能性もあります。この場合、強度関数は時間と、その時点までのイベントの履歴 $H(t)$ に基づいて決定されます。過去のイベントに依存する強度関数 $\lambda(t|H(t))$ を特に条件付き強度関数と呼びます。

Hawkes過程の強度関数

Hawkes過程は、一般化された点過程の特殊なケースとなります。Hawkes過程の条件付き強度関数は、ベースとなる発生率 $\mu$ と過去の全イベント $t_i$ に対する励起関数 $g(t - t_i)$ の総和で表されます。

$$
\lambda(t|H(t)) = \mu + \sum_{t_i \in H(t)} g(t - t_i)
$$

励起関数は、過去のイベントが未来のイベント発生にどれだけの影響を与えるかを示します。具体的には、この関数は過去のイベントからの時間経過に応じて減衰することが一般的で、これによりイベントの影響が時間と共に減少する様子を捉えます。このようにして自己励起性を持つように条件付き強度関数がデザインされています。

具体的な励起関数として、指数減衰する以下のような関数がよく用いられます。

$$
g(t) = \alpha \exp(-\beta t)
$$

この指数型の関数がよく用いられる理由については、後述します。

💡 Hawkes過程の対数尤度

Hawkes過程の対数尤度を求めるために、まずは一般化された点過程の対数尤度を確認します。

一般的な点過程の対数尤度

条件付き強度関数 $\lambda(t|H(t))$ の点過程を考えます。観察期間 $[0, T]$ 内で、イベント $\boldsymbol t = (t_1, t_2 ,..., t_n)$ が発生したときの対数尤度は、

$$
l(\theta | \boldsymbol t) = \sum_{i=1}^n \log\lambda(t_i|H(t_i)) - \Lambda(T)
$$

と表せます。¹ ここで、 $\Lambda(t)$ はcompensatorと呼ばれ、条件付き強度関数を $0$ から $t$ まで積分した値を表します。

$$
\Lambda(t) = \int_0^t \lambda(s|H(s)) ds
$$

Hawkes過程の対数尤度

一般化された点過程の対数尤度に、Hawkes過程の条件付き強度関数の形を代入して、Hawkes過程の対数尤度を求めます。

まず、観察期間 $[0, T]$ を

$$
[0, T] = [0, t_1) \cup (t_1, t_2] \cup \cdots \cup (t_{n-1}, t_n] \cup (t_n, T]
$$

と分割します。便宜上、 $t_{n+1} = T$ とします。このとき期間 $[0, t_1)$ の内の時刻 $s$ ではイベント履歴 $H(s)$ は空期間 $(t_i, t_{i+1}]$ 内の時刻 $s$ ではイベントの履歴 $H(s)$ は ${t_1, t_2, ...., t_i}$ となることに注意すると、 $M(t) := \int_0^t \mu(s)ds$ として、

\begin{align}
\Lambda(T) &= \int_0^T \lambda(s|H(s)) ds \\
&= \int_0^{t_1}\lambda(s|H(s))ds +  \sum_{i=1}^n \int_{t_i}^{t_{i+1}} \lambda(s|H(s))ds \\
&= \int_0^{t_1} \mu ds + \sum_{i=1}^n \int_{t_i}^{t_{i+1}} \left[ \mu + \sum_{j=1}^i g(s-t_j) \right] ds \\
&= \mu t_1 + \sum_{i=1}^n \left[ \mu (t_{i+1} - t_i) + \sum_{j=1}^i \left[ G(t_{i+1} - t_j) - G(t_i - t_j) \right] \right]\\
&= \mu T + \sum_{i=1}^n \sum_{j=1}^i \left[ G(t_{i+1} - t_j) - G(t_i - t_j) \right] \\
&= \mu T + \sum_{i=1}^n G(T-t_i)
\end{align}

ここで、 $G(t) := \int_0^t g(s)ds$ です。また、最後から２式目の二重和の多くの項が相殺され、最終式に帰着します。したがって、Hawkes過程の対数尤度は

$$
l(\boldsymbol \theta | \boldsymbol t) = \sum_{i=1}^n \left[ \log \left( \mu + \sum_{j=1}^{i-1} g(t_i - t_j) \right) \right] - \mu T - \sum_{i=1}^n G(T-t_i)
$$

となります。

🚨 計算量オーダーの問題

Hawkes過程の各イベント発生時刻 $t_i$ における条件付き強度

$$
\lambda(t_i | H(t_i)) = \mu + \sum_{j=1}^{i-1} g(t_i - t_j)
$$

を求めるには、 $O(n)$ の計算量が必要になります。そのため、対数尤度の第一項は $O(n^2)$ の計算量になります。イベント数 $n$ が増加すると、この方法では計算が現実的ではなくなります。

しかし、特定の励起関数を使用することで、この計算量オーダーを $O(n)$ まで削減できます。その一例が先に紹介した指数減衰する励起関数 $g(t) = \alpha \exp(-\beta t)$ です。以下のように計算します:

\begin{align*}
\lambda( t_{ i + 1 } | H( t_{ i + 1 } ) ) &= \mu + \sum_{ j = 1 }^i \alpha e^{ - \beta ( t_{ i + 1 } - t_j ) } \\
&= \mu + \left[  \sum_{ j = 1 }^i  \alpha e^{- \beta ( t_i - t_j ) } \right] e^{- \beta ( t_{ i + 1 } - t_i )} \\
&= \mu + \left[  \alpha + \sum_{ j = 1 }^{ i - 1 } \alpha e^{- \beta ( t_i - t_j ) } \right] e^{- \beta ( t_{ i + 1 } - t_i )} \\
&= \mu + \left[  \alpha - \mu + \mu + \sum_{ j = 1 }^{ i - 1 } \alpha e^{- \beta ( t_i - t_j ) } \right] e^{- \beta ( t_{ i + 1 } - t_i )} \\
&= \mu + \left[  \alpha - \mu + \lambda( t_i | H( t_i ) ) \right] e^{- \beta ( t_{ i + 1 } - t_i )}
\end{align*}

この式により、各イベントの発生時刻における条件付き強度は再帰的に計算可能です。これによって、対数尤度を $O(n)$ のオーダーで算出することが可能になります。

このような計算効率の向上は、励起関数の選択において重要な基準です。本例のように、直前のイベント発生時刻の条件付き強度を元に、次のイベント発生時刻の条件付き強度を求めることができる、ある種のマルコフ性を持つ関数が、励起関数として適しています。

💡 互いに励起するHawkes過程: 多変量化

金融市場では、特定の株式や商品の取引が他の関連市場や商品に影響を与えることがあります。例えば、ある企業の株価の急落が関連する業界や市場セクターの株価に影響を与える場合、これらの相互作用は多変量Hawkes過程を用いてモデル化できます。

多変量Hawkes過程の対数尤度を考えるために、まずは一般化された点過程に再び焦点を当てて、多変量バージョンに拡張していきます。

一般化された多変量点過程の対数尤度

$K$ 個の一般化された点過程を考えます。これらの点過程の条件付き強度関数を $\lambda_1(t), ..., \lambda_K(t)$ 、またcompensatorを $\Lambda_1(t), ..., \Lambda_K(t)$ とします。期間 $[0, T]$ に発生したすべてのイベントを $(t_1, k_1), ...., (t_n, k_n)$ と表すと、

$$
l(\boldsymbol \theta | \boldsymbol t) = \sum_{i=1}^n \log(\lambda_{k_i}(t_i | H(t_i)) - \sum_{k=1}^K \Lambda_k(T)
$$

となります。²

多変量Hawkes過程の対数尤度

述の通り、多変量Hawkes過程では、自己励起性（自己の履歴が将来のイベント発生に影響を与える）と相互励起性（他のイベントタイプがイベント発生に影響を与える）の両方の特徴を表現する必要があります。そこで、下記のような形の条件付き強度関数を想定します:

$$
\lambda_k(t | H(t)) = \mu_k + \sum_{(t_i, k_i) \in H(t)} g_{k_i, k}(t - t_i)
$$

ここで、 $\lambda_k(t|H(t))$ は時刻 $t$ におけるイベントタイプ $k$ の条件付き強度関数 $g_{j,k}(t-s)$ はタイプ $j$ の過去のイベントがタイプ $k$ に与える影響を表しています。

ここで、上記の条件付き強度関数を一般化された多変量点過程の対数尤度に適用すると、１変量の場合と同様に計算することで以下のようになります:

$$
l(\boldsymbol \theta | \boldsymbol t) = \sum_{i=1}^n \log \left[ \mu_{k_i} + \sum_{j=1}^{i-1} g_{k_j, k_i}(t_i - t_j) \right] - \sum_{k=1}^K \left[ \mu_k T + \sum_{j=1}^n G_{k_j, k}(T - t_i) \right]
$$

ここで、 $G_{j,k}(t) := \displaystyle \int_0^t g_{j,k}(s) ds$ です。

指数減衰する励起関数の場合

１変量の場合と同様に、指数減衰する励起関数を考慮すると、対数尤度の計算量オーダーを削減することができます。そこで次のような指数型の励起関数を考えます。

$$
\lambda_k(t) = \mu_k + \sum_{(t_i, k_i) \in H(t)} \alpha_{k_i, k} e^{-\beta_k (t-t_i)}, ~~k=1,...K
$$

ここでの目標は各 $\lambda_{d_i}(t_i), ~~i=1,...,n$ の計算量を減らすことです。そこで、少し余分な値を求めることになりますが、道筋をスッキリさせるために、各 $\lambda_k(t_i), ~~k=1,...,K, ~i=1,...,n$ を求めることにします。

まず以下の記法を用意します:

$\boldsymbol \lambda(t|H(t)) := (\lambda_1(t), ..., \lambda_K(t))^T$
$\boldsymbol \mu := (\mu_1, ..., \mu_K)^T$
$\boldsymbol \alpha_i := (\alpha_{i, 1}, ..., \alpha_{i, m})^T$
$\boldsymbol \beta := (\beta_1, ..., \beta_K)$

$\boldsymbol \alpha_i$ はタイプ $i$ が（自身を含む）各タイプに与える影響力をまとめたベクトルです。この記法を用いると、条件付き強度関数は次のように表せます:

$$
\boldsymbol \lambda(t|H(t)) = \boldsymbol \mu + \sum_{(t_i, k_i) \in H(t)} \boldsymbol \alpha_{d_i} e^{- \boldsymbol \beta (t-t_i) }
$$

この形式により、条件付き強度を以下のように再帰的に計算できます:

\begin{align*}
\boldsymbol{\lambda}(t_{i+1} | H(t_{i+1})) &= \boldsymbol{\mu} + \sum_{j=1}^i \boldsymbol{\alpha}_{d_j} e^{- \boldsymbol{\beta}(t_{i+1} - t_j)} \\
&= \boldsymbol{\mu} + \left[ \sum_{j=1}^i \boldsymbol{\alpha}_{d_j} e^{- \boldsymbol{\beta}(t_i - t_j)} \right] e^{- \boldsymbol{\beta}(t_{i+1} - t_i)} \\
&= \boldsymbol{\mu} + \left[ \boldsymbol{\alpha}_i e^{-\boldsymbol{\beta}(t_i - t_i)} + \sum_{j=1}^{i-1} \boldsymbol{\alpha}_{d_j} e^{- \boldsymbol{\beta}(t_i - t_j)} \right] e^{- \boldsymbol{\beta}(t_{i+1} - t_i)} \\
&= \boldsymbol{\mu} + \left[ \boldsymbol{\alpha}_i - \boldsymbol{\mu} + \boldsymbol{\mu} + \sum_{j=1}^{i-1} \boldsymbol{\alpha}_{d_j} e^{- \boldsymbol{\beta}(t_i - t_j)} \right] e^{- \boldsymbol{\beta}(t_{i+1} - t_i)} \\
&= \boldsymbol{\mu} + \left[ \boldsymbol{\alpha}_i - \boldsymbol{\mu} + \boldsymbol{\lambda}(t_i | H(t_i)) \right] e^{- \boldsymbol{\beta}(t_{i+1} - t_i)} \\
\end{align*}

また、Compensator $\boldsymbol{\Lambda}(t) := (\Lambda_1(t), ..., \Lambda_K(t))^T$ についても次のように計算できます:

$$
\boldsymbol{\Lambda}(t) = \boldsymbol{\mu} t + \sum_{(t_i, k_i) \in H(t)} \frac{1}{\boldsymbol{\beta}}(1-e^{-\boldsymbol{\beta}(t - t_i)})
$$

この方法により、１変量の場合と同様に、各イベント発生時刻での各タイプの条件付き強度を効率的に計算し、対数尤度の計算量を大幅に削減できます。

🚨 減衰パラメータ推定の不確実性

これまでの議論から、対数尤度を現実的な計算時間で求めることが可能であり、最尤推定を実行することができそうです。しかし、もう一つの問題が存在します。それは、負の対数尤度が、パラメータ $\boldsymbol{\beta}$ に対して凸性を持たない点です。この非凸性の問題は、最適化の過程で複数の局所的最適解が存在する可能性を示唆しています。つまり、初期値に依存して異なる解に収束する可能性があり、グローバルな最適解を見つけるのが困難になります。

そのため、最尤推定を行う際には、パラメータ $\boldsymbol{\beta}$ の推定を避け、経験的に得られた値を固定し、代わりにパラメータ $\boldsymbol{\mu}, A = (\alpha_{ij})$ のみを推定することが一般的です。パラメータ $\boldsymbol{\beta}$ を推定するためには、最尤推定ではなく、ベイジアンアプローチなどの別の方法を採用する必要があります。³

💡 Pythonで最尤推定

🙇‍♂ WIP: 記事執筆中です。もうしばらくお待ちください。

詳細な導出が知りたい方は近江崇宏・野村俊一『点過程の時系列解析』2021，共立出版を参照してください。とてもわかりやすい良書です！最高です！ ↩
1.同様 ↩
減衰パラメータの推定の不確実性については、以下の論文を参照してください:
Santos, T., Lemmerich, F., & Helic, D. (2021). Surfacing Estimation Uncertainty in the Decay Parameters of Hawkes Processes with Exponential Kernels. arXiv preprint arXiv:2104.01029. ↩