統計数理（2017）
第 65 巻 第 2 号 323–339©2017 統計数理研究所［研究詳解］

整数値自己回帰モデルの最近の発展
中嶋 雅彦1 ・酒折 文武2 ・川崎 能典3,4（受付 2017 年 3 月 31 日；改訂 7 月 4 日；採択 7 月 4 日）

要

旨は，自分自身の整数値自己回帰モデル（Integer-valued Autoregressive Models，INAR モデル）
過去（整数値）に直接依存させて現時点での値（整数値）を説明するものであり，潜在過程に基づく動的一般化線形モデルや一般化状態空間モデルとは，データへの接近法は大きく異なる．特に，自己回帰モデルの左辺と右辺において確率分布が整合的であるためには，周辺分布，自己回帰項，イノベーションのそれぞれの定式化が重要である．1980 年代後半から 90 年代初期に
いくつかの研究が行われたあとやや停滞していたこの分野は，2000 年代後半から活性化してきている．本稿は，最近の整数値自己回帰モデルの発展を，ポアソン分布に基づく INAR モデル
から出発して，定式化を変えた各種モデルを順に取り上げることで概観する．その過程で，必ずしも既存の文献で詳述ないし証明の与えられていない結果に関しては，著者らなりに補った
結果を付録に収めた．ポアソン分布の差から生成される分布に基づく INAR モデルについて
は，若干新しい提案と結果を加えることができた．その枠組みを利用した実データ解析例を，最後に紹介する．
キーワード：整数値時系列データ，間引き演算子，確率変数の分解，INAR(1) モデル，
INAR(p) モデル，モーメント法．

1.

はじめに

センサー技術など情報通信技術の発達に伴うデータ取得環境の変化の中で，様々な時系列
データが収集・活用されるようになり，時系列データに対する統計的モデリングの重要性が高まっている．時系列解析においては，自己回帰移動平均（ARMA）モデルに代表される，イノベーションに正規分布を仮定した線形モデルが，現在でも最も頻繁に用いられるモデルであり，その対象は理学・工学のみならず経済学等の社会科学にも広がって久しい．
しかしながら，このような枠組みが，観測される時系列が整数値である場合でも妥当かどうかは注意を要する．例えば，マクロ経済時系列データの多くは，その取りうる値が非負で整数値とはいっても，概ね連続量と見なして差し支えないであろう．従って，イノベーションならびにデータに対する連続分布の仮定は近似として許される範囲として，あとは分布形とモデルの定式化の問題となる．
一方，観測データが計数時系列（count time series）で，とりわけ稀なイベントの計測から生じ
1

中央大学大学院 理工学研究科：〒 112–8851 東京都文京区春日 1–13–27

2

中央大学 理工学部：〒 112–8851 東京都文京区春日 1–13–27

3

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

4

総合研究大学院大学 統計科学専攻：〒 190–8562 東京都立川市緑町 10–3

324

統計数理

第 65 巻 第 2 号 2017ているような場合，すなわち 0, 1, 2 からたかだか 10 程度の計数の連なりである場合はもはや正規近似は妥当性を欠き，データ分布に何らかの離散分布を仮定することが必要である．計数
時系列に対するモデリング法として多くの成功事例を集めながら拡張されてきたのは，一般化
線形モデル（generalized linear model, GLM; Nelder and Wedderburn, 1972）の時系列版である動的一般化線形モデル（Dynamic GLM，以下 DGLM と略）であった．このような問題設定では，一般化された状態空間表現に基づく更新アルゴリズムを実現するための新たな計算法の提案がによる数値積分フィルタが挙げられる．
必要とされた．例としては Kitagawa（1987）
一般化状態空間モデルないし階層構造モデルで少数計測データ（small count data）を分析す
ることの大きなメリットは，背後にある生起度（intensity）が仮に非定常に近い振る舞いをするような場合でも，生起度を表す潜在時系列モデルの側で柔軟性に富んだ適切なモデルを置いてやれば，少なくとも現象の記述には支障がない点であろう．DGLM に関しては Westを，一般状態空間モデルの推定法に関してはand Harrison（1997），Fahrmeir and Tutz（2001）
Kitagawa（2010）の第 14 章と第 15 章を参照されたい．
さて，本稿で取り上げる整数値自己回帰モデル（Integer-valued Autoregressive Models，INAR
モデル）
は，整数値時系列のモデリングのための別種のアプローチである．DGLM や一般化状
態空間モデルでは背後の潜在時系列に系列の履歴情報が縮約されるわけであるが，INAR モデ
ルは整数値の履歴自体にモデルを依存させる．そうすると，例えばイノベーションに連続分布を仮定した線形 AR(1) モデル Xt = αXt−1 + t からの類推で考えれば，自己回帰項 αXt−1 の部分が整数値となるよう定義を与え直し，かつ両辺の分布が整合する形でイノベーション t の分布を定式化しなければならない．INAR モデルにおける自己回帰項の定義は，2.1 節で説明
で与えられる．間引きの構造は，自己回帰項が 0 から Xt−1 の値
する通り「間引き（thinning）」
を取りうる形で定式化されており，それが 0 に近いのか Xt−1 に近いのか，換言すれば時系列が希薄なのか強いのかは，間引き項を生成する確率分布のパラメーとして持続性（persistency）
タ設定次第となっている．
以下では，INAR モデルとその拡張に関する文献をサーベイしておこう．著者らが文献上確認できた INAR に関する最も古い論文は McKenzie（1985）である．INAR(1) モデルを p 次に拡，Du and Li（1991）がある．Du and Li（1991）で張した研究として，Al-Osh and Alzaid（1990）
は INAR(p) モデルのパラメータ α1 , α2 , . . . , αp をユールウォーカー法により推定できること
では，INAR(p) モデルの周辺分布がポアソン分布のを示しており，Al-Osh and Alzaid（1990）
p
とき， i=1 αi ≤ 1 ならば 1 次のときと同様に確率変数が退化せずに分解可能であることが述べられている．
こうして 80 年代後半から 90 年代初頭に散発的に行われていた研究は，2000 年代後半に入ってようやく活性化する．Leonenko et al.（2007）は INAR(1) モデルの周辺分布をポアソン分布
ではなく負の二項分布とした NBD INAR(1) モデルを提案し，過分散なデータに対する INAR
モデルを提案している．Savani and Zhigljavsky（2007a, 2007b）では，NBD INAR(1) モデルのパラメータのモーメント法による推定について言及しており，3 種類の推定量を提案し，比較している．
Fokianos（2011）
は周辺分布をポアソン分布とした PINAR(1) モデルを中心に解説を加えているが，ハンドブックの一章ということもありサーベイ的な側面も強いことから，現時点でこの分野の研究を始める際にまず目を通しておくべき文献と言える．Rajarshi（2012）も INAR を含
む離散確率過程のモノグラフであり，間引き項を二項分布とし，周辺分布を一般的によく知られた離散型分布としたモデルについて広く述べている点に特徴がある．
間引き項を二項分布ではなく幾何分布や超幾何分布とした研究もある．Rinstić et al.（2009）は，間引き演算子を幾何分布，周辺分布を幾何分布とした NGINAR(1) モデルを導入し，パラ整数値自己回帰モデルの最近の発展

