R勉強会第六回は時系列分析だよ!!

Rで時系列データってどう扱うの?
データを定常にする
- ラグ処理
- 自己相関係数
シミュレーションでデータを生成
- シミュレーションでデータを生成してみよう
- sarimaのシミュレーション
モデルの同定
arimaモデル
- ARIMAモデル
モデルの診断
モデルを使った予測
参考文献

Rで時系列データってどう扱うの?

時系列オブジェクトを生成

4半期ごとのデータ

> ts(1:12, frequency = 4, start = c(2000, 1)) 
     Qtr1 Qtr2 Qtr3 Qtr4
2000    1    2    3    4
2001    5    6    7    8
2002    9   10   11   12

時系列オブジェクトの合併

2つの異なる系統のデータを一緒にしたい
例えば同じ年度の
- GDP
- 所得
ts.union関数

> a<-ts(1:12, frequency = 4, start = c(2000, 1)) 
> b<-ts(13:24, frequency = 4, start = c(2000, 1)) 
> ts.union(a,b)
         a  b
2000 Q1  1 13
2000 Q2  2 14
...
2002 Q3 11 23
2002 Q4 12 24

lhデータ

時系列データで使うデータの説明
ある女性の血液中の黄体形成ホルモンについて10分の感覚で測定したデータ

> class(lh)
[1] "ts"

データを覗いてみる

> lh
Time Series:
Start = 1 
End = 48 
Frequency = 1 
 [1] 2.4 2.4 2.4 2.2 2.1 1.5 2.3 2.3 2.5 2.0 1.9 1.7 2.2 1.8 3.2 3.2 2.7 2.2 2.2
[20] 1.9 1.9 1.8 2.7 3.0 2.3 2.0 2.0 2.9 2.9 2.7 2.7 2.3 2.6 2.4 1.8 1.7 1.5 1.4
[39] 2.1 3.3 3.5 3.5 3.1 2.6 2.1 3.4 3.0 2.9

UKgasデータ

もうひとつデータ
1960年から1986年までのイギリスのガス消費量を四半期ごとに観測した時系列データ

> UKgas
       Qtr1   Qtr2   Qtr3   Qtr4
1960  160.1  129.7   84.8  120.1
1961  160.1  124.9   84.8  116.9
...
1985 1087.0  534.7  281.8  787.6
1986 1163.9  613.1  347.4  782.8

出力が長いので要約した統計量を知りたい

開始時刻、終了時刻などなどを知りたい

開始時刻

> start(UKgas)
[1] 1960    1

終了時刻

> end(UKgas)
[1] 1986    4

時間ごとの観測回数

> frequency(UKgas)
[1] 4

時系列データの一部を取り出したい

window関数
第一引数:データ
第二引数:開始時刻を表すベクトル
第三引数:終了時刻を表すベクトル

> window(UKgas,c(1975,2),c(1978,3))
      Qtr1  Qtr2  Qtr3  Qtr4
1975       321.8 177.7 409.8
1976 593.9 329.8 176.1 483.5
1977 584.3 395.4 187.3 485.1
1978 669.2 421.0 216.1

時系列データをプロット

> plot(lh)

ガス消費量のプロット

> plot(UKgas)

複数種類の時系列データを一度に表示する

ldeaths
1974年から1979年までのイギリスで喘息、気管支炎、肺気腫による死亡数を月ごとに記録したもの
さらにこれを男女別に分けたmdeaths、fdeathsがある

ts.plot(ldeaths,mdeaths,fdeaths,lty=c(1:3),col=c(1:3))

データを定常にする

ラグ処理

時系列データでは、時間をシフトして比較することも少なくない
これをラグと呼ぶ
lag関数
- lag(x,k=)
  - kだけ戻る
  - 負の値も取れる
これ使うと差分も簡単

実際にやってみる

> ldeaths
      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1974 3035 2552 2704 2554 2014 1655 1721 1524 1596 2074 2199 2512
...
1979 3084 2605 2573 2143 1693 1504 1461 1354 1333 1492 1781 1915
> lag(ldeaths,k=-3)
      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1974                3035 2552 2704 2554 2014 1655 1721 1524 1596
...
1980 1492 1781 1915                                             
> lag(ldeaths,k=3)
      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1973                                              3035 2552 2704
1974 2554 2014 1655 1721 1524 1596 2074 2199 2512 2933 2889 2938
...
1978 1969 1870 1633 1529 1366 1357 1570 1535 2491 3084 2605 2573
1979 2143 1693 1504 1461 1354 1333 1492 1781 1915
> ldeaths-lag(ldeaths,k=12)
       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
