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

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

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用
茅根 脩司1 ・白石 博2（受付 2020 年 12 月 14 日；改訂 2021 年 3 月 31 日；採択 4 月 6 日）

要旨

多変量点過程の 1 つのクラスである多変量 Hawkes 過程は，自己励起性および相互励起性という特徴をもつことが知られている．本論文では，観測データが多変量 Hawkes 過程にしたがうと想定し，Hawkes 過程の特徴を表すカーネル関数を統計的に推定するための 2 つの手法を紹介する．1 つ目は最尤法を用いたパラメトリックな手法で，強度の従属性が指数型カーネル関数にしたがうと仮定してカーネル関数のパラメータ推定を行う．2 つ目はカーネル関数を特定化しないノンパラメトリックな手法で，連続時間確率過程の離散観測に基づく確率過程の近似を通じてカーネル関数を非母数的に推定する．シミュレーションにより 2 つの手法による推定結果を比較した後，パラメトリックな手法を用いて仮想通貨市場の価格変動と出来高変動をHawkes 過程に当てはめ，カーネル関数の特徴を可視化する手法の 1 つである Hawkes グラフ表
現を通じて伝播構造を可視化する．また，ノンパラメトリックな手法を用いて COVID-19 の関
東・関西の近隣の都府県の新規感染者数を Hawkes 過程に当てはめ，Hawkes グラフ表現により伝播構造を可視化する．
キーワード：Hawkes 過程，多変量点過程，INAR 過程，Hawkes グラフ，グラフィカ
ルモデリング，ノンパラメトリック推定．
1.

はじめに

社会で観測される様々な事象の裏には影響構造が存在する．特に自己励起性と相互励起性を持ったイベントは身近に多い．自己励起とはイベント発生が同種のイベント発生を引き起こしやすくする影響であり，相互励起は発生したイベントと異なるタイプのイベントの発生
確率が上がる影響である．これらの性質を持った点過程が Hawkes 過程である．Hawkes 過程
によって伝染病の感染構造のモデル化などの応用を見据えて提案された．
は，Hawkes（1971）
Hawkes 過程は様々な分野で応用され，疫学統計の分野では Rizoiu et al.（2018）
によって伝統
的に疫学統計で利用される SIR モデルとの関連付けが行われた．Park et al.（2020）は Hawkes
過程を 2014 年に西アフリカで流行したエボラウイルスのデータに適用し，SIR モデルの拡張である SEIR モデルに勝るとも劣らない結果を報告している．
疫学統計の分野以外にも，自己励起性や相互励起性を持つ身近なイベントとして地震が挙げられる．発生後に余震を引き起こす地震は Hawkes 過程への当てはまりがよく，Ogata（1988）1

慶應義塾大学大学院 基礎理工学研究科：〒 223–8522 神奈川県横浜市港北区日吉 3-14-1

2

慶應義塾大学 理工学部：〒 223–8522 神奈川県横浜市港北区日吉 3-14-1

182

統計数理

第 69 巻 第 2 号 2021

によって地震発生メカニズムのモデル化への応用研究が行われ，その後も様々なモデルの拡張が提案されている．また，ファイナンス関連の応用もされている．特に近年は高頻度取引
データを活用して，株式取引や価格変動をイベントとした Hawkes 過程の応用が増えている．
Bowsher（2007）
はニューヨーク証券取引所の “取引の発生” と “mid-quote の変化した時刻” に一般化した 2 変量 Hawkes 過程を当てはめて解析を行っている．Yang et al.（2018）は，株式市場のデータから得られた株価の増加・減少の 2 つの点過程と Thomson Reuters News sentiment（投資意欲）の増加・減少の 2 つの点過程を合わせた 4 変量
data から定義されたセンチメント
Hawkes 過程を当てはめて影響構造を解析しており，センチメントと株価の変動の複雑な関係性についての研究を行った．また，監物・中川（2019）は倒産発生の伝播構造のモデリングに多変量 Hawkes 過程を用いている．
Hawkes 過程を特徴付ける，時間的な影響変化を表す関数をカーネル関数という．従来，Hawkes 過程の推定はカーネル関数をパラメトリックに与えてパラメータの最尤推定を行うことが多かった．イベントが発生した瞬間からの影響の減衰の様子が経験的に観測されている，
または自然に考えられる地震分野やファイナンス分野の高頻度取引データの解析にはこの方法は適切である．ファイナンス分野における Hawkes 過程の応用研究をまとめた Bacry（2015）に
よると，指数型カーネル関数を用いた研究が多く行われている．しかし，疫学統計の分野においては，影響の時間変化をパラメトリックに特定することは困難な場合がある．感染症に対して．このような問題に対し，は経験的知識が十分でないことなどが理由である（Park et al., 2020）
Kirchner（2017）はカーネル関数をパラメトリックに特定せずに推定する方法を提案している．
Kirchner の手法では，Hawkes 過程は整数値をとる離散時間時系列に AR 過程（AutoRegressiveのような構造を組み込んだ INAR 過程（INteger valued AutoRegressive process）に関process）
連付けられ，INAR 過程の推定を通して Hawkes 過程の推定が行われる．Hawkes 過程の推定
結果を可視化する方法の 1 つとして Hawkes グラフがある．Hawkes グラフは Embrechts andが提案した重み付き有向グラフで各頂点がイベントタイプに対応し，頂点と各Kirchner（2018）
有向辺の重みがそれぞれイベント本来の発生しやすさと 1 回のイベント発生が及ぼす総合的な影響度合いを表している．
本論文では，最尤法と Kirchner の手法の比較を行うために（1）最尤法が想定する関数族に最尤法が想定する関数族にカーネル関数が属していない
カーネル関数が属している場合と（2）
場合の 2 つの状況でシミュレーションを行う．その後，それぞれの方法を用いて実データ解析を行う．まず，最尤法で解析するデータは日本における大手の仮想通貨取引所である bitFlyer
における仮想通貨の 15 分足データである．株価市場には取引できない時間があるのに対して，
仮想通貨市場は 24 時間開いているなどの違いがあるが，影響の時間変化は株価市場と類似しと同様，指数型カーネル関数を想定して解析を行う．イベンていると考え，Yang et al.（2018）
トとして定義されるのは価格と出来高の上昇・低下の中で変化の大きさがある閾値を超えた時刻から構成される 4 変量計数データを解析する．次に Kirchner の手法を用いて，2020 年に日本を含めた世界中で流行している COVID-19（COronaVIrus Disease 2019：新型コロナウイル
ス）の日本国内の新規感染者についてのデータを解析する．Hawkes 過程を用いた COVID-19 の解析事例については，Garetto et al.（2021），Lesage（2020），Gong et al.（2021）など幾つかの報告が既になされているが，COVID-19 は新しいウイルスであり，感染力等についての事前知識が十分でないなどの理由により，本論文では Kirchner の手法を用いてカーネル関数をノンパ
ラメトリックに推定する．関東および関西の都府県間の影響構造の推定を通して，日本国内における COVID-19 の感染拡大の第 1 波と第 2 波での影響構造の差異を解析する．また，人の移動と感染構造の関係についても解析を行う．
本論文は 2 節で Hawkes 過程と Hawkes グラフおよび 2 つの推定方法を紹介し，3 節でシミュ

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

183

レーションを通して 2 つの手法の比較を行う．4 節では実データ解析として，仮想通貨市場の影響構造を最尤法を用いて解析するとともに，COVID-19 の感染構造の解析を Kirchner の手法で行う．5 節でまとめを行い，今後の展望について述べる．
2.

多変量 Hawkes 過程と Hawkes グラフ

本節では，多変量 Hawkes 過程および Hawkes グラフの定義を紹介した上で，Hawkes 過程のパラメータを推定する 2 つの方法である最尤法と Kirchner の手法を紹介する．
2.1 多変量 Hawkes 過程
適当なフィルター付き完備確率空間 (Ω, F , (Ft )t≥0 , P ) を考える．正の値をとる Ft - 停止時刻の増大列 {T (j) }j∈N を Ft - 適合な点過程と言う．T (j) を j 回目のイベントの発生時刻とみなし，T (j) < T (j+1) P -a.s. とする．さらに点過程 {T (j) }j∈N から定義される計数過程 N = (Nt )t≥0 をNt := j∈N 1{T (j) ≤t} とする．今後 N も点過程と呼ぶ．
Hawkes 過程は点過程のモデルであり，次の強度を用いて特徴付けられる．