325

表 1．INAR モデルの概要．

メータ推定にも言及している．Weiß（2008）
は，間引き演算子を超幾何分布，周辺分布を二項分布とした BINAR(1) モデルを考え，確率変数の取りうる値が有限な場合の INAR モデルを提案した．
整数値全体をとる時系列データに対する INAR モデルとして Freeland（2010）は，ポアソン（以下スケラム分布と記す）に着目し，2 つの独立な分布の差が従う分布としての Skellam 分布
PINAR(1) モデルの差による SINAR(1) モデル（Skellam INAR(1) モデル）を提案した．なお，
Freeland（2010）
では，同じパラメータを持つポアソン分布の差を取った研究しか述べられていない．また，Baretto-Souza and Bourguignon（2013）は，2 つの独立な NGINAR(1) 過程の差を考えた非対称離散ラプラス分布による STINAR(1) モデルを提案し，さらに負の相関をもつ場
合への拡張を与えた．ここまで名前を挙げたモデル群を，その特徴と拡張性の現状という観点から整理したのが，表 1 である．
以上を踏まえ本稿では，まず PINAR モデルから始めて，NBD INAR モデル，BINAR モデ
ル，SINAR モデルの定義と性質について紹介する．文献上必ずしも詳述ないし証明が与えられていない結果に関しては，著者らなりに補った結果を付録に収めている．本稿は基本的には研究詳解であるが，SINAR(1) モデルの拡張とパラメータ推定法に関しては，若干新しい提案を加えることができたので 4.2 節に記し，5 節で実データへの応用を試みた結果を報告する．
2.

平均と分散が等しい場合：PINAR モデル

非負整数値全体をサポートとする基本的な分布であるポアソン分布を Xt の周辺分布とする
PINAR モデル
（Rajarshi, 2012）
は，INAR モデルの中で最も代表的なモデルである．
2.1 PINAR(1) モデル
PINAR(1) モデルは以下で表される：(2.1)

Xt = α ◦ Xt−1 + εt ,Xt ∼ P o(λ).

いま，Xt を平均が λ であるポアソン分布 P o(λ) に従う確率変数とする．また，α ∈ (0, 1], Xを非負整数値とし，{Yi } を独立にベルヌーイ分布 B(1, α) に従う確率変数とする．このとき，
X
α ◦ X を確率変数として α ◦ X = i=1 Yi と定義する．すなわち，α ◦ X のとりうる値は 0 から
X までの整数値であり，X を “間引く” という意味で ◦ を間引き演算子（thinning operator）と呼ぶ．このとき，Xt−1 ∼ P o(λ), α ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , α) より，α ◦ Xt−1 が平均 λα のポア
ソン分布に従うことが示せる．したがって，ポアソン分布の再生性より， {εt } ∼ P o(λ(1 − α))が成り立つこともわかり，PINAR(1) のパラメータは α と λ のみで表されることが確認できる．
PINAR(1) モデルに関する基本的な性質を見ておこう．条件付き期待値，条件付き分散はE(Xt |Xt−1 ) = E(α ◦ Xt−1 + εt |Xt−1 ) = αXt−1 + (1 − α)λ,

326

統計数理

第 65 巻 第 2 号 2017

V (Xt |Xt−1 ) = V (α ◦ Xt−1 + εt |Xt−1 ) = α(1 − α)Xt−1 + (1 − α)λとなる．また，E(Xt ) は条件付き期待値の性質を用いるとE(Xt ) = E(E(Xt |Xt−1 )) = E(αXt−1 + (1 − α)λ) = λ,V (Xt ) は条件付き分散の性質を用いるとV (Xt ) = V (E(Xt |Xt−1 )) + E(V (Xt |Xt−1 )) = λとなり，周辺分布であるポアソン分布の期待値と分散に一致していることが確認できる．
PINAR(1) モデルは，1 時点前の値に依存しているため 1 次のマルコフ性をもち，状態 i から jに推移する確率は i zmin(i,j)pij = P (Xt+1 = j|Xt = i) =zz=0α (1 − α)i−z

e−λ λj−z(j − z)!

