概要

前回大数の法則の視覚化から理想の推定量を考える - ill-identified diary の最後に上げた具体例の, 時系列分析の場合についても, 推定量の違いから生じる結果を視覚化してみた.
時系列はあまり詳しくないので操作変数編より内容が薄い.
安定な自己回帰 (AR) モデルと, 自己回帰移動平均 (ARMA) モデルの場合のみ.
時系列分析の話なのでそれなりに数式が出てくる. AR・ARMAモデルの基本的な話は, 日本語なら沖本 (2010, 経済・ファイナンスデータの計量時系列分析) とか Rで計量時系列分析：AR, MA, ARMA, ARIMAモデル, 予測 - 六本木で働くデータサイエンティストのブログとかを参考に. もう少し詳細厳密な話が知りたい場合は, Hamilton (2006, 時系列解析 (上) 定常過程編), Hayashi (2000, Econometrics) とかを.
文章量は PDF換算で 9 ページ (やはり大半が画像とサンプルコード)

自己回帰モデル

まずは, 通常最小二乗法 (OLS) で一致性のみが保証されるケースとして, 自己回帰モデルを簡単に説明する.

$\begin{align} y_{t}= & \alpha+\rho y_{t-1}+\varepsilon_{t}\label{eq:ar1}\tag{1}\end{align}$

のように, 被説明変数の過去の値が説明変数となるようなものを自己回帰 (Auto-Regressive) モデルという. 1期前の値を参照するので, これは 1次の自己回帰モデルであり, AR (1) と表現する. 説明変数に $y_{t-2}$ , $y_{t-3}$ , というふうに過去の値が連なるような, AR (p) モデルというのも考えられるが, ここでは簡単のため, 1次の場合だけを扱う. AR (1) の係数 $\rho$ がゼロ以上1未満なら, AR (1) は安定的 (または定常的) と呼ばれる*1. 安定的でない場合は, 例えば「単位根過程」とか「和分過程」とか呼ばれる時系列データになるのだが, これの説明まで含めると長いので, 今回は安定的な　AR (1) についてのみ考える.

ところで, 時系列というのは, 日々の積み重ねであり, 今日の誤差項 $\varepsilon_{t}$ と昨日の誤差項 $\varepsilon_{t-1}$ が互いに独立しているとは考えにくい. そこで, 時系列においては, $\varepsilon_{t}$ は互いに独立ではなく, ホワイトノイズと呼ばれる確率変数の数列, 確率過程として扱う. といっても, 平均がゼロ, 分散が $\mathrm{V}(\varepsilon_{t})=\sigma^{2}$ で, 自分の過去の値との共分散 (自己共分散) が, どの過去の値と参照してもゼロ $\mathrm{Cov}(\varepsilon_{t},\,\varepsilon_{t-k})=0,\,\forall k=1,\,2,\,\cdots$ というものなので, これまでの仮定とそう大して変わらない. なお, ホワイトノイズでも, 特に異なる時点で互いに独立なものを独立ホワイトノイズというが, 誤差項が独立ホワイトノイズなら, もう今までの仮定と全く変わらない*2.

しかし, 大数の法則 (と中心極限定理) は, それぞれの確率分布が互いに独立である, という前提条件があることを覚えている人は, 大数の法則が適用できないのではないか, と考えるかもしれない. 実は大数の法則はこの前提条件を緩和しても成り立つ*3ので, 今回のケースではこの問題をそれほど気にする必要はない. しかし, 次のような場合は要注意である.

誤差項が自己相関 (系列相関) する: つまり, (1) と同様に, $\varepsilon_{t}=\phi\varepsilon_{t-1}+u_{t}$ のような関係が成り立っている
誤差項 $\varepsilon_{t}$ と説明変数 $y_{t-1}$ が相関する

1. のケースは, これはいわゆる不均一分散なので, OLS が有効性を失う原因になる. なお, これは自己回帰モデルでなくても ( $y_{t}=\alpha+\beta x_{t}+\varepsilon_{t}$ のような) 時系列分析全般で起こりうる*4. 次に 2. は, 操作変数の場合と同じで OLS の一致性そのものを失うことになる. さらに, 特に時系列分析の場合は, 1. に誘発するかたちで 2. の問題も発生することがある. つまり, 自己相関が存在し,

$\begin{align} \varepsilon_{t}= & \phi\varepsilon_{t-1}+u_{t}\label{eq:ma}\tag{2}\end{align}$

とする. $u_{t}$ は, 今度こそ他の一切の変数と独立で, 自身の過去の値とも独立なホワイトノイズと考える (これで他の変数との相互関係を考えなくて済む).