定義 1.（点過程の強度）[0, ∞) 上の Ft - 適合な点過程 N を考える．正の値をとる Ft - 可予測な確率過程 λ = (λt )t≥0 が，任意の 0 ≤ s ≤ t を満たす s, t に対して

 t

E [Nt − Ns |Fs ] = Eλu du|Fs

P -a.s.

sを満たすとき，λ を N の (P, Ft ) - 強度という．

t[0, T ] 上の点過程 N が (P, Ft ) - 強度 λ を持つとする．任意の t に対して 0 λs ds < ∞ P -a.s. が成り立ち，λ のサンプルパスが左連続かつ右極限を持ち，任意の t ∈ [0, T ] について λt ≤ Y P -a.s.
を満たす可積分確率変数 Y が存在するとき，強度は次の表現を持つ（Aalen, 1978, Lemma 3.3）
．
lim

1

Δ→+0 Δ

E[Nt+Δ − Nt |Ft ] = λt+ P -a.s.

ただし，λt+ = limΔ→+0 λt+Δ である．この式において，N の強度 λt の右極限 λt+ は時刻 t ま
での情報が与えられた下での瞬間的なイベント発生確率と解釈することができる．
Hawkes 過程は過去のイベント発生が現在の強度に影響を及ぼす点過程モデルである．
定義 2.（単変量 Hawkes 過程）点過程 N が Hawkes 過程であるとは，強度 λ が th(t − s)N (ds)λt = η +−∞で与えられるものをいう．ただし，η ≥ 0 は基底強度，h(t) ≥ 0 はカーネル関数と呼ばれ t < 0に対しては h(t) = 0 とする．
定義 2 からわかるように，カーネル関数に h を持つ Hawkes 過程において，時刻 s で発生したイベントは時刻 t(> s) において h(t − s) だけ強度を押し上げる．これはイベントの自己励起性を表し，イベントが続いて起こりやすい地震の解析に Hawkes 過程が広く使われている理由である．
次に，この単変量 Hawkes 過程を多変量に拡張する．自然数 M を与え，各 i ∈ {1, 2, . . . , M }に対して，正の値をとる Ft - 適合な点過程 {Ti(j) }j∈N を考える．各 i ∈ {1, 2, . . . , M } に対して
(j){Ti }j∈N から定義される点過程 N i = (Nti )t≥0 を Nti := j∈N 1{T (j) ≤t} とし，これを要素にも
i

統計数理

184

第 69 巻 第 2 号 2021

つ N = (Nt = (Nt1 , . . . , NtM ))t≥0 を M - 変量点過程と言い，i ∈ {1, 2, . . . , M } をイベントのタイプと呼ぶ．多変量点過程 N の各要素 N i (i = 1, . . . , M ) の強度 λi = (λit )t≥0 は N i の (P, FtN ) - 強度である．ここで FtN は点過程 N の自然なフィルトレーションを表している．多変量 Hawkes
過程は自己励起性だけでなく相互励起性も考慮したモデルであり，ある時刻の強度が他の種類のイベントも含む過去の発生履歴に依存して定義される．
定義 3.（多変量 Hawkes 過程）M - 変量点過程 N が M - 変量 Hawkes 過程であるとは，多変量強度過程 Λ = (Λt = (λ1t , λ2t , . . . , λMt ))t≥0 が

 t

Λt = η +H(t − s)N (ds)−∞で与えられるものをいう．ただし，各タイプの基底強度 ηi を並べた η = (ηi )M
i=1 は基底強度ベ
クトル，タイプ j のイベント発生がタイプ i の強度に及ぼす影響を表すカーネル関数 hij (t) を要素に持つ H(t) = (hij (t))M
i,j=1 はカーネル関数行列と呼ばれる．また，各 i = 1, . . . , M に対
して

 t



M

H(t − s)N (ds)−∞ thij (t − s)N j (ds)

=
i

j=1

−∞

である．
カーネル関数は 1 回のイベント発生が引き起こす影響の時間変化を定義する関数で，Hawkes
過程を特徴付ける重要な関数である．従来の研究ではパラメトリックにカーネル関数が与えられることが多かった．特に指数型カーネル
h(t) = αe−βt , t ≥ 0が仮定されることが多い．このカーネル関数は 2 つのパラメータ (α, β) ∈ R2>0 = {(x1 , x2 ) ∈R2 |x1 > 0, x2 > 0} で決定される．ファイナンス分野においても指数型カーネルは広く利用さや Yang et al.（2018）は多変量 Hawkes 過程に対して次のカーネル関れており，Bowsher（2007）
数行列を用いた解析を行っている．
(2.1)(H(t))i,j = αij e−βij t , t ≥ 0

このカーネル関数行列を持つ M - 変量 Hawkes 過程は M (2M +1) 個のパラメータ (ηi , αij , βij ; i,j = 1, . . . , M ) で特徴付けられる．
これに対して，疫学統計の分野では影響の減衰の仕方についての知見が乏しく，カーネル関
数をパラメトリックに特定できないことがある．そのようなときに，カーネル関数を誤って特定すると間違った解析結果を与える危険がある．それを回避する方法である Kirchner の手法
では，カーネル関数行列 (hij (t))M
i,j=1 はパラメトリックに特定されず，代わりに十分小さな正数 Δ と十分大きな自然数 p に対して h(kΔ), k = 1, . . . , p の推定量 ĥ(kΔ) を非母数的に求める．
2.2 Hawkes グラフ
M - 変量 Hawkes 過程は，自己励起性と相互励起性の構造を含んでいる．これらの関係は複雑であり，どのタイプのイベントがどのタイプのイベントに重大な影響力を持っているのかを比較することは困難である．3 章で行うシミュレーションのようにカーネル関数行列の要素がそれぞれ異なるパラメトリック関数族に含まれる場合に比較が困難であることは当然であるが，（α と β がカーネル関数が全て指数型カーネルで式（2.1）の表現を持つときも，大きく一時的な影響同士の比較は容易ではない．そこで，各影
大きい）影響と弱く持続的な（α と β が小さい）

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

響の重大さを表す次の指標を用いる．
(2.2)

185

 ∞

aj,i =hij (t)dt
0aj,i はタイプ j の一度のイベント発生がタイプ i の強度に与える影響の総和である．また，式（2.2）から，aj,i はタイプ j の一度のイベント発生が誘発するタイプ i のイベント発生回数の期待値であると解釈することもできる．Embrechts and Kirchner（2018）が提案した Hawkes グラ
フは式（2.2）を用いて次のように定義される．

定義 4.（Hawkes グラフ）定義 3 で定義される M - 変量 Hawkes 過程 N を考える．このとき，
N の Hawkes グラフとは頂点に基底強度，辺に式（2.2）で定義される値を重みとして持つ重み付
き有向グラフ G := (V, E) であり，頂点集合 V と辺集合 E は次式で与えられる．



V = {(i; ηi ), i = 1, . . . , M },

E=



 ∞

i, j; ai,j =hji (t)dt , i, j = 1, . . . , M
0
α

ji
特に，式
（2.1）
で定義されるカーネル関数行列を持つ M - 変量 Hawkes 過程の場合は ai,j = βjiとなる．また，Kirchner の手法を用いるときは，カーネル関数行列の t = Δ, . . . , pΔ の時点で
pの値 hi,j (kΔ) , k = 1, . . . , p を用いて積分を離散近似することにより，ai,j = k=1 Δhji (kΔ) とする．

2.3 最尤法による推定と Hawkes グラフ表現[0, T ](0 < T < ∞) 上の M - 変量点過程 N における多変量強度過程 Λ がある有限次元空間
Θ ⊂ Rq (q ∈ N) 上のパラメータ θ ∈ Θ によって特定化され Λ(θ) = (Λt (θ) = (λ1t (θ), λ2t (θ), . . . ,λMt (θ)))t≥0 と表されるとする．このとき，任意の ω ∈ Ω に対する尤度関数は

 M 

(2.3)(1 − λis (θ|ω))ds +L(θ|ω) = exp
0

i=1





T

(0,T ]log λis (θ|ω)N i (ds|ω)で与えられる（Daley and Vere-Jones, 2003, Proposition 7.3. III）．特にカーネル関数行列が式
M
M
M（2.1）で与えられるとき，θ = (θi )M
i=1 = (ηi , (αij )j=1 , (βij )j=1 )i=1 の対数尤度関数は

 M 

log L(θ) = log

exp
i=1