1974   102  -337  -234    57   144   -71   114   -21   200   287   123  -325
...
1978  -269   532   106  -174   177   129    68    12    24    78  -246   576

データの定常化

差分を取る
季節性を考える

diff関数

差分を取るための関数
differences
- 差分の回数を指定
- ディフォルトでは1回
lag
- 差分を取るための間隔を指定
- ディフォルトでは前の期

x <- 1:50
x
plot(x)

diff(x)
plot(diff(x))

diff(x,differences=2)
plot(diff(x,differences=2))

x <- c()
#何か適当に季節性を持たせてみた
for(i in 1:50){
  if(i%%7==0){
    x[i] <- rnorm(1)+10
  }else{
    x[i] <- rnorm(1)
  }pp
}
x 
plot(as.ts(x))

plot(diff(as.ts(x),lag=7))

plot(diff(UKgas,differences=1,lag=4))

自己相関係数

時系列の相関係数版
ただただ、相関を見るだけじゃなくって、データの定常性を見ることができる
acf
データが定常じゃないと直線になってしまう
データが定常だと
- 指数的に減衰
- sinカーブを描きながら、減衰

acfを使って、データが定常になるまで差分を取る

co2のデータで遊びます
diffを取ったり、lagを変えてみたりして定常にしてみよう

co2
tsp(co2)
acf(co2)

答え

acf(diff(co2))

acf(diff(co2,lag=12))

acf(diff(co2,lag=12,differences=2))

acf(diff(co2,lag=12,differences=3))

acf(diff(diff(co2,lag=12)))

もう一個練習

以下のデータを
- plot
- acf
を使って、どのくらい差分や季節性を追えば定常になるか実験してください
ts_example.R · GitHub

答え

季節性を考えた差分を取って、その差分の差分を取る

plot(my.data)

acf(my.data)

acf(diff(my.data))

acf(diff(my.data,lag=12))

acf(diff(my.data,lag=12,differences=2))

acf(diff(my.data,lag=12,differences=3))

acf(diff(diff(my.data,lag=12)))

ちなみに

自己相関係数のプロットで主観的に決めるんじゃなくって、統計量に基づいて決める方法もある
単位根の検定(Unit Root)
- 紹介程度で

library(fseries)
adf.test(lh)
PP.test(lh)

シミュレーションでデータを生成

シミュレーションでデータを生成してみよう

さっきのデータは人工的に作ったもの
Rではモデルに従う時系列データを生成できる
arima.sim
差分を適当に変化させてacfを取ったりして遊べる
- 定常性を見るための練習に使える

x <-arima.sim(n=1000, list(ar=.6, ma=.3, order=c(1,1,1)))
plot(x)
acf(x)
acf(diff(x))

x <-arima.sim(n=1000, list(ar=.6, ma=.3, order=c(1,2,1)))
plot(x)
acf(x)
acf(diff(x))
acf(diff(diff(x)))
acf(diff(x,differences=2))

sarimaのシミュレーション

arima.simはあるけど季節性を持つデータが作れない
- じゃあ、作っちゃえ
さっき定常になるように差分取っててもらったデータはこれから生成したもの

my.sarima.sim <- function (
   n = 20,
   period = 12,
   model,
   seasonal
) {
 x <- arima.sim( model, n**period )
 x <- x[1:(n**period)]
 for (i in 1:period) {
   xx <- arima.sim( seasonal, n )
   xx <- xx[1:n]
   x[i + period ** 0:(n-1)] <-
     x[i + period ** 0:(n-1)] + xx
 }
 x <- ts(x, frequency=period)
 x
}

x <-my.sarima.sim(
 20,
 12,
 list(ar=.6, ma=.3, order=c(1,1,1)),
 list(ar=c(.5), ma=c(1,2), order=c(1,1,2))
)
#実はさっきやっていたこと
x
plot(x)

plot(diff(diff(x,lag=12)))

acf(x)
acf(diff(x,lag=12))
acf(diff(x,lag=12,differences=2))
acf(diff(diff(x,lag=12)))

モデルの同定

偏自己相関係数の直感的理解

自己相関はなんとなく分かる
- じゃあ偏自己相関って何だ?
計量時系列分析
- ユールウォーカ方程式の解が…
- 厳密にはそうなんですけど、(直感的な)意味が分かりません
直感的理解をしよう
- arモデルを考える
今のデータはずーと昔のデータと(ちょっとだけど)関係がある
- 自己相関係数が減衰しながらも残っていることからも分かる
  - y_tとy_(t-k)のデータの相関を見たいときに、その間のデータが関連している
