統計数理（2021）
第 69 巻 第 2 号 165–180©2021 統計数理研究所

特集「Hawkes 過程の新展開と応用」［原著論文］

Hawkes 型計数時系列モデル
小山 慎介†（受付 2020 年 12 月 25 日；改訂 2021 年 4 月 17 日；採択 5 月 6 日）

要旨

本論文では，離散時間の計数時系列に対して Hawkes 過程に類似する時系列モデルを提案する．特に Hawkes 過程と分岐過程の対応に着目し，同様な性質を計数時系列モデルで構築する．
このときに，計数時系列モデルの構成を Hawkes 過程とパラレルに展開し，両者を対比する．
最後に，提案するモデルを感染症のモデリングに応用する．Wallinga と Teunis による実効再生産数の推定量を Hawkes 型計数時系列モデルから導出し，過分散の場合に一般化する．
キーワード：Hawkes 過程，分岐過程，計数時系列モデル，負の二項分布，実効再生
産数．

1.

はじめに

イベント（事象）の発生が引き金となって更なるイベントが引き起こされる，そのような現象は
至るところで観察される．地震や感染症は典型的な例である．Alan G. Hawkes は，不規則に発生するイベントの続発性や相互誘発作用を記述する点過程を導入した（Hawkes, 1971a, 1971b）．
Hawkes 過程はベースラインの発生率に過去のイベントの影響を足し合わせた強度を持つ点過程である．線形性のため，Hawkes 過程は分岐過程からなるクラスターの重ね合わせで表現．イベントの連鎖的な発生やクラスター性は様々することができる（Hawkes and Oakes, 1974）
な現象で共通にみられる特徴であり，地震（Ogata, 1988），神経活動（Chornoboy et al., 1988;
Pernice et al., 2011），感染症（Meyer et al., 2012; Chiang et al., 2020; Koyama et al., 2021），
ファイナンス（Bacry et al., 2015; Hawkes, 2018），ソーシャル・ネットワーキング・サービス（Fox et al., 2016; Kobayashi and Lambiotte, 2016; Koyama and Shinomoto, 2020），犯罪や紛争（Mohler et al., 2011; Lewis et al., 2012; Zhuang and Mateu, 2019），交通事故（Kalair et al.,2021）など多くの分野のモデリングに Hawkes 過程は応用されている．
Hawkes 過程の普及が急速に進んでいる背景のひとつに，計測データの精密化を挙げることができる．個々のイベントが発生する場所や時刻が正確に記録され，発生のタイミングや相関に興味があるとき，連続時間上で定義された Hawkes 過程は適当なモデルである．
一方で，個々のイベント発生時刻は計測されず，イベント発生件数が区間毎に集計された
データも多く存在する．例えば，本稿で取り上げる新型コロナウィルス感染症（COVID-19）では，1 日の新規陽性者数が毎日報告されている．感染もしくは発病した時点をイベントとす
れば，感染症の伝播を Hawkes 過程でモデル化することはできる．しかし，実際には 1 日の発生件数だけがデータとして記録されるので，個々のイベントの発生時刻は観測されない潜在変数とみなされる．このようなデータに Hawkes 過程を当てはめるアプローチも考えられる†

統計数理研究所：〒 190–8562 東京都立川市緑町 10–3

統計数理

166

第 69 巻 第 2 号 2021

が（Cheysson and Lang, 2020），ここでは Hawkes 過程をデータ生成モデルとして想定するのではなく，離散時間の計数時系列 “そのもの” をモデリングするアプローチをとる．
本論文では，Hawkes 過程に類似する離散時間の計数時系列モデルを提案する．計数時系列
に対するよく知られたモデリング法として，動的一般化線形モデル（West and Harrison, 1997;（あるいは一般状態空間モデル，Kitagawa, 2010）や整数値自己回帰モ
Fahrmeir and Tutz, 2001）
が挙げられるが，これらとは異なるアプローチを展
デル（Kirchner, 2016, 2017; 中嶋 他, 2017）
開する．ここでは特に Hawkes 過程と分岐過程の対応に着目し，計数時系列モデルで同様の性質を構成する．第 2 節では Hawkes 過程と分岐過程の対応関係をまとめ，第 3 節で Hawkes 型
計数時系列モデルを構築する．第 4 節では Hawkes 型計数時系列モデルから実効再生産数の推定量を導き，COVID-19 のデータに応用する．第 5 節で他の方法との関連を議論する．
2.

Hawkes 過程

時刻 t ∈ R までに生起したイベント数を N (t) とし，対応するイベント発生時刻を ti(i = 1, 2, . . .) とする．時刻 t までのイベントの発生履歴 Ht = {ti |ti < t} が与えられた下で次の瞬間にイベントが発生する確率が，条件付き強度関数 λ(t) を用いてP {N (t + Δt) − N (t) = 1|Ht } = λ(t)Δt + o(Δt)
P {N (t + Δt) − N (t) > 1|Ht } = o(Δt)

で与えられるとする．Hawkes 過程は条件付き強度関数が

 t

(2.1)

