Davis, R., Fokianos, K., Holan, S.H., Joe, H., Livsey, J., Lund, R. ,Pipiras, V. Ravishanker, N. (2021) Count Time Series: A Methodological Review. Journal of the American Statistical Association, 116(535), 1533-1547.
仕事の都合でカウント時系列について考えることになり、とりあえずRのtscountパッケージの解説に目を通したら、知らない話題がたくさんあることに気が付き、あわてて手に取った。
1. イントロダクション
本論文はカウント時系列モデルについて概観する。
もともとカウント時系列はGLMで記述されることが多かった。分野の発展に伴い、相関のあるカウント時系列のモデリングのためのいろいろなアプローチが登場した。ガウシアン時系列ではARMAが主流だが、カウント時系列では整数ARMA, 離散ARMA, 一般化ARMAなどが登場し、統一理論無しで発展した。
2. 間引き演算子ベースのカウント時系列モデル
2.1 概観
初期の研究では、非負整数値をとり、いろんなACFと周辺分布(ポアソンとか)を持つ強定常系列を構築するというのが目標とされることが多かった。
現在の研究では、ひとつの方向として、整数自己回帰モデルINAR(p)、ないし一般化INAR(p)がある。\(p\)個の過去の系列値とイノベーション変量(ポアソン分布やNBD)に対して間引き演算子thinning operatorを適用するというアプローチである。
基本的なアイデアを示すために、INAR(1)について考える。$$ X_t = \alpha * X_{t-1} + \epsilon_t$$ \(*\)は自己回帰における乗算を模倣する演算子で、成功率\(\alpha\)のiidなベルヌーイ試行を\(\{B_i\}^\infty_{i=1}\)として\(\alpha * W = \sum_{i=1}^W B_i\)である。\(\{\epsilon_t\}\)はiidなカウント確率変数。周辺分布を(たとえば)ポアソン分布にしたければ、平均\((1-\alpha)\lambda\)のiidなポアソン確率変数にすればよい。
[ああ、なるほどね… たとえば\(X_1 = 3\)だったら、\(X_2\)は「成功率\(\alpha\)のベルヌーイ試行を3回行って成功した回数 + パラメータ\((1-\alpha)\lambda\)のポアソン確率変数」にするわけだ。確かに、AR(1)なACFを持つ非負整数系列になるね。間引き演算子ってそういう意味か]
いくつかの細かい仮定の下で、\(\{X_t\}\)は一意な定常分布を持つマルコフ過程となる。
[…]
平均や他のパラメータを共変量の関数として非定常へと拡張するとか、他の間引き演算子を用いるという拡張とか、いろいろある。
間引き演算子ベースのモデルの推定には、最尤法が用いられることが多いが、他の方法もある。ベイジアンアプローチもある。
2.2 一般化間引き演算子
以下では特記しない限り、間引き演算子は非負整数の集合に対して働き、パラメータ\(\alpha \in [0,1]\)を持ち、カウント確率変数\(X\)を別のカウント確率変数\(X_\alpha\)にマップし、\(X\)が有限の平均を持つならば\(E(X_\alpha) = \alpha E(X)\)である。ふつうは\(X_\alpha \leq X\)だがそうでなくてもよい。
\(X\)を非負確率変数とし、そのCDFを\(F_X\)とする。いま\(X\)が、すべての\(0 \lt \alpha \lt 1\)について、$$ X = \alpha * X + \epsilon(\alpha)$$ を満たすような、\(X\)とは独立なイノベーション確率変数\(\epsilon(\alpha)\)を持つとしよう[数式の等号の上に\(\mathcal{D}\)がついている。どういう意味だろう…?]。このとき、\(X\)は離散自己分解可能(DSD)であるという。定常カウントモデル\(X_t = \alpha * X_{t-1} + \epsilon_t(\alpha)\)はすべての\(0 \lt \alpha \lt 1\)について周辺分布\(F_X\)を持つ。なお、DSD確率変数は無限分解可能であり[任意の数の独立同一分布の確率変数の合計の確率分布として表現できるという意味らしい]、従って混合ポアソン分布であり、分散は少なくとも1であることが知られている。
\(\alpha \in (0,1)\)であり、\(\epsilon_t\)のCDF\(F_\epsilon\)を固定したとき、\(X_t = \alpha * X_{t-1} + \epsilon_t(\alpha)\)を再帰的に適用してもやはり正当なカウント時系列が定義される。つまり、ペア\((\alpha, F_\epsilon)\)で\(F_X\)の定常分布が決まるわけだ。
[このくだり、いまいち主旨を理解できなかった…]
間引き演算子のうち、非負整数確率変数\(K(\alpha), \alpha \in [0,1]\)の合成に基づく間引きを期待間引きという。その例として二値間引きがある。\(\{K_i(\alpha)\}^infty_{i=1}\)を\(K(\alpha)\)のiidな反復として、$$ K(\alpha) * l = \sum_{i=1}^l K_i(\alpha)$$ [数式のなかの\(*\)は、原文では丸の中のアスタリスク]。
ガウシアンAR(p)の整数版を作りたかったら、合成に基づく間引き演算子を用いればよい。以下を一般化整数自己回帰GINAR(p)という。$$ X_t = \sum_{j=1}^p K(\alpha_j) * X_{t-j} + \epsilon_t = \sum_{j=1}^p \sum_{i=1}^{X_{t-j}} K_{ijt}(\alpha_j) + \epsilon_t$$
もっとも広く用いられているINAR(p)は二値間引きである。
なお、これらのモデルは基本的には負の系列相関を持てない。
3. 新しい展開
以下では長期記憶という言葉が出てくるが、これは\( \sum_{h=0}^\infty | cov(X_t, X_{t+h}) |= \infty\)という意味である。[わ、わからん…]
3.1 二値系列からの構築
任意の離散分布は独立コイン投げから生成できることが知られている。同様に、相関のあるベルヌーイ過程を組み合わせて、目指す周辺分布を持つ系列を構築しようという試みがある。
\(\{B_t\}^\infty_{t=1}\)を二値確率変数の系列とし、そのラグ\(h\)の共分散を\(\gamma_B(h) = cov(B_t, B_{t+h})\)とする。ここからカウント系列を作りたいんだけど、その前に、定常で相関のあるベルヌーイ系列を作る方法をふたつ示そう。
- renewal-based \(\{B_t\}\)。ライフタイム\(L \in \{1, 2, \ldots, \}\)と初期ディレイ\(L_0 \in \{1, 2, \ldots, \}\)を定義しておく。\(L\)と同じ分布を持つiidな変数\(L_1, L_2, \ldots\)をつくり、ランダムウォーク\(S_n = L_0 + L_1 + \cdots + L_n\)をつくる。で、そのどこかで\(S_n = t\)になることがあったら\(B_t = 1\)とする。こうしてできる \(\{B_t\}\)が定常になる条件は…[めんどくさい話になってきたので読み飛ばした]
- Gaussian clipped \(\{B_t\}\)。[だめだ、もう私の体力はゼロだ… スキップ]
話をカウント系列の構築に戻すと…[パス。いつか元気のあるときに読みなおそう]
3.2 潜在ガウシアン変換
周辺分布\(F_X\)を持つ定常カウント系列は、probability integral transformと定常ガウシアン系列から作ることもできる。[ガウシアン・コピュラとかいうやつね。残念ながら、この節はたとえ元気なときでも理解できないだろうと思われる。パス]
4. 一般化状態空間モデル
[ああ、なんとか理解できそうな話になった。気を取り直して頑張りましょー]
一般化状態空間モデリングは時系列分析のもっとも一般的で柔軟なフレームワークのひとつである。
ごく簡単に説明すると、状態空間モデルというのは、状態過程\(\alpha_t\)のもとでの観察時系列\(X_t\)のモデルと、\(\{\alpha_t\}\)の時間発展のモデルからなる。
感圧時系列のモデル指定には既知の確率質量関数を使うことが多い。たとえばポアソンモデルでは$$ X_t | \alpha_t \sim Poisson(\exp(\alpha_t)) $$ となる。確率資料関数としては他にNBDや一般化ポアソン分布などが使われる。
一般化すると、1パラメータ指数関数族から条件付き確率質量関数を得るとき、$$ P(X_t = x_t | \alpha_t) = \exp \left( \phi(x_t) + \alpha_t x_t – A(\alpha_t) \right) $$ \(A(\alpha)\)は確率質量関数の合計が1になるように決める基準化定数である。\(\phi(x)\)はふつう既知の関数で、ポアソンならば\(-\log(x!)\)である[ポアソン分布の確率質量関数は\( \frac{\lambda^x}{x!} \exp(-\lambda) = \exp \left(-\log(x!) + x \log(\lambda) – \lambda \right) \)だからね]。しかし\(\phi(\cdot)\)をパラメータとみてもよい。
状態過程\(\{\alpha_t\}\)のほうについてはパラメータ駆動モデリングと観察駆動モデリングがある。
- パラメータ駆動モデリングについて。
\(\{\alpha_t\}\)がそれ自体の確率的メカニズム(ふつうはマルコフ過程)に従うと仮定する。ポアソン・モデルではガウシアンAR(1)を仮定することが多い。$$ \alpha_t = \phi_0 + \phi_1 \alpha_{t-1} + \epsilon_t, \ \{\epsilon_t\} \sim \ iid \ N(0, \sigma^2)$$ っていう感じね。
もっと一般的に、1パラメータ指数関数族についていうと、リンク関数を使って\(E(X_t|\alpha_t) = B(\alpha_t)\)とする。GLMと同じである。
\(p\)次元共変量\(\mathbf{c}_t\)を組み込む場合は $$ \alpha = \mathbf{c}^\top_t \mathbf{\beta} + Z_t $$ とする。\(\{Z_t\}\)はたとえば定常ガウシアン過程とする。これは一般化線形混合モデルだといえる。
ポアソンの場合だと、無条件平均は$$ E(X_t) = E(\exp(\mathbf{c}^\top_t \beta + Z_t)) = \exp(\mathbf{c}^\top_t \beta) E(\exp(Z_t)) = \exp(\mathbf{c}^\top_t \beta + \frac{var(Z_t)}{2}) $$つまり、共変量は伝統的なGLMと同じように解釈することができる。
パラメータ駆動モデリングは階層ベイズモデルと捉えることができる。しかし、完全にベイジアンな定式化の場合は、データ・モデルのパラメータとパラメータ・モデルのパラメータの両方に事前分布が設定されるという点に注意。
パラメータ駆動モデルの推定はなかなか難しい。MCMCとか重点サンプリングとか、尤度の近似とかペアワイズ合成尤度[?]とかが使われるが、あまり信頼できない。 - 観察駆動モデリングについて。
\(\alpha_t\)を過去の観察\(X_s\)の関数と仮定する。1パラメータ指数関数族における一般的なチョイスは、条件付き平均\(\lambda_t = E(X_t | \alpha_t) = B(\alpha_t)\)について$$ \lambda_t + d + \alpha_1 \lambda_{t-1} + b_1 X_{t-1}$$ というGARCH過程を仮定するというもの。これが有限の平均の強定常な解を持つ必要十分条件は\(a_1 + b_1 \lt 1\)である。[なるほど、ポアソン・パラメータについてGARCH(1,1)を仮定し、状態方程式側には撹乱項を入れない、というのが一番シンプルな定式化だってわけか]
観察駆動モデリングの難点は、状態過程に定常性のような安定性を持たせるのが難しいという点にある。上の定式化は特殊ケースで、\(a_1 + b_1 \lt 1\)でありさえすれば一意な強定常解が存在し、\(\{X_t\}\)は幾何レートの\(\beta-\)混合になる。[???]
4.1 観察駆動モデル
カウント系列に対するARモデルの単純な一般化として線形・対数線形ARモデルがある。
- 線形ARポアソンモデルは、過程の履歴を\(\mathcal{F}_t = \sigma(X_t, X_{t-1}, \ldots) \)とし[難しい記号を使うのやめてもらえませんかね？]、\(d\)と\(\{b_j\}^p_{j=1}\)を非負として、$$ X_t | \mathcal{F}_{t-1} \sim Poisson(\lambda_t)$$ $$ \lambda_t = d + \sum_{j=1}^p b_j X_{t-j}$$ というモデルである。\(X_t\)の周辺分布がポアソン分布だという意味ではない点に注意。ランダム要素がポアソン、システム要素が\(d + \sum_{j=1}^p b_j X_{t-j}\)、リンク関数自体はidentityであるようなGLMである。[なるほど、たしかに、リンク関数はidentityなのね。混乱しないようにしないと]
- 対数線形モデルは、\(v_t = \log (\lambda_t)\)に対して$$ v_t = \sum_{j=1}^p b_j \log(X_{t-j} +1) $$ というモデル。加法的に共変量を入れることもできる。
線形であれ対数線形であれ、短期記憶的な自己共分散しか持てないという点でAR(p)やARCH(p)に近い。記憶減衰をもっと遅くするためにGARCHモデルを使おうという提案もある。
- 線形のほうだと、たとえば\( X_t | \mathcal{F}_{t-1} \sim Poisson(\lambda_t) \)として $$ \lambda_t = d + \sum_{i=1}^p a_i \lambda_{t-i} + \sum_{j=1}^q b_j X_{t-j} $$ というモデルである。\(p\)を大きくすれば減衰を遅くできる。\(\lambda_t\)とは\(E(X_t|\mathcal{F}_{t-1})\)であるとともに\(Var(X_t|\mathcal{F}_{t-1})\)でもあるので、これを整数GARCH(INGARCH)モデルという。
- 対数線形モデルのほうでも同じように、$$ v_t = d + \sum_{i=1}^p a_i v_{t-1} + \sum_{j=1}^q b_j \log(X_{t-j}-1)$$ というモデルが広く使われている。
[閾値モデル、一般化線形ARMAモデル、一般化ARスコアモデルの紹介。後ろめたいけど時間がないので丸ごとパス]
4.2 マルコフモデルと隠れマルコフモデル
[半頁くらい。パス]
5. 多変量モデル
[ここからが最も関心あるパート。気を取り直して頑張ろう… でも間引き演算子とコピュラのパートは飛ばそう…]
5.1 多変量INARモデル
[パス]
5.2 コピュラ・ベースのモデル
[パス]
5.3 パラメータ駆動モデル
すでに述べたように、パラメータ駆動モデルは(過去の過程の値ではなくて)観察されていない過程によって駆動される。
Jorgensen et al.(1999 Biometrika)はこういうモデルを考えた。\(X_{i,t}\)の条件付き分布が平均\(a_{i,t} \lambda{t}\)のポアソンだとする。\(\lambda_0 = 1\)とし、\(\lambda_t\)はガンマ分布に従い、その平均を\(b_t \lambda_{t-1}\)、変動係数の平方を\(\sigma^2/\lambda_{t-1}\)とする。\(b_t = \exp ( \mathbf{c}_t — \mathbf{c}_{t-1})^\top \beta \) とする。Jorgensenらによれば、このとき\(\mathbf{X}_t\)の分散は「ポアソン分散」とoverdispersionに分解される。このモデルはカルマンフィルタで推定できる。
[自分の関心にジャストフィットなんだけど、この簡潔な紹介では到底理解できない。元論文にあたるか… 辛いなあ…]
[挙げられている研究のうち入手できそうなのは: Jorgensen, Lundbye-Christenten, Song, & Sun (1999 Biometrika), Jung, Liesenfeld, Richard (2011 J.Business & Enon.Stat.), Ravishanker, Serhiyenko, Willig(2014 Stat.&ItsInterface), Ravishanker, Venkatesan, Hu(2016 in “Handbook of Discrete-Valued Time Series”)]
5.4 観察駆動モデル
単変量の場合の線形モデルも対数線形モデルも、多変量への拡張がある。\( X_{it} | \mathcal{F}_{t-1} \sim Poisson(\lambda_{it})\)、\(\mathbf{v}_t = \log(\mathbf{\lambda}_t) \) (要素ごとに対数をとっている)、\(\{\mathbf{X}_t\}\)は\(d\)次元とする。
- 線形モデルの場合、 $$ \mathbf{\lambda_{t}} = \mathbf{r} + \mathbf{A} \mathbf{\lambda}_{t-1} + \mathbf{B} \mathbf{X}_{t-1} $$ という定式化[しれっとGARCHになっておる…]。\(\mathbf{r}\)は\(d\)次元(要素はすべて正)、\(\mathbf{A}, \mathbf{B}\)は\(d \times d\)の未知の行列(要素はすべて正)。
- 対数線形モデルの場合、$$ \mathbf{v}_t = \mathbf{r} + \mathbf{A} \mathbf{v}_{t-1} + \mathbf{B} \log (\mathbf{X}_{t-1} + \mathbf{1}_d) $$ という定式化。要素の符号制約はいらない。こっちのほうが共変量を入れやすい。
[挙げられている研究のうち入手できそうなのは: Heinen & Rengifo(2007 J.EmpiricalFinance), Lee, Lee, & Tjostheim (2018 TEST), Fokianos, Stove, Tjostheim, & Doukhan (2019 Bernoulli)]
6. ベイジアン動的GLM
まずはベイジアン階層モデリングについて。一般化状態空間モデルより柔軟にモデル化できる。
データ\(\mathbf{X} = \{\mathbf{X}_t\}\)は過程(ふつうは潜在ガウシアン過程)\(\mathbf{Z} = \{\mathbf{Z}_t\}\)とパラメータ\(\mathbf{\theta}_X\)に条件づけられ、\(\mathbf{Z}\)はパラメータ\(\mathbf{\theta}_Z\)に条件づけられていると考える。二種類のパラメータが条件付き独立だとすると、$$ [\mathbf{X}, \mathbf{Z} | \mathbf{\theta}_X, \mathbf{\theta}_Z] = [\mathbf{X} | \mathbf{Z}, \mathbf{\theta}_X] [\mathbf{Z} | \mathbf{\theta}_Z] $$ と書ける。条件付き尤度\([\mathbf{X} | \mathbf{Z}, \mathbf{\theta}_X]\)が「データモデル」、その後ろの\([\mathbf{Z} | \mathbf{\theta}_Z]\)が「プロセスモデル」である。
潜在ガウシアン過程とパラメータの事後分布は$$ [\mathbf{Z}, \mathbf{\theta}_Z, \mathbf{\theta}_X | \mathbf{X}] \propto [\mathbf{X}| \mathbf{Z}, \mathbf{\theta}_X][[\mathbf{Z} | \mathbf{\theta}_Z][\mathbf{\theta}_X, \mathbf{\theta}_Z]$$ と書ける。右辺はさらに下位要素へと分割できる。
6.1 単変量カウント系列のモデリング
データモデルにはふつう指数分布族をつかう。$$ [X_t | \alpha_t] = \exp(\phi_t(X_t) + \alpha_t X_t – A(\alpha_t)) $$ \(\mu_t = E(X_t|\alpha_t)\)とし、リンク関数を\(g(\cdot)\)として、次の形式が使われることが多い。共変量を\(\mathbf{c}_t\)として $$ g(\mu_t) = \mathbf{c}^\top_t \beta + Z_t$$ \(\{Z_t\}\)はiidな平均0のガウシアン誤差。モデルとしては以下がある。
- ポアソン。リンク関数は\(\log\)。この形式でもoverdispersionが表現できている。
- 負の二項分布。リンク関数はやはり\(\log\)。
- Conway-Maxwell-Poisson。underdispersionも表現できる。
- ほかにも、ゼロ過剰ポアソン(ZIP)とか、混合ポアソンとか。
基準化定数の計算が大変になることもあるが、漸近近似とか、OpenMPやCUDAをつかった並列計算が行われている。
[Ravishanker, Venkatesan, Hu(2016 前掲)の紹介。興味深いが、これは原文を当たった方が良さそうなので、メモ省略]
6.2 多変量モデリング
3つのアプローチを紹介しよう。
- level correlated models (LCM)。Aikchison & Ho(1989 Biometrika), Chib & Winklemann(2001 J.Business&Econ.Stat.), Ma, Kockelman, Damian(2008 Accident Analysis & Prevention), Serhiyenko, Ravishanker, Venkatesan(2017 Applied Stochastic Models in Business & Industory), Rue, Martino, Chopin (2009 JRSS)をみよ。
カウントベクトル \(\mathbf{X}_{it} = (X_{1it}, \ldots, X_{dit}\) [あれえええ？ 添字\(i\)ってなに？] の要素間の依存性を組み込むために、まず、観察方程式を$$ X_{jit}|\theta_{jit} \sim UDC_j (\theta_{jit})$$ とする。UDCは「カウントの単変量分布」の略。要素によって異なる分布を使っても良い。で、\( \mu_{jit} = E(X_{jit}) \)について$$ \log(\mu_{jit}) = \beta_{ji0} + \gamma_{jt} + \mathbf{c}^\top_{jit} \beta_{ji} + \alpha_{jit}$$ $$ \gamma_{jt} = \gamma_{j,t-1} + w_{jt} $$ $$ w_{jt} \sim N(0, 1/W)$$ $$ \alpha_{it} \sim MVN(\mathbf{0}, \mathbf{\Sigma}_i)$$ とする。ここではリンク関数はすべての要素について\(\log\)だが、変えても良い。また状態誤差の分散\(1/W\)を\(1/W_j\)としてもよい。
RueらはこれをR-INLAで推定している。[へえー] - 階層多変量動的モデル(HMDM)。
以下ではポアソンモデルについてのみ説明する。\(d\)変量のポアソン確率ベクトル\(\mathbf{X}\), \(q\)次元の未観察の独立ポアソン変量を\(\mathbf{Z}\)として\(\mathbf{X}=\mathbf{A}\mathbf{Z}\)と考える。行列\(\mathbf{A}\)の要素は0か1とする。\(q\)次元パラメータ行列を\(\mathbf{\mu}\)として、確率質量関数は…[式省略]。
[1pくらい説明が続く。いまちょっと時間がないのでパスするけれど、なるほどね、個人ごとにカウント時系列があって、強度の系列が未知の群ごとに異なると考えられるような場合には実にもっともなモデルだ]
推定についてはRavishanker, Venkatesan, Hu(2016 前掲)をみよ。標本サイズが大きすぎると現実的でない。LCMのほうがよいだろう。 - 潜在多変量ログガンマ潜在過程を持つ多変量ポアソンカウントモデル。Bradley, Holan, Wikle(2018 Bayesian Analysis)による。
[この節、良く理解できないのでほぼ逐語訳する]
彼らのモデルは一般的な多変量時空間過程のために開発された。時系列バージョンは場所がひとつであるような特殊ケースとなる。この枠組みの主な利点として、多変量ログガンマ分布がポアソン分布の共役分布なので、モデリングが柔軟で計算が速くなるという点がある。
互いに独立な\(m\)本のログガンマ確率変数を\(\mathbf{w} = (w_1, \ldots, w_m)^\top\)とし、\(i=1, \ldots, m\)について\(w_i \sim LG(\alpha_i, \kappa_i)\)とする。\(\alpha_i, \kappa_i\)は正のパラメータである。
次に\(\mathbf{q} = \mathbf{b} + \mathbf{Vw}\)を定義する。\(\mathbf{V}\)は\(m \times m\)行列である。\(\mathbf{q}\)を多変量ログガンマ(MLG)確率ベクトルと呼ぶ。確率密度関数を導出できる。
ポアソンMLGモデルは以下となる。データモデルをポアソンとし、$$ g(\mathbf{\mu}_t = \mathbf{c}^\top_t \beta + \phi^\top_t \eta_t + \xi_t$$ \(\eta_t\)はVAR構造を持ちMLG分布に従う。\(\{\phi^t\}\)はたとえば時間基底関数の決定論的な集まりである。\(\{\xi_t\}\)は独立ノイズで、これもMLG分布に従う。
ハイパーパラメータについて便利な事前分布を置き、完全に共役なモデルとすることができる。潜在過程の周辺分布があきらかに非ガウシアンである状況において性能を向上できる。
[なんだかよくわかんないんだけど、データモデルをポアソンとしたとき、パラメータの事前分布をうまいこと決めて共役にする、という話だろうか…]
7. 将来の方向
[いろいろ書いてある。パス]
——-
tscountの解説よりこっちを先に読めばよかったな。それにしても、いろんな話題があるものだ…