となる．推移確率が時間 t に依存しないので，この時系列過程は斉時性が成り立つことが分かる．また，このマルコフ系列は任意の i, j に対して pij > 0 となることから，既約で非周期的であることが分かる．
2.2 PINAR(p) モデル
PINAR(1) モデルは p 次へ自然に拡張できる（Al-Osh and Alzaid, 1990 参照）．PINAR(p) モ
デルは
Xt = α1 ◦ Xt−1 + α2 ◦ Xt−2 + · · · + αp ◦ Xt−p + εtpで定義される．ここで，Xt ∼ P o(λ), i=1 αi ≤ 1, αi ◦ Xt−i ∼ B(Xt−1 , αi ) である．このとき，
pα ◦ Xt−j ∼ P o(αj λ) となり，εt ∼ P o((1 − i=1 αi )λ) であることが分かる．
PINAR(p) モデルにおいて，Xt の条件付き確率分布は
P (Xt = k|Xt−1 = x1 , . . . , Xt−p = xp ) =p min(xi ,k)  


xi
z
i=1

zi =0

zi

αi (1 − αi )条件付き期待値，条件付き分散は(2.2)

E(Xt |Xt−1 , Xt−2 , . . . , Xt−p ) =pk−li −zi λ
αi Xt−i +

1−

i=1

(2.3)

V (Xt |Xt−1 , Xt−2 , . . . , Xt−p ) =

p

i=1

p


(j −

zi −λ

e
,
z
)!
i=1 i

i=1

p


αi λ,

i=1


αi (1 − αi )Xt−i +

p

1−

p



αi λ

i=1

で与えられる．この証明は Al-Osh and Alzaid（1990）では示されていないので，付録 A.2 に記しておく．
2.3 PINAR モデルにおけるパラメータ推定
PINAR モデルにおけるパラメータ推定はモーメント法を用いるのが最も単純である．ここでは PINAR(1) モデルにおけるパラメータ α, λ のモーメント法での推定量について述べておく．まず，1 次の自己共分散 Cov(Xt , Xt−1 ) は(2.4)

Cov(Xt , Xt−1 ) = E(Xt−1 (α ◦ Xt−1 + εt )) − E(Xt )E(Xt−1 ) = αV (Xt−1 )整数値自己回帰モデルの最近の発展327である．定常性より V (Xt ) = V (Xt−1 ) であるので，1 次の自己相関 ρ(1) はCov(Xt , Xt−1 )ρ(1) = =αV (Xt )V (Xt−1 )となることが分かる．よって，α のモーメント法による推定量は，1 次の標本自己相関 ρ̂(1) を用いて α̂ = ρ̂(1) とすればよい．ただし，α は間引きの確率を表すため非負であるが，標本自己相関は負の値をとることもありうるため注意が必要である．さらに，Xt ∼ P o(λ) なので，E(Xt ) = λ であることから，λ のモーメント法による推定量は標本平均 X̄ を用いて λ̂ = X̄ と1であることから，不偏性と有効性をもつこすればよい．なお，λ̂ は，E(λ̂) = λ, V (λ̂) = nI(λ)とがわかる．
もちろん最尤推定量などを考えることも可能であるが，推定量は陽関数としては得られず数値的に解くことになるため，最尤推定量のもつ漸近的な性質以上のことを言うことが困難となる．
3.

その他の非負整数値自己回帰モデル

本節では，まず周辺分布を負の二項分布とした NBD INAR(1) モデルについて述べる．次に周辺分布を二項分布とした BINAR(1) モデルについて述べる．
3.1 過分散の場合：NBD INAR(1) モデル
これまで，INAR(1) の周辺分布をポアソン分布とした PINAR モデルについて述べた．ポア
ソン分布には平均と分散が等しいという特徴があるので，データの平均と分散がほぼ等しけれ
ば PINAR(1) は妥当なモデルであると考えられる．しかしながら，平均よりも分散が大きい，
すなわち過分散の場合には PINAR(1) モデルは適切でない．そこで，過分散なモデルとして，
INAR(1) モデルの周辺分布を負の二項分布（NBD）とした NBD INAR(1) モデルについて説明する．なお，これらにおけるパラメータ推定も PINAR(1) モデルと同様にモーメント法などで行うことができる．詳しくは Leonenko et al. （2007），Savani and Zhigljavsky（2007a, 2007b）を参照されたい．
3.1.1 NBD INAR(1) モデル
NBD INAR(1) モデルの式を次のように定義する：
Xt = α ◦ Xt−1 + εt ,α ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , α), Xt ∼ N BD(γ,β) .

ここで N BD(γ,β) はパラメータ γ, β の負の二項分布であり，その確率関数は



f (x) =

γ+x−1
x



β
β+1

γ 

1
β+1

x

(x = 0, 1, 2, . . . )βである．このとき，α ◦ Xt−1 ∼ N BD(γ, α ) である．証明は Leonenko et al.（2007）を参照されたい．
この NBD INAR(1) モデルにおいては，Xt と α ◦ Xt−1 の従う負の二項分布の第 2 パラメー
タが異なるため，εt の分布を再生性から求めることはできない．そこで，εt の確率母関数を求めると，

−1 −γGXt (s)βαsγGεt (s) =1 − (1 − α)1−=αGα◦Xt−1 (s)
α+β
α+β
β
β
となる．これは，負の二項幾何分布 NBG(γ, α+β, α) の確率母関数なので，εt ∼ NBG(γ, α+β
, α)

統計数理

328

第 65 巻 第 2 号 2017であることが分かる．なお，パラメータ m, p, θ の負の二項幾何分布 NBG(m, p, θ) は，以下の確率関数をもつ確率分布である：
f (x) =

∞ 
k+x−1 kk=0p (1 − p)xx

m+k−1 mθ (1 − θ)kk(x = 0, 1, 2, . . . ).

この確率母関数はG(s) = θm (1 − (1 − θ)p(1 − (1 − p)s)−1 )−m



=

