[連載]フリーソフトによるデータ解析・マイニング第35回

Rと時系列(3)

 1.GARCHモデル

 時系列データが、条件付き平均gt、条件付き分散htの正規分布N(gt,htに従うとき、htの変動を

自己回帰条件付き分散不均一モデル
で表現するモデルを自己回帰条件付き分散不均一(ARCH: AutoRegressive Conditional Heteroscedastic) モデルと呼び、通常ARCHモデルと呼ぶ。式の中のωは定数である。
 さらに、ARCHモデルを次のように拡張したものをGARCH(Generalized ARCH)モデルと呼ぶ。
GARCHモデル
また、条件付き分散の非対称性を表すTGARCH(Threshold GARCH)モデル、GARCHモデルとTGARCHモデルを統合したAPARCH(Asymmetric Power ARCH)モデルも提案されている[1][2][3]。
 パッケージtseries、fSeriesにはこれらのモデルに関連する関数が用意されている。これらのパッケージは、いずれもCRANミラサイトからダウンロードできる。パッケージtseriesには、GARCHモデルを当てはめる関数garchがあるが、ここではパッケージfSeriesの中の関数garchFitを用いたARMA-GARCHモデルを当てはめる例を示す。データはUKgasを用いるが、トレンドを除くため関数diffを用いて差分を取っている。
 関数garchFitでは、条件付き平均gtのARMAモデルを引数formula.mean で指定し、条件付き分散htのGARCHモデルやAPARCHモデルは引数formula.varで指定する。次に条件付き平均をARMA(2,0)、条件付き分散をGARCH(1,1)にしたdiff(UKgas)を当てはめる例を示す。関数garchFitは詳細な計算結果を返すので、データによっては夥しい結果が画面上に出力される。関数summaryは、係数の推定値および検定統計量、残差の検定統計量、情報量規準(AICやBICなど)などを返す。

> library(fSeries)
> UKg.d<-diff(UKgas)
> UKg.m <-garchFit(formula.mean = ~arma(1, 1),formula.var=~garch(1, 1),series = UKg.d)
> summary(UKg.m)

 返された結果を小数点以下第3位まで丸めたモデルを次に示す。

小数点以下第3位まで丸めたモデル

 ARIMA(p,d,q)モデルの次数p,d,qを決めるときと同じくプログラムを用いて最適なモデルを求めることができる。ただし、関数garchFitは指定するパラメータに対し計算ができない場合も少なくない。
 情報量規準は@fit$ics、残差は@residuals、当てはめ値は@fitted.valuesで返すことができる。

> UKg.m@fit$ics

AIC BIC SIC HQIC
-12.20320 -12.05332 -12.20906 -12.14244

 関数garchFitの結果を関数plotに代入し、図1に示す13種類のグラフを作成することができる。関数plotを実行すると作成可能なグラフのリストを返す。「選択:」の右にグラフの種類の番号を入力し、Enterキーを押すと指定されたグラフが作成される。数値0を入力し、Enterキーを押すと終了する。

> plot(UKg.m)

garchFitの結果
関数garchFitの診断図

図1 関数garchFitの診断図

 次に、用いた値(diff(UKgas))と予測値のプロットを作成するコマンドの例とその結果を図2に示す。

> UKg.fit <-ts(UKg.m@fitted.values, start=c(1960,2),frequency=4)
> ts.plot(UKg.d, UKg.fit,lty=c(1,2), col=c(1,2))
> legend(locator(1),c("差分値","予測値"), lty=c(1,2),col=c(1,2))

用いたデータと予測値

図2 用いたデータと予測値

2.成分の分解

 時系列データには、周期的に変動する周期性質を持つものも少なくない。一定の期間で周期的に変動する時系列データについて、それを幾つかの成分に分解することによって、より詳細に分析を行うことができる。
 典型的な例として次のモデルがある。

観測値=トレンド+周期変動+残差

 ここの場合のトレンドは、比較的に滑らかな長期変動を示すものであり、周期変動は、週・月・四半期など時間あるいは時節とともに一定のリズムで周期的に変動する成分で、季節変動(あるいは季節成分)とも呼ばれている。
 Rには時系列データを、トレンド、周期変動、残差に分解する関数stlがある。関数stlを用いてデータUKgasをトレンド、季節変動、残差に分解したプロットを図3に示す。関数stlでは、引数s.windowに季節成分を抽出する条件として、文字列"periodic"(略してperでもよい)を指定するか、あるいはラグの値を指定する必要がある。

> plot(stl(UKgas, s.window ="per"))

成分ごとのプロット

図3 成分ごとのプロット

3. 多変量時系列

 多変量時系列解析に関しては、最も広く知られているのはベクトル自己回帰 (VAR: Vector AutoRegressive) モデルである。
 多変量時系列の標本データを

多変量時系列の標本データ

とするとVARモデルは次のように定義される。

多変量時系列の標本データ

 モデルの中の係数A1,A2,…,Apは行列、残差Etはベクトルである。
 VARモデルの当てはめには、関数arを用いることが可能である。ここではパッケージtseriesの中に用意されている、1954年の第1四半期から1987年第4四半期までのアメリカの4 変量の経済時系列データUSeconomicの中の先頭の2変量を用いることにする。第1変量はlog(GNP)=実質GNP の対数値(季調済み)、第2変量はlog(M1)=実質M1 の対数値(季調済み)である。データにはトレンドがあるので、1階差分を取って用いる。

> library(tseries);data(USeconomic)
> USe.d<-diff(USeconomic[,1:2])
> ts.plot(USe.d,lty=c(1,2),col=c(1,2))
> legend(locator(1),c(colnames(USe.d)[1], colnames(USe.d)[2]),lty=c(1,2),col=c(1,2))

diff(USeconomic[,1:2])のプロット

図4 diff(USeconomic[,1:2])のプロット

 自己相関は、単変量の場合と同じく次のように関数acfで求めることができる。

> acf(USe.d, na.action=na.pass)

自己相関プロット

図5 自己相関プロット

 また、相互相関は、次のように関数ccfで求めることができる。

>ccf(USe.d[,1],USe.d[,2],main="d.log(M1) & d.log(GNP)")

相互相関プロット

図6 相互相関プロット

 関数arを用いたVARモデルを当てはめる例を次に示す。関数arでは、引数aic=TRUEで計算すると自動的に最適な次数の推定値に基づいて結果を返す。ここでは、返す結果とモデルの対応関係を説明するため、次のように強制的にVAR(2)モデル(必ずしも最適なケースではない)の推定値を求める。

> (USe.ar<-ar(USe.d,order.max=2,aic=F))

VAR(2)モデルの推定値

 返された結果を用いた2変量のモデルを次に示す。

2変量のモデル

 式の中のMはlog(M1)の1階差分値で、Mはlog(GNP)の1階差分値である。残差のプロットを次に示す。

> plot(USe.ar$res)

VAR(2)モデルの残差のプロット

図7 VAR(2)モデルの残差のプロット

 作成したモデルの予測値は次のように求めることができる。ただし、現在のバージョンでは多変量時系列の予測値の推定標準誤差は計算できないので、se.fit=Fにする。

> USe.f<-predict(USe.ar,n.ahead=12,se.fit=F)

4.カオス時系列

 不規則に変動する時系列データを非線形的に解析する手法としてカオス理論に基づいた方法がある。
 カオス時系列解析には、多次元状態空間に埋め込む(embedding)方法、リカレンスプロット(recurrence plots)法がある。
  Rにはカオス時系列解析のためのパッケージtseriesChaosがある。パッケージtseriesChaosはCRANミラーサイトからダウンロードできる。誌面上の都合により、ここではパッケージの中の関数の使用例およびその結果を示すにとどめる。
 まず時系列データrossler.tsを3次元に埋め込み、3次元空間で図示するコマンドを次に示し、その結果を図8に示す。

> library(tseriesChaos)
> library(scatterplot3d)
> ro.em<-embedd(rossler.ts,m=3,d=8)
> scatterplot3d(ro.em,type="l",color=4)

埋め込み3次元プロット

図8 埋め込み3次元プロット

 埋め込み図示は、時系列データの振る舞い(平衡、周期、カオス、ノイズ)を考察するのに有効である。返された図からデータrossler.tsはカオス的時系列であることが分かる。
 次に、同じデータの一部を用いたリカレンスプロットを示す。リカレンスプロットは時系列の比較分析に有効である。近年、脈拍データのリカレンスプロットを用いて、健常な人とそうでない人を判別する研究事例も報告されている。

> recurr(lorenz.ts, m=3, d=0, start.time=10, end.time=15)

リカレンスプロットの例

図9 リカレンスプロットの例

5.その他のパッケージと関数

 以上で紹介したパッケージ以外にも、Rには幾つかの時系列データ解析の専用パッケージ、あるいは時系列データ解析の関数が含まれているパッケージがある。例えば、パッケージade、itsには豊富な時系列データの操作やモデル作成に関する関数が用意されている。また、長期記憶時系列モデルに関するパッケージfracdiff、時系列因子分析パッケージtsfa、マイクロアレイ時系列解析パッケージGeneTSなどがある。なお、パッケージbootには、時系列のブートストラップ関数tsboot、パッケージdynlmにはダイナミックな時系列データの回帰関数dynlm、時空間データ解析パッケージpastecsには時系列データの成分を分解する関数tsd 、decdiff、deccensusなどがある。
 最後に、統計数理研究所が1972年に発表し、1984まで開発し続けたTIMSACシリーズのRパッケージtimsacについて触れておく。
 パッケージtimsacは、現時点では、次のサイトからダウンロードでき、またインストールおよびマニュアルなどに関する情報を入手することができる。パッケージtimsacには、日本語のマニュアルが付いている。
 http://jasp.ism.ac.jp/timsac/
 パッケージtimsacのバージョン1.0には、1972に発表されたTIMSAC-72を利用する15種類の関数、1974に年発表されたTIMSAC-74を利用する11種類の関数、1978年に発表されたTIMSAC-78を利用する12種類の関数、1984年に発表されたTIMSAC-84を利用する関数decompと幾つかのデータセットが用意されている。

 本連載「Rと時系列」(2006年4号〜6号)について、旭川大学の姜興起教授から協力とコメントを頂いた。またパッケージtimsacについては、統計数理研究所の中野純司教授から情報を頂いた。誌面を借りて謝意を表す。

[1] 小暮厚之(2005):時系列解析法、講義資料。http://web.sfc.keio.ac.jp/~kogure/
[2] 刈屋武昭・田中勝人・矢島美寛・竹内啓(2003):経済時系列の統計―その数理的基礎、岩波書店。
[3] Diethelm Wurtz, Yohan Chalabi, Ladislav Luksan: Parameter Estimation of ARMA Models with GARCH/APARCH Errors- An R and SPlus Software Implementation
(http://www.itp.phys.ethz.ch/econophysics/ R/pdf/garch.pdf)