統計数理（2019）
第 67 巻 第 2 号 181–192©2019 統計数理研究所

創立 75 周年記念号［総合報告］

時系列解析における状態空間モデルの利用
北川 源四郎†（受付 2019 年 1 月 7 日；改訂 3 月 18 日；採択 3 月 26 日）

要旨

時系列解析において状態空間モデルは，統計的制御における最適制御則の導出や ARMA モ
デルの最尤推定など陽な形で解が得られない問題に対する解を求めるための逐次計算を実現するための計算手段として導入されたが，1980 年頃からは非定常時系列モデリング，非線形モデ
リング，信号分離，異常値処理，自己組織型モデリング，データ同化など，様々な時系列モデ
リングを統一的に取り扱うためのプラットフォームとしての役割りを果たしてきた．その一方
で，状態空間モデルに関連した状態推定のために，カルマンフィルタ，ガウス和フィルタ，非ガウス型フィルタ，粒子フィルタなどの様々なアルゴリズムが開発されてきた．
本項では，統計数理研究所の研究を中心に，状態空間モデルと関連する計算法およびその応用について概観する．
キーワード：状態推定，非定常モデル，非線形モデル，カルマンフィルタ，非ガウス型フィルタ，粒子フィルタ．
1.

状態空間モデルとカルマンフィルタによる状態の推定

線形・ガウス型状態空間モデルは，dm 次元時系列 yn の変動を下記のように 2 つのモデルで
表現する．
(1.1)

xn = Fn xn−1 + Gn vn （システムモデル）(1.2)

yn = Hn xn + wn

（観測モデル）

ただし，xn は dk 次元状態ベクトル，Fn , Gn , Hn はそれぞれ dk × dk , dk × d , dm × dk の行列である．vn ∼ N (0, Qn ) と wn ∼ N (0, Rn ) はそれぞれ d 次元および dm 次元の正規分布に従う
白色雑音でシステムノイズ，観測ノイズと呼ばれる．通常はこの線形・ガウス型状態空間モデ
ルのことを単に状態空間モデルと呼ぶことが多い．物理的なシステムでは，過去の情報は現在を経由して将来に伝えられるが，状態空間モデルにおいてはこの状態 xn だけが分かれば，時系列の将来の動きは予測できることになる．
状態空間モデルに関して，時系列 Yj = {y1 , . . . , yj } が与えられたとき状態 xn を推定する問題が状態推定である．観測値の最終時点 j と推定する状態の時刻 n の大小関係により，予測，フィルタ（n = j の場合）
，平滑化（n < j の場合）と区別される．時系列の予測，（n > j の場合）
制御，補間，パラメータ推定，成分分解などの問題のほとんどがこの状態推定を利用することによって統一的に解決できる．
†

東京大学 数理・情報教育研究センター：〒113–8556 東京都文京区弥生 2–11–16

182

統計数理

第 67 巻 第 2 号 2019

初期状態 x0 の分布が正規分布 N (x0 , V0 ) で与えられるとき，状態 xn の条件付き分布 p(xn |Yj )は正規分布となる．したがって，線形・ガウス型状態空間モデルに対しては，カルマンフィル
タによって効率的に条件付き平均 xn|j と分散共分散行列 Vn|j を求めることによって状態推定が実現できる（Anderson and Moore, 1979; 片山, 2005; 北川, 2005）．
このように，状態空間モデルを利用すると，時系列解析に用いられる様々なモデルを統一的に取り扱うことができ，しかも予測，平滑化，成分分解，パラメータ推定などに必要となる状
態推定をデータ数に比例する計算量で効率的に計算できる．当初は，統計的制御における最適制御則の導出や ARMA モデルの最尤推定など陽な形で最適解が得られない問題に対して，逐次計算によって解を求めるための計算手段として導入されたが，1980 年頃からは非定常時系
列モデリング，非線形モデリング，信号分離，異常値処理，自己組織型モデリング，データ同化など，様々な時系列モデリングの問題を統一的に取り扱うためのプラットフォームとしての役割りを果たしてきた．その一方で，状態空間モデルに関連した状態推定のために，カルマン
フィルタ，ガウス和フィルタ，非ガウス型フィルタ，粒子フィルタなどの様々なアルゴリズムが開発されてきた．
本項では，統計数理研究所の研究を中心に，状態空間モデリングの発展と関連する計算法およびその応用について概観する．
2.