だとすると, $\varepsilon_{t}$ は過去の $\varepsilon_{t-1}$ に影響して決まっていることになる. 次に, (1) を1期前に戻すと,

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

となる. よって, $y_{t-1}$ が $\varepsilon_{t-1}$ に影響して決まっているのは明らかなので, 結局, 誤差項の自己相関を仮定するだけで, $\varepsilon_{t-1}$ が $y_{t-1}$ と相関するという結論が導ける.

つまり, 誤差項の系列相関がある場合, OLS は一致性からして保証されない. よって, 重み付き最小二乗法 (WLS) やコクラン=オーカット法などもこの場合は使えない*5.

しかし, まずは, (1) の誤差項 $\varepsilon_{t}$ が, 自己相関のない, ホワイトノイズであると仮定する. このとき, OLS は一致性のみ保証される.

自己回帰移動平均モデル

実際には, 自己回帰モデルを仮定して, 残った誤差項 $\varepsilon_{t}$ がホワイトノイズの条件を満たすという場合は少ない. 自己回帰に加えて, 過去の誤差項の移動平均で表現される

$\begin{align} y_{t}= & \alpha+\rho y_{t-1}+\phi\varepsilon_{t-1}+\varepsilon_{t}\label{eq:arma}\tag{3}\end{align}$

ような時系列モデルは, (1, 1)-次の自己回帰移動平均モデルといい, ARMA (1, 1) と表す. これは先に説明した (1) の誤差項が自己相関するような場合と似ている. (2) の $u_{t}$ が $\varepsilon_{t}$ に置き換わっているだけが異なる点である. 誤差項の値は観測できないので, 実際にパラメータを推定する場合は, 移動平均 $\phi\varepsilon_{t-1}+\varepsilon_{t}$ の部分をひとまとめにして1つの誤差項として扱うことになる (このような1期前の誤差項で表される時系列モデルを, 1次の移動平均 (Moving Average; MA) モデルといい, MA(1) と表す. AR(1) と MA(1) の合成だから, (3) は ARMA(1, 1) である). すると, さきほどに説明したように (3) は誤差項と説明変数が相関する. ここから, ARMA モデルの係数 $\rho$ を推定するとき, OLS では, 一致性・不偏性・有効性はどれも保証されないことがわかる.

誤差項と説明変数の相関の問題なので, ARMA モデルでも操作変数法で対処することは可能である. $y_{t}$ が安定な ARMA (1, 1) モデルならば, 無限次の移動平均であり, つまり $y_{t}$ は誤差項 $\varepsilon_{t}$ の過去のすべての値の影響を受ける. これを反転可能性という (反転可能性については詳しくは先に上げた参考文献を参照). つまり,

$\begin{align} y_{t}= & \varepsilon_{t}+\theta_{1}\varepsilon_{t-1}+\theta_{2}\varepsilon_{t-2}+\cdots\end{align}$

のように, $y_{t}$ は1期前の誤差項 $\varepsilon_{t-1}$ とそれ以前のすべての誤差項に影響される.

一方で, $\varepsilon_{t}$ はホワイトノイズであり, 自らの過去の値とは相関しないので, $k\geq1$ に対して, $\mathrm{E}\varepsilon_{t}\varepsilon_{t-k}=0$ である. よって, $y_{t-1}$ と $\varepsilon_{t}$ の共分散は (ホワイトノイズの平均はゼロなので, 交差モーメント $\mathrm{E}y_{t-1}\varepsilon_{t}$ が共分散に等しくなるから), 誤差項の部分を $e_{t}=\varepsilon_{t}+\phi\varepsilon_{t-1}$ とおくと,

$\begin{align} \mathrm{E}y_{t-1}e_{t}= & \mathrm{E}(\varepsilon_{t-1}+\theta_{1}\varepsilon_{t-2}+\theta_{2}\varepsilon_{t-3}+\cdots)(\varepsilon_{t}+\phi\varepsilon_{t-1})\\ = & \mathrm{E}\phi\varepsilon_{t-1}^{2}=\phi\sigma^{2}>0\end{align}$

ということで, 積を展開すると, $t-1$ 期の誤差項の二乗が残るので, $y_{t}$ の説明変数である $y_{t-1}$ が誤差項と相関していることが再確認された. 次に, $y_{t-2}$ と $\varepsilon_{t}$ の共分散を見ると,

$\begin{align} \mathrm{E}y_{t-2}e_{t}= & \mathrm{E}(\varepsilon_{t-2}+\theta_{1}\varepsilon_{t-3}+\theta_{2}\varepsilon_{t-4}+\cdots)(\varepsilon_{t}+\phi\varepsilon_{t-1})\\ = & 0\end{align}$

