[連載] フリーソフトによるデータ解析・マイニング 第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)

Phillips-Perron Unit Root Test
data: lh
Dickey-Fuller = -3.6999, Truncation lag parameter = 3,  p-value = 0.03465

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

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

Augmented Dickey-Fuller Test
data: UKgas
Dickey-Fuller = -1.6079, Lag order = 4, p-value = 0.7393
alternative hypothesis: stationary

  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までの関係式

自己回帰モデル

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

(1) 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))

Call:
ar(x = lh)
Coefficients:

1 2 3
0.6534 -0.0636 -0.2269

Order selected 3 sigma^2 estimated as 0.1959

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

> summary(lh.ar)

Length Class Mode
order
1 -none- numeric
ar
3 -none- numeric
var.pred
1 -none- numeric
x.mean
1 -none- numeric
aic
17 -none- numeric
n.used
1 -none- numeric
order.max
1 -none- numeric
partialacf
16 -none- numeric
resid
48 ts numeric
method
1 -none- character
series
1 -none- character
frequency
1 -none- numeric
call
2 -none- call
asy.var.coef
9 -none- numeric

  自己回帰モデルの次数は $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) モデルを次に示す。

yt = 0.653yt-1 - 0.064yt-2 - 0.227yt-3 + et
 

(2) モデルの診断

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

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

Box-Ljung test
data: lh.ar$res
X-squared = 0.025652, df = 1, p-value = 0.8728

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

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

Jarque Bera Test
data: temp
X-squared = 7.913, df = 2, p-value = 0.01913

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

(3) 予測

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

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

$pred
Time Series:
Start = 49
End = 58
Frequency = 1
[1] 2.461588 2.272267 2.199151 2.262914 2.352194 2.423066
[7] 2.449223 2.441544 2.418779 2.398456

$se
Time Series:
Start = 49
End = 58
Frequency = 1
[1] 0.4425687 0.5286675 0.5525786 0.5527502 0.5592254
[6] 0.5665903 0.5688786 0.5689385 0.5692396 0.5697534

  次にデータ 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)))

Call:
arima(x = lh, order = c(2, 0, 1))
Coefficients:

ar1 ar2 ma1 intercept
1.1765 -0.5044 -0.5080 2.3946
s.e. 0.3990 0.2190 0.4517 0.0944

sigma^2 estimated as 0.1827: log likelihood = -27.6, aic = 65.2

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

yt = 2.3946 + 1.1765yt-1 - 0.5044yt-2 - 0.508et-1 + et

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

> summary(lh.ari)

Length Class Mode
coef
4 -none- numeric
sigma2
1 -none- numeric
var.coef
16 -none- numeric
mask
4 -none- numeric
loglik
1 -none- numeric
aic
1 -none- numeric
arma
7 -none- numeric
residuals
48 ts numeric
call
3 -none- call
series
1 -none- character
code
1 -none- numeric
n.cond
1 -none- numeric
nobs
2 -none- call
model
10 -none- list

  モデル評価に関連する情報として、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)))

Call:
arima(x = lh, order = c(3, 0, 0))
Coefficients:

ar1 ar2 ar3 intercept
0.6448 -0.0634 -0.0634 2.3931
s.e. 0.1394 0.1668 0.1421 0.0963

sigma^2 estimated as 0.1787: log likelihood = -27.09, aic = 64.18

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

yt = 2.3946 + 0.6448yt-1 - 0.0634yt-2 - 0.2198yt-3 + et

  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)

Call:
fracdiff(x = AirPassengers, nar = 3, nma = 1, dtol = AP.d)
Coefficients:

Estimate Std. Error z value Pr(>|z|)
d
4.604e-01 1.442e-01 3.192 0.00141
**
ar1
-1.159e-02 5.996e-05 -193.222 < 2e-16
***
ar2
3.295e-01 3.103e-01 1.062 0.28826
ar3
-1.266e-01 7.667e-05 -1650.804 < 2e-16
***
ma
-9.009e-01 1.197e-01 -7.527 5.2e-14
***

Signif. codes: 0 ‘***’  0.001 ‘**’  0.01 ‘*’  0.05 ‘.’  0.1 ‘ ’  1
sigma[eps] = 29.69193
[d.tol = 0.1, M = 100, h = 7.322e-06]
Log likelihood: -693.9 ==> AIC = 1399.788 [6 deg.freedom]

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

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