線形定常時系列解析における状態空間モデルの利用

2.1 統計的制御
化学プラント，火力発電所や大型船舶のように巨大で強い外乱を受けるシステムの制御では，制御対象となるシステムのモデルを求めることが困難なために，最適制御理論の適用が困
は多変量 AR モデルを制御入力 rn を持つ状態空間モデル難であった．Akaike（1970）
xn = F xn−1 + Grn + vn

(2.1)

yn = Hxnで表現し，2 次評価基準の期待値(2.2)

J(x0 , r) = E

 L
 T

T
xn Qxn + rn−1
Rrn−1





n=1

を最小化する制御入力を決定する方式を提案した．この方法は，セメントの焼成炉や火力発電所ボイラーの制御や船舶のオートパイロットの設計に適用され実用化されたが（赤池・北川,，これが時系列解析における，状態空間モデル利用のさきがけ
1994, 1995; Ohtsu et al., 2015）となった．
2.2 時系列モデルの尤度計算とパラメータ推定
状態空間モデルに含まれる未知のパラメータを θ とし，長さ N の時系列 YN = {y1 , . . . , yN }が与えられるとき，時系列モデルの尤度は同時確率密度関数を用いて L(θ) = fN (YN |θ) と定義
される．ここで fn (Yn |θ) = fn−1 (Yn−1 |θ)p(yn |Yn−1 , θ) という分解を繰り返し適用して，dm 次
元の密度関数に分解していくと最終的には(2.3)

L(θ) =N

p(yn |Yn−1 , θ)n=1と表現できる．ただし εn = yn − Hn xn|n−1 , Σn = HnT Vn|n−1 Hn とするとき，p(yn |Yn−1 , θ) =

時系列解析における状態空間モデルの利用

183(2π)−m/2 |Σn |−1/2 exp{−εTn Σ−1n εn /2} である．したがって，一般に時系列モデルの対数尤度は(2.4) log L(θ) =Nlog pn (yn |Yn−1 , θ) = −

n=1

dm N
1  T −11
log |Σn | −εn Σn εnlog 2π −
2
2
2
N

N

n=1

n=1

によって与えられる．
ここで，xn|n−1 と Vn|n−1 はカルマンフィルタで得られるので，状態空間モデルで表現できる時系列モデルに対してはカルマンフィルタの副産物として，対数尤度が統一的かつ自動的に計算できることになる．時系列モデルのパラメータの最尤推定値を求めるためには，数値的最適化によりこの対数尤度を最大とするパラメータを求めればよい．
2.2.1 ARMA モデルの最尤推定
AR（自己回帰）モデルに対しては Yule-Walker 法のほか，Parcor 法（Burg, 1967），Householder最小二乗法（Akaike et al., 1975, 1979; Kitagawa and Akaike, 1978, 1981）などの多くの実用的な推定法が提案されている（北川, 2005）．これらは近似的な最尤推定量を与えるものと解釈でき
るが，ARMA モデルに対しては残差のフィードバックがあるために精度の良い同様の近似解を得ることができなかった．Box and Jenkins（1970）では，backcasting によってイノベーショ
ン系列を生成し，近似尤度を計算する方法が利用されている．
Akaike（1974, 1978）
は ARMA（m, ）
モデル
yn −

(2.5)

m


aj yn−j = vn −

j=1

(2.6)

0
⎢ ..
⎢
F =⎢ .
⎣ 0
ak

1
..
.
0
ak−1

bj vn−j ,j=1に対して H = [100 · · · 0]，Q = σ 2 ，R = 0，

⎡




···
..
.
···
···

⎤

⎡

⎤

0
1
⎢ r1 ⎥.. ⎥⎢⎥. ⎥⎥, G = ⎢ . ⎥
.
⎦
⎣
. ⎦
1
a1
rk−1

