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

Rと時系列(2)

 1.ランダムウォークと単位根

 回帰分析と同様に時系列データ解析の主要な目的は、収集したデータを用いてモデルを作成し、将来の予測やシステムの制御などを行うことである。
 時系列データのモデルは

yt=ayt-1+et

で表現でき、かつ|a|=1(単位根)である場合、ランダムウォークと呼ばれる。ランダムウォークで表現される時系列データは非定常である。時系列データを分析する際には、まずデータの変動がランダムウォークで表現できるか、そうでないかを調べることが重要である。 時系列がランダムウォークで表現できるか否かを検定することを単位根検定と呼ぶ。単位根検定は「単|a|=1が存在する」という帰無仮説検定で、1970年代後半にDickey-Fullerによって初めて考案され、その後Phillips-Perron検定、McKinnons's検定などが提案されている。
 Rには、Phillips-Perron検定に関する関数PP.testがある。データlhのPP.testの使用例およびその結果を示す。

> PP.test(lh)

PP.testの結果

 p値は約0.03であることから単位根があると判断しがたい。
 パッケージtseriesには、Dickey-Fullerの単位根検定関数adf.testがあり、パッケージfSeries には単位根検定関数adfTest 、unitrootTestなどがある。これらのパッケージはCRANサイトからダウンロードできる。
 データUKgasを用いた関数adf.testの使用例を示す。

> library(tseries)
> adf.test(UKgas)

関数adf.testの結果

 p値は約0.7になり、仮説が採択される。これはデータUKgasのプロットでトレンドが見られる視覚情報と一致する。
 データによっては、非定常データの対数や差分を取ることで定常化することができる。例えば、データUKgasの1階差分のadf.testを用いた単位根検定のp値は0.01まで下がる。

> adf.test(diff(UKgas))$p.value

[1] 0.01

2.ARモデル

 時系列データの時間t-pからtまでの関係式

自己回帰モデル
を自己回帰モデルと呼ぶ。ai,i=1,2,・・・,pを自己回帰係数と呼び、pを次数(order)と呼ぶ。次数pであるARモデルを通常AR(p)で表す。eiは平均0、分散σ2の正規分布に従う残差である。自己回帰分析は、最も適切な次数pと自己回帰係数。aiを推定することであり、このプロセスをモデルの当てはめ、あるいはモデルの推定と呼ぶ。
 ARモデルを当てはめる方法としては、ユールウォーカー(Yule-Walker)法、最小2乗法、最尤法、Burg法など幾つかの方法が提案されている。モデルの評価にはAIC、BICなどの情報量規準が多く用いられている。

◇ AR関数による当てはめ
 Rには自己回帰モデルを求める関数arがある。関数arの書き式を次に示す。

ar(x, aic = TRUE, method= "",order.max = NULL, ...)

 引数xは時系列データオブジェクト、aicはモデルを評価する情報量規準AICを用いるか否かを指定する引数で、デフォルトはTRUEになっているので、AIC情報量規準に基づいてモデルの次数が決められる。
 引数methodには、自己回帰モデルを推定する方法、"yule-walker"、"ols"、"mle"、"burg"などの中から1つ選択して指定する。デフォルトはユールウォーカー(Yule-Walker)法が指定されている。olsは最小2乗法、mleは最尤法、burgは文字どおりBurg法である。引数order.maxには、次数の最大値を指定することができる。
 次に関数arにおけるユールウォーカー法によるデータlhの自己回帰モデルを推定する例を示す。

> (lh.ar <-ar(lh))

関数arにおけるユールウォーカー法によるデータlhの自己回帰モデルの推定結果

 関数arで作成した自己回帰モデルに関連する項目のリストは、関数summaryを用いて確認することができる。

> summary(lh.ar)

関数arで作成した自己回帰モデルの確認

 自己回帰モデルの次数は$oder、係数は$ar、aicの値は$aic、残差は$residで返すことができる。

> lh.ar$order

[1] 3

 この結果から自動的に求めたデータlhのARモデルの次数は3であることが分かる。ARモデルの係数を次に示す。