となり, 展開すると異なる時期どうしの誤差項の積しかないので, 相関がゼロであることがわかる. 同様に考えれば $y_{t-1}$ と $y_{t-2}$ が相関していることは明らかだから, これを GMM のモーメント条件とみなせば, $y_{t-2}$ を $y_{t-1}$ の操作変数として扱える. さらに, $y_{t-3},\, y_{t-4},\,\cdots$ も同様に誤差項と無相関なはずなので, これらも理論上は操作変数として使える. しかし, 一般に安定的な時系列モデルは, 過去の値ほど現在の値との相関が弱くなるので, 過剰識別の問題が発生するから, 操作変数の数は抑えておくのが良い (Hamilton (2006, 時系列解析 (上) 定常過程編) には適切な操作変数の選び方が書かれているらしいが, 現在所持してないので確認できない). 今回は $y_{t-2}$ のみを操作変数として, 2段階最小二乗法で ARMA モデルのパラメータの一致推定量を求める. または, 最尤法でもいい.

コードと実行結果

まずは下準備. ここは前回までとほとんど同じなので特に解説しない.

# 本体 ver. 3.2.0
library(ggplot2)  # ver. 1.0.1
library(ggthemes) # ver. 2.1.2 
library(reshape2) # ver. 1.4.1
library(MASS)     # ver. 7.3-39 
library(AER)      # ver. 1.2-3 
library(tseries)  # ver. 0.10-34
size     <- 10^3 # 最大サンプルサイズ
size.min <- 10   # 最小値
ser      <- 10   # 系列の数
# rho は |rho| < 1 の範囲にすること
rho      <- .5   # X の自己相関 (AR) 
e.sigma  <- 5    # ホワイトノイズの分散
phi      <- .4   # MAの係数

# グラフ描画関数
ggLLN.ts <- function(data, width=1){
  ggplot(data) + geom_path(aes(x=N, y=value, group=variable, linetype=series, color=type), size=.7) + scale_linetype_discrete(guide=F) + 
    geom_hline(yintercept = rho, linetype="dashed", color=5, size=1.2) +
    labs(x="N (対数)", y=expression(beta) ) +
     theme(axis.title.y =element_text(angle = 0), legend.title=element_blank()) +
    scale_color_colorblind() +
    scale_x_log10() + coord_cartesian(xlim=c(size.min,1000), ylim=c(rho-width/2, rho+width/2))
}

R で時系列データを扱う場合, 少なくとも tseries パッケージの各オブジェクトの扱い方を知る必要がある*6.

まず, arima.sim() は, AR や ARMA モデルにしたがう時系列データの乱数列を作成できるので, 今回はこれを使って作成している. この関数で作成されるのは, 行列でもデータフレームでもなく, ts という時系列データ用のオブジェクトである. 次に, 説明変数の値のベクトルとして, $y_{t}$ (コードでは x というオブジェクトを使っているのに注意) の1期前の値を用意する必要がある. lag(x,k) 関数はラグを取ってくれるのだが, この関数は ts オブジェクトにしか効果がなく, 通常のデータフレームやベクトルに適用しても, エラーも警告もなく元の値がそのまま返されるので注意が必要である. 最小二乗法をする lm() 関数は, データとして ts オブジェクトを与えることができるが, 関数内で行列に変換するため,

lm(x~lag(x, -1), x)

と書くのは誤りである (x で x を回帰するという全く無意味なことになる). lm() 関数で自己回帰モデルのパラメータを計算するには, まず ts.union() 関数を使い, 元の時系列のベクトルと, ラグをとった時系列のベクトルを並べた多変量時系列データのオブジェクトを作成する. 下記コードでは,

df <- ts.union(x,x.l=lag(x,-1))

と書いているので, ラグを取った変数は x.l という名前になる.

こうすれば lm() 関数に入れても正しく認識される. しかし, 実際には,同じく tseries パッケージで提供される ar() 関数や arima() 関数で計算したほうが楽である. いずれも method= 引数で計算方法を指定できるが, arima() の方は AR(I)MA モデル*7一般を推定するための関数なので, 最小二乗法は使えない. (1) のモデルを当てはめて OLS で計算したいなら,

ar(x, order=1, method="ols")

とする. ar() なら method=”mle” に, arima() なら method=”ML” にすれば最尤法で計算する.

また, ts オブジェクトから一部の期間だけを抜き出すときに, 行列のように x[i:j] とすると, ts オブジェクトが行列に変換されたものが返されてしまう. このような場合は, window() 関数で取り出したほうが安全である.