と定義することによって線形・ガウス型状態空間モデルで表現できることを利用して最尤推定を実現した．ここで，rj は ARMA モデルから定まるインパルス応答関数，k = max(m,  + 1) でではカルマンフィルタは用いられておらず，また Akaike（1978）
ある．ただし，Akaike（1974）
では，より高速な Morf et al.（1974）のアルゴリズムが用いられておよび Akaike et al.（1979）
いる．
2.3 欠測値の処理
欠測値の処理は統計学において重要な課題であり，EM アルゴリズムが利用されることが多い．しかしながら，状態空間表現とカルマンフィルタを用いると欠測値の処理は極めて明快か
つ簡単である（Jones, 1980）．観測値 yn が実際には観測されなかった場合には，Yn = Yn−1 となるので，p(xn |Yn ) = p(xn |Yn−1 ) が成り立つ．これは予測分布とフィルタ分布が等しいことを意味するので，ある観測値が欠測となった場合には，対応するフィルタのステップを省略するだけでよいことになる．
特に，欠測値を推定したい場合には，フィルタリングを行った後に，平滑化を行えば，その平滑化分布が欠測値の推定値を与えることになる．ただし，パラメータ推定などでは，欠測値の補間をする必要はなく，実際に観測されたデータだけを用いて尤度を計算することができ

統計数理

184

第 67 巻 第 2 号 2019

る．Kitagawa and Akaike（1981），北川（2005）には，顕著なスペクトルピークを持つ時系列の場合には，数 10 パーセントの欠測値があってもよいパラメータ推定や欠測値の補間ができる事が示されている．
2.4 不等間隔時系列
時系列が不等間隔に観測されている場合には，推定すべき信号が確率微分方程式で表される
連続時間のモデル
(k)

(2.7)

ft

(k−1)

= a1 ft+ · · · + ak ft + vt

によって生成されているものと仮定する．ただし，ft(k) は ft の k 階微分，vt はウィーナー過程 Wt の独立増分で分散 σ 2 とする．このとき，この不等間隔時系列は係数 a1 , . . . , ak から求められる遷移行列 Ft を用いて
xt = Ft−s xs + Gvt,s(2.8)

yt = Hxt + wtと状態空間表現が得られる．ただし，G = Ik , H = [10 · · · 0], xt = [ft , ft(1) , . . . , ft(k−1) ]T

 t

(2.9)

vt,s =

 t

Ft−u BdWu ,Cov(vt,s ) =ssT
Ft−u BB T Ft−uduである．これによって，確率微分方程式モデルを不等間隔データにあてはめることができる（Kitagawa, 1984）．特に，ft(k) = dWt というモデルを想定すると，不等間隔データの平滑化を行うことができる．
3.

非定常時系列解析への利用

においてベイズモデリングに基づく経済時系列の新しい季節
赤池は Akaike（1980a, 1980b）
調整法を提案した．季節調整法では時系列を(3.1)

yn = Tn + Sn + wnとトレンド成分 Tn ，季節成分 Sn ，不規則成分 wn の 3 つ以上の成分に分解するが，未知数がデータ数より多いために，通常の最小二乗法や最尤法では意味のある解が得られない．赤池は対数尤度関数にペナルティ項(3.2)

d2 (Tn − 2Tn−1 + Tn−2 )2 + r2 (Sn − Sn−12 )2 + z 2 (Sn + · · · + Sn−11 )2を付加した最小化問題を考え，更にベイズ型情報量規準 ABIC を導入することによって，新し．
い季節調整法を提案した（Akaike, 1980b）
この方法は，単に新しい季節調整法を開発したということに止まらず，データ数を超える大規模パラメトリックモデリングに適用可能な汎用的な方法を提供するものであり，現代的に言えば L2 正則化パラメータを客観的に決定する方法を提案したという意味で画期的であった．
この方法を季節成分がない最も簡単な場合について，モデルの観点から見直してみるとトレンド成分に関して(3.3)

Tn = 2Tn−1 − Tn−2 + vn ,vn ∼ N (0, τ 2 )yn = Tn + wn ,wn ∼ N (0, σ 2 )という 2 つのモデルを想定することと同等であることがわかる．この二つのモデルは 2 次元の

時系列解析における状態空間モデルの利用

185

状態空間モデルとして表現できることから，トレンド推定や季節調整の問題が状態空間モデル
によって解決できること，および平滑化パラメータ d2 が分散比 σ 2 /τ 2 に対応することを示唆している．
3.1 非定常時系列モデリングによる成分分解
一般に，状態空間モデルによる成分分解は次のように定式化できる．時系列 yn を次のように J 個の成分 x(j)n , j = 1, . . . , J と観測ノイズ wn に分解するものとする．
(J)yn = x(1)n + . . . + xn + wn

(3.4)

さらに，それぞれの成分 x(j)n は状態空間モデルで
(j)

