時系列編の続き: サンプルサイズが小さいときの情報量基準

概要

前回の[[R] 回帰分析で適切な方法を使わないとどうなるか (時系列編) - ill-identified diary]で, 「時系列分析の場合は線型過程のラグ項の次数が分からないことが普通ではないか (なので実用性に欠ける用例でないか)」という指摘をいただいたので, 過去3回*1で推定量の漸近的性質 (大数の法則) を説明してきた流れで, 次数の決定の目安によく使われる AIC や BIC などの情報量基準の漸近的性質について調べてみた. PDF換算で7ページ程度.

情報量基準の性質

詳しくは他の文献を当たって欲しいが, そもそも情報量基準は与えられたデータに対して, 当てはまりと予測性能のバランスのよいモデルを選ぶ指標であるため, 必ずしも「正しいモデル」を選ぶとは限らない. さらに情報量基準といっても, AIC と BIC では式が違うので, 両者の基準で選んだとして, 必ずしも同じモデルが選ばれるとは限らない. 例えば沖本 (2010, 経済・ファイナンスデータの計量時系列分析) では, 一般に BIC のほうが過剰適合に対するペナルティが厳しいため, AIC より次数が少ないモデルを選ぶ傾向があると述べている.

2016/10/09 加筆: AIC と　BIC の結果が異なるのは, 両者が何を目的として計算されているかが異なっているため. 派生型も含め, AIC 系は「予測精度」を, BIC 系は「真のモデルを選ぶ確率」を最適化する目的で作られている. この辺の話は, 小西 & 北川 (2004, 情報量規準), 赤池, 甘利, 北川, 樺島, & 下平, (2007, 赤池情報量規準 AIC – モデリング・予測・知識発見) などを読むと良い.

少し古いが, Hurvich & Tsai (1989, Regression and time series model selection in small samples) ではタイトルの通りサンプルサイズの小さい時に情報量基準で変数選択した場合の問題について言及されている. 特に時系列の場合, サンプルサイズが10~20と非常に小さい場合, AR(2) の変数選択は有限修正AIC (AICc) が最も適しているという結果を導いている. そこで今回は, 真のモデルを ARMA(1, 1) として, サンプルサイズが10~100の時に AIC, BIC, AICc でモデル選択を行った場合の結果をシミュレーションしてみた.

まず, ARMA モデルの特性について簡単に言っておくと, ARMA(1, 1) であるには,

$\begin{align} y_{t}= & \rho y_{t-1}+\theta\varepsilon_{t-1}+\varepsilon_{t}\end{align}$

の自己回帰過程の係数 $\rho$ について, $|\rho|<1$ である必要がある. 絶対値がこれを超えると, 定常性を失い, 和分過程となるため, ARIMA(1, 1, 1) モデルになってしまう. なぜそうなるか, そして他の次数の場合はどうなるかの詳しい話は, 沖本 (2010, 経済・ファイナンスデータの計量時系列分析) の2章や Hayashi (2000, Econometrics) の6章などを見るとよい. 一方で, 移動平均過程の係数 $\theta$ は特に制約はない. そこで今回は, $\rho$ , $\theta$ についてそれぞれ 0.1, 0.5, 0.9 の場合の組み合わせでシミュレーションを実施した.

結果

サンプルサイズを, 10から10きざみで100までの10通りについて100回繰り返し推定し, それぞれのケースで選ばれた次数をヒストグラムに表した. まずは $\rho$ , $\theta$ いずれも 0.1 の場合

f:id:ill-identified:20150630194116p:plain

f:id:ill-identified:20150630194120p:plain

f:id:ill-identified:20150630194125p:plain

各行は上から順に AR, MA, 和分過程の次数であり, 上2つは1が, 3行目は 0 になるのが正しい. 上のケースでは, どの情報量基準でも, どのサンプルサイズでも定常過程 (ゼロ次の和分過程) と判定されることが最も多いようだ. しかし, AR と MAの次数もゼロが最も選ばれやすいようだ. 特に BIC ではそれが顕著である. これが, 係数が小さいので次数ゼロと判断されたのだと考えると, 係数がより大きい場合なら, 正しいモデルが選ばれるようになるかもしれない. そこで, $\rho=0.5$ , $\theta=0.5$ で実施してみる.

f:id:ill-identified:20150630194351p:plain

f:id:ill-identified:20150630194359p:plain

f:id:ill-identified:20150630194403p:plain

サンプルサイズが大きくなるほど, AR と MA の正しい次数が選ばれる頻度が増えていくことがわかる. AIC と AICc ではあまり違いがないが, BIC はより次数の少ない側に分布が偏っているようだ. また, いずれの情報量基準でも言えることだが, サンプルサイズ100程度では, 正しく次数を特定できたのは50~60%程度しかない (サンプルサイズ100においてすべての次数が同時に一致するのも全体の50%程度だった).

では, さらに $\rho=.9$ とした場合の結果はどうだろうか.

f:id:ill-identified:20150630194436p:plain

f:id:ill-identified:20150630194441p:plain

f:id:ill-identified:20150630194448p:plain

今度は相対的に, AR の次数が正しく選ばれるケースが減っている. 3行目の和分の次数に着目すると, 1が選ばれることが増えている. $\rho$ が 1に近いため, 非定常と判定されることが多くなっている. その場合, 差分をとってから ARMA を当てはめるので, 次数が正しく選ばれない場合が増えるのだと考えられる. 一方, 移動平均過程の部分はほとんど変わらない. この傾向は $\theta=0.9$ の場合でも同じである.