θ
p
1 − (1 − θ) 1−(1−p)smであり，負の二項分布の確率母関数の中に幾何分布の確率母関数が含まれている．これが負の二項幾何分布とよばれる理由である．
3.1.2 Stochastic thinning による NBD INAR(1) モデル
先の NBD INAR(1) モデルでは，Xt の周辺分布および間引き項 α ◦ Xt−1 は負の二項分布であったが，誤差項 εt は NBG 分布となった．誤差分布も同じ負の二項分布とするために，
Leonenko et al.（2007）
では，NBD INAR(1) モデルにおけるパラメータ α を定数ではなくベーとするモデルを提案している．
タ分布に従って変動する Stochastic thinning（確率的間引き）
いま，
Xt = At ◦ Xt−1 + εt , Xt ∼ N BD(γ,β)とし，間引きの確率を表すパラメータ At が定数ではなく，独立にベータ分布に従うとする．
すなわち，At ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , At ), At ∼ Beta(a, γ − a), Xt ∼ N BD(γ,β) である．このとき，At ◦ Xt−1 ∼ N BD(a,β) が成り立つことが示される．このことより，Xt と At ◦ Xt−1 の
第 2 パラメータがともに β のため，εt の分布が負の二項分布の再生性より εt ∼ N BD(γ−a,β)となることがわかる．このモデルを，stochastic thinning による NBD INAR(1) モデルという．
3.1.3 Mixed thinning による NBD INAR(1) モデル
3.1.1 節の間引きパラメータが定数 α のモデルと，3.1.2 節の確率変数 At のモデルを混合し，αAt を間引きパラメータとする一般的なモデル
Xt = αAt ◦ Xt−1 + εt

も Leonenko et al. （2007）
で 提 案 さ れ て い る ．こ こ で ，αAt ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , αAt ),βAt ∼ Beta(a, γ − a), Xt ∼ N BD(γ,β) で あ る ．こ の と き ，αAt ◦ Xt−1 ∼ N BD(a, α ) , εt ∼ββ, α) となることが示される．このモデルを，mixed thinning（混合NBD (γ, α ) · NBG(γ, α+β間引き）による NBD INAR(1) モデルという．
3.2 サポートが有限の場合：BINAR(1) モデル
ここまで，Xt のとりうる値が非負の整数全体のモデルについて考えてきた．しかし実際には，ある有限の値 n までしか値をとらない場合も考えられる．ここでは，周辺分布を最も代表的な二項分布としたモデルについて説明する．
単純に INAR(1) モデルの周辺分布を二項分布とした場合，誤差分布が一般的に知られた分布にならない．その改善策として 2 つの方法が提案されている．誤差分布を用いず，2 つの二項演算子を用いる方法と，超幾何演算子を用いる方法である．詳細は，Weiß（2008, 2009），
Rajarshi（2012）を参照されたい．原論文では両方とも BINAR(1) モデルとして言及されているが，ここでは区別のためにそれぞれ BINARb (1) と BINARh (1) と呼ぶこととする．

整数値自己回帰モデルの最近の発展

329

3.2.1 BINARb (1)：2 つの二項演算子を用いたモデルp, − 1−p), 1],
まず，2 つの二項演算子を用いたモデルを説明する．Xt ∼ B(n, p), ρ ∈ [max(− 1−ppα ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , α) とする．また，β = p(1 − ρ), α = β + ρ とする．このとき，
Xt = α ◦ Xt−1 + β ◦ (n − Xt−1 )(t ≥ 1)を BINARb (1) モデルという．
右辺により与えられる分布が Xt の分布，すなわち二項分布 B(n, p) であることを確認しておこう．右辺の確率母関数はE sα◦Xt−1 +β◦(n−Xt−1 ) = E((1 − α + αs)Xt−1 )E((1 − β + βs)n−Xt−1 )n= (1 − β + βs)



1−p+p

1 − α + αs
1 − β + βsnとなる．ここで，β = p(1 − ρ), α = β + ρ を代入すると，
E(sα◦Xt−1 +β◦(n−Xt−1 ) ) = (1 − p + ps)nとなることが示される．
BINARb (1) モデルにおける Xt の 1 次の自己相関 ρ(1) は ρ である．証明は付録 A.3 に記した．
3.2.2 BINARh (1)：超幾何演算子を用いたモデル
次に，間引き演算子に超幾何演算子を用いたモデルについて説明する．いま，Xt ∼ B(N, p)とし，n/N  Xt−1 |Xt−1 ∼ HG(N, Xt−1 , n) とする．ここで HG(N, m, n) は超幾何分布であり，その確率関数は
f (x) =

m
x

N −m
n−x
N
n

(x = 0, 1, 2, . . . , min(m, n))である．このとき，次のモデルを BINARh (1) モデルという：
Xt = n/N  Xt−1 + εt(t ≥ 1).

このモデルにおいて，n/N  Xt−1 ∼ B(n, p) となることが示される（付録 A.4 参照）．したがって，二項分布の再生性より {εt } ∼ B(N − n, p) となることが分かる．
また，このモデルは PINAR(p) モデルと同様，二項分布の再生性より p 次への拡張が可能で（2008, 2009）
では提案されていない．
あるが，Weiß4.

Z 上の整数値自己回帰モデル

これまでは，Xt が非負整数値をとる場合のモデルについて言及してきた．本節では，負の整数値も含む整数全体 Z に拡張したモデルとして，周辺分布をスケラム分布とするモデルについて述べる．
4.1 SINAR(1) モデル
ここでは，Freeland（2010）で提案されている 1 次の SINAR モデル（Skellam INAR モデル）について述べる．
まずはスケラム分布について確認しておく．X, Y が独立でそれぞれパラメータ λ1 , λ2 のポアソン分布に従うとする．このとき，Z = X − Y の従う分布を，パラメータ λ1 , λ2 のスケラ
ム分布という．ここでは，パラメータが λ1 , λ2 のスケラム分布を S(λ1 , λ2 ) と表記する．スケ

統計数理

330

第 65 巻 第 2 号 2017

ラム分布 S(λ1 , λ2 ) の確率関数は
f (z) = e−(λ1 +λ2 ) λz1∞

x=max(0,−z)(λ1 λ2 )x(x + z)!x!

(x ∈ Z)である．スケラム分布の性質については，例えば Alzaid and Omair（2010）を参照されたい．
スケラム分布を用いて，SINAR(1) モデルは次のように定義される．Xt , Yt をそれぞれ独立でパラメータの等しい PINAR(1) モデルとする．つまり，Xt = α ◦ Xt−1 + δt ,Xt ∼ P o(λ),α ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , α),Yt = α ◦ Yt−1 + ηt ,Yt ∼ P o(λ),α ◦ Yt−1 |Yt−1 ∼ B(Yt−1 , α)と す る ．こ の と き ，2 つ の PINAR(1) モ デ ル の 差 を 取 り ， Zt = Xt − Yt , αα ◦ Xt−1 − α ◦ Yt−1 , εt = δt − ηt とおくとZt = α Zt−1 + εt