(3.5)

(j)
x(j)
xn−1 + G(j) vn （システムモデル）
n = F

(3.6)

yn = H (j) xn + wn(j)と表現できるものとする．このとき，

⎡
⎢

⎤

F (1)

F =⎣

..

（観測モデル）⎡⎥⎢⎦, G = ⎣

.

⎤

G(1)
..

F (J)
(3.7)

⎥
⎦,

.
G(J)

⎡

(1)

⎤

xn
⎢ .. ⎥
xn = ⎣ . ⎦(J)xn

H = [H (1) , . . . , H (J) ]とおくと，時系列の変動を J 個の成分の和として表現する状態空間モデルが得られる．それぞれの成分モデルを具体的に定めることによって，成分分解のための具体的な状態空間モデルが(J)定まる．このモデルを用いて，時系列を J 個の成分 x(1)n , . . . , xn と観測ノイズ wn に分解する
ことができる．以下では，この方法にもとづく，季節調整法，時変分散の推定法，時変係数モデルの推定法について順次紹介する．
3.2 状態空間モデルによる季節調整
観測した時系列をトレンドと呼ばれる傾向成分と季節成分および偶然変動に成分分解する方法が季節調整法である．
トレンドの代表的モデルとしては，2 次のトレンドモデル tn = 2tn−1 − tn−2 + vn が用いられることが多い．一方，時系列の変動のうち毎年繰り返して現れる変動パターンが季節成分である．1 年周期の季節成分モデルとしては，sn = sn−12 + un が最も自然であるが，トレンドのモデルと共通の固有根を持つため，識別性がないことが知られている．そのため，このトレンドモデルから共通の固有根を除去することによって得られるモデル，すなわち 1 周期の合計がほぼ 0 になるとしたモデル
sn + · · · + sn−p+1 = vn , vn ∼ N (0, τ 2 )(3.8)

が利用される（北川, 2005）．この季節成分モデルの状態空間表現は

⎡

⎤

sn
⎢ sn−1 ⎥
⎢
⎥
x(2)
⎥,
n = ⎢ .
.
⎣ .
⎦sn−p+1(3.9)

H (2) = [1 0 · · · 0 ]

⎡

−1
⎢ 1
⎢

F (2) = ⎢

⎣

−1
..

···

.
1

−1

⎤

⎥
⎥
⎥,
⎦

⎡

⎤

1
⎢ 0 ⎥
⎢ ⎥
G(2) = ⎢ . ⎥ ,⎣ .. ⎦
0

統計数理

186

第 67 巻 第 2 号 2019

で与えられる．p は周期の長さで，月次データの
p = 12 である．

 1 年周期をモデル化する場合はF (1) OO F (2)

この標準的季節調整モデルは F =, G=

G(1) OO G(2), H = [H (1) , H (2) ] とすること

により状態空間モデルで表現できる．ただし，O は適当な次元の 0 行列である．状態空間モ
デルを利用する方法の特長は，この標準的な季節調整法を簡単に拡張や改良できることであり，曜日効果項や定常 AR 成分を含むモデルに容易に拡張することができる（Kitagawa, 1981;．
Gersch and Kitagawa, 1983; Kitagawa and Gersch, 1984, 1996; 北川, 2005）
3.3 時変係数 AR モデル
自己回帰モデル(3.10)yn =m

aj yn−j + wn , wn ∼ N (0, σ 2 )j=1

は定常時系列解析の標準的モデルとして利用されている（赤池・中川, 1972; 北川, 2005）．このAR モデルにおいて自己回帰係数 aj を時間とともに変化する係数 an,j で置き換えたものが時変係数 AR モデルである．ここで係数変化のモデルとしてランダムウォークモデル(3.11)

an,j = an−1,j + vn,j , vn,j ∼ N (0, τ 2 )を導入し，m 個の時変 AR 係数を縦に並べて状態ベクトルを xn = [an,1 , . . . , an,m ]T と定義すると，自己回帰モデルと係数の変化を表すランダムウォークモデルは状態空間モデルによって表現することができる．ただし，vn ∼ N (0, τ 2 Im )，（Im は m × m の単位行列）(3.12)