##### ts オブジェクトの扱い方 #####
# AR(1)
x <- arima.sim(list(ar=rho), n=size, sd=e.sigma)
df   <- ts.union(x,x.l=lag(x,-1))
ols  <- lm(x~x.l,df)
# arima() か ar() を使うほうが楽
arima(x, order=c(1,0,0), method="ML") # arima では最小二乗法不可
ar(x,order=1, method = "ols") # ar なら ols 可能
ar(window(x,1, 100), order=1, method="ols") # ts から一部を抜き出すときは window() を使う
# ols  <- lm(x~lag(x, -1), x) # 内部でオブジェクト変換されるのでこれは間違い

実際には, 時系列データの分析で AR・ARMA の次数があらかじめ分かっているということは少ないので, 情報量規準を元に変数選択したり, Box-Pierce 検定量を使うなどして次数を決める (詳しくは沖本 (2010, 経済・ファイナンスデータの計量時系列分析) や Rで計量時系列分析：AR, MA, ARMA, ARIMAモデル, 予測 - 六本木で働くデータサイエンティストのブログ, あるいは Hayashi (2000, Econometrics) などを参照). ar(), arima() では, order= を省略すると情報量規準で自動的に次数を決定してくれるが, 今回は正確な次数があらかじめ分かっているものとして計算する*8.

以上を踏まえて, 実際に計算したコードと結果が以下である. なお, AR, ARMA 両方の計算で 30分程度かかる. 時間がもったいない場合, 適宜 ser や size を減らすとよい (とくに size <- 500 くらいにするとかなり速くなる).

##### AR(1) の乱数生成と推定 ######
set.seed(1)
x <- list()
res.ar1.ols <- matrix(nrow=size, ncol=ser)
res.ar1.ml <-  matrix(nrow=size, ncol=ser)
for (i in 1:ser){
  # 乱数生成
  x[[i]]  <- arima.sim(list(ar=rho), n=size, sd=e.sigma)
  for (j in size.min:size){
    # 推定
    print(paste("series:", i, "in", ser, ", size:", j, "/", size))
    # OLS
    tryCatch(expr = {
        res.ar1.ols[j,i] <- ar(window(x[[i]], 1, j), order=1, method = "ols")$ar
      },
      error=function(x) NA )
    # ARMA として MLで推定
    tryCatch(expr = {
        res.ar1.ml[j,i]  <- ar(window(x[[i]], 1, j) , order=1, method = "mle")$ar
      },
      error = function(x) NA)
  }  
}

# 結果を ggplot 用に加工
colnames(res.ar1.ols) <- paste("OLS", 1:ser, sep="")
colnames(res.ar1.ml) <- paste("ML", 1:ser, sep="")
res.all <- cbind(res.ar1.ols, res.ar1.ml)
res.all <- data.frame(N =  1:size, res.all)
res.melt.ar1 <- melt(res.all, measure.vars = colnames(res.all)[-1])
res.melt.ar1$type <- "OLS"
res.melt.ar1$type[grep("ML", res.melt.ar1$variable)] <- "ML"
res.melt.ar1$series <- sub(x= res.melt.ar1$variable, pattern = ".+([0-9]+)$", replacement= "\\1")

ggLLN.ts(data = res.melt.ar1)

図1: AR の推定

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

AR (1) の場合, 予想通り OLS と最尤法 (ML) の結果はいずれも真の値 $\rho=0.5$ に収束している. 分散も似たようなものだろう. 本題は ARMA (1, 1) の結果である.

#### ARMA の推定 ####
# OLS, 2SLS, ML で推定して比較
set.seed(1)
x.arma <- list()
res.arma11.ols <- matrix(nrow=size, ncol=ser)
res.arma11.ml <-  matrix(nrow=size, ncol=ser)
res.arma11.iv <-  matrix(nrow=size, ncol=ser)
for (i in 1:ser){
  # ARMA(1,1) 乱数生成
  x.arma[[i]]  <- arima.sim(list(ar=rho, ma=phi), n=size, sd=e.sigma)
  # OLS/IV 用に data.frame で作成
  x.df <- ts.union(x=x.arma[[i]], x.l = lag(x.arma[[i]],-1), x.l2 = lag(x.arma[[i]], -2), dframe = T )

  for (j in size.min:size){
    # 推定
    print(paste("series:", i, "in", ser, ", size:", j, "/", size))  
    
    # OLS
    tryCatch(expr = {
      res.arma11.ols[j,i] <- ar(window(x.arma[[i]], 1, j), order=1)$ar },
    error=function(x) NA )
    
    # ML
    tryCatch(expr = {
      res.arma11.ml[j,i]  <- arima(window(x.arma[[i]], 1, j) , order=c(1,0,1), method = "ML")$coef[1]
    },
    error = function(x) NA)
    # IV
    tryCatch(expr={
      res.arma11.iv[j,i] <- ivreg(formula = x~x.l|x.l2, data = x.df[1:j,])$coefficients[2]},
      error=function(x) NA)
    
  }
}