> lh.ar$ar

[1] 0.65340168 -0.06362084 -0.22694020

 上記の値を小数点以下3第位までに丸めた、データlhのAR(3)モデルを次に示す。

AR(3)モデル

◇ モデルの診断
 作成したモデルの適切さを判断するためには、残差分析が必要である。ARモデルにおける残差は平均0の正規分布に従い、残差の間には相関がないことを期待する。
 Rには、時系列データの独立性を検定する関数Box.testがある。関数Box.testは、Box-Pierce検定とLjung-Box検定を行う。この関数を用いて残差の独立性を検定することが可能である。
 次にAR(3)モデルの残差のLjung-Box検定の例を示す。結果から分かるように、残差が独立である仮説が採択される。

> Box.test(lh.ar$res, type="Ljung")

AR(3)モデルの残差のLjung-Box検定結果

 残差の正規性に関しては、パッケージtseriesに時系列データの正規性を検定する関数jarque.bera.testがある。関数jarque.bera.testは欠損値を含んでいるとエラーメッセージを返す。関数arが返す残差の中の先頭のp個(モデルの次数)は欠損値であるので、次のように切り取って用いる。

> temp<- window(lh.ar$res,start=4)
> jarque.bera.test(temp)

関数jarque.bera.testの結果

 これらの検定関数が返す結果はあくまでも参考までにしてほしい。その理由は、乱数を生成させ検定を行ってみると理解できるであろう。
 また、残差の正規性に関してはqqnormを用いて考察を行うことも1つの方法である。

◇ 予測
 モデルを求める目的の1つは予測である。自己回帰モデル関数arで求めたモデルによる予測値は関数predictを用いて求めることができる。
 作成したモデルを用いた、関数predictの使用例を次に示す。引数n.aheadには予測の期間を指定する。関数predictは、予測値$pred、標準残差$seを時系列オブジェクトとして返す。

> (lh.pr<-predict(lh.ar,n.ahead=10))

自己回帰モデル関数arで求めたモデルによる予測値

 次にデータlhを用いた関数predict による10期の自己回帰予測値および2倍の標準残差のプロットを作成する例を示す。

> SE1<-lh.pr$pred+2*lh.pr$se
> SE2<-lh.pr$pred-2*lh.pr$se
> ts.plot(lh,lh.pr$pred,SE1,SE2,gpars=list(lt=c(1,2,3,3),col=c(1,2,4,4)))
> legend(locator(1),c("実測値","予測値","2*SE"), lty=c(1,2,3),col=c(1,2,4))

データlhの自己回帰予測値

図1 データlhの自己回帰予測値

3.ARMAとARIMAモデル

 自己回帰モデル自己回帰モデルに移動平均モデル移動平均モデルを結合したモデル

自己回帰移動平均モデル
を自己回帰移動平均(ARMA; AutoRegressive Moving Average model)モデルと呼び、通常ARMA(p,q)で表す。定義式から分かるようにARモデルはARMAモデルの特殊なケースである。
 yid階の差分演算子差分演算子のARMAモデルを自己回帰和分移動平均 (ARIMA; AutoRegressive Integrated Moving Average ) モデルと呼び、ARIMA(p,d,q)で表す。ARIMAモデルは、1960年代末にBoxとJenkinsにより提案されていることから、彼らの名前を冠してbox-Jenkins法とも呼ばれている。
 Rには、単変量時系列データを当てはめるARIMAモデル関数arimaがある。次に関数arimaの書き式を示す。

arima(x, order = c(0,0,0),…)

 引数xは、時系列データオブジェクトであり、orderは自己回帰の次数 、差分の階数 、過去の残差の移動平均の次数 を指定する引数である。次にデータlhを用いたarima(2,0,1)の例を示す。

> (lh.ari<-arima(lh,order = c(2,0,1)))

データlhを用いたarima(2,0,1)の例

 返された係数を用いたARIMA(2,0,1)モデルを次に示す。

ARIMA(2,0,1)モデル

 関数arima結果の項目は、関数summaryで確認できる。