F = G = Im , Hn = (yn−1 , . . . , yn−m )である．したがって，時間変化する AR 係数をデータから推定することができる．また，推定された時変係数 aj,n から時間とともに変化するスペクトルを推定することができる（Kitagawa,．
1983）
3.4 時変分散モデル
時系列 yn は正規分布に従うが，その分散は時々刻々と変化し，yn = σn εn ，εn ∼ N (0, 1) となる場合を考える．このとき(3.13)

log yn2 = log σn2 + log ε2nとなるので，時変分散とノイズが分離される．ここで分散 σn2 の時間変化に，例えば，ランダ
ムウォークモデル(3.14)

2log σn2 = log σn−1
+ vnを仮定すると，この変換したデータのトレンドから時変分散を推定することができる．ただし，観測ノイズ log ε2n は二重指数分布に従うが，カルマンフィルタによる推定のためにはこれを正規分布で近似する必要がある（北川, 2005）．この方法は，地震波の分散の推定や，金融時系列の確率的ボラティリティの推定に用いられる．
4.

非線形・非ガウス型状態空間モデリング

状態空間モデルは様々なモデルを統一的に取り扱うことができ，しかも効率的な計算法が利用できる便利なモデリング・プラットフォームであったが，1980 年代の後半になると線形・ガ時系列解析における状態空間モデルの利用

187

ウス型モデルでは対応できない問題が多くなり，一般化した非線形・非ガウス型状態空間モ
デル
(4.1)

xn = Fn (xn−1 , vn ) （システムモデル）(4.2)

yn = Hn (xn ) + wn （観測モデル）が利用されるようになった．ここで，vn と wn はそれぞれ，確率密度関数 q(v)，r(w) に従う d
次元と dm 次元の白色雑音，Fn は dk 次元，Hn は dm 次元の非線形関数である．
ただし，モデルが非線形あるいはノイズの分布が非ガウス型の場合には，状態の分布は非ガウス型の分布となるので，直接カルマンフィルタを適用することができず，フィルタリングのためにはこれらの分布を何らかの形で近似することが必要である．非ガウス型の分布の近似法の違いによって様々なフィルタのアルゴリズムが得られている．
（1）正規近似による方法：非ガウス型の分布をひとつの正規分布で近似する方法で，非線形時系列モデルに適用できる拡張カルマンフィルタ（Anderson and Moore, 1979; 片山, 2011）が知られている．ただし，真の分布が 2 峰以上の場合には，よい近似は得られない．
m（2）混合正規近似による方法：非ガウス型の分布を p(xn ) = i=1 αin ϕi (xn ) と複数の正規分
布 ϕi の加重和で近似する方法である．この場合，モデルが線形であれば予測分布およびフィルタ分布もガウス分布の和で表現でき，しかもそれらのパラメータを m 個のカルマンフィル
タで計算できるという大きな利点がある（Kitagawa, 1994; Kitagawa and Gersch, 1996）．ただし，実際には，正規分布の項数が時間の進行とともに爆発的に増大するので，この方法の実装のためには各時間ステップで，項数の削減を行うことが不可欠である．
（3）数値的近似による方法：非ガウス型分布を多数の分点を持つ階段関数で近似するもので，分点数を増やせば任意の分布を高い精度で近似できる．非ガウス型フィルタ（Kitagawa, 1987;
北川, 2005) では，予測分布，フィルタ分布および平滑化分布を階段関数近似などを用いて数値的に表現し，状態推定の公式を数値積分などの数値計算によって実現している．計算機の高速化によって，この方法では，低次元の状態空間モデルに対しては，極めて精緻な結果を得ることができる．ただし，数値積分が必要なので 4 次元以上の連続状態の状態空間モデルに対しては適用が困難である．
（4）粒子近似による方法：粒子近似では，他の密度関数を近似する方法とは異なり，独立に生成されたとみなせる多数の粒子を用いて状態の分布を表現する．この近似に基づくフィルタとして粒子フィルタ（Gordon et al., 1993; Kitagawa, 1993, 1996）が提案され，複雑なモデルの場合でも実装が容易なために様々な問題に適用されている（Doucet et al., 2001）．ただし，必然的にサンプリング誤差を伴うので，少数の粒子を用いる場合には注意が必要である．粒子数の変化に伴う，粒子フィルタの精度の違いに関しては Kitagawa（2014）に詳しい．
4.1 非ガウス型分布モデル
状態空間モデルのノイズ項を適切な非ガウス型分布に変更することによって，正規分布モデ
ルとは異なる面白い結果が得られる．
4.1.1 構造変化の検出
システムノイズをコーシー分布や，ピアソン分布族などの裾の重い分布に変更することによって，トレンドや信号のジャンプが自動的に検出できる．混合ガウス分布や正規分布と一様分布の和なども利用される（Kitagawa, 1987, 1996; Kitagawa and Gersch, 1996）．