g(t − u)dN (u)λ(t) = μ +0で与えられる点過程である（Hawkes, 1971b）．式（2.1）の右辺第 1 項の μ はベースラインの発生率を表し，第 2 項は過去のイベントの影響を表す．過去のイベントの影響の時間変化を表す
カーネル関数 g(τ ) は g(τ ) ≥ 0 および g(τ ) = 0 (τ ≤ 0) を満たすとする．期間 (0, T ] に n 個のイ（2.1）を用いてベントが時刻 {t1 , . . . , tn } に発生する確率密度関数は，条件付き強度関数(2.2)

p(0,T ] (t1 , . . . , tn ) =

 n






=

i=1


μ+


j<i


λ(t)dt

0

i=1
n

  T

λ(ti ) exp −g(ti − tj ) exp −μT −

n  T

j=1


g(t − tj )dttjで与えられる（Daley and Vere-Jones, 2003, Chap.7）．
は過去のすべてのイベントの影響を受けるので，各々のイベントが過条件付き強度関数（2.1）
去のどのイベントによって引き起こされたのかは定まらない．ここで，各イベントに対してそれを引き起こした “親イベント” を割り当てることで，Hawkes 過程に対応する分岐過程を構成することができる（Hawkes and Oakes, 1974; 近江・野村, 2019, 第 5 章）．i 番目のイベントを引き起こした親イベントの番号を zi ∈ {0, 1, . . . , i − 1} とする．zi = 0 の場合は親イベントを持たないとする．イベント発生時刻と親イベントの番号を合わせた時系列 {(ti , zi )|i = 1, . . . , n}は以下のルールに従って生成されるとする．
（a）親を持たないイベント {(ti , zi )|zi = 0} は強度 μ の Poisson 過程に従って発生する．
（b）j (≥ 1) 番目のイベントを親に持つイベント {(ti , zi )|zi = j} は強度 g(ti − tj ) の Poisson
過程にしたがって発生する．

Hawkes 型計数時系列モデル

167

図 1．上：分岐過程の系統図．黒丸と白丸はそれぞれ親を持たないイベント（a）と子イベント（b）を表し，矢印は親子関係を表す．下：イベント発生時刻だけを見ると Hawkes 過程に従う．

このようにして発生したイベントの親子関係を繋げていくと，親を持たないイベントを祖先とする分岐過程の系統図を描くことができる（図 1 上）．
と（b）の重ね合わせで与えられるとすると，期間 (0, T ] におけるイベント時系列全体は（a）{(ti , zi )|i = 1, . . . , n} の確率密度関数は(2.3)

p(0,T ] ({ti , zi |i = 1, . . . , n}) = p({ti , zi |zi = 0})np({ti , zi |zi = j}|tj )j=1

で与えられる．ここでp({ti , zi |zi = 0}) =

(2.4)





μ exp(−μT )i:zi =0

および
(2.5)





p({ti , zi |zi = j}|tj ) =



 
T

g(ti − tj ) exp −g(t − tj )dt

tj

i:zi =j

はそれぞれ
（a）と
（b）
のイベントの寄与を表す．式（2.4）と式（2.5）を式（2.3）に代入して整理すると確率密度関数は(2.6)

p(0,T ] ({ti , zi |i = 1, . . . , n}) =

 n






ψizi exp −μT −

i=1

n  T

j=1



g(t − tj )dt

tjと求められる．ここでμ,g(ti − tj ),

ψij =

j=0
1≤j<i

とおいた．
式（2.6）はイベント発生時刻と親イベントの情報を含む分岐過程のモデルである．これを親イベント {z1 , . . . , zn } について周辺化すると(2.7)

p(0,T ] (t1 , . . . , tn ) =01
z1 =0 z2 =0=

 n i−1
i=1 zi =0

···

n−1

zn =0



p(0,T ] ({ti , zi |i = 1, . . . , n})ψizi exp −μT −

n  T

j=1

tj


g(t − tj )dt

統計数理

168

=

n



μ+

i=1

第 69 巻 第 2 号 2021







g(ti − tj ) exp −μT −

j<i

n  T

j=1


g(t − tj )dt

tjとなり Hawkes 過程の確率密度関数（2.2）
が導かれる．すなわちイベント間の因果関係は観測されず，イベント発生時刻だけ観測されたものが Hawkes 過程であると解釈することができる（図 1 下）．
（2.6）
か ら ，イ ベ ン ト 時 刻 列 {t1 , . . . , tn } が 与 え ら れ た 下 で の 親 イ ベ ン ト
式（2.2）と 式{z1 , . . . , zn } の条件付き確率は
p(0,T ] (z1 , . . . , zn |t1 , . . . , tn ) ==p(0,T ] ({ti , zi |i = 1, . . . , n})p(0,T ] (t1 , . . . , tn )
n

ψiz

i

i=1

λ(ti )=np(zi |t1 , . . . , ti )i=1と求められる．ここでp(zi |t1 , . . . , ti ) =

(2.8)

ψizi
λ(ti )

は {t1 , . . . , ti } が与えられた下での zi の条件付き確率である．
親イベントの条件付き確率は様々な用途に応用されている．Zhuang et al.（2002）はこの確率に基づいて親イベントをサンプリングするアルゴリズムを提案した．また，親イベントを潜在変数と見なして，パラメータ推定のための EM アルゴリズムを構成することもできる（Veen．
and Schoenberg, 2008）
3.