Zt−1 ={Zt } ∼ S(λ, λ)である．このモデルを SINAR(1) モデルという．周辺分布がスケラム分布であり，Zt のとりうる値は整数全体となっていることがわかる．
Zt が正負の値をとりうるとなると，その自己相関が正の場合のみならず負の場合についても積極的に検討が必要となる．Freeland（2010）では，奇数と偶数で Xt と Yt を入れ替える方法が提案されている．すなわち，先ほどと同様に
Xt = α ◦ Xt−1 + δt , Yt = α ◦ Yt−1 + ηt ,Xt ∼ P o(λ1 ), Yt ∼ P o(λ2 )

とし，
Zt =

εt =

Xt − Yt(t = 0, 2, 4, . . . )Yt − Xt(t = 1, 3, 5, . . . )δt − ηt(t = 0, 2, 4, . . . )ηt − δt(t = 1, 3, 5, . . . )のように偶数と奇数で Xt と Yt の差の順序を入れ替えると，ρ(k) = (−α)k が成り立ち，1 次の自己相関が負の場合を表現することができる．
4.2 SINAR(1) モデルの拡張
では Xt , Yt の 2 つの PINAR(1) モデルのパラメータが等しい場合のみ
さて，Freeland（2010）を論じている．ここでは， λ が異なるモデルへの拡張を提案する．Xt , Yt をそれぞれ独立な
PINAR(1) モデルとし，そのポアソン分布の平均は異なるとする，すなわちXt = α ◦ Xt−1 + δt

{Xt } ∼ P o(λ1 ), α ◦ Xt−1 |Xt−1 ∼ B(Xt−1 , α),Yt = α ◦ Yt−1 + ηt

{Yt } ∼ P o(λ2 ), α ◦ Yt−1 |Yt−1 ∼ B(Yt−1 , α)とする．このとき，この 2 つの PINAR(1) モデルの差を取ったモデル
Zt = α Zt−1 + εt

{Zt } ∼ S(λ1 , λ2 )を考える．ここで，Zt = Xt − Yt , α Zt−1 = α ◦ Xt−1 − α ◦ Yt−1 , εt = δt − ηt である．λ1 = λ2のときは 4.1 節で述べた SINAR(1) モデルであり，このモデルはそれを拡張したものになっている．以降では，この拡張したモデルのことを改めて SINAR(1) モデルと呼ぶことにする．

整数値自己回帰モデルの最近の発展

331

この（拡張した）SINAR モデルの性質を述べる．εt のモーメント母関数，期待値，分散はそれぞれMεt (s) = exp −((1 − α)λ1 + (1 − α)λ2 ) + (1 − α)λ1 es − (1 − α)λ2 e−sE(εt ) = (1 − α)(λ1 − λ2 )(4.1)

V (εt ) = (1 − α)(λ1 + λ2 )であり
（付録 A.5 を参照）
，このモーメント母関数より εt ∼ S((1 − α)λ1 , (1 − α)λ2 ) であることがわかる．そして，
MαZt−1 (s) = exp −(αλ1 + αλ2 ) + αλ1 es + αλ2 e−sであり，α Zt−1 ∼ S(αλ1 , αλ2 ) であることもわかる．条件付き期待値，条件付きモーメント
母関数はそれぞれE(Zt |Zt−1 = z) = αz + (1 − α)(λ1 − λ2 )
(4.2)

 (λ1 (1 − α + αes ))z+y (λ2 (1 − α + αe−s ))yMεt (s)P (Zt−1 = z)(z + y)!
y!
∞E(exp {sZt }|Zt−1 = z) =

y=0

である
（付録 A.6 を参照）．
パラメータ推定は，PINAR(1) モデルと同様にモーメント法などで行うことができ，推定量⎧⎪⎪⎨α̂ = ρ̂(1)λ̂ = 1 (S 2 + z̄)

1
Z
2
⎪
⎪
⎩λ̂2 = 1 (S 2 − z̄)2Zが得られる．ここで，SZ2 は Z の不偏標本分散である．また，λ̂1 , λ̂2 はそれぞれ，パラメータ
λ1 , λ2 の不偏推定量である．
5.

実例