> summary(lh.ari)

関数arima結果確認

 モデル評価に関連する情報として、sigma2(残差の分散σ2)、対数尤度$loglik、AIC情報量規準$aic、残差$residualsなどがある。
 関数arimaを用いる際、最も基本的な問題は、引数order のpdqを幾らにするかである。pdqを求める1つの方法は、ある範囲内のpdqのすべての組み合わせの中から、情報量規準(AICやBICなど)値が最も小さい組み合わせを用いる方法である。
 次にデータlhを用いてpdqを求める例として、pは1から4まで、dは0から1まで、qは0から4までのすべての組み合わせで最適なpdqを求める簡単なプログラムを示す。情報量規準はAICを用いる。

> data<-lh; T<-0
> for(p in 1:4)
  for(d in 0:1)
   for(q in 0:4){
    fit<-arima(data,order=c(p,d,q))
    T<-T+1
    if(T==1){
     minaic<-fit$aic
     orderP<-p; orderD<-d;orderQ<-q;
    }else{
     if (fit$aic<-minaic){
     minaic<- fit$aic;
     orderP<-p; orderD<-d;orderQ<-q;
    }
   }
  }
>cat("結果: p=",orderP,"d=",orderD,"q=",orderQ,"AIC=",minaic,"\n");

結果: p= 3 d= 0 q= 0 AIC= 64.18482

 求めたpdqはそれぞれ、3、0、0である。この結果は、AR(3)モデルと同じである。

> fit<-arima(lh,order = c(3,0,0))

(p,d,q)=(3,0,0)のfitの結果

 返される結果を用いたARIMA(3,0,0)モデルを次に示す。

ARIMA(3,0,0)モデル
 ARIMAモデルを診断(残差分析)するツールとして関数tsdiagがある。データlhのARIMA(3,0,0)モデルの診断図を次に示す。診断図の上部は、残差のプロット、中部は残差の自己相関のプロット、下部は引数gof.lagに対応するLjung-Box検定のp値のプロットである。

> tsdiag(fit,gof.lag=12)

ARIMA(3,0,0)診断図

図2 lhのARIMA(3,0,0)診断図

 関数arimaで求めたモデルの予測値は、関数predictを用いて求めることができる。

4.ARFIMAモデル

 ARIMAモデルは、過剰の差分が起こり得ることがあると指摘されている。その短所を克服するため、差分の階数dが整数に限らず、任意の実数に一般化した方法が自己回帰実数和分移動平均(ARFIMA; AutoRegressive Fractionally Integrated Moving Average)モデルである。ARFIMA モデルは差分要素を2項級数展開することで実数dを取ることを可能にしている[1]。
 パッケージfracdiffは、ARFIMA分析の専用パッケージであり、CRANからダウンロードすることができる。パッケージの中には、データの適切な差分階数を推定する関数fdGPH、モデルの係数を推定する関数fracdiffなどがある。次に関数fracdiffの書き式を示す。

fracdiff(x, nar = 0,dtol = NULL,nma = 0,…)

 引数xは時系列データで、引数narは自己回帰の次数p、引数dtolは差分の階数d、引数nmaは移動平均の次数qを指定する引数である。
 ここではデータAirPassengersを用いて、まず差分階数を推定し、その結果を用いてARFIMA(3,d,1)モデルの係数を推定する例を示す。データAirPassengersは、ある航空会社の1949年から1960年までの国際旅客数に関する時系列データである。

> library(fracdiff)
> (AP.d<-fdGPH(AirPassengers)$d)

[1] 0.6513916 #差分パラメータd

> AP.fra<-fracdiff(AirPassengers,nar=3, dtol= AP.d,nma=1)
> summary(AP.fra)

ある航空会社の1949年から1960年までの国際旅客数に関する時系列データの要約

 最大対数尤度は$log.likelihood、自己回帰係数は$ar、移動平均係数は$maで返すことができる。

 [1] 刈屋武昭・田中勝人・矢島美寛・竹内啓(2003):経済時系列の統計―その数理的基礎、岩波書店。