Hawkes 型計数時系列モデル

ここからは，イベント発生件数が区間毎に集計された離散時間の計数時系列を考える．このようなデータに対して，Hawkes 過程と同様の性質を持つ計数時系列モデルを構成する．
3.1 モデルの定義
時点 i ∈ {1, 2, . . .} におけるイベント数を ni とする．期待値 λ = E(n) をパラメータに持つイベント数の確率分布 p(n, λ) を用いて，過去のイベント数 {n1 , . . . , ni−1 } が与えられた下での niの条件付き確率がP (ni |n1 , . . . , ni−1 ) = p(ni , λi )で与えられる計数時系列モデルを考える．Hawkes 過程の条件付き強度関数（2.1）からの類推で期待値パラメータを(3.1)

λi = μ +i−1
gi−j njj=1で与える．右辺第 1 項の μ はベースラインの頻度であり，第 2 項は過去のイベントの影響を表す．gτ は過去のイベントの影響の時間変化を表し，gτ ≥ 0 および gτ = 0 (τ ≤ 0) を満たすとする．L 時点までの時系列 {n1 , . . . , nL } の結合確率分布は(3.2)

P (n1 , . . . , nL ) = P (n1 )

L

i=2

P (ni |n1 , . . . , ni−1 ) =

L

i=1

p(ni , λi )

Hawkes 型計数時系列モデル

169

で与えられる．
イベント数の確率分布 p(n, λ) には Poisson 分布もしくは負の二項分布を用いることにする．
Poisson 分布の確率質量関数は(3.3)

p(n, λ) =

λn −λ
e
n!

で与えられる．平均と分散は等しく E(n) = Var(n) = λ で与えられる．集計区間内の各々のイベント発生が互いに独立であるときイベント数は Poisson 分布に従うので，このような状況を近似的にも想定できる場合は Poisson 分布は妥当な選択である．
一方，集計区間内でイベントの自励効果を無視できないとき，イベント数の分散は平均よりも大きくなる．これを過分散（over-dispersion）という．このような場合には過分散を持つ確率分布を用いることが望ましい．負の二項分布はそのような分布のひとつである．ここでは以下の確率質量関数を持つ負の二項分布を考える（Koyama and Fujiwara, 2019）：(3.4)

p(n, λ, ρ) =Γ(n + λρ )Γ(n + 1)Γ( λρ )

ρ
1+ρ

n

1
1+ρ

λ
ρ

平均と分散はそれぞれ E(n) = λ および Var(n) = (1 + ρ)λ で与えられる．ρ(> 0) は過分散の度（3.4）は Poisson 分布（3.3）に収束する．負の二項分布合いを表すパラメータであり，ρ → 0 で式については補足 A にまとめた．
（3.3）
もしくは負の二項分布（3.4）で与えられる時系列
イベント数の確率分布が Poisson 分布（3.2）–
を “Hawkes 型” 計数時系列モデルと呼ぶことにする．他の確率分布でも時系モデル（3.1）
列モデルは定義されるが，特にこれら二つの分布を用いる理由は第 3.3 節で明らかになる．
3.2 親イベントの割り当て
Hawkes 過程では，各イベントに親イベントをひとつ割り当てることで分岐過程と対応付けることができた．同様に Hawkes 型計数時系列モデルに対しても親イベントを割り当ててみよう．時点 i のイベント数 ni のうち，過去の時点 j(< i) のイベントに引き起こされたイベント
数を yij とし，親を持たないイベント数を yii とする．つまり {yi1 , . . . , yii } は ni の発生要因についての内訳であり
(3.5)

ni =

i


yij

j=1

を満たす．イベント数の内訳 yij の確率分布について以下を仮定する．
（a’）親を持たないイベント数 yii は期待値が E(yii ) = μ の確率分布 P (yii ) に従う．
（b’）過去の時点 j(< i) のイベントに引き起こされたイベント数 yij は期待値が E(yij |nj ) =gi−j nj の確率分布 P (yij |nj ) に従う．
時点 i に発生するイベント数の内訳 {yi1 , . . . , yii } は，過去のイベント数 {n1 , . . . , ni−1 } が与えられた下で互いに独立であるとする：(3.6)

P (yi1 , . . . , yii |n1 , . . . , ni−1 ) = P (yii )i−1P (yij |nj )j=1

すると，時点 1 から L までのイベント数の全内訳 Y1:L := {yij |i = 1, . . . , L, j = 1, . . . , i} の結合確率分布は

170

統計数理

(3.7)

P (Y1:L ) = P (y11 )第 69 巻 第 2 号 2021
L

P (yi1 , . . . , yii |n1 , . . . , ni−1 )

i=2

=

L


P (yii )

i=1

i−1


P (yij |nj )j=1で与えられる．式（3.7）
は各時点のイベント数の内訳情報を含むモデルであり，Hawkes 過程にに対応している．
対する分岐過程（2.6）（3.5）
を満たすように周辺化すると計数時系列 {n1 , . . . , nL } の確率式（3.7）を Y1:L について式（3.5）を満たす Yi の値の集合を Yi とす
分布が得られる．すなわち Yi = {yi1 , . . . , yii } とし，式
ると
(3.8)