M  T(1 − λis (θ))ds +(1 − λis (θ))ds +

=
0

i=1

0





T

(0,T ]log λis (θ)N i (ds)

(0,T ]log λis (θ)N i (ds)M=li (θi )i=1という表現を得る．ここで Tli (θi ) =0(1−λis (θ))ds+


(0,T ]M= (1−ηi )T −j=1 T (k) ∈(0,T ]j

であり，Rij (k) は再帰的にlog λis (θ)N i (ds)(k)αij
−β (T −Tj )(1−e ij
)+
βij

M

log
(k)

Ti

∈(0,T ]

ηi +αij Rij (k)
j=1

統計数理

186

第 69 巻 第 2 号 2021

⎧
(k)
(l)

(k)
(k−1)
−β (T
−Tj )
−βij (Ti −Ti)⎪
Rij (k − 1) + T (l) ∈[T (k−1), T (k) ) e ij i(i = j, k ≥ 1)⎪⎨e
Rij (k) =(k)−βii (Ti
⎪e

⎪
⎩

(k−1)

−Ti

j

)

i

i

(1 + Rii (k − 1))(i = j, k ≥ 1)

0

(k = 0)と定義する．
本論文においては準ニュートン法の一種である BFGS 法を用いて，この関数を最大化するパ
ラメータを求める．
M
M
上記の対数尤度関数から得られる最尤推定量 θ̂ = (η̂i , (α̂ij )M
j=1 , (β̂ij )j=1 )i=1 を ai,j の式に代入することで Hawkes グラフの辺の重みの推定量が次のように定まる．
âi,j =

α̂ji
β̂ji

本論文において Hawkes グラフの辺 (i, j), i, j = 1, . . . , M の有無は，辺の重みの推定量 âi,j の漸近分布の 95% 信頼区間に負の値が含まれるか否かで判断をする．前述の通り α̂ji , β̂ji は最尤推定量であるため，漸近的に正規分布にしたがう．これらの比の信頼区間を求めるために
Fieller（1932）
が提案した Fieller’s method を用いる．Fieller’s method は正規分布にしたがう 2
つの確率変数の比の信頼区間を求める方法で，(α, β) が 2 変量正規分布にしたがっているとし，その比を α/β = k とする．このとき，k が与えられたときの α − βk の条件付き分布はN(0, Var(α) − 2Cov(α, β)k + Var(β)k2 ) となるから，標準正規分布の下側 x% 点を Zx% とすると
P (|α − βk| ≤ |Z97.5%

Var(β)k2 − 2Cov(α, β)k + Var(α)|) = 0.95と表すことができる．この式の ( ) 内の両辺を 2 乗して k について整理すると，222{β 2 − Z97.5%
Var(β)}k2 − 2{αβ − Z97.5%
Cov(α, β)}k + {α2 − Z97.5%
Var(α)} ≤ 0となる．この k についての 2 次関数が下に凸で判別式が正の値をとる時，k についての範囲で書き直すことができるため，k の信頼区間の下限は2{αβ−Z97.5%
Cov(α, β)}−