以上から, ARMAモデルの次数を情報量基準で特定する場合,

サンプルサイズが小さい場合, ARMA モデルが正しく特定される可能性はかなり低い
係数の大小によっても, 次数を誤ったり和分過程と判定されたりする可能性が変動する.

ということがわかる.

Rのコード

以下がシミュレーションに使用した Rのコードである. 使用しているパッケージが多いが, ggplot2, ggthemes, grid はグラフの体裁を整えるためだけで, compiler は実行速度を上げるためにある. rho や theta は真のモデルの係数で, e.sigma は誤差項 (独立ホワイトノイズ) の分散, size は最大のサンプルサイズで, times は反復回数である. ある. これらの数値をいじればほかの場合のシミュレーションができる. ただし, このままでも実行に数分かかるので, サンプルサイズを増やすとかなり時間がかかる可能性があるので注意.

#################################
# Asymptotic properties of 
# Information Criteria
#################################
# 本体: ver. 3.2.1
library(ggplot2)  # ver. 1.0.1
library(ggthemes) # ver. 2.1.2 
library(tseries)  # ver. 0.10-34
library(forecast) # ver. 6.1
library(tidyr)    # ver. 0.2.0
library(dplyr)    # ver. 0.4.2
library(grid)     # ver. 3.2.1
library(compiler) # ver. 3.2.1

# 乱数データ作成とARIMAパラメータ推定
arima.simulation <- cmpfun(function(rho, theta, size, times, e.sigma){
  sizes <-seq(from=10, to=size, by=10)
  res.order <- data.frame(rep=rep(1:times, each=length(sizes)), size=rep(sizes,times), 
                          ar.aic=NA, ma.aic=NA, ar.bic=NA, ma.bic=NA,
                          ar.aicc=NA, ma.aicc=NA,
                          aic=NA, bic=NA,aicc=NA,
                          d.aic=NA, d.bic=NA, d.aicc=NA)
  for (j in 1:times){
    x <- arima.sim(model = list(order=c(1,0,1), ar=rho, ma=theta), n=size, sd=e.sigma)
    #x <- arima.sim(model = list(ar=rho, ma=theta)), n=size, sd=e.sigma)
    for ( i in 1:length(sizes)){
      k <- i+(j-1)*length(sizes)
      res.order[k,"rep"] <- j
      temp.aic  <- auto.arima(x[1:sizes[i]], stepwise=T, max.p=10, max.q=10, max.d=2, ic="aic")
      temp.aicc <- auto.arima(x[1:sizes[i]], stepwise=T, max.p=10, max.q=10, max.d=2, ic="aicc")
      temp.bic  <- auto.arima(x[1:sizes[i]], stepwise=T, max.p=10, max.q=10, max.d=2, ic="bic")
      res.order$ar.aic[k]  <- temp.aic$arma[1]
      res.order$ma.aic[k]  <- temp.aic$arma[2]
      res.order$ar.bic[k]  <- temp.bic$arma[1]
      res.order$ma.bic[k]  <- temp.bic$arma[2]
      res.order$ar.aicc[k] <- temp.aicc$arma[1]
      res.order$ma.aicc[k] <- temp.aicc$arma[2]
      res.order$d.aic[k]   <- temp.aic$arma[6]
      res.order$d.bic[k]   <- temp.bic$arma[6]
      res.order$d.aicc[k]  <- temp.aicc$arma[6]
      res.order$aic[k]  <- temp.aic$aic
      res.order$bic[k]  <- temp.bic$bic
      res.order$aicc[k] <- temp.aicc$aicc
    }
  }
  return(res.order)
})

# ヒストグラム作成関数
draw.hist <- function(data, selection="AIC"){
  data <- data %>% select(size, starts_with("ar"), starts_with("ma"), starts_with("d") ) %>% gather(key = param, value=order, -size)
  if (selection=="AIC")      data <- data%>% filter(grepl(pattern ="aic$" , x=param) )
  else if(selection=="BIC")  data <- data%>% filter(grepl(pattern ="bic" , x=param) )
  else if(selection=="AICc") data <- data%>% filter(grepl(pattern ="aicc" , x=param))
  data$param <- with(data, factor(toupper(substr(param,start = 1, stop = 2)), levels = c("AR","MA","D.")))
  print(ggplot(data) + geom_histogram(aes(x=order), binwidth=1) + 
          labs(x="次数", y ="カウント", title=selection )+
          coord_cartesian(xlim = c(0,max(data$order)+1))+
          facet_grid(param~size) + theme_fivethirtyeight(base_size = 12) +
          theme(title=element_text(size=20),
                axis.title.y=element_text(angle=0, vjust=1), 
                #text=element_text(size = 8),
                panel.margin=unit(2,"char") )
        )
}

# -1< rho < 1 の範囲で設定すること
rho   <- .1 # AR
theta <- .9 # MA
size  <- 10^2
times <- 10^2
e.sigma <- 10

#実際の計算
output <- arima.simulation(rho = rho, theta = theta, size = size, times = times, e.sigma = e.sigma)
png(filename = paste0("rho",rho,".theta",theta,".aic.png"), width = 1000)
draw.hist(output)
dev.off()
png(filename = paste0("rho",rho,".theta",theta,".bic.png"), width = 1000)
draw.hist(output, "BIC")
dev.off()
png(filename = paste0("rho",rho,".theta",theta,".aicc.png"), width = 1000)
draw.hist(output, "AICc")
dev.off()