P (n1 , . . . , nL ) =

Y1 ∈Y1=

P (Y1:L )YL ∈YLL ···

Y1 ∈Y1=···YL ∈YL i=1

L 

P (yii )i=1 Yi ∈Yi= P (n1 )P (yii )i−1P (yij |nj )

j=1

i−1


P (yij |nj )

j=1

L


P (ni |n1 , . . . , ni−1 )i=2となる．ここで(3.9)

P (ni |n1 , . . . , ni−1 ) =P (yii )Yi ∈Yii−1P (yij |nj )j=1の期待値は，仮定
（a’）
と（b’）
より
E(ni |n1 , . . . , ni−1 ) = E(yii ) +i−1E(yij |nj )

j=1

=μ+

i−1


gi−j njj=1となり式（3.1）に一致する．すなわち式（3.7）
をイベント数の内訳について周辺化することでHawkes 型計数時系列モデルが得られた．これは，分岐過程（2.6）を親イベントについて周辺化することで Hawkes 過程が得られることに対応している．
3.3 加法性とイベント内訳の条件付き確率分布
Hawkes 過程に対する親イベントの条件付き確率（2.8）と同様に，各時点のイベント数 ni が与えられた下での内訳 yij の条件付き確率を導こう．そのために，yij に対して以下の “加法性”の条件を加える．

定義.

n 個の独立な確率変数の和y = y1 + · · · + yn ,yi ∼ p(yi , λi )

Hawkes 型計数時系列モデル

171が同一の確率分布族 p(y, λ1 + · · · + λn ) に従うとき加法的であるという．
イベント数の内訳 yij が加法的な確率分布 p(yij , ψij ) に従うとする．ここでμ,gi−j nj ,ψij = E(yij ) =

j=i
1≤j<i

とおいた．すると，式（3.5）
および式（3.9）より，ni は互いに独立な yij (j = 1, . . . , i) の和で与えられるから，加法性より同一の確率分布族
P (ni |n1 , . . . , ni−1 ) = p(ni , λi )(3.10)に従う．ここで λi は式（3.1）
で与えられる．つまり加法性を仮定すると，各時点のイベント数
ni もその内訳 yij (j = 1, . . . , i) も同一の確率分布族に従うのである．
（3.10）
より，時点 i までのイベント数 {n1 , . . . , ni } が与えられた下での時点 i の式（3.6）と式
イベント数の内訳 Yi = {yi1 , . . . , yii } の条件付き確率分布はp(Yi |n1 , . . . , ni ) =p(Yi |n1 , . . . , ni−1 )p(ni |n1 , . . . , ni−1 )ip(yij , ψij )j=1=(3.11)

p(ni , λi )と求められる．ここまで来ると，イベント数の分布に Poisson 分布（3.3）と負の二項分布（3.4）
を採用した理由が明らかになる．すなわち，これらはともに加法的である．それぞれの確率質に代入することで，Yi の条件付き確率を具体的に導くことができる．
量関数を式（3.11）
Poisson 分布の場合 式（3.3）
を式（3.11）に代入すると，Yi の条件付き確率分布は多項分布
p(Yi |n1 , . . . , ni ) =(3.12)ni !

i
ψij yijiy !
j=1 ij j=1λiとして求められる．平均と分散はそれぞれni ψijλi
ni ψijVar(yij |n1 , . . . , ni ) =λiE(yij |n1 , . . . , ni ) =(3.13)(3.14)

で与えられる．
負の二項分布の場合(3.15)

1−

ψij
λi



を式
（3.11）に代入すると，Yi の条件付き分布確率分布は式（3.4）

p(Yi |n1 , . . . , ni , ρ) =i
Γ(ni + 1)Γ( λρi ) Γ(ni + λρi )

ψij
)
ρ
ψij
Γ(yij + 1)Γ( ρ )j=1