整数値自己回帰モデルの実データへの適用として，サッカーにおける各試合の得失点差のデータを，拡張した SINAR(1) モデルにあてはめた例を紹介する．データとしては，2015 年
サッカー J1 リーグ戦における，鹿島アントラーズの各試合の得失点差を用いる．データは J
リーグウェブサイトより取得した．図 1 は全 34 試合（節）における得失点差の分布である．標
本平均は 0.47，標本分散は 1.5 であった．また，図 2 はその推移を表す折れ線グラフである．
正負の値をとりながら，比較的上昇と下降を交互に繰り返している傾向があり，負の相関がありそうに見える．この時系列データのコレログラムは図 3 のようになり，1 次の自己相関係数
は −0.16 と，想定されたように負の値である．すなわち，−0.16 と値こそ大きくはないが相関があり，各試合での得失点差は独立ではないといえる．
そこで，SINAR(1) モデル
Zt = α Zt−1 + εtによりモデル化することを考える．1 次の自己相関が負であることから，
Zt =Xt − Yt(t = 0, 2, 4, . . . )Yt − Xt(t = 1, 3, 5, . . . )

332

統計数理

第 65 巻 第 2 号 2017

図 1．2015 年鹿島アントラーズの得失点差のヒストグラム．

図 2．2015 年鹿島アントラーズの得失点差のグラフ．

εt =δt − ηt(t = 0, 2, 4, . . . )ηt − δt(t = 1, 3, 5, . . . )のように偶数と奇数で得点数 Xt と失点数 Yt の差の順序を入れ替える，すなわち得失点差 Ztの正負を交互に入れ替える．パラメータをモーメント法で推定すると

⎧
⎪
⎪
⎨α̂
λ̂

= −ρ̂(1) = 0.16
= 1.331⎪⎪⎩λ̂2 = 0.86が得られる．しかし，この結果自体の解釈は難しい．

整数値自己回帰モデルの最近の発展

333

図 3．2015 年鹿島アントラーズの得失点差のコレログラム．

1 次の自己相関が小さく，このモデルを用いての予測の精度はあまり高くはないが，1 期先の得失点差は 1 であり，これ以予測などを行うことは可能である．例えば，第 34 節（最終節）
上の試合はなかったが，もし次の試合があったとした場合の得失点差を予測してみよう．例えば Z35 = 0 となる確率は，Z34 = 1 である条件の下で，P (Z35 = 0|Z34 = 1) = · · · + P (α̂ Z34 = −3|Z34 = 1)P (ε35 = 3)
+ P (α̂ Z34 = −2|Z34 = 1)P (ε35 = 2)+ ···+ P (α̂ Z34 = 2|Z34 = 1)P (ε35 = −2)+ P (α̂ Z34 = 3|Z34 = 1)P (ε35 = −3) + · · ·のように計算することができる．ここで，α̂ = −0.16 < 0 であり，t = 35 が奇数であることから，
P (α Z34 = v|Z34 = z34 ) =P (α ◦ X34 − α ◦ Y34 = v, Z34 = z34 )P (Z34 = z34 )= P (Z341=z34 )xP (α ◦ X34 = w|X34 = x)P (Y34 = x − z34 )w· P (X34 = x)P (α ◦ Y34 = w + v|Y34 = x − z34 )により求まる．この結果をまとめたものが表 2 である．表 2 には，得失点差が独立なスケラム
分布 S(1.33, 0.86) に従うとした場合の確率も併記している．いずれのモデルを用いた場合においても得失点差が 0 の確率が最も大きく，続いて得失点差 1，得失点差 −1 の確率が大きいことがわかる．しかし，その確率自体を比較してみると，相関を考慮した SINAR(1) モデルのほうが，全体的に得失点差が負の方向に分布が寄っていることがわかる．実際，得失点差の推移をみると，得失点差が 1 であった場合の次の節においては，ほとんどのケースで得失点差が負（とくに −1）になっており，SINAR(1) モデルではこの結果を反映していると解釈することがで
きる．

統計数理

334

第 65 巻 第 2 号 2017

表 2．SINAR モデルによる推定．

ところで，同じデータに対して DGLM で分析することは可能だろうか．少なくとも，得失点差系列自体に DGLM をあてはめるために，ポアソン過程の生起度（intensity）として単一の系列をどうモデル化して与えるかは自明ではないように思われる．この点に鑑みれば，正負双方を取り得る整数値時系列データに対しては，DGLM のようなパラメータ駆動型のモデルよりは，INAR モデルのような観測値駆動型のモデルのほうが構築しやすさからメリットがあると言えよう．
ここでは，得点系列と失点系列にそれぞれ別々の非定常ポアソンモデルをあてはめた結果とその含意を述べる．まずポアソンパラメータは時変（λt ）とし，非負制約を課して log λt に対して log λt = log λt−1 + vt というランダムウォーク型のシステム方程式を仮定する．シス
テムノイズには正規分布を仮定し，vt ∼ N(0, τ 2 ) とする．観測時系列 Yt は，ポアソン分布
P (Yt |λt ) = e−λt (λt )Yt /Yt ! に従うものとする．ここでは状態変数の分布を短冊状に近似してKitagawa（1987）の数値積分フィルタ・平滑化で，時変生起度 λt|T (t = 1, . . . , T ) を求める．状態の初期分布は，各系列の標本平均を λ̄ と書くとすると，N(λ̄, λ̄2 ) で与えた．
は，得点系列に関しては 1.12 × 10−4 と，失点系列に関しては
システムノイズの分散（τ 2 ）4.01 × 10−3 と推定された．このとき最終の第 34 節でのフィルタ値はそれぞれ 1.657 と 1.067と得られたので，ランダムウォーク予測で仮想的第 35 節もこの生起度で得失点が発生すると仮定する．10 万回のシミュレーションから得られた得失点差の相対頻度を，表 2 の最終行に示している．傾向としては，独立スケラム分布の結果に近いが，今回の DGLM では得失点差系
列に見られる負の相関（高周波成分の強い平均回帰性）
がモデルに取り込めていないことを考えれば自然である．2 本の対数生起度系列（潜在時系列）に関する多変量モデリングが現象説明能力の改善に資するかどうかは，今後の研究課題である．
独立であっても，得点と失点に DGLM をあてはめることの効用はある．時変生起度の平滑化分布を描画してみれば，これらの整数値時系列が定常という仮定が妥当であったかどうかのチェックを，少なくとも視覚的に行うことが可能である．図 4 が今回の推定結果である．ポア
ソンの生起度に非定常な振る舞いを許す枠組みで分析したわけだが，平滑化分布の中央値に多少の凹凸は見られるものの，得点，失点ともに定常ポアソン過程とみなして問題ないであろう．
6.

