AR(1)モデルのメモ

R 統計科学

今更感あるが、諸事情によりメモ。やりたいことは(定常性を仮定した)AR(1)モデルの自己相関係数がちゃんと理論値と一致するのかを見たいって話。

AR(1)モデルとは

AR(1)は、適当な確率過程 x=\{x_t; t\ge 1 \}が、
x_t=d+\rho x_{t-1} + \eta_t
に従うようなモデル。ここで、

  • \eta_tは独立な平均0、標準偏差\sigmaの正規分布、dは適当な定数
  • \{x_t; t\ge 1 \}に(弱)定常性を仮定(AR(1)に関しては|\rho| < 1が条件)

としておく。この時、その平均値 E(x_t)と分散Var(x_t)は時間tに依らずに一定となりそれぞれ

  • E(x_t)=d/(1-\rho)
  • Var(x_t) = \sigma^2/(1-\rho^2)

となる。

標本自己相関関数 \gamma_hは、標本自己共分散C_h

  • C_h=\frac{1}{T}\sum^{T}_{t=h+1}\left(x_t-\mu\right)\left(x_{t-h}-\mu\right)
  • \mu=\frac{1}{T}\sum_{t=1}^{T} x_t (標本平均)

を用いて

  • \gamma_h = \frac{C_h}{C_0

で定義される。特に、AR(1)モデルの場合、 \gamma_hの関数系をちゃんと計算することが出来て、

  • \gamma_h = \rho^{h}

になる。ちゃんとこうなるのかをシミュレーションを通して確かめることが目的。

Rでやりたい

arima.sim関数を使ってもいいが、AR(1)のシミュレーションぐらいならサクっと

AR1 <- function(size, d, rho, sigma)
{
  Reduce(function(x, noise){d+rho*x+noise}, rnorm(size, mean=0, sd=sigma), init=d/(1-rho), accumulate=TRUE)
}

こんな感じで関数を用意できるので、こいつを使うことにする。
まずは、適当なパラメーターに対して100万個のデータxを生成。この時x_1に相当するもの(init引数)は、平均値になるようにAR1関数内で調整している。

rho <- 0.8
sigma <- 0.5
d <- 1
x <- AR1(10^6, d, rho, sigma)

生成したデータxの標本平均が理論値と合うかのチェック。ほぼ合ってますよと。

> mean(x)
[1] 5.001267
> d/(1-rho)
[1] 5

同じくデータの標本分散が理論値と合うかチェック。こちらもほぼ合ってますよと。

> var(x)
[1] 0.6930294
> sigma^2/(1-rho^2)
[1] 0.6944444

最大のラグを30とした場合の自己相関を計算し、合わせて描画も行う。ちなみにacf関数のplot引数をfalseと設定すると描画はされない。

lag <- 30
x.acf <- acf(x, lag.max = lag)$acf

これが理論的な自己相関と一致するかを描画して比較してみる。結果はたしかにそれっぽくは見える。

matplot(cbind(x.acf, rho^(0:lag)), type="l")

こういう指数関数ぽい形になっている関数はlogを取って線形に変換( h\log(\rho))した方が比較しやすいので、logを取ってもう一度描画。
ラグ30あたりの端以外はまぁいい感じか。

matplot(cbind(log(x.acf), (0:lag)*log(rho)), type="l")


参考