Γ(yij +と導かれる．これはディリクレ多項分布と呼ばれ，平均と分散はそれぞれ以下で与えられる：(3.16)(3.17)

ni ψijλi
ni ψijVar(yij |n1 , . . . , ni ) = κiλi
E(yij |n1 , . . . , ni ) =ここで分散に掛かる係数 κi は
(3.18)

κi =

λi + ρni(≥ 1)
λi + ρ

1−

ψij
λi



統計数理

172

第 69 巻 第 2 号 2021

で与えられる．多項分布の分散（3.14）よりも大きいことは，元になる負の二項分布が過分散であることに対応している．
4.

感染症モデルへの応用

前節で構成した Hawkes 型計数時系列モデルを感染症の実効再生産数 (eﬀective reproduction
number) の推定に応用する．再生産数は 1 人の感染者が引き起こす 2 次感染者数の平均である．
実効再生産数はすでに感染が広がっている状況における再生産数の総称で，幾つかの異なる定義がある（Fraser, 2007; Nishiura and Chowell, 2009）．ここでは Wallinga and Teunis（2004）による実効再生産数を Hawkes 型計数時系列モデルから導き，過分散の場合に一般化する．
4.1 Wallinga and Teunis（2004）の方法
感染者 k が感染者 l から感染したとし，それぞれの感染時刻を tk , tl とする．世代時間が確率密度関数 φ(τ ) に従うとすると，感染源候補の中で k が l から感染した確率はp(k,l) = φ(tk − tl )φ(tk − tm )m=kで与えられる．全ての感染者が独立に 2 次感染者を感染させていたと仮定すると，感染者 l の2 次感染者数は
Rl ∼

(4.1)



Bernoulli(p(k,l) )k



に従う．ここで k は時刻 tl 以降の感染者すべてについての和である．t 日目の実効再生算数

Rt を，t 日目に感染した患者の 2 次感染者数の平均 Rt = n1t
Rl で定義する．ここで nt
tl =t
は t 日目の感染者数である．式（4.1）
から t 日目の実効再生算数の平均と分散はそれぞれ
1 R̂t = E(Rt ) =(4.2)
p(k,l)
nt
tl =t

k



(4.3)


1  
s = Var(Rt ) = 2
p(k,l) (1 − p(k,l) ) −
nt
2

k

tl =t




p(k,l) p(k,m)

tl =t tm =t,m=lで与えられる（Cowling et al., 2008）．これが Wallinga and Teunis（2004）によって与えられた
実効再生産数の推定量である．
4.2 Hawkes 型計数時系列モデルからの導出
感染症の伝播を Hawkes 型計数時系列モデルで表す．ni を i 日目の新規感染者数とし，頻
度が
(4.4)

λi = rii−1
j=1

φi−j nj∞

で与えられるとする．ここで φτ は世代時間の分布であり τ =0 φτ = 1 を満たす．ri は i 日目の感染力を表すパラメータとする．このとき，イベント数の内訳 yij は “i 日目の感染者のうち
j 日目の感染者から感染した人数” となる．Wallinga and Teunis（2004）に従うと，j 日目の実
効再生産数はこの日の感染者 1 人当たりの 2 次感染者数であるから，yij を用いて

Hawkes 型計数時系列モデル

Rj =

173

1 
yij
nj
i>j

で与えられる．
L 日目までの新規感染者数の時系列データ {n1 , . . . , nL } が与えられているとする．感染者
数が Poisson 分布に従うとすると，データが与えられた下での yij の条件付き確率は多項分布（3.13）と式（3.14）を用いて（3.12）に従うので，Rj の条件付き期待値と分散は式
R̂j = E(Rj |n1 , . . . , nL )

(4.5)

=

L
1 
E(yij |n1 , . . . , ni )
nj
i=j+1


L

=

ni φi−j

i−1

k=1

i=j+1

φi−k nkおよび(4.6)

s2 = Var(Rj |n1 , . . . , nL )=L
1 
Var(yij |n1 , . . . , ni )
n2j
i=j+1

L
1 
ni φi−j=i−1nj
φi−k nki=j+1


φi−j nj

1 − i−1k=1

k=1

φi−k nkと求められ，それぞれ式（4.2）と式（4.3）に一致する（証明は付録 B を参照）
．すなわち，Poisson
分布を用いた Hawkes 型計数時系列モデルから Wallinga and Teunis（2004）による実効再生産数の推定量と分散が導かれた．
に従うとすると，実効再生産数の条件付き期待値は式
一方，感染者数が負の二項分布（3.4）（3.17）より（4.5）と同じであるが，条件付き分散は式(4.7)

Tni φi−j1 s =κi i−1nj
φi−k nk
2

i=j+1





φi−j nj

1 − i−1k=1

k=1

φi−k nkとなる．過分散のパラメータ κi が掛かることに着目しよう．κi = 1 のときは式（4.6）に帰着さは Wallinga and Teunis（2004）の方法を過分散の場合に一般化した
れる．したがって，式（4.7）
ものとみなせる．
4.3 COVID-19 への応用
（4.7）を新型コロナウィルス感染症（COVID-19）の–上で導いた実効再生産数の推定量（4.5）
データに応用する．世界各国の 1 日の新規陽性者数は毎日報告され，オープンデータとして公開されている．ここでは，厚生労働省が公開するデータ（https://www.mhlw.go.jp/stf/covidを用いた．
19/open-data.html）
図 2 上に日本の新規陽性者数を示す．陽性者数にみられる 1 週間の周期変動は，平日と週末の検査態勢の違いに起因する．この周期変動の実効再生産数への影響を抑えるため，曜日毎の新規陽性者数の平均を取り，
1
7

Saturday

Sundayβi = 1

統計数理

174

第 69 巻 第 2 号 2021

図 2．上：日本の新規陽性者数．下：推定した実効再生産数（平均値 ±SD）．薄い灰色はデー
タの過分散を考慮に入れた場合を表し，濃い灰色は考慮に入れない場合を表す．

と規格化して求めた重み βi を用いて，もとの新規陽性者数 ni (i = 1, . . . , L) を ñi = ni /βi(i = 1, . . . , L) と変換した（小数部は最も近い整数に丸め込んだ）．
分布 φτ には，Nishiura et al.（2020）の報告に従い，平均 4.7 日，標発症間隔（serial interval）
準偏差 2.9 日の対数正規分布を採用し，累積分布関数
(4.8)

Φτ =

1
log τ − μerfc − √22σ 2,を用いて φτ = Φτ − Φτ −1 とした．対数正規分布のパラメータ (μ, σ 2 ) は，平均 m = 4.7，標準√偏差 s = 2.9 として，μ = log(m2 / s2 + m2 ) および σ 2 = log(1 + s2 /m2 ) で与えた．
のパラメータ ρ の値はデータから
負の二項分布（3.4）
1  (ñi − λ̄i )2
−1
L
λ̄i
L

ρ̂ =

3

i=1

と推定した．ここで λ̄i =ñ /7 は 1 週間の新規陽性者数の平均である．これと式
j=−3 i+jのパラメータ κi は
（3.18）より式（4.7）κ̂i =λ̄i + ρ̂ñi
λ̄i + ρ̂と求められる．
（平均値 ±SD）を示す．薄い灰色はデータの過分散を考慮に
図 2 下に推定した実効再生産数を表入れた場合を表し，濃い灰色は考慮に入れない場合（Wallinga and Teunis, 2004 の方法）の方法では推定量の分散がす．実際のデータは過分散であるため，Wallinga and Teunis（2004）
過少に評価されている．

Hawkes 型計数時系列モデル

175

実効再生産数は 4 月に下がり始め 1 を下回る．緊急事態宣言が発出されたのは 4 月 7 日である．しばらく実効再生産数は 1 を下回るが，5 月 25 日に緊急事態宣言が全国で解除され，夏にかけて再び 1 を超える．8 月に実効再生産数は再び 1 を下回るが，秋から冬にかけて陽性者数は急増している．
5.

まとめ

本稿では，Hawkes 過程の性質，特に分岐過程との対応に着目して，同様な性質を持つ計数時系列モデルを提案した．第 3 節の Hawkes 型計数時系列モデルの構成を，第 2 節の Hawkes 過
程とパラレルに展開させ，両者の対応をみた．ただし計数時系列モデルに対して，確率分布の加法性という条件を課した．これは Hawkes 過程にはない要請である．第 4 節では，Wallinga
and Teunis（2004）による実効再生産数の推定量を Hawkes 型計数時系列モデルをから導出した．さらに感染者数に負の二項分布を仮定することで，実効再生産数の推定量が過分散の場合に一般化された．
Hawkes 過程と計数時系列を関連付けるアプローチとして，Kirchner（2016, 2017）は整数値自己回帰モデル（Integer-valued Autoregressive Models，INAR モデル）に基づく方法を提案している．この方法では，離散観測された Hawkes 過程を INAR(p) モデルで近似し，条件付き最小 2乗法で求めた推定量を Hawkes 過程のカーネル関数に対応させている．比較的計算が容易な行列演算でカーネル関数のノンパラメトリックな推定量を与えていることが利点である．パラメ
トリックな関数を事前に想定できないとき，Kirchner の方法でノンパラメトリックに推定し，それをここで提案した時系列モデルのカーネル関数に使うこともできるであろう．
本論文では，もう一つのアプローチである一般状態空間モデル（あるいは動的一般化線形モ（3.1）にシデル）
を扱わなかったが，Hawkes 型計数時系列モデルを観測モデルとし，パラメータ
ステムモデルを組み込むことで，容易に一般状態空間モデルとして定式化することもできる．
Koyama et al.（2021）は，このアプローチで非ガウス平滑化アルゴリズムを用いた実効再生産数の推定方法を提案している．
Hawkes 過程の分岐過程による表現は，デクラスタリング（Zhuang et al., 2002），シミュレーション（Moller and Rasmussen, 2005），EM アルゴリズム（Veen and Schoenberg, 2008）など
様々なアルゴリズムや応用の基礎になる．Hawkes 型計数時系列モデルに対する同様の展開は今後の研究課題である．
付
A.

録

負の二項分布負の二項分布は通常，成功率 p のベルヌーイ試行で r 回成功する前に失敗する回数 y の確率
分布
(A.1)

p(y, p, r) =y+r−1 rp (1 − p)yr−1として表される．平均と分散はそれぞれE(y) =

である．パラメータをr(1 − p),p

Var(y) =r(1 − p)
p2

統計数理

176

第 69 巻 第 2 号 2021

p=

1
,
1+ρ

r=

λ
ρ

と変換して階乗をガンマ関数 Γ(x + 1) = x! で表すと式（3.4）の形が得られる．
式（3.4）のキュムラント母関数はK(s) = E(esy ) = −λlog[1 − (es − 1)ρ]
ρ

で与えられる．ρ で展開して主要項を取り出すとK(s) =λ s[(e − 1)ρ + o(ρ)]ρとなるので，ρ → 0 で Poisson 分布のキュムラント母関数 K(s) = λ(es − 1) に収束する．すなは ρ → 0 で Poisson 分布（3.3）に収束する．
わち負の二項分布（3.4）
yi (i = 1, . . . , n) が互いに独立な負の二項分布（3.4）に従うとすると，y = y1 + · · · + yn の分布のキュムラント母関数は

nKn (s) = −

i=1

λi

ρ

log[1 − (es − 1)ρ]となるので，y は期待値パラメータ λ1 + · · · + λn を持つ同一の分布に従う．したがって負の二項分布（3.4）は加法的である．
式（A.1）のパラメータを
1
1
p=
, r=
1 + ρλρと変換すると一般化線形モデルでよく使われるもう一つのパラメータ表示(A.2)

p(y, λ, ρ) =Γ(y + ρ−1 )Γ(y + 1)Γ(ρ−1 )yρλ1 + ρλ1
1 + ρλ1ρが得られる（Hilbe, 2011）．平均と分散はそれぞれ E(y) = λ および Var(y) = λ + ρλ2 で与えられる．式（3.4）の分散と異なることに注意しよう．Cameron and Trivedi（1986）は式（3.4）と式（A.2）をそれぞれ “Negbin I” および “Negbin II” と呼んで区別している．Negbin I は加法性を持つが，Negbin II は加法的ではない．このため Hawkes 型計数時系列モデルには Negbin I を用いる．
B.

Wallinga and Teunis（2004）の方法との等価性の証明

感染者 k と l の感染した日をそれぞれ tk = s, tl = t とすると φ(tk − tl ) = φs−t である．
φ(tk − tl ) = s ns φs−t に注意すると式（4.2）はkφ(tk − tl )ns φs−t
1 

=
R̂t =
s−1
nt
φ(t
−
t
)
m
k
φ
n
m=k
u=1 s−u ut =t ks=t+1
Llとなり式（4.5）に一致する．
の右辺第 1 項と第 2 項はそれぞれ同様に式（4.3）
φ(tk − tl )
1 
1 

p
(1
−
p
)
=
(k,l)
(k,l)
n2t
n2t
φ(tk − tm )
m=k
k

tl =t

tl =t

k

φ(tk − tl )φ(tk − tm )
m=k

1− 

Hawkes 型計数時系列モデル=

L
ns φs−t1 s−1nt
φs−u nu
u=1

s=t+1

1 
n2t
k



177

φs−t
1 − s−1φnu=1 s−u up(k,l) p(k,m)

tl =t tm =t,m=lnt − 1 
=
nt



k

φ(tk − tl )1 nt
φ(tk − tj )
j=k


1
nt − 1tl =t

tm =t,m=lφ(t − tm ) kφ(tk − tj )j=k

ns φs−tφs−tnt − 1 
s−1
s−1
nt
φs−u nu
φs−u nu
L

=

s=t+1

u=1

u=1

となるので，これらを式（4.3）の右辺に戻すとs2 =L
ns φs−t1 s−1nt
φs−u nu

φs−t
1 − s−1φnu=1 s−u uns φs−t1 s−1nt
φs−u nuφs−t nt
1 − s−1φnu=1 s−u u

s=t+1

u=1

L

=

s=t+1

u=1

ns φs−tφs−tnt − 1 
s−1
s−1
nt
φ
n
φ
n
u=1 s−u u
u=1 s−u u
s=t+1
L

−

となり式（4.6）に一致する．

参

考

文

献

Bacry, E., Mastromatteo, I. and Muzy, J. F. (2015). Hawkes processes in ﬁnance, Market Microstructureand Liquidity, 1, 1550005.
Cameron, A. C. and Trivedi, P. K. (1986). Econometric models based on count data: Comparisons andapplications of some estimators and tests, Journal of Applied Econometrics, 1, 29–53.
Cheysson, F. and Lang, G. (2020). Strong mixing condition for Hawkes processes and application to
Whittle estimation from count data (preprint), arXiv:2003.04314.
Chiang, W. H., Liu, X. and Mohler, G. (2020). Hawkes process modeling of COVID-19 with mobility
leading indicators and spatial covariates (preprint), medRxiv:2020.06.06.20124149.
Chornoboy, E. S., Schramm, L. P. and Karr, A. F. (1988). Maximum likelihood identiﬁcation of neuralpoint process systems, Biological Cybernetics, 59, 265–275.
Cowling, B. J., Ho, L. M. and Leung, G. M. (2008). Eﬀectiveness of control measures during the SARS
epidemic in Beijing: A comparison of the Rt curve and the epidemic curve, Epidemiology and
Infection, 136, 562–566.
Daley, D. and Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes Vol.1 : Elementary Theory and Methods, 2nd ed., Springer-Verlag, New York.
Fahrmeir, L. and Tutz, G. (2001). Multivariate Statistical Modeling Based on Generalized Linear Models, 2nd ed., Springer-Verlag, New York.
Fox, E. W., Short, M. B., Schoenberg, F. P., Coronges, K. D. and Bertozzi, A. L. (2016). Modeling e-mail networks and inferring leadership using self-exciting point processes, Journal of theAmerican Statistical Association, 111, 564–584.
Fraser, C. (2007). Estimating individual and household reproduction numbers in an emerging epidemic,
PLoS ONE, 2, e758.
Hawkes, A. G. (1971a). Point spectra of some mutually exciting point processes, Journal of the Royal
Statistical Society. Series B (Methodological), 33, 438–443.

178

統計数理

第 69 巻 第 2 号 2021

Hawkes, A. G. (1971b). Spectra of some self-exciting and mutually exciting point processes, Biometrika,
58, 83–90.
Hawkes, A. G. (2018). Hawkes processes and their applications to ﬁnance: A review, Quantitative
Finance, 18, 193–198.
Hawkes, A. G. and Oakes, D. (1974). A cluster process representation of a self-exciting process, Journalof Applied Probability, 11, 493–503.
Hilbe, J. M. (2011). Negative Binomial Regression, 2nd ed., Cambridge University Press, Cambridge.
Kalair, K., Connaughton, C. and Loro, P. A. D. (2021). A non-parametric Hawkes process model of primary and secondary accidents on a UK smart motorway, Journal of the Royal Statistical Society
Series. C (Applied Statistics), 70, 80–97.
Kirchner, M. (2016). Hawkes and INAR(∞) processes, Stochastic Processes and Their Applications,
126, 2494–2525.
Kirchner, M. (2017). An estimation procedure for the Hawkes process, Quantitative Finance, 17, 571–
595.
Kitagawa, G. (2010). Introduction to Time Series Modeling, Chapman & Hall/CRC, Boca Raton.
Kobayashi, R. and Lambiotte, R. (2016). TiDeH: Time-dependent Hawkes process for predicting retweet
dynamics, ICWSM 2016, 191–200.
Koyama, S. and Fujiwara, Y. (2019). Modeling event cascades using networks of additive count sequences, Journal of Statistical Mechanics: Theory and Experiment, 2019, 023402.
Koyama, S. and Shinomoto, S. (2020). Statistical physics of discovering exogenous and endogenous
factors in a chain of events, Physical Review Research, 2, 043358.
Koyama, S., Horie, T. and Shinomoto, S. (2021). Estimating the time-varying reproduction number of
COVID-19 with a state-space method, PLOS Computational Biology, 17, e1008679.
Lewis, E., Mohler, G., Brantingham, P. J. and Bertozzi, A. L. (2012). Self-exciting point process modelsof civilian deaths in Iraq, Security Journal, 25, 244–264.
Meyer, S., Elias, J. and Hohle, M. (2012). A space-time conditional intensity model for invasive
meningococcal disease occurrence, Biometrics, 68, 607–616.
Mohler, G., Short, M. B., Brantingham, P. J., Schoenberg, F. P. and Tita, G. E. (2011). Self-exciting
point process modeling of crime, Journal of the American Statistical Association, 106, 100–108.
Moller, J. and Rasmussen, J. G. (2005). Perfect simulation of Hawkes processes, Advances in Applied
Probability, 37, 629–646.
中嶋雅彦, 酒折文武, 川崎能典 (2017). 整数値自己回帰モデルの最近の発展, 統計数理, 65, 323–339.
Nishiura, H. and Chowell, G. (2009). The eﬀective reproduction number as a prelude to statistical estimation of time-dependent epidemic trends, Mathematical and Statistical Estimation Approaches
in Epidemiology, 103–121, Springer, Dordrecht.
Nishiura, H., Lintona, N. M. and Akhmetzhanov, A. R. (2020). Serial interval of novel coronavirus(COVID-19) infections, International Journal of Infectious Diseases, 93, 284–286.
近江崇宏, 野村俊一 (2019).『点過程の時系列解析』, 共立出版, 東京.
Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes, Journal of the American Statistical Association, 83(401), 9–27.
Pernice, V., Staude, B., Cardanobile, S. and Rotter, S. (2011). How structure determines correlations
in neuronal networks, PLoS Computational Biology, 7(5), e1002059.
Veen, A. and Schoenberg, F. P. (2008). Estimation of space-time branching process models in seismology using an EM-type algorithm, Journal of the American Statistical Association, 103, 614–624.
Wallinga, J. and Teunis, P. (2004). Diﬀerent epidemic curves for severe acute respiratory syndrome
reveal similar impacts of control measures, American Journal of Epidemiology, 160, 509–516.
West, M. and Harrison, J. (1997). Bayesian Forecasting and Dynamic Models, 2nd ed., Springer-Verlag,
New York.
Zhuang, J. and Mateu, J. (2019). A semiparametric spatiotemporal Hawkes — Type point process

Hawkes 型計数時系列モデル

179

model with periodic background for crime data, Journal of the Royal Statistical Society. Series
A (Statistics in Society), 182, 919–942.
Zhuang, J., Ogata, Y. and Vere-Jones, D. (2002). Stochastic declustering of space-time earthquake
occurrences, Journal of the American Statistical Association, 97, 369–380.

180

Proceedings of the Institute of Statistical Mathematics Vol. 69, No. 2, 165–180 (2021)

Hawkes-type Count Time Series Models
Shinsuke Koyama
The Institute of Statistical Mathematics

We propose a “Hawkes-type” count time series model. By emphasizing the branching process representation of the Hawkes process on a real number line, our model has asimilar representation. The applicability of the proposed model is demonstrated using an
epidemiological example where the eﬀective reproduction number proposed by Wallinga
and Teunis (2004) is generalized for over-dispersion.

Key words: Hawkes process, branching process, count time series model, negative binomial distribution,eﬀective reproduction number.