時変係数自己回帰モデル

統計数理第34巻第2号（1986）
時変係数自己回帰モデル
プログラムTVCARの紹介
統計数理研究所北川源四郎
（1986年9月受付）
1．はじめに
この解説の目的は，時系列解析のためのプログラムパッケージTIMSAC−84（Akaikeet a1．，
1985）中のプログラムTVCARの原理を簡単に説明し，実例を用いたがらその使用法や出力の
意味などを示すことにある．
TVCARは共分散に関し非定常な時系列に対して時間とともに係数が変化する自己回帰モ
デルをあてはめる．そのために，Kitagawa（1983）で用いられたスペクトル構造の時間変化に
関するモデルとKitagawa−Gersch（1985b）で用いたスペクトルの滑らかさに関するモデルの
2つを制約モデルとして採用している．TVCARでは，これらのモデルが状態空間表現により，
統一的に記述できることを利用して，効率的な時変係数の推定および尤度計算を実現している．
モデル中の未知のパラメータは，非線形最適化を用いて最尤法により求めている．
第2節ではTVCARのモデルおよび推定法を簡単に説明する．第3節ではTVCARの入力
の方法と出力の見方を示し，いくつかの注意を与えている．
第2節のモデルおよび計算法は，Kitagawa（1983）およびKitagawa−Gersch（1985a，b）と
極めて関連が深いが，スペクトルの滑らかさと時問的変化の両方にベイズ的モデルを導入した
TVCARのモデルを紹介した文献は本稿が初めてである．たお，TIMSAC−84中のLOCCAR
およびNONSPAも非定常スペクトルの推定を行なうものでTVCARと関連が深い．LOC−
CARの原理はAkaike（ユ979）に示されている．
2．TVCARのモデル
2．1時変分散モデル
プ（m）は正規白色雑音でその分散は時間とともに変化するものとする．すたわち，
7（m）∼jV（O，σ2（m））
とし，簡単のために，σ2（2m）＝σ2（2m−1），m＝1，2，… がたりたつものとする．ここで，
（1）
1
∫（m）＝一（プ（2m）2＋プ（2m−！）2）
・2
という変数を考えると，
62（2m）。
∫（m）∼
λ。
2
274
統計数理第34巻第2号 1986
という形をしていることがわかる．これより，さらに対数変換
（2）
∠（m）＝Io9∫（m）
を適用すると
！
才（m）∼1o9山λξ十109σ2（2m）
（3）
2
1
となり，1ogσ2（2m）とそのまわりの確率的な変動という形に分解できる．ここで1og一λξの密
2
度関数はexp（卜e‘）で与えられるが，その平均および分散は，それぞれγ一1og2，π2／6であ
る．（ただし，γはオイラー定数，d．57721）以下では，この分布を正規近似することにする．
ここで，σ2（m）の推定のために，次のようだモデルを導入することにする．
▽㌦（m）＝ε（m）
（4）
才（m）＝m（m）十δ（m）
ただし，ε（m）∼N（O，ソ2），δ（m）∼N（0，π2／6）．このモデルは容易に状態空間表現の形に変形
できる．したがって，カルマンフィルターおよび固定区間スムーザーを用いてオ（m）のトレン
ド，ρ（m）を推定できる．このときn（m）十γ一1og2が1ogσ2（2m）の，したがってexp／売（m）
十γ一1og2｝が時変分散σ2（2m）の推定値とたる．
2．2時変係数ARモデル
ッ（m）に対しAR係数が時間的に変化するモデル
（5）
ツ（m）＝Σα（ゴ，m）ツ（m一ゴ）十〃（n）
ゴ＝1
を考える．ただし，〃（m）は平均O，分散σ2（m）の正規白色雑音とする．
ここで，このARモデノレに対して時変（逆）周波数応答関数
（6）
m
1
1
λ（∫，・）＝1一君α（プ・・）…（一2π〃）■丁≦！≦丁
を考えることにする．このとき，y（m）の瞬間的なスペクトルが
σ2（m）
（7）
1
1
州・）r〃、）1・一丁≦！≦T
によって定義できる．また，力（！，m）を時刻mの関数と考えると時変スペクトルが得られる．
TVCARでは時変AR係数α（ノ，m）の推定のために，λ（！，m）すたわちスペクトルの滑らか
さに関連した2種類の制約モデルを導入している．
スペクトルの時間的変化
λ（∫，m）の々階差分の2乗積分値を考えると
（・）
1
∫：1・1州・）12〃一茗1・W，・）12
2
となる．したがって，AR係数の々階差分の2乗和を小さくすることにより，（左辺の意味で）ス
ペクトノレの時間変化を小さくすることができる．これは
時変係数自己回帰モデル
プログラムTVCARの紹介
275
▽莞α（ノ，n）＝o（ノ，n）（ノ＝！，…，m）
（9）
o（ノ，n）∼N（0，τ2）
というモデルを用いることに相当している．モデル（9）はKitagawa（1983）ではAR係数の時
問的変化を小さくするという目的で導入されたものである．
スペクトルの滑らかさ
各時刻nに対してはKitagawa−Gersch（1985b）で用いたスペクトルの滑らかさの制約を用
いることにする．すなわち
（1・）
見一∫引㌣m）r・
＝（2π）2足Σノ2尾。（ノ，m）2
．ゴ＝1
は周波数応答関数の后隠微分の意味で，スペクトルの滑らかさを表わしているので，右辺の項
を小さくすることにより，スペクトルの滑らかな推定値が得られる．これは
α（ノ，m）＝z4点（ノ，n）ノ：1，…，m
m庄（ノ，m）∼N（O，（〃庄）一2）
という制約モデルを導入することにより実現できる．后としては0，1，2，… などが考えられる
が，TVCARでは々＝0と2の二つを同時に考えて
（11）
δ（プ，m）＝m（ノ，m）
m（ノ，m）∼W（O，（λ看十ブ4λ…）一］）
というモデルを用いている．
時変係数ARモデル（5）と2つの制約モデル（9），（11）は状態空間表現の形で表現できる．簡単
のために以下では々二1の場合だけを示す．
［二（二1））」十に二））」・「1（ン1））l
／ザ∵∴川∴∴」
ただし
（∵・（（1∴））
（∵（け∴∵））
276
統計数理第34巻第2号 1986
である．このように，m次の状態ベクトル，（m＋1）次の観測ベクトルを用いて，3つのモデル
が表現できる．したがって，カルマンフィルター，固定区間スムーザーを用いて，α（ノ，n）（ノ＝
1，…，m；m＝1，…，N）が推定できる．状態空間表現に含まれるパラメータ，τ2，λ看，λ姜はカル
マンフィルターで求められる尤度の値を最大化することによって定められる．
2．3TVCARの推定手続き
TVCARでは次の順にモデル推定を行なう．
（1）原データの変動分散を2．1の方法で推定する．この平方根をプログラムではENVE−
LOPEFUNCTIONと呼んでいる．原データをENVELOPEFUNCTIONで割ることに
より，NORMALIZED DATAが得られる．このデータは近似的に分散1の分散定常時系
列と見たせる．
（2）NORMALIZED DATAに2．2の方法で変動係数ARモデルをあてはめる．
（3）変動係数ARモデルの残差系列にさらに2．1を適用して，その変動分散を推定する．
3．TVCARの使用法
3．1TVCARへの入力
TVCARを実行させるためには，NAMELIST文による制御パラメータと解析の対象となる
原データの二種類の入力が必要である．
＜制御パラメータの入力〉
プログラムの制御パラメータとしての次のようたものがある．制御パラメータには標準値が設
定されているので，変更する必要があるものだけをNAMELIST文を用いて入力すればよい．
制御パラメータ
MT
10
MAR
KAR
NOBS
標準値
8
AR次数
2
AR係数の変化を表わす差分モデルの階差
10
IDIF
説
明
原データの入力機番
＝5：カードリーダ
1
何ステップ毎に係数を変化させるか
Gradient推定のための数値差分
＝1片側差分
＝2両側差分
MESH
1
グラフの座標上にmeshを描く
＝0 描かない
NOUT
TAU1
TAU2
TAU3
0
10■4
0，3
10−5
システムの急激な変化の回数
τぞ
τ多
τξ
ただし，NOUTに正の値を入れたときは，NAMELIST文のあとに，急激た変化がおこった場
所を示すNOUT個の整数をブリーフォーマットで入力しなければたらたい．
後で示す計算例の場合には
時変係数自己回帰モデループログラムTVCARの紹介一
277
△＆PARAM△NOUT＝2，＆END
635△1030
の2行を入力した．ただし△は空白を表わす．
＜原データの入力＞
原データは機番10の入力装置（NAMELIST文により変更された場合はその機番）から入力
される．原データはDECOMPと同様次のフォーマットに従っているものと仮定する．
レコード番号
内
容
1
データのタイトル
2
データ数
フォーマット
（20A4）
（I5）
3
データのフォーマット
（20A4）
4以降データY（1），…，Y（N）
3で指定されたフォーマット
以下の計算例では次の様なデータが用いられた．
1…＾RTHQ∪＾1く1≡ D＾T＾：同YElF
2600
（20F4 0）
・1・1−1−1
0−1
0 −10−1
0−2−2−2−1
＾2 −2 −1 −1 −2 −1
−1 −2
−1 −30−1
−1 −1 0−1−1−2
0 −1 −1
1
1 −1
0 1 0−1−1−1 0 0−2−2−3−1−2−2−3−3−2−1−2−2
−2
−3−11−1−1−21010−3−122−1−1030−1
20001−101000−200−1−30−1−2−4
−2
−3
−4
−2
−4
−5 −4 −2 −2 −4 −7 −3 −3 −3 −4
−5
−6
−1
−2
−3
−5−400−2−20001−1032213322
112201−1001−1−1−2−1−2−3−3−3−3−3
−5−4−3−3−2−3−4−3−1−1−3−1
0−1 0−1−1 1 0−1
0100102010012200120−1
以下略
3．2TVCARの出力の説明
TVCARの出力には数値と図がある．数値はラインプリンタ（又は端末），図はレーザープリ
ンタ（又はプロッタ）から出力される．
＜数値出力〉
TVCARが実行され，NAMELIST文が入力されると，まず次の様な出力が得られる．
・一・
PR06R＾”
MT
g 8
1D I F 昌
TVC＾R
一・・
1
T＾u1 昌 O．0001000
丁＾u2 ， 0．3000000
T＾∪3 ＝ 0．0000100
N08S 昌 10
NO∪T g 2
Ml≡SH g 1
■○■
・（O∪丁（I，・ I，1．．．．，∼OU1I ？
業■■
これによってNAMELIST父および標準値によって設定されたパラメータを確認できる．
NOUTに正の値を設定した場合には，ここでその個数だけ構造変化の場所KOUT（ゴ）を入力
する必要がある．
次に，機番MTから原データが入力されると次のようだ出力が得られる．
〈。㎝IGIN阯O＾T＾X（I）（ト1．”）・・ N・2600 F㎝…い（Z昨仙）
同Y l…1ド
278
統書十数理第34巻第2号 1986
これにより，入力されたデータのデータ長，入力フォーマット，タイトル等が確認できる．
次にサブルーチンEPARA2が非線形最適化法によりτ子，τ妻およびτξの最尤推定値を求め
る．この過程で出力される情報は次のようなものである．
PARAMETER：非線形最適化中のパラメータベクトル．ただし，適当た変数変換がほどこ
されてし・る．
GRADIENT ：関数値のパラメータに関する偏微分
RAM
E
：線形探索のステップ幅
：そのときの関数値
最尤推定値が求まると改ぺ一ジされて，次のように，データ名，プログラム名，モデル，計
算日時，初期推定値，最尤推定値がまとめられた出力が得られる．
一一一
@ D＾T＾
一一一
一一一
@ MODEL
一一＿
一一一
@ PROGR＾M
MY E1F
M＾R＝8
N
I（＾R二2
TV＾R
DATE： 85−03−04
＜〈＜
一一一
14 48 06
TIME：
INITIAL ESTIM＾TES
一一
@ PAR＾METER VI≡CTOR
一一
@ GR＾D工ENT
＞＞〉
一一
一0．921034D＋01 0．120397D＋01
0．220927D＋00
I
1
2
3
＝
0．330481D＋02 −0．898171D−03
一2727，768
5463，536
TAU2〈I）
0．10000000D−03
0．30000000D＋00
0．10000000D−04
＜＜〈
FIN＾L l…STIM＾TES
＞〉〉
@ PAR＾METER VECTOR 一一
一〇．114116D＋02 0．937851D＋00
一一
一一
@ GR＾DIENT
＾IC
＝
0．115162D＋02
一一
0・226783D−02
LI1〈ELIHOOD ≡
I
0．115129D＋02
一一
LIKELIトlOOD ＝
＾IC
2600
0．828236D−02 −0．315199D−03
一2723，253
5454，506
T＾u2（工）
1
0’11066335D−04
3
0．99674242D−05
20．39146833D＋00
＜図による出力〉
図1−1は原データである．この地震データの場合にはm＝630．1030の付近で急激に分散が増
大し，その後，徐々にもとの分散にもどっていく様子がわかる．急激た分散の増大はP波とS
波の到着に対応している．
図↓一2は（1），（2）により変換したデータと分散の対数の推定値である．（2）の変換により分散
の増大は平均値の上昇に対応していることがわかる．
図1−3は推定された時変分散により規格化されたデータである．この規格化により，非定常
時系列がほぼ分散一様の時系列に変換されていることがわかる．時変係数ARモデルはこの
データから推定される．
時変係数自己回帰モデループログラムTVCARの紹介一
60
279
．0RlC川＾L口＾T＾
40
20
0
ト三鮒1 榊紬砂〆
，20
．40
−60
200 400 600
800 −000 12口O 1400 1600 1800 2000 220口 2400 2600
図1−1．地震波データ
ENVELOPE I＝UNCT10N
l l
1．
2 ’ 一■r I■
一1
榊榊≡脇」I■
1ミ
．．
@ 岬
｝
一2
200 ’O0 600 800 −0口0 1200 1400 】600 1800 2000 2200 2400 2600
図ユー2、変換したデータの分散の対数
∼0RH＾LIZE口 1〕＾1一＾
榊
榊，1
榊
榊
榊榊
一2
一4
200 400 600 8口0 i000 1200 i400 1600 1800 2000 2200 2400 2600
図1−3 規格化したデータ
図2は推定された偏自己回帰係数の時間的変化を図示したものである．P波および∫波が到
着したn＝630．1030の付近では，偏自己回帰係数が急激に変化し，変動の分散だけでなく共分
散構造自体も大きく変化したことがわかる．またm＝2600では，ほ惇P波到着前の状態にも
どっていることがわかる．
図3は，（7）によって求めた時変スペクトルである．P波およびS波の到着による分散の変
化，スペクトルの変化がよくわかる．とくに，5波到着後，二つのスペクトルのピークが高周
波側へ移動しながら，徐々に常時微動のスペクトルにもどっていく様子がよくわかる．このよ
うな，周波数による到着時刻の変化は表面波については知られているが，S波に関してはよく
判っていない．
280
統計数理第34巻第2号 1986
CH＾I，1ClNG P＾RCOR
1．O
“’⊥ u．■一1ドーI「1’■．l L⊥．1．一一I−1一’■．1’■．一一F＾｛n．1■．一一1一一■・Tσ丁1
O．5
ノ！：
I
、
；
l
l
l
l
l
l
l
0．O
−O．5
一．O
0
200
1，O
“uI
400
60‘］
800
】OO口 1200 i400 1600 i800 2000 2200 2400 2600
u■一■■「一．一I「■一…1’＾一一1一一一一皿「■■■1一一一一■F■爪引■1−II−1二丁呵η
O．5
0．O
−O．5
−1．O
200
1 ．O
0 ．5
5
一一■一
400
600
8口0
1000 12口0 14D0 160口 18口0 2口O0 2200 2400 2600
u一一一■「一一“「．一■1一一一一「I一■T“一■1一■一一一1I’一’河百i−i一一■I1一“■TO百丁1
一
￨
1
P
1
I
I 1
I 』
− 1
− 1
■ 1
v
I
■
l 1 l
?1 l
l
? 一l 一l
一
I
h l
l l
? ，
1 1
00
、O
@■
一
1
lC 1l
．
1
一＿＿＿＿＿＿＿＿■＝一＿
1
II一■一
l
1
−1 ．O
0
1．O
I
h
■
一
，
−
I
P
■
■
．
I
一一■一一■一■■■’
h
’
I
h
l
■
一
■l
一■一一一一一．’’■一一一■■一一’．■■
’ 1
I ．
’ I
一 I ■
P 1
｡ I
v 1
｡一， 1
I 一
一 1
I 一
gー’一一ト■一一一トー一■一I皿一．・’1’’■一一1■一一■一1一一■一一一■’一．．一一■一I−I一’一一『一一一一H一一’■『
l l
一 1
h l
h 一
P 1 1
P 一 ’
P l 一
， 1
P −l l
l 1 1
I
．5
−05
一
＿＿＿＿＿ v‘凹’’L＿一＿＿L＿＿＿＿＿1＿＿＿＿＿一一一一’一■＿’＿＿＿I＿＿＿■＿1−I一一’一一＿＿＿＿一＿＿■I一’’＿＿■＿＿＿一■
l 1
−
1
ll ll I1
一一
1
I
，
2ロロ
40D 600
60口 800
800 1000
10Dロ−200
一2口0
14001600
1600 1800
18D0 2000
2000．2200
240D 2600
260［
200 400
！400
2200 2400
■一「一一一一「一一■TI“■1．一一■一1一■I一■「“1一．．一■FI＾σi−1一■一一1〔丁匝τ1
O．5
0．0
−O，5
−1．O
0200400500800100D12001400160018002000220024002600
図2、偏自己回帰係数の時間変化
3．3パラメータの急激な変化
TVCARで用いたモデルは本来，パラメータのゆっくりした変化に追従するためのものであ
る．したがって，本稿で示した地震波の様に急激た変化をともたう場合には特別な取扱いが必
要である．TVCARでは，この様な場合を想定して，システム雑音に異常値の存在を許してい
る．この異常値により，前後のパラメータ間の制約条件は無視されることにたり，したがって，
パラメータの急激な変化が可能とたる．この異常値の処理を全く行なわずに推定した図4と図
3を比較するとこの効果がはっきりと削る．ただし，現在のプログラムでは，異常値の数および
その場所はユーザーが指定しなければたらない．異常値の場所が簡単には判らない場合には，原
データおよび残差系列の異常値の自動検出を行なってみるのも一つの方法である．浪花（1986）
は異常値の個数と場所に対応する多くのモデルを求めAICによって，その場所を特定すること
を試みている．
時変係数自己回帰モデル
281
プログラムTVCARの紹介
2600
400
一1
一2
一2
0．O
図4．
謝
O．1 0．2
O．3
0・4
0．5
異常値処理しないで求めた時変スペクトル
辞
地震波のデータを提供していただいた北海道大学理学部の高波鉄夫氏に感謝いたします．
参考文献
Akaike，H．（1979）、On the Construction of Composite Time Series Mode1s，〃。α42ma∫e∬わm∫〃、
∫6α比ム∫m∫乙
Akaike，H．et a1．（1985）．
TIMSAC−84Part1，Com〃〃∫cゴmce Momogm助∫，22，The Institute of
282
統計数理第34巻第2号 1986
Statistica1Mathematics，Tokyo．
Kitagawa，G．（1983）．Changing Spectrum Estimation，∫ommZげ∫oma m∂γ必m‘ゴ。m，89，No．3，
433−445．
Kitagawa，G．and Gersch，W．（1985a）．A Smoothness PriorsTime−VaryingAR Coe伍。ient Mcde1i㎎
of Nonstationary Covariance Time Series，旭E亙rmm∫．λ．σ，30，No，1，48−56．
Kitagawa，G．and Gersch，W．（1985b）．A Smoothness Priors Long今R Mode1Method for Spectra1
Estimation，皿五五rmm∫．λ、C。，30，No．1，57−65．
浪花貞夫（1986）トレントを除去した経済時系列の非定常性について
金融研究，5，第4号，97−139．
構造的変化の統計的検討
，
Proceedings of the Institute of Statistica1Mathematics Vo1．34，No，2（！986）
Time−Varying Autoregressive ModeI
An Introduction to the Program TVCAR
Genshiro Kitagawa
（The Institute of Statistica1Mathematics）
A procedure for the identiication and spectra1estimation of nonstationary time series
is described here．The basic mode1is an autoregressive mode1with time−varying
coe田。ients．Two types of constraint mode1s are introduced to get reasonab1e estimates of
the coe箭。ients．The irst is the mode1for the time evo1ution of the coe冊。ients and the
second for the smoothness of the spectrum at each time．The basic mode1and these two
constraint mode1s are incorporated into a state space mode1form．Computationary
e節。ient Ka1man趾er and the丘xed interva1smoother are used to obtain the time−varying
coe伍。ients and the1ike1ihood of the mode1．The va1ues of the parameters which contro1
the smoothness of the estimated spectrum are determined so as to maximize the like1ihood
of the mode1．Time varying variance of the process is a1so estimated by a simi1ar state
space mode1．
The usage of the computer program TVCAR which rea1izes the procedure is
exemp1iied by exp1aining the samp三e inputs and the samp1e outputs for an actual seismic
data．Some remarks on the abrupt changes of the parameters are a1so given．
283

Download Report