# 結果を ggplot 用に加工
colnames(res.arma11.ols) <- paste("OLS", 1:ser, sep="")
colnames(res.arma11.ml)  <- paste("ML", 1:ser, sep="")
colnames(res.arma11.iv)  <- paste("IV", 1:ser, sep="")
res.all <- cbind(res.arma11.ols, res.arma11.ml, res.arma11.iv)
res.all <- data.frame(N =  1:size, res.all)
res.melt.arma11 <- melt(res.all, measure.vars = colnames(res.all)[-1])
res.melt.arma11$type <- "OLS"
res.melt.arma11$type[grep("ML", res.melt.arma11$variable)] <- "ML"
res.melt.arma11$type[grep("IV", res.melt.arma11$variable)] <- "IV"
res.melt.arma11$series <- sub(x= res.melt.arma11$variable, pattern = ".+([0-9]+)$", replacement= "\\1")

ggLLN.ts(data = res.melt.arma11)

図2: ARMA の推定結果

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

理論通り, OLS は全く異なる値に収束しているため, 一致推定量ではないとわかる. 一方, ML と 2SLS (IV) は一致推定量だが, 分散が明らかに AR のときより大きくなっている (サンプルサイズが小さいときは, 若干 ML のほうが分散が小さい?).

参考文献

Farnworth, G. V. (2008). Econometrics in R. Retrieved from http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf

Hamilton, J. D. (2006). 時系列解析 (上) 定常過程編. シー・エー・ピー出版. 沖本竜義・井上智夫訳, 原題: ``Time Series Analysis''

時系列解析〈上〉定常過程編

作者: J.D.ハミルトン,James D. Hamilton,沖本竜義,井上智夫
出版社/メーカー: シーエーピー出版
発売日: 2006/03
メディア: 単行本
クリック: 59回
この商品を含むブログ (2件) を見る

Time Series Analysis

作者: James D. Hamilton
出版社/メーカー: Princeton Univ Pr
発売日: 1994/01/11
メディア: ハードカバー
購入: 1人クリック: 5回
この商品を含むブログ (8件) を見る

Hayashi, F. (2000). Econometrics. Princeton University Press.

Econometrics

作者: Fumio Hayashi
出版社/メーカー: Princeton University Press
発売日: 2011/12/12
メディア: Kindle版
この商品を含むブログを見る

沖本竜義 (2010). 経済・ファイナンスデータの計量時系列分析. 朝倉書店.

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

作者: 沖本竜義
出版社/メーカー: 朝倉書店
発売日: 2010/02/01
メディア: 単行本
購入: 4人クリック: 101回
この商品を含むブログ (6件) を見る

*1:沖本本では「平均回帰的」とも表現される. つまり時間が進む毎に一定の値 (平均値) に収束するような時系列データなので, 「安定的」なのである.

*2:一応注意しておくが, 独立ならば相関はゼロだが, 相関がゼロだからといって独立とは限らない. ホワイトノイズと今までの仮定との違いはここである. 「相関がゼロでも独立でない」場合をイメージできないなら, "相関"の話＆そのついでに"21世紀の相関（MIC）"の話（ややマニア向け） - Take a Risk：林岳彦の研究メモなどが参考になる.

*3:ここでは証明は書かない.

*4:むしろ最初から疑ってかかったほうが良いくらい.

*5:蛇足だが, 自己回帰モデルにおいて系列相関があるかの仮説検定としてダービン・ワトソン検定は使えない.

*6:これに限らず, 計量経済学でよく使う計算を R でやるときの解説として, Farnworth (2008, Econometrics in R) などが参考になる.

*7:ARIMA とは ARMA モデルの拡張版であるが, ここでは詳しく知らなくても問題ない.

*8:ARIMA とは ARMA モデルを拡張したもので, ARMA が2つのパラメータ (p, q) だったのに対し, (p, d, q) という3つのパラメータで表現されるが, d=0 のときは ARMA (p, q) と同じである. よって, ARMA (1, 1) として推定したいなら order=c(1, 0, 1) と書けば良い.