お気軽にAR(1)の構造転換点を区間推定してみる

学術的な意義は無いのですが、時系列データの構造転換点を求めたいときがあります。記述的に転換点を仮定して推定することが多いわけですが、主観的になりやすく言い合いになることがあるからです。統計学者に殺されそうな荒業ですが、尤度を使って構造転換点を探し、さらに区間推定してみましょう。
時系列データであれば何でも良いのですが、AR(1)を考えます。

1. モデル

$y_t = C + \beta y_{t-1} + \gamma (t \gt T) D + \epsilon$

$y_t$ は $t$ 期の観測値、 $C$ は切片項、 $T$ は推定する構造転換点、 $t \gt T$ はtがTより大のときに1、他は0をとる構造転換ダミー、 $\beta$ と $\gamma$ は推定する係数、 $\epsilon$ は誤差項です。
なお、AR(1)にダミーを入れていいのかは謎です。

2. データセット

モデルにそって、データセットを作成します。

set.seed(20150214)
n <- 100 # 観測数
y <- numeric(n) # 観測値
C <- 3 # 切片項
gamma <- 4 # 構造転換ダミー
T <- round(n/2) # 構造転換点
beta <- 0.8 # β
epsilon <- rnorm(n) # ε
y[1] <- 15 + epsilon[1]
for(i in 2:n){
y[i] <- C + beta * y[i-1] + epsilon[i] + (i>T)*gamma
}
# 生成データを確認する
plot(y, type="l", main="データセット", xlab="t", ylab="y")

3. 構造転換点が分かっている場合

構造転換点Tが固定で分かっているのであれば、簡単に推定できます。

d <- (2:n>50)*1 # tが50以下は0、51以上は1になる
r <- lm(y[-1] ~ y[-n] + d)
summary(r)

しかし、実際にはT=50は分かりません。40〜60のどれを仮定しても、それらしい結果になります。

4. 対数尤度を最大にする点を探す

Tを変えて説明変数の数は増えたりしないので、対数尤度を最大化する点を探しましょう。対数尤度関数を書いてニュートン・ラフソン法などで探したくなりますが、構造転換点は不連続なため上手く推定できません。ループして力技で最大値を求める方が確実です。

logLik <- numeric(n) # 対数尤度の保存用ベクトル
logLik[1] <- logLik[n] <- NA
ml_scp <- 0 # 尤度が最大の点
r <- lm(y[-1] ~ y[-n]) # 尤度が最大の推定結果
for(scp in 2:(n-1)){
  d <- (2:n>scp)*1
  tmp <- lm(y[-1] ~ y[-n] + d)
  logLik[scp] <- logLik(tmp) # 対数尤度は蓄えておく
  if(logLik(r) < logLik(tmp)){
    r <- tmp
    ml_scp <- scp
  }
}
summary(r)
sprintf("最尤推定された構造転換点: %d", ml_scp)

5. 構造転換点の区間推定を行なう

構造転換点が分かれば十分なときも多いわけですが、信頼性を疑われるので区間推定が行なえるように標準誤差を求めたいです。しかし、対数尤度の集合は離散データなので、このままでは最尤推定値の標準誤差は求まりません。

そこでスプラインで補間して、連続な対数尤度関数モドキを作ってしまいましょう。

t <- 2:(n-1)
sp <- smooth.spline(t, logLik[t])
f <- function(p){
predict(sp, p)$y
}

プロットするとそこそこの精度で近似できていることが分かります。全体としては丸くなるので、区間推定の範囲は広くなりそうですが、狭くなるよりは良いでしょう。

plot(logLik, type="p", main="対数尤度", xlab="t", ylab="y")
lines(f(1:n))

作った対数尤度関数モドキから、構造転換点を再推定します。

r_ml <- nlm(function(p){
-f(p)
}, ml_scp, hessian=TRUE)
SEs <- sqrt(diag(solve(r_ml$hessian))) # 標準誤差

これで構造転換点の標準誤差が出ました。区間推定もしてみましょう。

interval <- function(beta,se,range,nof){
a <- 1 - range/2
sprintf("%.3f(%d%%信頼区間%.3f〜%.3f)",beta,as.integer((1-range)*100),beta-se*qt(a, nof),beta+se*qt(a, nof))
}
interval(r_ml$estimate[1], SEs[1], 0.05, summary(r)$df[2])

厳密な方法ではありませんが、tは49から51ぐらいを見とけば良いとなります。