188

統計数理

第 67 巻 第 2 号 2019

4.1.2 異常値の検出
観測ノイズの分布に，コーシー分布や混合正規分布などを用いると異常値の影響を軽減で
き，ロバスト化が自動的に実現できる（Kitagawa and Matsumoto, 1996）．Kitagawa（1994）にはガウス和フィルタによって，季節調整などの高次元状態空間モデルに対しても混合正規分布を用いた異常値処理が可能であることが示されている．
4.1.3 非対称分布の利用：変動分散の推定
3.4 節で示した分散変動の推定では，二重指数分布 r(w) = exp(x − ew ) をこれと同じ平均，分散をもつ正規分布で近似していたが，非ガウス型フィルタあるいは粒子フィルタを用いると，二重指数分布をそのまま用いて，より正確に変動分散の推定値を求めることができる（北川,．
2005）
さらに分散 σn2 の時間変化が(4.3)

2
log σn2 = log σn−1
+ vnと表されるとき，状態を xn = log σn2 と定義すると，時系列の状態空間モデルは(4.4)

xn = xn−1 + vn
yn = exp(xn /2)wnと 表 さ れ る ．こ の よ う に ，時 変 分 散 モ デ ル は 一 般 型 の 非 線 形 状 態 空 間 モ デ ル に お い てF (x, v) = x + v ，H(x, w) = exp(x/2)w とおけばよい．
4.2 非線形平滑化
システムモデルが非線形の場合，従来は拡張カルマンフィルタによって状態を推定することが多かった．しかしながら，真の予測分布が多峰になることがある場合には，拡張カルマン
フィルタは全く見当違いの結果を与える場合があることが知られている．しかしこのような場
合でも，非線形モデルをそのまま利用して，非ガウス型フィルタあるいは粒子フィルタにより平滑化分布を求めると，真の状態の変化をよく再現できることが知られている（Kitagawa,．
1991; 北川, 2005）
4.3 自己組織型の状態空間モデル：状態とパラメータの同時推定，（1.2）
が未知パラメータ θ を含む場合，拡大した状態を zn = [xTn , θT ]T
状態空間モデル（1.1）と定義し，θ に対する時間変化のモデルを例えば θn = θn−1 + un と仮定すると，状態 zn に対する状態空間モデルが得られる．このモデルを用いた非線形フィルタ・平滑化によって，状態とパラメータの同時推定が実現できる（Kitagawa, 1998）．上記のモデルのノイズ項 un の分散を 0 とすると固定パラメータの推定，分散に正値を与えると時間的に変化するパラメータが推定できる．
4.4 離散過程
非ガウス型フィルタの導出過程を検討すると，状態 xn+1 および時系列 yn の条件付き分布が状態 xn だけで定まることが本質的であることがわかる．したがって，非線形・非ガウス型状
態空間モデルを更に一般化した一般化状態空間モデル(4.5)

xn ∼ qn (·|xn−1 ) （システムモデル）(4.6)

yn ∼ rn (·|xn )（観測モデル）時系列解析における状態空間モデルの利用189に対しても全く同様にフィルタを導出することができる．
この拡張したモデルを用いると，二項系列やポアソン系列などに従う離散系列も通常の時系
列と全く同じように取り扱うことができる．
4.5 データ同化
データ同化は観測データとの統合によりシミュレーションモデルを現実に合わせる方法であり，海洋学，気象学に止まらず多くの分野で使われるようになっている．逐次型のデータ同化
は，状態空間モデル観点からは極めて自然な方法であり，アンサンブル・カルマンフィルタあ
るいは粒子フィルタが利用されている（樋口 他, 2011）．
5.

まとめ

