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

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


グラフィックスの復習

  • ltyって何でしたっけ?
  • colって何でしたっけ?
  • legendって何でしたっけ?
  • locatorって何でしたっけ?
    • 教えてません

ちなみに

  • lhは50個のデータなのに500期先を予測するとかバカなことをやってみた
    • 計量時系列分析の授業の最後で説明なかったところ
  • 予測値の分散
    • 直観的は期が離れるほど、予測値もぶれそう
    • ある程度期が離れると、分散は一定の値に収束してくる
  • だから、筒っぽい(?)形になってくる

モデル作成にいたるまでのフローチャート

  • 同定
  • 推定
  • 診断
  • 予測

参考文献

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

時系列解析入門

時系列解析入門

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

基礎からわかる時系列分析 ―Rで実践するカルマンフィルタ・MCMC・粒子フィルター (Data Science Library)

基礎からわかる時系列分析 ―Rで実践するカルマンフィルタ・MCMC・粒子フィルター (Data Science Library)