他の影響を取り除いた、y_(t-k)がy_tに及ぼしている影響を自己相関でははかることができない＞＜
- じゃあ、他の影響を取り除いた上での相関を見てみようよ!!
- →偏自己相関係数

自己相関係数、偏自己相関係数から見るAR、MA、ARMAの特徴

プリントの2枚目を見てみよう
実際に偏自己相関係数を計算してみよう
- すごく簡単

pacf(x)

ARモデル

Auto Regressive Model
ar
モデルの次数
- 偏自己相関係数から大体判断できる
  - どこで切断されているか
- AICで(勝手に?)決めれられる
推定の方法
- yule-walker
- mle
- ols(最小二乗法)
- などなど
- ディフォルトではユールウォーカー法

使うデータ

lh
?lh
ある女性の血液中の黄体形成ホルモンについて10分ごとで計測したデータ

モデルを推定する

> par(mfrow=c(3,1))
> plot(lh)
> acf(lh)
> pacf(lh)
> 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
> summary(ar(lh))
             Length Class  Mode
order         1     -none- numeric
ar            3     -none- numeric
var.pred      1     -none- numeric
...
asy.var.coef  9     -none- numeric

arimaモデル

ARIMAモデル

ARMAモデルの差分も考えたもの
(p,d,q)がパラメータ
- pはARの次数
- dは差分の次数
- qはMAの次数
次数の指定が必要
- 定数項に気を付けろ!!

> arima(lh,order=c(1,1,1))

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

Coefficients:
         ar1      ma1
      0.6060  -0.9919
s.e.  0.1383   0.3126

sigma^2 estimated as 0.2033:  log likelihood = -30.34,  aic = 66.68

> arima(lh,order=c(1,0,1))

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

Coefficients:
         ar1     ma1  intercept
      0.4522  0.1982     2.4101
s.e.  0.1769  0.1705     0.1358

sigma^2 estimated as 0.1923:  log likelihood = -28.76,  aic = 65.52

モデルの診断

モデルの同定

Idetification
どの次数が適切か考えないといけない
必要になるもの
- 差分
- 自己相関係数
- 偏自己相関係数

par(mfrow=c(2,1))
acf(lh)
pacf(lh)

差分、自己相関係数、偏自己相関係数すべて見る

par(mfrow=c(3,1))
plot(lh)
acf(lh)
pacf(lh)

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,fill=TRUE)

arima(lh,order=c(2,0,1)
arima(lh,order=c(3,0,0))

季節性も入ったarimaモデル→sarimaモデル

配布資料を見てみよう
diffを何回もするのとdiff(diff)の違いを押さえる
- dとd'の次数を決めるときにちゃんと分かってないとだめ

arima(lh,order = c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12))

説明変数も入っているarimaモデル→arimaxモデル

arimaxモデルはarimaモデルにt期の何かの説明変数を付け加えたもの
arimaモデルが分かっていれば、もう説明いらないよね?

arima(y,order=c(p,d,q),xreg=cbind(x1,x6)

モデルの診断

多変量回帰のモデルに診断が必要だったように、時系列分析にもモデルの診断が必要
なぜ?
- モデルを作る上でいろいろ仮定を置いていたから
  - 正規性
- 分散の不均一性
  - 誤差項が相関を持たない

qqnorm(fit$res)
qqline(fit$res)
abline(0,0.35)

モデルの診断2

残差プロット
- arimaモデルとかの約束事
- 誤差項の平均や分散が時間とともに変化するような傾向がないかチェック
残差の自己相関
- 今の残差が前の残差と相関関係を持っていないかチェック
- 相関関係は自己相関係数に出てくる

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

モデルを使った予測

予測

Rで予測と言えば?
- predict関数
n.aheadでどれだけ先を予測したいかを書く
予測値とその標準誤差を返す

> lh.pr <- predict(fit,n.ahead=6)
> lh.pr
$pred
Time Series:
Start = 49
End = 54
Frequency = 1
[1] 2.889903 2.850933 2.805867 2.764254 2.723736 2.683554

$se
Time Series:
Start = 49
End = 54
Frequency = 1
[1] 0.4984075 0.6868654 0.8065043 0.9091524 1.0002672 1.0812370

予測の情報を使って、プロットしてみよう

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(lty=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))