2
2
2
{αβ−Z97.5%
Cov(α, β)}2 −{β 2 −Z97.5%
Var(β)}{α2 −Z97.5%
Var(α)}22
β −Z97.5% Var(β)

 ji ), Var(β ji ),となる．これに α̂ji , β̂ji および Var(αji ), Var(βji ), Cov(αji , βji ) の一致推定量 Var(αCov(αji , βji ) を代入することで信頼区間の一致推定量を得ることができる．ここで，最尤推定量の漸近分散共分散行列はフィッシャー情報行列であり，推定値の計算では，対数尤度関数の22
 ji ) < 0．ただし，β̂ji− Z97.5%Var(βヘッセ行列の逆行列を用いる（Bowsher, 2007, Theorem 2.2）のとき，k についての 2 次関数が上に凸となるため，信頼区間は [−∞, ∞] とする．
2.4 Kirchner の手法による推定と Hawkes グラフ表現
多変量 Hawkes 過程における未知の影響構造をノンパラメトリックに推定する Kirchner の手法を紹介する．なお，最尤法は連続時間観測されたサンプルパスに基づく推定手法であるが，Kirchner の手法は離散時間観測データでも可能な推定手法である．
基底強度 η ，カーネル関数行列 H を持つ多変量 Hawkes 過程にしたがう [0, T ] 上の点過程 Nを考える．N から生成される Δ 時間毎のイベント発生件数の列 X̃ (Δ) = {X̃k(Δ) }k=1,...,n を
(Δ)

X̃k

とする．ここで n := T /Δ である．

:= NkΔ − N(k−1)ΔHawkes 過程における 2 つの推定手法の比較と実データ解析への応用
(Δ)

(Δ)

187
(Δ)

X (Δ) = {Xk }k=1,...,n を生起パラメータ a0 := Δη ，間引き
（thinning）係数行列 Ak :=ΔH(kΔ), k ∈ N を持つ M 変量 INAR(∞) 過程にしたがう確率変数列とする．また，X (Δ,p) =
(Δ,p)
(Δ,p)
(Δ)
(Δ,p)
(Δ)
}k=1,...,n を生起パラメータ a0:= a0 ，間引き係数行列 Ak:= 1{k≤p} Ak , k ∈ N{Xkを持つ M - 変量 INAR(p) 過程にしたがう確率変数列とする．このとき，十分小さい Δ > 0 と十分大きい pΔ > 0 に対して
ddX̃ (Δ) ≈ X (Δ) ≈ X (Δ,p)という近似が成り立つことを Kirchner（2017）は主張している．この主張の下で，Δ 毎に離散観測された Hawkes 過程を INAR(p) 過程で近似し，CLS 法（Conditional Least Square method,
条件付き最小 2 乗法）で求めた推定量を用いて Hawkes 過程のカーネル関数行列の推定に帰着させることができる．以下では点過程 N から誘導される整数値確率過程を X (Δ) と表す．
定義 5.（Kirchner の手法による推定量）点過程 N が，基底強度ベクトル η ，カーネル関数
行列 H = (hij )M
i,j=1 を持つ [0, T ] 上の M - 変量 Hawkes 過程にしたがうとし，ある Δ > 0 に対し(Δ)
N が時間間隔 Δ で離散的に観測されているとする．この観測列を N (Δ) = {Nk }k=1,...,n と(Δ)(Δ)= {Xk }k=1,...,n をし，各時間間隔における N の増分過程として X(Δ)

Xk:= NkΔ − N(k−1)Δとする．ここで n := T /Δ である．このとき Δ < s < T を満たすあるサポート s を適当に選び p = s/Δ として，興味のあるパラメータを H (Δ,s) = (H(Δ), . . . , H(pΔ), η) ∈ RM ×(M p+1)と定める．この H (Δ,s) に対する推定量 Ĥ (Δ,s) = (Ĥ1(Δ,s) , . . . , Ĥp(Δ,s) , η̂ (Δ,s) ) をĤ (Δ,s) :=と定義する．ここで，

⎛

(Δ)

Xp
(Δ)
⎜ Xp−1⎜⎜ ..
Z := ⎜
.
⎜ (Δ)
⎝ X
1
1
(Δ)

(Δ)

Xp+1
(Δ)
Xp
..
.
(Δ)
X2
1

1
Y Z  (ZZ  )−1
Δ
···
···
..
.
···
···

(Δ)

⎞

Xn−1
(Δ)
Xn−2 ⎟⎟⎟..
(M p+1)×(n−p)
⎟∈R
.
⎟
(Δ)
Xn−p ⎠
1(Δ)

Y := (Xp+1 , Xp+2 , . . . , Xn(Δ) ) ∈ RM ×(n−p)

である．
(Δ,s)

Ĥk

(k = 1, . . . , p) は，カーネル関数行列の Δ 区間毎の推定量となり，⎛

ĥ11 (kΔ)⎜ ĥ21 (kΔ)
⎜
(Δ,s)
Ĥk
:= ⎜
..
⎝
.
ĥM 1 (kΔ)

···
···
..
.

ĥ12 (kΔ)
ĥ22 (kΔ)..
.
ĥM 2 (kΔ)···⎞

ĥ1M (kΔ)
ĥ2M (kΔ) ⎟
⎟
⎟
..
⎠
.
ĥM M (kΔ)

と表現される．Kirchner（2017）は p を固定した上で，n → ∞ としたときに vec(Ĥ) ∈ RM (M p+1)
が漸近正規性をもつことを示した．

d(n − p){vec(Ĥ) − vec(H)} → NM 2 p+M (0, V )approx.

11
この結果から，vec(Ĥ) ∼ NM 2 p+M (vec(H), n−pV ) と表すことができる．この S 2 = n−pV
2の経験バージョン Ŝ として，Kirchner（2017）は以下を提案している．

統計数理

188

Ŝ 2 :=

第 69 巻 第 2 号 2021


1((ZZ  )−1 IM )W ((ZZ  )−1 IM ) ∈ RM (M p+1)×M (M p+1)
2
Δ

ここで，
は ク ロ ネ ッ カ ー 積 ，IM は M 次 の 単 位 行 列 ，W :=RM (M p+1)×M (M p+1) ,
(Δ)

(Δ)

(Δ)

wk := (((Xk−1 ) , (Xk−2 ) , . . . , (Xk−p ) , 1)



n

k=p+1

wk wk ∈

p

IM )

(Δ)

Xk

(Δ)

− Δη̂ −

ΔĤ(lΔ)Xk−ll=1である．なお，本節の冒頭で述べた通り，Kirchner の手法は Hawkes 過程からの離散観測デー
タを用いて推定することができるが，上記の理論は観測時間間隔 Δ が十分小さい場合に成立する手法であり，Δ が十分小さいかどうかは実用の際には注意が必要である．
3.

シミュレーションによる精度の比較

本節では 2 つの状況下でのシミュレーションを通して，最尤法と Kirchner の手法による推定
精度の比較を行う．イベントタイプの数を M = 2 とし，基底強度ベクトル η およびカーネル
関数行列 H を定めた上で，長さ T のサンプルパスを 1000 回生成する．その上で Hawkes グラ
フのパラメータである (ηi )i=1,2 ，(ai,j )i,j=1,2 の推定値 (η̂i )i=1,2 ，(âi,j )i,j=1,2 を計算し，その推定精度を比較する．最尤法の場合は N のサンプルパスが [0, T ] の全区間で観測されているとし，カーネル関数は指数型カーネルを仮定する．一方，Kirchner の手法の場合では N のサンs カーネル関数行列の少なくプルパスが観測幅 Δ で離散的に観測されているとし，サポート （とも 1 つの成分が正となる t の上限）をある有限の値で固定した上で Δ を変化させ，その影響を見る．まず，Case 1 として多変量 Hawkes 過程の真のカーネル関数が指数型カーネルである
場合を想定してデータを生成し，最尤法と Kirchner の手法で推定を行う．次に Case 2 として真のカーネル関数が指数型カーネルでない場合のデータを生成し，最尤法では指数型カーネルを想定した上で 2 つの推定結果を比較する．
単調減少なカーネル関数を持つ Hawkes 過程の生成アルゴリズムは近江・野村（2019）に掲載されている．本節においては単調減少でないカーネル関数を持つ Hawkes 過程を考えるた
め，変更を加えた Algorithm 1 に基づいてサンプルパスの生成を行った．Hawkes 過程は，Ftまで条件付けたもとで，次のイベント発生までは強度が決定論的に定まる．また，ある時刻で発生したイベントのタイプが i である確率は発生時刻における強度の比で定義される．すなわち，タイプを区別しない k 回目のイベント発生時刻が与えられたとき，k + 1 回目のイベント発生時刻は非一様ポアソン過程にしたがい，そのイベントタイプの確率分布は強度の比に基づく．Algorithm 1 の Step. 3–10 は非一様ポアソン過程のサンプルパスの生成手法である
thinning theorem（Lewis and Shedler, 1979）を表しており，Step. 11–13 で得られたイベントに対してタイプを決定している．
Algorithm 1 Simulation of Multivariate Hawkes Processes
Ensure: T > 0 ∧ M ∈ N
1: t = 0, k = 0, ki = 0 (i = 1, . . . , M ), λi0 = ηi (i = 1, . . . , M ), T0 = 02: while t < T do
3:C=
λit + maxj=1,...,M maxt≥0 hij (t)i4:
以下の乱数を生成する．
τ ∼ Exp(C)x ∼ Uniform([0, 1])
5:

if

i

λit+τ

C

> x thenHawkes 過程における 2 つの推定手法の比較と実データ解析への応用

6:
7:
8:
9:
10:
11:

189

t←t+τ
else
k ←k+1t←t+τ
Ti ← t以下の乱数を生成する．
i ∼ Multi 1,

λit

ki ← ki + 1(k )13:
Ti i ← t
14:
end if
15: end while

i=1,...,M

12:

3.1 Case 1
基底強度ベクトル η ，カーネル関数行列 H を以下のように定めた 2 変量 Hawkes 過程のサンプルパスを生成する．
0.2

η=

(3.1)

0.1

,

H(t) =

0.3 exp(−2t)

0.5 exp(−1t)

0.4 exp(−2t)

0.2 exp(−1t)

生成する時系列の長さを T = 100, 500, 1000 とし，それぞれ Δ = 1.0, 0.5, 0.1 に設定した離散時間観測に基づく Kirchner の手法による推定と区間 [0, T ] での連続時間観測に基づく最尤法による推定を行い，Hawkes グラフの各辺および頂点の重みの推定値についてまとめたものが表 1である．式（3.1）のように H(t) を定めた場合，(ai,j )i,j=1,2 はそれぞれ ∞a1,1 =

 ∞0.3 exp(−2t)dt = 0.15,

a1,2 =0 ∞a2,1 =

0.4 exp(−2t)dt = 0.20 ∞0.5 exp(−1t)dt = 0.5,

0a2,2 =

0.2 exp(−1t)dt = 0.2
0となる．第 k 番目のパスに対するこれらの推定値を (â(k)i,j )i,j=1,2 (k = 1, . . . , 1000) とする．表 1
1000 (k)¯i,j = 1の値は 1000 回のシミュレーションによる平均値 ââ を表し，( ) 内の数値は1000
k=1 i,j1
標準誤差 SEi,j = 1000



1000 (k)¯i,j )2 を表している．Kirchner の手法におけるサポート(âi,j − â
k=1

は s = 6 に固定した．
表 1 から，T が大きく，Δ が小さくなるほど，Kirchner の手法の場合の標準誤差が小さくなる傾向にあり，T = 1000, Δ = 0.1 のときの推定値が真値の近くに分布する傾向にあることがわかる．Kirchner の手法による推定は離散時間観測であり，最尤法に比べて情報量が少ないため，最尤法に比べると標準誤差が大きくなる傾向にある．それに加えて，Kirchner の手法によ（バイアス）るカーネル関数の積分値 (ai,j )i,j=1,2 の推定結果は真の値に対して負の方向の偏りを持っているように見える．これは，カーネル関数のサポートが本来は s = ∞ であるところを
s = 6 に仮定して推定しているためと考えられる．これに対して，最尤法による (ai,j )i,j=1,2 の推定量は不偏性があるように見える．これらの結果から，カーネル関数がパラメトリックな関数として想定できる場合は，最尤法を用いたほうが推定精度が高い傾向にあることがわかる．
表 2 は T = 1000, Δ = 0.1 の Kirchner の 手 法 と 最 尤 法 の そ れ ぞ れ の MSE (MSEi,j =1000 (k)1(âi,j − ai,j )2 ) とその比である Relative Eﬃciency（ 最尤法の MSEi,j /Kirchner の手
1000
k=1
法の MSEi,j ）をまとめた表である．
MSE の意味においても最尤法の方が精度が高い傾向にあることがわかる．

190

統計数理

第 69 巻 第 2 号 2021

表 1．Case 1. 真のカーネル関数が指数型カーネルである場合のシミュレーション結果．1000回のシミュレーションによる ηi (i = 1, 2) および ai,j (i, j = 1, 2) の推定値の平均（標準偏差）を記載している．上 3 段は Kirchner の手法でそれぞれ Δ = 1.0, 0.5, 0.1 の場合の結果であり，下段は最尤法の結果である．平均が真値に最も近いものを太字で表し，標準誤差が最も小さかったものに下線を引いた．

表 2．T = 1000, Δ = 0.1 の Kirchner の手法と最尤法の MSE と Relative Eﬃciency（最尤法の MSE/Kirchner の手法の MSE）．

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

191

3.2 Case 2
次に，Kirchner（2017）と同様にして，基底強度ベクトル η ，カーネル関数行列 H を以下のように定めた 2 変量 Hawkes 過程のサンプルパスを生成する．
(3.2)

0.5

η=

0.25

,

H(t) =

0

0.5(1+t)21{1<t<3} 0.251{t<π} 0.25 sin(t)生成する時系列の長さを T = 100, 500, 1000 とし，それぞれ Δ = 1.0, 0.5, 0.1 に設定した離散時間観測に基づく Kirchner の手法による推定と区間 [0, T ] での連続時間観測に基づく最尤法による推定を行い，Hawkes グラフの各辺および頂点の重みの推定値についてまとめたのが表 3 である．式（3.2）のように H(t) を定めた場合，(ai,j )i,j=1,2 はそれぞれ ∞a1,1 =

 ∞0dt = 0,a1,2 =0 ∞a2,1 =
01{1<t<3} 0.25dt = 0.5
00.5
dt = 0.5,(1 + t)2 ∞a2,2 =

1{t<π} 0.25 sin(t)dt = 0.50¯i,j を表し，( ) 内の数値は標準となる．表 3 の値は 1000 回のシミュレーションによる平均値 â
表 3．Case 2. 真のカーネル関数が指数型カーネルでない場合のシミュレーション結果．1000回のシミュレーションによる ηi (i = 1, 2) および ai,j (i, j = 1, 2) の推定値の平均（標準偏差）を記載している．上 3 段は Kirchner の手法でそれぞれ Δ = 1.0, 0.5, 0.1 の場合の結果であり，下段は最尤法の結果である．平均が真値に最も近いものを太字で表し，標準誤差が最も小さかったものに下線を引いた．

統計数理

192

第 69 巻 第 2 号 2021

図 1．Kirchner の手法で推定したカーネル関数 (ĥij (t))i,j=1,2 を区間 [0, 6] で表示した結果．
ただし，T = 1000, Δ = 0.1 としている．

誤差 SEi,j を表している．Kirchner の手法におけるサポートは s = 6 に固定した．
表 1 と同様，表 3 からも T が大きく，Δ が小さくなるほど，Kirchner の手法の標準誤差が小さくなる傾向にあり，T = 1000, Δ = 0.1 のときの推定値が真値の近くに分布する傾向にある
ことがわかる．理論的には Kirchner の手法における各推定量は，T → ∞, Δ → 0 のもとで一致性を持つため，この結果がシミュレーションに表れていることがわかる．一方，最尤法を考えると，a1,1 の推定に大きな偏り（バイアス）が生じているように見える．このことの詳細を調べるため，T = 1000, Δ = 0.1 の下で，Kirchner の手法で推定したカーネル関数を区間 [0, 6] で表示した結果が図 1 であり，連続時間観測のもとで最尤法で推定したカーネル関数を区間 [0, 6]で表示した結果が図 2 である．
図 1 は，推定結果が真のカーネル関数の周りに分布しているのに対し，図 2 では，h11 だけでなく，h21 や h22 も正しく推定できていないことがわかる．これは真のカーネル関数とは異なる指数型の関数を当てはめているためである．この結果から，誤ったパラメトリック関数族を仮定することで推定結果に大きなバイアスが生じる恐れがあることがわかった．一方，Kirchnerの手法は，カーネル関数を Δ 間隔毎に推定しており，推定値は図 1 にあるように不安定となっているが積分の近似値である ai,j の推定値の標準誤差は比較的安定している．これは，スペク√
トル密度関数 f (λ) のピリオドグラムを用いた推定量 fˆn (λ) の一致性のオーダーは一般には n√より低いが，その積分汎関数は n - 一致性をもつ（谷口, 2005）ことと類似している．このことから，カーネル関数の推定に関しては推定精度が悪く，オーバーフィッティングも懸念されるが，積分値である ai,j の推定については比較的精度の良い推定ができていると思われる．
4.

実データ解析例

本節では，2 つの多変量計数データに対して Hawkes 過程への当てはめを行う．
（https://bitﬂyer.com/ja-jp/（ 最終アクセス 2019 年 8 月 15
まず 1 つ目のデータは，bitFlyer
日））の 15 分足データである．このデータは離散時間観測データであるが，15 分毎の高頻度

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

193

図 2．最尤法で推定したカーネル関数 (ĥij (t))i,j=1,2 を区間 [0, 6] で表示した結果．ただし，
T = 1000 としている．

データであり，ほとんど連続時間観測データと考えられること，Bacry（2015）や Yang et al.
（2018）などの先行文献で扱っているデータと類似していることから，指数型カーネルを使って最尤法によりパラメータ推定を行う．
2 つ 目 の デ ー タ は ，日 本 国 内 に お け る COVID-19 の デ ー タ（https://gis.jag-japan.com/「都道府県別新型コロナウイルス感染者数マップ（ジャッグジャパン株式会社
covid19jp/，提供）」）である．入手できたデータは日毎に観測された離散時間観測データであること，影響
構造に関して少なくても現時点ではカーネル関数族を特定化するほどの共通認識を持っていないことから，Kirchner の手法によりノンパラメトリックにカーネル関数を推定する．
4.1 仮想通貨の値動きについての最尤法による推定と Hawkes グラフ表現
ここでは，bitFlyer のデータを解析する．解析対象とする期間は 2019 年 6 月 16 日 00:00 から 2019 年 7 月 25 日 00:00 の 39 日間である．この解析ではイベント発生を観測値がある閾値
を超えた場合とする．このイベント発生を定義する閾値の設定のために解析期間より広い 2019年 6 月 6 日 00:00 から 2019 年 7 月 25 日 00:00 の 49 日間を観測期間とした．図 3 は観測期間に（対数変換し前期との差を取ったもの）の推移をプ
おける 15 分毎の価格およびその対数成長率
ロットしている．この対数成長率の絶対値の上側 80% 点を超過した時刻を価格の変化のイベ
ント発生時刻とし，その中でも正の変化を “価格の上昇” を表すタイプ 1 のイベント発生時刻と定義し，負の変化を “価格の低下” を表すタイプ 2 のイベント発生時刻と定義する．
（前期との差を取ったもの）の推
図 4 は観測期間における 15 分毎の出来高およびその変化量移をプロットしている．図 3 と同様にして，変化量の絶対値の上側 80% 点を超過した時刻を
出来高の変化のイベント発生時刻とし，その中でも正の変化を “出来高の上昇” を表すタイプ 3のイベント発生時刻と定義し，負の変化を “出来高の低下” を表すタイプ 4 のイベント発生時
刻と定義する．
解析対象となる期間内の 4 つのタイプのイベント発生時刻で線を引いたものが図 5 であり，（タイプ 2；i = 2），左下が出来高の上左上が価格の上昇（タイプ 1；i = 1），右上が価格の低下

194

統計数理

第 69 巻 第 2 号 2021

図 3．観測期間における 15 分毎の価格およびその対数成長率の推移のプロット．下図の閾値は対数成長率の絶対値の上側 80% 点であり，これを超過した時刻を価格の変化のイベ
ント発生時刻とし，その中でも正の変化を “価格の上昇” を表すタイプ 1 のイベント発
生時刻，負の変化を “価格の低下” を表すタイプ 2 のイベント発生時刻とする．

図 4．観測期間における 15 分毎の出来高およびその変化量の推移のプロット．下図の閾値は変化量の絶対値の上側 80% 点であり，これを超過した時刻を出来高の変化のイベント
発生時刻とし，その中でも正の変化を “出来高の上昇” を表すタイプ 3 のイベント発生
時刻，負の変化を “出来高の低下” を表すタイプ 4 のイベント発生時刻とする．

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

195

図 5．タイプ毎のイベント発生時刻のプロット．イベント発生時刻に線が引かれている（左，右上：価格の低下（タイプ 2；i = 2）
，左下：出来高
上：価格の上昇（タイプ 1；i = 1）（タイプ 4；i = 4））．の上昇（タイプ 3；i = 3），右下：出来高の低下

昇（タイプ 3；i = 3），右下が出来高の低下（タイプ 4；i = 4）のイベントを表している．それぞれポアソン過程にしたがう点過程には通常見られない “ムラ” が確認できる．これはあるイベ
ント発生が同種のイベントを誘発することによって生じるイベント発生クラスターであり，自己励起性を持っていることが推察される．また，タイプ 1 とタイプ 2 は特に近い時刻に “ムラ”が発生しているように見えるが，これはこの 2 つのイベントの間に相互励起性が存在することを示唆している．
4.1.1 パラメータの推定結果
図 5 のデータが N = (Nt = (Nt1 , . . . , NtM ))t∈[0,T ] から観測されたものとする．ただし，Nti
は，i = 1 のときタイプ 1（価格の上昇），i = 2 のときタイプ 2（価格の低下），i = 3 のときタ
（ 出来高の上昇），i = 4 のときタイプ 4（ 出来高の低下）の各イベント発生回数を表し，
イプ 3
M = 4, T = 39 とする．このとき，式（2.3）の尤度関数 L(θ) を最大にするようにして求めたパラメータ θ の最尤推定値 θ̂M L の推定結果は表 4 のようになった．
表の値は最尤推定値と推定量の標準誤差を表している．この標準誤差に用いる漸近分散の推 2 log L(θ)−1
|θ=θ̂M Lの対角成分で計算している．
定値は L(θ) のヘッセ行列を用いて T1 ∂ ∂θ∂θ

表 4 左上において α13 の推定値が相対的に大きいことから，“出来高の上昇” が “価格の上昇”の強度を大きく押し上げる傾向にあることがわかる．しかし対応する β13 の推定値がかなり大きいことから，その影響は限定的であることが推測される．一方，α12 の推定値は α13 ほど大きくないが，β12 の推定値が小さいため，押し上げられた強度が緩やかに減衰する．したがってその影響は持続する傾向にあることが推察される．同様にして，表 4 右上において α23 の推定値が相対的に大きいが，β23 の推定値がかなり大きいため，その影響は限定的と思われる．
一方で，α22 の推定値はそれほど大きくないが，β22 の推定値も小さいことから，その影響は持続することが推測される．表 4 左下の場合，α31 , α32 と β31 , β32 の推定値がそれぞれ大きいことから，“出来高の上昇” が受ける影響は，価格の変化によるものが大きいがその影響は瞬間的であることがわかる．表 4 右下において α43 の推定値が相対的に大きく，β43 の推定値はそれほど大きくないため，“出来高の上昇” が “出来高の低下” に与える影響が顕著に現れることが

196

統計数理

第 69 巻 第 2 号 2021

表 4．bitFlyer のデータに Hawkes 過程を最尤法で推定したときのパラメータの推定値と標（i = 2）が受ける影
準誤差（左上：価格の上昇（i = 1）が受ける影響，右上：価格の低下（i = 4）が受ける
響，左下：出来高の上昇（i = 3）が受ける影響，右下：出来高の低下影響）．

推測される．また，η については，η3 の推定値が η1 , η2 , η4 の推定値に比べて大きい値になっており，“出来高の上昇” は過去のイベントに関係なく発生しやすいイベントであることが推測される．
4.1.2 Hawkes グラフ
4.1.1 節で考察したように，パラメータの推定結果から，イベント発生時に他のイベント発生
へ与える影響の特徴を大まかに把握することができた．しかし，異なる特徴を持つ影響同士を総合的に比較することは表からは困難なことが多い．Hawkes グラフを描くことで，影響構造
を視覚的に表現することができ，前述の問題点も解消することができる．各有向辺 (i, j) の重みの推定値は 2.3 節にある通り，α̂ji /β̂ji で与えられ，表 5 にまとめられている．それぞれの信頼区間は Fieller’s method を用いて算出した．
これらの値から描いた Hawkes グラフが図 6 である．なお，有向辺の重みの 95% 信頼区間
「ai,j > 0」に負の値を含んだ辺は存在しないものとした．これは帰無仮説「ai,j = 0」，対立仮説の片側検定を有意水準 2.5% で行っていることと同等であり，帰無仮説が棄却されたものがグ
ラフ上に有向辺として描かれている．図 6 において頂点と有向辺の重みは値で記されるとともに，それぞれ円の面積，有向辺の太さで表現されている．Hawkes グラフから視覚的にわかる
ように，“出来高の上昇” が “出来高の低下” の強度に与える影響が非常に重大であることがわかる．同様に，“価格の低下” が “価格の上昇” を引き起こしやすいことを表している．また，“価格の低下” と “出来高の上昇” のみが自己励起性を持っていると推定された．
4.1.3 考察
Hawkes グラフ表現から，“出来高の上昇” は独立に発生しやすいイベントであり，これに対して “出来高の低下” は独立には発生しづらいイベントであるが，一度 “出来高の上昇” が発生すると “出来高の低下” を誘発しやすくなることが推察される．逆に，“出来高の低下” が “出来Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

表 5．有向辺の重みの 95% 信頼区間．信頼区間は Fieller’s method を用いて算出した．

図 6．bitFlyer のデータに関する Hawkes グラフ表現．頂点の円の面積は η̂i (i = 1, 2, 3, 4)の大きさに連動しており，有向辺の太さは âi,j (i, j = 1, 2, 3, 4) の大きさに連動している．有向辺の有無は 95% 信頼区間に負の値を含むか否かで判定している．

197

198

統計数理

第 69 巻 第 2 号 2021

高の上昇” を誘発することは考えにくく，“出来高の上昇” がさらなる “出来高の上昇” を誘発する自己励起性の傾向が見られた．一方，価格の上昇と低下はお互いに影響を及ぼしあうが，低下が上昇を誘発する影響の方が総合的に重大であることが推察される．“価格の低下” については，独立には発生しづらいイベントであるが，低下がさらなる低下を誘発する自己励起性が見られた．また，出来高と価格の関係性については，出来高の上昇と価格の変化の間に相互励起構造があると推定されたが，その程度は限定的であり，顕著な関係があるとは言えないと思われる．
4.2 COVID-19 の感染拡大についての Kirchner の手法による推定と Hawkes グラフ表現
次に，Kirchner の手法による解析を日本国内における COVID-19 のデータを用いて行う．
COVID-19 は 2019 年末に初めて観測された新しいウイルスであり，報道によれば，中国武漢から世界中に広がり大きな影響を与えている．主な感染経路は人から人であると考えられ，日本
国内でも感染が広がっている．2020 年 3 月 24 日には東京オリンピックの延期が決定され，同年 4 月 5 日に東京の累積感染者が 1,000 人を突破し，その 2 日後には埼玉，千葉，東京，神奈
川，大阪，兵庫，福岡を対象に緊急事態宣言が発令された．その後 4 月 16 日に緊急事態宣言の対象は全国に拡大され，5 月 25 日までに全ての地域で解除されたが，2020 年 12 月現在第 3波と呼ばれる流行により東京は 1 日 500 人以上の新規感染者を報告する日もあり，急速な再拡大が懸念されている．このような感染症の解析例として Park et al.（2020）は Hawkes 過程と伝統的な疫学統計モデルである SEIR モデルの比較を西アフリカにおけるエボラウイルスの感染
データを用いて行っている．この論文では，ある地域内の感染者数の時系列データを，単変量
Hawkes 過程に当てはめている．すなわち，この論文での主たる興味は，“エボラウイルスへの感染” というイベントがその後の感染者の発現に影響を及ぼす自己励起性の有無である．また，その影響の減衰の仕方についての知見が十分でないなどの理由から，カーネル関数はノンパラ
メトリック法で推定されている．本解析においても，COVID-19 の感染波及構造についても事
前知識が十分でなく，データも離散観測であるため，カーネル関数はノンパラメトリック法の1 つである Kirhcner の手法で推定する．さらに，地域間の影響構造の有無に興味があるため，
本解析においては地域毎の日毎の新規感染者数（PCR 検査の偽陽性・偽陰性の問題や無症状者の問題があるため，感染者数と定義することは適切とは言えないが，本研究ではデータに基づ
く解析を行いたいため，便宜上感染者数と呼ぶこととする）を 2.4 節で定義した多変量 INAR 過と仮定する．特に，データから定義した国内
程からの観測列 X (Δ)（ただし，Δ = 1 日とする）
における流行の第 1 波と第 2 波の影響構造の変化の有無に着目する．
観測期間は 2020 年 1 月 1 日から 2020 年 10 月 10 日の 284 日間で，解析対象とする地域は東京，埼玉，神奈川，大阪，兵庫，京都の 1 都 2 府 3 県とする．これらの各地域を関東の東京，と関西の大阪，兵庫，京都（以後，第 2 グループと埼玉，神奈川（以後，第 1 グループと呼ぶ）
呼ぶ）の 2 つのグループに分け，各グループ内の自己・相互励起性の有無を調べる．観測期間内の国内新規感染者数（確定日で集計）をプロットしたものが図 7 である．図 7 には 7 日間の移
動平均
（7 項移動平均，以降 7 日移動平均と呼ぶ）が太い線で描かれている．本解析では，国内新規感染者数の 7 日移動平均が初めて 100 を超えた 2020 年 3 月 26 日から，次に 100 を下回った 2020 年 5 月 7 日までの期間を流行の第 1 波とする．同様に，その後 7 日移動平均が再び 100を超えた 2020 年 6 月 27 日を第 2 波の始まりとするが，7 日移動平均が算出できる 2020 年 10
月 7 日まで 100 を下回っていないため，第 2 波の終わりを 2020 年 10 月 7 日と定義する．したがって，本解析において第 1 波とする期間は 43 日間で，第 2 波とする期間は 104 日間である．
第 1 グループ，第 2 グループの各都府県における新規感染者数をプロットしたものが図 8，9である．図 10 は，第 1 波における，第 2 グループの各地域間の相互相関関数（cross-correlation

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

199

図 7．2020 年 1 月 1 日から 2020 年 10 月 10 日の日本の新規感染者数．太線は 7 日移動平均を表し，これが 100 を上回った 2020 年 3 月 26 日から 2020 年 5 月 7 日までの期間と，2020 年 6 月 27 日以降をそれぞれ流行の第 1 波，第 2 波と定義する．

図 8．第 1 グループの新規感染者数．粗い破線が東京，実線が埼玉，細かい点線が神奈川の新規感染者数を表す．

function; ccf）
を描いたものである．左図は大阪と京都の新規感染者数の時刻のズレに関する影響度合いを表す．同様にして，中図は兵庫と大阪の影響度合い，右図は兵庫と京都の影響度合
いを表す．
図 10 の左図および右図から非対称性が見られる．特に，左図からは大阪と過去の京都との相関関係が見られ，右図からは兵庫と過去の京都との相関関係が存在することが示唆される．
解析するデータは 1 日毎の新規感染者数であるから，観測幅 Δ = 1 の離散観測と考える．

200

統計数理

第 69 巻 第 2 号 2021

図 9．第 2 グループの新規感染者数．粗い破線が大阪，実線が京都，細かい点線が兵庫の新規感染者数を表す．

図 10．第 1 波における第 2 グループの各地域間の相互相関関数（cross-correlation function;ccf）．左図は大阪と京都の新規感染者数の時刻のズレに関する影響度合いを表す．同様にして，中図は兵庫と大阪の影響度合い，右図は兵庫と京都の影響度合いを表す．

Kirchner の手法を適用するにあたり，影響の残る期間を定める p を与える必要があるが，本解に則って，以下で与えられる AIC が最小となる p を適用する．
析においては Kirchner（2017）
AIC(p) := log(det Σ̂(p)) +ここで Σ̂(p) は 2.4 節で準備したベクトル wk を用いて

2pM 2
T −p

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用Σ̂(p) :=

1
T −p

201

T

wk wk
k=p+1とし，T は第 1 波の場合は 43，第 2 波の場合は 104 とする．
4.2.1 推定結果
第 1 グループ（東京，埼玉，神奈川）における第 1 波と第 2 波の推定結果の Hawkes グラフ表
現は，それぞれ図 11，12 のようになった．図 11，12 においても，4.1.2 節と同様に帰無仮説，対立仮説「ai,j > 0」の仮説検定を行い，有意水準 2.5% で帰無仮説が棄却されたもの「ai,j = 0」のみ有向辺として描いている．また，頂点，有向辺の重みの表現方法は図 6 と同様であるが，重みの推定値に加えて 95% 信頼区間も併記している．図 11 が第 1 波の Hawkes グラフ表現，
図 12 が第 2 波の Hawkes グラフ表現である．
各都県における頂点の大きさから独立な感染者の発生確率を表す η の推定値は，どちらの期間においても，東京が神奈川と埼玉に比べて大きいことがわかる．ただし，η の推定値の解釈は本来であれば独立な感染者の発生確率を表すべきであるが，本解析においては陽性と判明していない感染者からの感染やグループ外の地域からの感染者の流入による影響が反映されていると推察される．また，2 つのグラフを比較すると，第 1 波の期間は相互励起性が大きく推定されているように見えるが，第 2 波の期間は相互励起性の影響はほとんど見られないように見える．一方で，第 1 波の期間では自己励起性は埼玉のみに見られたが，第 2 波の期間では全ての地域で自己励起性があると推定されている．以上のことから，第 1 グループの第 2 波における影響構造は，第 1 波の期間に存在した相互の影響構造から，各地域の自己励起構造に変化したことが推察される．
の推定結果の Hawkes グラフ表現は図 13，14 である．
第 2 グループ（大阪，兵庫，京都）各府県における頂点の大きさから独立な感染者の発生確率を表す η の推定値は，3 府県の間に第 1 グループのような大きな違いは見られなかった．また，2 つのグラフを比較すると，相互励起性に関しては第 1 波では京都から兵庫への弱い相互励起性が見られたのみで大阪とは相互励起性が見られなかったが，第 2 波では京都から大阪への強い相互励起性と大阪から兵庫への弱い相互励起性が見られた．自己励起性については第 1 波では兵庫と大阪に自己励起性が見られたが，第 2 波では京都に自己励起性が見られ，第 1 波と第 2 波での影響構造が大きく変化したことが推察される．
なお，AIC 基準で選択された p は，第 1 波の第 1 グループと第 2 グループがそれぞれ p = 2, 1
で，第 2 波においてはそれぞれ p = 5, 9 であった．
第 1 波と第 2 波で影響構造の変化があることが確認できたが，次に影響構造と人の移動構造の間の関係を調べるためにさらに解析を行う．ここでは日常的な人の移動の流れを把握するた
めに通勤・通学を考える．対象とするデータは東京 23 区の区毎の感染者発生時系列である．
人の移動についてのデータは平成 27 年国勢調査における居住区毎の通勤・通学地別人口を用いた．解析は 23 区を 5 つの地区に分けて行った．地区と該当する区は以下の通りである．
• 都心地区：千代田区，中央区，港区，新宿区，渋谷区，文京区• 城西地区：中野区，杉並区• 城南地区：品川区，目黒区，大田区，世田谷区• 城北地区：豊島区，北区，板橋区，練馬区• 城東地区：台東区，墨田区，江東区，荒川区，足立区，葛飾区，江戸川区

解析対象とした期間は 2020 年 4 月 1 日から 2020 年 5 月 14 日の 44 日間である．平成 27 年国

202

統計数理

第 69 巻 第 2 号 2021

図 11．第 1 グループの第 1 波での Hawkes グラフ表現．頂点の円の面積は η̂i (i = 1, 2, 3) の大きさに連動しており，有向辺の太さは âi,j (i, j = 1, 2, 3) の大きさに連動している．
有向辺の有無は 95% 信頼区間に負の値を含むか否かで判定している．

図 12．第 1 グループの第 2 波での Hawkes グラフ表現．頂点の円の面積は η̂i (i = 1, 2, 3) の大きさに連動しており，有向辺の太さは âi,j (i, j = 1, 2, 3) の大きさに連動している．
有向辺の有無は 95% 信頼区間に負の値を含むか否かで判定している．

勢調査による地区から地区への人の移動を表す有向グラフが図 15 に示されている．ここで，人の移動を表すグラフの有向辺の重みは “終点に通勤・通学している人の中で，始点の地区に住んでいる人の割合” である．なお，煩雑さを避けるため，居住地区と通勤・通学先の地区が一致している人口に対応する辺は描いていない．
図 16 は都内を上記の地区毎に分けたデータから得られた推定結果の Hawkes グラフ表現を，対立仮説「ai,j > 0」の仮説検定を行い，表している．なお，このグラフのみ帰無仮説「ai,j = 0」

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

203

図 13．第 2 グループの第 1 波での Hawkes グラフ表現．頂点の円の面積は η̂i (i = 1, 2, 3) の大きさに連動しており，有向辺の太さは âi,j (i, j = 1, 2, 3) の大きさに連動している．
有向辺の有無は 95% 信頼区間に負の値を含むか否かで判定している．

図 14．第 2 グループの第 2 波での Hawkes グラフ表現．頂点の円の面積は η̂i (i = 1, 2, 3) の大きさに連動しており，有向辺の太さは âi,j (i, j = 1, 2, 3) の大きさに連動している．
有向辺の有無は 95% 信頼区間に負の値を含むか否かで判定している．

有意水準 15% で帰無仮説が棄却されたもののみ有向辺として描いている．2 つの図を見比べると，人の移動が COVID-19 の感染構造と傾向が類似しているように見える．例えば，都心部は城東地区および城南地区からの流入が多いが，都心部への影響を持つと推定されたのはこの 2地区である．また，城西地区への流入は城北地区からが最も多いが，影響構造も同様な特徴を持つと推定された．一方，城東地区への流入は比較的少ない（自地区内での移動が約 85%）が，

204

統計数理

第 69 巻 第 2 号 2021

図 15．平成 27 年国勢調査による人の移動．有向辺の重みは “終点に通勤・通学している人の中で，始点の地区に住んでいる人の割合” であり，居住地区内での通勤・通学割合を表す有向辺は省略している．

図 16．都内の推定結果の Hawkes グラフ表現．都内を 5 つの地区毎に分けた新規感染者数
データから Kirchner の手法により推定した結果を Hawkes グラフで表現している．

Hawkes 過程における 2 つの推定手法の比較と実データ解析への応用

205

他地区からの影響は受けないと推定された．
4.2.2 考察第 1 グループの場合，第 1 波で見られた相互励起性を持った影響構造が，第 2 波では各都県
で自己励起性を持った影響構造に変化したことが推察される．逆に，第 2 グループの場合，第1 波では相互励起性は弱かったが，第 2 波では大阪を中心とした影響構造に変化したことが推察される．また，東京 23 区についての解析により，人の移動と影響構造が類似していることがわかった．これらの結果から，人の移動を制御することによる感染状況の制御の有効性が示
唆される．
5.

おわりに

本研究では，従来 Hawkes 過程の推定に用いられてきた最尤法による推定手法と，Kirchner（2017）が提案したノンパラメトリックな推定手法を比較し，応用例として実データを用いた解析を行った．2 つの状況を想定したシミュレーションを通してそれぞれの手法の使い分けが重要であることが確認された．真のカーネル関数が予想できる場合には従来の最尤法による推定の方が MSE の観点で精度が高い傾向にあることがわかった．一方で，真のカーネル関数と異なる関数族を想定した場合，最尤法による推定は推定精度が落ちてしまうため，モデルの誤特定のリスクが大きい状況では Kirchner の手法を用いるべきであると思われる．
実データを用いた解析として，仮想通貨市場の影響構造の解析と COVID-19 の感染構造の解析を行った．仮想通貨市場の解析においては，Bowsher（2007）と同様に指数型カーネルを想定した最尤法を用いて推定を行い，出来高の変化と価格の変化の間の影響構造を推定した．推定結果の Hawkes グラフ表現から，価格の変化と出来高の変化に関して，自己励起性および相互励起性が存在することが示唆された．一方，COVID-19 のデータに対して Kirchner の手法を用いて Hawkes グラフを構築した．近隣の都府県で構成した 2 つのグループにおける COVID-19の感染構造を推定した Hawkes グラフ表現からは，第 1 波と第 2 波で影響構造の変化が観察された．また，東京 23 区を地区別にグループ分けして得られた Hawkes グラフと平成 27 年国勢調査が示す通勤・通学の流れの間に類似する点があることが確認された．
本研究においてはイベントの発生時刻にのみ注目して解析を行った．しかし，どちらの解析においてもイベントの規模を考慮することはより的確な影響構造の推定を与えるかもしれない．例えば，本研究では同種のイベントとして認識した “価格の上昇” だが，閾値と同程度の上昇よりも閾値を大きく上回る上昇が与える影響の方が投資家に与える影響が大きいことが予想される．また，COVID-19 の感染者の中でも，軽症者のグループと重症者のグループは影響
構造が異なることが予想される．今後は，このような “閾値を越えた量” や “感染発覚時点での重症度” をマークとして持った Hawkes 過程の解析を行いたい．
謝辞

本稿に対して多くの有益なコメントを下さった編集委員および査読者の皆様に深く感謝いた
します．また，本研究は JSPS 科研費 JP16K00036 の助成を受けたものです．

206

統計数理

参

第 69 巻 第 2 号 2021

考

文

献

Aalen, O. (1978). Nonparametric inference for a family of counting processes, The Annals of Statistics,
6, 701–726.
Bacry, E., Mastromatteo, I. and Muzy, J. F. (2015). Hawkes processes in ﬁnance, Market Microstructureand Liquidity, 1(01), 1550005.
Bowsher, C. G. (2007). Modelling security market events in continuous time: Intensity based, multivariate point process models, Journal of Econometrics, 141, 876–912.
Daley, D. J. and Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes: Volume I :Elementary Theory and Methods, 2nd ed., Springer-Verlag, New York.
Embrechts, P. and Kirchner, M. (2018). Hawkes graphs, Theory of Probability and Its Applications,62, 132–156.
Fieller, E. C. (1932). The distribution of the index in a normal bivariate population, Biometrika, 24,
428–440.
Garetto, M., Leonardi, E. and Torrisi, G. L. (2021). A time-modulated Hawkes process to model the
spread of COVID-19 and the impact of countermeasures, arXiv preprint arXiv:2101.00405.
Gong, T., Chen, Y. and Zhang, W. (2021). An environmentally-adaptive hawkes process with anapplication to COVID-19, arXiv preprint arXiv:2101.09942.
Hawkes, A. G. (1971). Spectra of some self-exciting and mutually exciting point processes, Biometrika,
58, 83–90.
監物輝夫, 中川秀敏 (2019). 多次元 Hawkes 過程を用いた倒産リスク伝播構造の推定 Hawkes グラフ表
現による可視化, ジャフィー・ジャーナル, 17, 15–44.
Kirchner, M. (2017). An estimation procedure for the Hawkes process, Quantitative Finance, 17,
571–595.
Lesage, L. (2020). A Hawkes process to make aware people of the severity of COVID-19 outbreak:Application to cases in France, Doctoral Dissertation, Université de Lorraine; University of
Luxembourg.
Lewis, P. W. and Shedler, G. S. (1979). Simulation of nonhomogeneous Poisson processes by thinning,
Naval Research Logistics Quarterly, 26(3), 403–413.
Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes, Journal of the American Statistical Association, 83, 9–27.
近江崇宏, 野村俊一 (2019).『点過程の時系列解析』, 共立出版, 東京.
Park, J., Chaﬀw, A. W., Harrigan, R. J. and Schoenberg, F. P. (2020). A non-parametric hawkes modelof the spread of ebola in west africa, Jounal of Applied Statistics, 1–17.
Rizoiu, M. A., Mishra, S., Kong, Q., Carman, M. and Xie, L. (2018). SIR-Hawkes: Linking epidemic
models and Hawkes processes to model diﬀusions in ﬁnite populations, Proceedings of the 2018
World Wide Web Conference, 419–428.
谷口正信 (2005).『数理統計・時系列・金融工学』, 朝倉書店, 東京.
Yang, S. Y., Liu, A., Chen, J. and Hawkes, A. (2018). Applications of a multivariate Hawkes process tojoint modeling of sentiment and market return events, Quantitative Finance, 18(2), 295–310.

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

207

Comparison of Two Estimation Methods for Hawkes Processes andApplication to Actual Data Analysis
Shuji Chinone1 and Hiroshi Shiraishi2
1 Graduate School of Science and Technology, Keio University
2 Faculty of Science and Technology, Keio University

We consider multivariate Hawkes processes that are a class of multivariate point
processes with self and mutual excitation properties. In this paper, we assume that theobserved data follow a multivariate Hawkes process, and consider two statistical estimation methods of the kernel function, which represents the characteristics of the Hawkes
process. The first method is a parametric approach using the maximum likelihood method.
It estimates the parameters of the kernel function on the assumption that the intensity dependency follows the exponential kernel function. The second method is a non-parametricapproach that does not specify the kernel function. The kernel function is estimated nonparametrically by approximating the continuous-time stochastic process based on discreteobservation. After comparing the estimation results of the two methods using simulation,
we fit the price and volume fluctuations in the cryptocurrency market to the Hawkes process based on the parametric approach, and visualize the propagation structure throughthe Hawkes graph representation, which is a method for visualizing the characteristics ofkernel functions. Furthermore, the number of people in the Kanto and Kansai regionsin Japan who were newly infected COVID-19 is fitted to the Hawkes process based on
the non-parametric approach, and the propagation structure is seen by the Hawkes graphrepresentation.

Key words: Hawkes process, multivariate point process, INAR process, Hawkes graph, graphical modeling,
non-parametric estimation.