おわりに

本研究では，整数値をとる時系列データに対するモデルとして近年研究が進んでいる PINAR
モデル，NBD INAR モデル，BINAR モデル，SINAR モデルについて，その定義や性質，そしてパラメータ推定法を説明した．
さらに本論文では，SINAR(1) モデルに関する拡張を提案し，それに伴いいくつかの性質を示すことができた．これらの結果により，今までよりもさらに複雑な構造をもつ整数値時系列デー
タに対するモデリングの可能性が広がったといえよう．しかしながら，拡張したモデルにおけるパラメータ構造の問題が残されており，これは今後の課題である．また，他のモデルに関する p 次への拡張への可能性，そしてそのパラメータ推定についての研究も残されている．最後に，パラメータ推定法が確立されていないモデルでの推定法の提案も今後の研究課題といえる．

整数値自己回帰モデルの最近の発展

335

図 4．得点（左パネル中 ◦）失点
（右パネル中 ×）の時系列と，DGLM で推定した時変生起度の平滑化系列．平滑化系列は太実線が 50% 点，細破線は下が 3% 点，上が 97% 点．横軸
は試合の開催節．

謝

辞

改稿につながるコメントを寄せてくださった匿名の査読者にこの場を借りて感謝申し上げた
い．本研究を進めるにあたり，中嶋は統計数理研究所の特別共同利用研究員制度（平成 28 年度）を利用した．
付録. 証明A.1 PINAR(p) モデルの誤差分布の導出
まず，αi ◦ Xt−i はMαi ◦Xt−i (s) = exp {λαi (s − 1)}pである．ここで，i=1

αi ◦ Xt−i の確率母関数は，ポアソン分布の再生性より
Gp

(A.1)

i=1

αi ◦Xt−i (s) = exp

Gp

GXt (s)

i=1

となる．



αi λ(s − 1)i=1となる．したがって，εt の確率母関数は，
Gεt (s) =p

αi ◦Xt−i (s)


= exp

1−

p

i=1


αi

λ(s − 1)

統計数理

336

第 65 巻 第 2 号 2017

A.2 （2.2），（2.3）の証明
まず，条件付き期待値は，E(Xt |Xt−1 , Xt−2 , . . . , Xt−p ) = E(α1 ◦ Xt−1 + · · · + αp ◦ Xt−p + εt |Xt−1 , Xt−2 , . . . , Xt−p )= E(α1 ◦ Xt−1 |Xt−1 ) + · · · + E(αp ◦ Xt−p |Xt−p ) + E(εt )
=

p




αi Xt−i +

1−

p


i=1

αi

λ

i=1

となる．続いて条件付き分散はV (Xt |Xt−1 , Xt−2 , . . . , Xt−p ) = V (α1 ◦ Xt−1 + · · · + αp ◦ Xt−p + εt |Xt−1 , Xt−2 , . . . , Xt−p )= V (α1 ◦ Xt−1 |Xt−1 ) + · · · + V (αp ◦ Xt−p |Xt−p ) + V (εt )
=

p




αi (1 − αi )Xt−i +

i=1

1−

p


αi

λ

i=1

となる．
A.3 BINARb (1) モデルにおける Xt の 1 次の自己相関 ρ(1) の導出
1 次の自己共分散はCov(Xt , Xt−1 ) = Cov(α ◦ Xt−1 + β ◦ (n − Xt−1 )) = ρV (Xt−1 )となるので，1 次の自己相関は ρ(1) = ρ となることが分かる．
A.4 n/N  Xt−1 ∼ B(n, p) の証明
X と n/N  X の同時分布を求め，そのあと n/N  X の周辺分布を求める．まず，X と
n/N  X の同時分布は
P (n/N  X = k, X = x) = P (n/N  X = k|X = x)P (X = x)
=

x
k

 

N −x
n−k
N
n

N x
p (1 − p)N −x
x

 
=





n k
N − n x−kp (1 − p)n−kp(1 − p)N −n−(x−k)
kx−kである．ここで，x の範囲が k ≤ x ≤ N − n + k であることに注意を払うと，n/N  X の確率関数は

 n k

N −n+kP (n/N  X = k) =x=kkp (1 − p)

 
=


n k
p (1 − p)n−k
k


 
=

n−k

N −n

x =0





N − n x−kp
(1 − p)N −n−(x−k)
x−k




N − n xp (1 − p)N −n−xxn kp (1 − p)n−kkとなるので，n/N  X ∼ B(n, p) となることが分かる．

(x = x − k とおく)

整数値自己回帰モデルの最近の発展

337

A.5 （4.1）
の証明
まず，
(A.2)