時系列解析は 1960 年代まではスペクトル解析などの周波数領域解析が中心であったが，1970
年代になるとフィードバックを持つ多変量システムの解析や制御のためには，時間順序を明示的に表現した時系列モデルの利用が必要であることが認識されるようになった．当初は，ARモデルや ARMA モデルなどの定常時系列モデルの同定や推定されたモデルを用いた分析や予測・制御の方法が個別に研究されていたが，状態空間モデルの導入により殆どの線形・定常時系列モデルが統一的かつ効率的に取り扱うことができるようになった．このような，モデリン
グの過程において情報量規準 AIC が指導原理として大きな役割を果たしたことは言うまでもない．
1980 年代に入ると，状態空間モデルの自由にモデルを合成できるというメリットを活かして様々なタイプの非定常時系列モデルが開発されたが，合成モデルは一種のベイズモデルであり，非定常時系列解析法の発展はベイズモデリングの発展と軌を一にしていたといえる．更に1990 年代になると，非線形性や非ガウス性を持った時系列の分析が重要になり，非線形・非ガウス型状態空間モデルが利用されるようになった．ただし，この場合には，従来のカルマンフィルタに替わって状態空間モデルが非線形の場合やノイズ分布が非ガウス型の場合でも適用できる新しいタイプの逐次フィルタおよび平滑化の方法が必要となった．その後，粒子フィル
タを始めとする統計計算法の開発によって，複雑なモデルを駆使した時系列解析が可能にな
り，状態空間モデルは広範な応用分野で利用されるようになっている．
ビッグデータが出現した現在では，データ同化などにおいて巨大な次元の状態空間モデルが必要になっている．今後は，そのような場合にも適用可能な計算法の開発だけでなく，データ
サイエンスの課題に即した新しいモデル形式自体の開発も期待される．

参

考

文

献

Akaike, H. (1970). Autoregressive model ﬁtting for control, Annals of the Institute of Statistical Mathematics, 23(2), 163-180.
Akaike, H. (1974). Markovian representation of stochastic processes and its application to the analysis
of autoregressive moving average processes, Annals of the Institute of Statistical Mathematics,
26(3), 363-387.
Akaike, H. (1978). Covariance matrix computation of the state variable of stationary Gaussian process,Annals of the Institute of Statistical Mathematics, 30, 499-504.
Akaike, H. (1980a). Likelihood and the Bayes procedure, Bayesian Statistics (eds. J. M. Bernardo etal.), 1-13, University Press, Valencia.

190

統計数理

第 67 巻 第 2 号 2019

Akaike, H. (1980b). Seasonal adjustment by a Bayesian modeling, Journal of Time Series Analysis,1(1), 1-13.
赤池弘次, 北川源四郎 編 (1994).『時系列解析の実際 I』, 朝倉書店, 東京.
赤池弘次, 北川源四郎 編 (1995).『時系列解析の実際 II』, 朝倉書店, 東京.
赤池弘次, 中川東一郎 (1972).『ダイナミックシステムの統計的解析と制御』, サイエンス社, 東京.
Akaike, H., Arahata, E. and Ozaki, T. (1975). TIMSAC-74, A time series analysis and control program
package-(1), Computer Science Monographs, No.5, The Institute of Statistical Mathematics,
Tokyo.
Akaike, H., Kitagawa, G., Arahata, E. and Tada, F. (1979). TIMSAC-78, Computer Science Monographs, No.11, The Institute of Statistical Mathematics, Tokyo.
Anderson, B. D. and Moore, J. B. (1979). Optimal Filtering, Prentice-Hall, Englewood Cliﬀs.
Box, G. E. and Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control, John Wileyand Sons, New Jersey.
Burg, J. P. (1967). Maximum entropy spectral analysis, Proceedings of 37th Annual International
Meeting, Society of Exploration Geophysics, Oklahoma City, Oklahoma, Oct. 31, 1967.
Doucet, A., De Freitas, N. and Gordon, N. (2001). Sequential Monte Carlo Methods in Practice,Springer, New York.
Gersch, W. and Kitagawa, G. (1983). The prediction of time series with trends and seasonalities,
Journal of Business and Economic Statistics, 1(3), 253-264.
Gordon, N. J., Salmond, D. J. and Smith, A. F. (1993). Novel approach to nonlinear/non-Gaussian
Bayesian state estimation, IEE Proceedings F (Radar and Signal Processing), 140(2), 107-113.
樋口知之, 上野玄太, 中野慎也, 中村和幸, 吉田 亮 (2011).『データ同化入門＝次世代のシミュレーション, 朝倉書店, 東京.
技術』
Jones, R. H. (1980). Maximum likelihood ﬁtting of ARMA models to time series with missing observations, Technometrics, 22(3), 389-395.
片山 徹 (2011).『非線形カルマンフィルタ』, 朝倉書店, 東京.
Kitagawa, G. (1981). A nonstationary time series model and its ﬁtting by a recursive ﬁlter, Journal ofTime Series Analysis, 2(2), 103-116.
Kitagawa, G. (1983). Changing spectrum estimation, Journal of Sound and Vibration, 89(3), 433-445.
Kitagawa, G. (1984). State space modeling of nonstationary time series and smoothing of unequally
spaced data, Time Series Analysis of Irregularly Observed Data, 189-210, Springer, New York.
Kitagawa, G. (1987). Non-Gaussian state space modeling of nonstationary time series, Journal of theAmerican Statistical Association, 82(400), 1032-1063.
Kitagawa, G. (1991). A nonlinear smoothing method for time series analysis, Statistica Sinica, 1(2),371-388.
Kitagawa, G. (1993). Monte Carlo ﬁltering and smoothing method for non-Gaussian nonlinear state
space model, Research Memo., No.462, The Institute of Statistical Mathematics, Tokyo.
Kitagawa, G. (1994). The two-ﬁlter formula for smoothing and an implementation of the Gaussian-sum
smoother, Annals of the Institute of Statistical Mathematics, 46(4), 605-623.
Kitagawa, G. (1996). Monte Carlo ﬁlter and smoother for non-Gaussian nonlinear state space models,
Journal of Computational and Graphical Statistics, 5(1), 1-25.
Kitagawa, G. (1998). A self-organizing state-space model, Journal of the American Statistical Association, 93(443), 1203-1215.
北川源四郎 (2005).『時系列解析入門』, 岩波書店, 東京.
Kitagawa, G. (2014). Computational aspects of sequential Monte Carlo ﬁlter and smoother, Annals of
the Institute of Statistical Mathematics, 66(3), 443-471.
Kitagawa, G. and Akaike, H. (1978). A procedure for the modeling of non-stationary time series, Annals

