[連載] フリーソフトによるデータ解析・マイニング 第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(1,1)、条件付き分散を GARCH (1,1) にした diff (UKgas) を当てはめる例を示す。関数 garchFit は詳細な計算結果を返すので、データによっては夥しい結果が画面上に出力される。関数 summary は、係数の推定値および検定統計量、残差の検定統計量、情報量規準 ( AIC や BIC など) などを返す。

> library(fGarch)
> UKg.d<-diff(UKgas)
> UKg.m <-garchFit(formula= ~ arma(1,1) + garch(1,1), data = UKg.d, trace=F)
> 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.42750 12.57738 12.42165 12.48826

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

> plot(UKg.m)

Make a plot selection (or 0 to exit):

1: Time Series 
2: Conditional SD 
3: Series with 2 Conditional SD Superimposed 
4: ACF of Observations 
5: ACF of Squared Observations 
6: Cross Correlation 
7: Residuals 
8: Conditional SDs 
9: Standardized Residuals 
10: ACF of Standardized Residuals 
11: ACF of Squared Standardized Residuals 
12: Cross Correlation between r^2 and r 
13: QQ-Plot of Standardized Residuals

Selection:

garchFitの結果

garchFitの結果

図1 関数 garchFit の診断図

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

> UKg.fit <-ts(UKg.m@fitted, 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))

garchFitの結果

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

2.成分の分解

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

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

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

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

成分ごとのプロット

図3 成分ごとのプロット

3.多変量時系列

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

Y1 , … , Yt-p , … , Yt , … , Yt+p , … , Yn

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

Yt = A1Yt-1 + … + ApYt-p + Et

  モデルの中の係数 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))

成分ごとのプロット

図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))

Call:
ar(x = USe.d, aic = F, order.max = 2)

$ar
, , 1
log(M1) log(GNP)
log(M1)
0.5446 -0.1888
log(GNP)
0.1981 0.1295
, , 2
log(M1) log(GNP)
log(M1)
0.1231 0.04513
log(GNP)
0.1371 0.07110
$var.pred
log(M1) log(GNP)
log(M1)
1.063e-04 1.506e-05
log(GNP)
1.506e-05 8.494e-05

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

Mt = 0.545M t-1 + 0.123M t-2 - 0.189G t-1 + 0.045G t-2 + E t
Gt = 0.198M t-1 + 0.137M t-2 + 0.130G t-1 + 0.071G t-2 + E t

  式の中の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、deccensu などがある。
  最後に、統計数理研究所が1972年に発表し、1984まで開発し続けた TIMSAC シリーズのRパッケージ timsac について触れておく。
  パッケージ timsac は、現時点では、次のサイトからダウンロードでき、またインストールおよびマニュアルなどに関する情報を入手することができる。パッケージ timsac には、日本語のマニュアルが付いている。

  http://jasp.ism.ac.jp/ism/timsac/

  パッケージ timsac のバージョン1.0には、1972に発表された TIMSAC-72 を利用する15種類の関数、1974に年発表されたT IMSAC-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-stat.wharton.upenn.edu/~steele/Courses/956/RResources/GarchAndR/WurtzEtAlGarch.pdf