MαZt−1 (s) = E(es(α◦Xt−1 −α◦Yt−1 ) )= exp −(αλ1 + αλ2 ) + αλ1 es + αλ2 e−sとなる．ここで，（A.2）の結果を用いると，Mεt (s) はMεt (s) =MZt (s)MαZt−1 (s)= exp −((1 − α)λ1 + (1 − α)λ2 ) + (1 − α)λ1 es − (1 − α)λ2 e−sと な る ．こ れ は S((1 − α)λ1 , (1 − α)λ2 ) の モ ー メ ン ト 母 関 数 で あ る こ と か ら ，εt ∼S((1 − α)λ1 , (1 − α)λ2 ) である．このことから，期待値，分散はそれぞれ，E(εt ) = (1 − α)(λ1 − λ2 )V (εt ) = (1 − α)(λ1 + λ2 )であることが分かる．
の証明
A.6 （4.2）
まず，条件付き期待値はE(Zt |Zt−1 = z) = E(E(α Zt−1 |Zt−1 = z)) + E(εt )= αz + (1 − α)(λ1 − λ2 )である．一方，条件付きモーメント母関数はE(exp {sZt }|Zt−1 = z) = E(exp {s(α Zt−1 + εt )}|Zt−1 = z)Mεt (s)P (Xt−1 = z + y, Yt−1 = y)Mα◦Xt−1 (s)Mα◦Yt−1 (−s)
P (Zt−1 = z)
∞

=

y=0

=

Mεt (s)P (Zt−1 = z)∞
(λ1 (1 − α + αes ))z+y (λ2 (1 − α + αe−s ))y(z + y)!

y=0y!

となる．これ以上は解析的に解くことができない．

参

考

文

献

Al-Osh, M. A. and Alzaid, A. A. (1990). An integer-valued pth-order autoregressive structure(INAR(p)) process, Journal of Applied Probability, 27(2), 314–324.
Alzaid, A. A. and Omair, M. A. (2010). On the poisson diﬀerence distribution inference and applications, Bulletin of the Malaysian Mathematical Sciences Society, 33(1), 17–45.
Baretto-Souza, W. and Bourguignon, M. (2013). A skew true INAR(1) process with application, arXiv:
1306.0156 [stat.ME].
Du, J. G. and Li, Y. (1991). The integer-valued autoregressive (INAR(p)) model, Journal of Time
Series Analysis, 12, 129–142.
Fahrmeir, L. and Tutz, G. (2001). Multivariate Statistical Modelling Based on Generalized Linear

338

統計数理

第 65 巻 第 2 号 2017

Models, 2nd ed., Springer-Verlag, New York.
Fokianos, K. (2011). Count time series models, Handbook of Statistics Vol.30 (eds. T. S. Rao, S. S.
Rao and C. R. Rao), 315–347, North-Holland, Oxford.
Freeland, R. K. (2010). True integer value time series, AStA Advances in Statistical Analysis, 94,
217–229.
Kitagawa, G. (1987). Non-Gaussian state-space modeling of nonstationary time series (with discussion),
Journal of the American Statistical Association, 76, 1032–1064.
Kitagawa, G. (2010). Introduction to Time Series Modeling, Chapman & Hall/CRC, Boca Raton.
Leonenko, N., Savani, V. and Zhigljasky, A. (2007). Autoregressive negative binomial processes, Annales de l’Institut de Statistique de l’Universite de Paris, 51(1), 25–47.
McKenzie, E. (1985). Some simple models for discrete variate time series, Journal of the American
Water Resources Association, 21(4), 645–650.
Nelder, J. A. and Wedderburn, R. W. M. (1972). Generalized linear models, Journal of the Royal
Statistical Society Series A, 135, 370–384.
Rajarshi, M. B. (2012). Statistical Inference for Discrete Time Stochastic Processes, Springer, New
Delhi.
Rinstić, M. M., Bakouch, H. S. and Nastić, A. S. (2009). A new geometric ﬁrst-order integer-valued
autoregressive (NGINAR(1)) process, Journal of Statistical Planning and Inference, 139, 2218–
2226.
Savani, V. and Zhigljavsky, A. (2007a). Eﬃcient parameter estimation for independent and INAR(1)negative binomial samples, Metrika, 65(2), 207–225.
Savani, V. and Zhigljavsky, A. (2007b). Eﬃcient estimation of parameters of the negative binomial
distribution, Communications in Statistics — Theory and Methods, 35(5), 1–17.
Weiß, C. H. (2008). Thinning operations for modeling time series of counts — A survey, AStA Advancesin Statistical Analysis, 92(3), 319–341.
Weiß, C. H. (2009). A new class of autoregressive models for time series of binomial counts, Communications in Statistics — Theory and Methods, 38(4), 447–460.
West, M. and Harrison, J. (1997). Bayesian Forecasting and Dynamic Models, 2nd ed., Springer-Verlag,
New York.

Proceedings of the Institute of Statistical Mathematics Vol. 65, No. 2, 323–339 (2017)

339

Recent Development of Integer-valued Autoregressive Models
Masahiko Nakajima1 , Fumitake Sakaori2 and Yoshinori Kawasaki3,4
1 Graduate School of Science and Engineering, Chuo University
2 Faculty of Science and Engineering, Chuo University
3 The Institute of Statistical Mathematics

4 Department of Statistical Science, SOKENDAI

Integer-valued Autoregressive Models (INAR models) express current observations,which are integers that depend on past integer-value observations. INAR models completely diﬀer from Dynamic Generalized Linear Models and General State-Space Models, which employ an unobservable/latent process to model integer-valued time series. In
INAR models, the choice of the marginal and/or innovation distribution and the deﬁnitionof the ‘autoregressive part’ are very important to ensure compatibility (in the distributional sense) of both sides of an autoregressive model. Although publications sporadicallyappeared from the mid 1980s to the early 1990s, notable contributions were not reported
for more than a decade. However in the late 2000s, new results began appearing on aregular basis. In this paper, we introduce recent developments in INAR models beginning with Poisson INAR. For some theorems and propositions without proofs in existing
literature, we compiled our notes in the appendix. Additionally, we include new results
on an INAR model based on the diﬀerence of two Poisson distributions and use real data
analysis as an illustration.

Key words: Integer-valued time series data, thinning operator, decomposition of random variables, INAR(1)model, INAR(p) model, method of moment.