時系列解析における状態空間モデルの利用

191

of the Institute of Statistical Mathematics, 30(2), 351-363.
Kitagawa, G. and Akaike, H. (1981). On TIMSAC-78, Applied Time Series Analysis II (ed. D.
Findley), 499-547, Academic Press, Massachusetts.
Kitagawa, G. and Gersch, W (1984). A smoothness prior state space modeling of time series with trend
and seasonality, Journal of American Statistical Association, 79(386), 378-389.
Kitagawa, G. and Gersch, W. (1996). Smoothness Priors Analysis of Time Series, Lecture Notes in
Statistics, No.116, Springer, New York.
Kitagawa, G. and Matsumoto, N. (1996). Detection of coseismic changes of underground water level,
Journal of the American Statistical Association, 91(434), 521-528.
Morf, M., Sidhu, G.S. and Kailath, T. (1974). Some new algorithms for recursive estimation in constant,linear, discrete-time systems, IEEE Transactions on Automatic Control, AC-19, 315-323.
Ohtsu, K., Peng, H. and Kitagawa, G. (2015). Time Series Modeling for Analysis and Control: Advanced Autopilot and Monitoring Systems, Springer, Tokyo.

192

Proceedings of the Institute of Statistical Mathematics Vol. 67, No. 2, 181–192 (2019)

Use of a State Space Model in Time Series Analysis
Genshiro Kitagawa
Mathematics and Informatics Center, University of Tokyo

In time series analysis, the state space model was initially used as means for realizing
the maximum likelihood estimation of the ARMA model and for obtaining an optimum
gain in statistical control. Since 1980, it has played a role as a platform for handling various kinds of time series modeling, such as nonstationary time series modeling, nonlinear
modeling, signal extraction, missing-value processing, self-organizing state-space modelingand data assimilation. In pararell to the development of the use of the state-space model,
various algorithms extending the conventional Kalman filter, such as Gaussian-sum filter,
non-Gaussian filter, particle filter, etc. have been developed for state estimation of thegeneralized state-space model. In this article, we focus on the research of the Institute of
Statistical Mathematics, and outline the development of state-space models and related
calculational methods and their applications.

Key words: State estimation, nonstationary model, nonlinear model, Kalman ﬁlter, non-Gaussian ﬁlter, particle ﬁlter.