自己回帰モデル（AR モデル）は、ある時刻 t の値を、時刻 t よりも古いデータを使って回帰するモデルである。自己回帰モデルは、自己相関の高い時系列データのモデリングに有効である。回帰式は次のように表せる。φi は自己回帰係数、p は次数である。また、εt は誤差項であり、平均 0 かつ分散 σ2 のホワイトノイズである。
\[ y_{t} = \phi_{0} + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \cdots + \phi_{p}y_{t-p} + \epsilon_{t} \]自己回帰 AP(1) モデル
平均
AP(1) の平均値は、yt = φ0 + φ1yt-1 + εt により、
\[ E[y_{t}] = \phi_{0} + \phi_{1}E[y_{t-1}] + E[\epsilon_{t}] \]が成り立つ。また、自己回帰モデルを適用できるデータは定常性を持つ。定常性が仮定されるとき、E[yt] = μ かつ E[εt] = 0 が成り立つので、AP(1) の平均値は次のように計算できる。
\[ \begin{eqnarray} && \mu = \phi_{0} + \phi_{1}\mu \\ & \Longleftrightarrow & \mu = \frac{\phi_{0}}{1 - \phi_{1}} \end{eqnarray} \]自己分散
次に自己分散を計算するために、AP(1) のモデル式の両辺に yt をかけて、期待値を計算する。
\[ E[y_{t}^{2}] = \phi_{0}E[y_{t}] + \phi_{1}E[y_{t}y_{t-1}] + E[y_{t}\epsilon_{t}] \]自己回帰モデルのとき、データの定常性が仮定されることから、上式は次のように置き換えられる。
\[ \gamma(0) + \mu^{2} = \phi_{0}\mu + \phi_{1}(\gamma(1) + \mu^{2}) + \sigma^{2} \]μ = φ0/(1-φ1) を代入して、式を整理すると、
\[ \gamma(0) =\phi_{1}\gamma(1) + \sigma^{2} \]同様にして、AP(1) のモデル式の両辺に yt-1 をかけて、期待値を計算することで次のような等式を得ることができる。
\[ \gamma(1) = \phi_{1}\gamma(0) \]2 式より、AP(1) の時差 j における自己共分散は次のように計算される。
\[ \gamma(j) = \phi_{1}^{j}\gamma(0) = \frac{\phi_{1}^{j}\sigma^{2}}{1 - \phi_{1}^{2}} \]また、AP(1) の時差 j における自己相関は次のように計算される。
\[ \rho(j) = \frac{\gamma(j)}{\gamma(0)} =\phi_{1}^{j} \]シミュレーションデータ
次に、乱数を使って自己回帰モデルに基づくシミュレーションデータを作成し、その時系列変化とコレログラムを描いてみる。
AP(1) モデルのコレログラム
まず、AP(1) モデルのデータを作成してみる。モデルを以下の式のように構築する。
\[ y_{t} = 2 + 0.8y_{t-1} \]set.seed(2019)
y <- rep(0, 100)
y[1] <- 9
for (i in 2:length(y)) {
y[i] <- 2 + 0.8 * y[i - 1] + rnorm(1, mean = 0, sd = 1)
}
plot(y, type = 'l')
acf(y)
シミュレーションデータを生成するモデルは、「ある時点のデータはその 1 つ前の時点のデータから生成される」モデルとなっている。つまり、コレログラムでは、タイムラグ 1 のところで相関が高く、タイムラグ 2 以降では相関が低くなるはずである。実際のコレログラムをみると、ある時点のデータはその 1 つ前の時点のデータとの相関がほぼ 1 で、それよりも前の時点のデータとの相関が低いことが確認できる。過去に遡るほど、相関が見られなくなる。
また、コレログラムをみると、2 時点前との相関が 0.6 でありやや高いことが確認できる。これは、以下のように yt が yt-2 として書き表されることによる見せかけの相関によるものだと考えられる。
\[ y_{t} = 2 + 0.8y_{t-1} = 2 + 0.8 (2 + 0.8y_{y-2}) = 3.6 + 0.64y_{t-2}\]AP(2) モデルのコレログラム
AP(2) モデルのデータを作成してみる。モデルは以下の式で表せるものとする。
\[ y_{t} = 2 + 0.4y_{t-1} + 0.3y_{t-2}\]set.seed(2019)
y <- rep(0, 100)
y[1] <- 4
y[2] <- 8
for (i in 3:length(y)) {
y[i] <- 2 + 0.4 * y[i - 1] + 0.3 * y[i - 2] + rnorm(1, mean = 0, sd = 1)
}
plot(y, type = 'l')
acf(y)
この AR(2) モデルは、ある時点のデータとその 1 つ前の時点のデータとの相関が高いことが見られる。モデルでは 2 次の回帰式を使用してデータを生成しているが、ある時点とその 2 つ前の時点のデータとの相関が必ずしも高くないことがわかる。
ナイル川流量
R の標準パッケージには、ナイル川の流量の時系列データが保存されている。このナイル川の流量に対してコレログラムを描くと以下のようになる。
data(Nile)
plot(Nile, type = 'l')
acf(Nile)
R を使ったモデリング方法
AR モデルによる時系列データのモデリングは、R の ar 関数で行うことができる。ar 関数の戻り値にモデリング結果が保存され、計算された次数は order 要素、各係数は ar 要素に保存される。
set.seed(2019)
y <- rep(0, 100)
y[1] <- 9
for (i in 2:length(y)) {
y[i] <- 2 + 0.8 * y[i - 1] + rnorm(1, mean = 0, sd = 1)
}
obj <- ar(y)
obj
## Call:
## ar(x = y)
## Coefficients:
## 1
## 0.6304
## Order selected 1 sigma^2 estimated as 0.8022
set.seed(2019)
y <- rep(0, 100)
y[1] <- 4
y[2] <- 8
for (i in 3:length(y)) {
y[i] <- 2 + 0.4 * y[i - 1] + 0.3 * y[i - 2] + rnorm(1, mean = 0, sd = 1)
}
obj <- ar(y)
obj
## Call:
## ar(x = y)
## Coefficients:
## 1 2
## 0.2104 0.1631
## Order selected 2 sigma^2 estimated as 0.8949次に、AR(1) モデルに従うデータを 120 個生成して、最初の 100 個で AR モデルを構築し、残り 20 個のデータを予測するサンプルを示す。
set.seed(2019)
y <- rep(0, 120)
y[1] <- 9
for (i in 2:length(y)) {
y[i] <- 2 + 0.8 * y[i - 1] + rnorm(1, mean = 0, sd = 1)
}
obj <- ar(y[1:100])
obj
## Call:
## ar(x = y)
## Coefficients:
## 1
## 0.6304
## Order selected 1 sigma^2 estimated as 0.8022
y.predicted <- predict(obj, n.ahead = 20)
y.predicted
## $pred
## Time Series:
## Start = 101
## End = 120
## Frequency = 1
## [1] 8.869750 9.168805 9.357327 9.476169 9.551086 9.598312 9.628084 9.646851 9.658682 9.666140 9.670842 9.673806 9.675674 9.676852 9.677594 9.678062 9.678357 9.678543 9.678661 9.678735
## $se
## Time Series:
## Start = 101
## End = 120
## Frequency = 1
## [1] 0.8956729 1.0587870 1.1170132 1.1393259 1.1480724 1.1515298 1.1529008 1.1534452 1.1536615 1.1537474 1.1537815 1.1537951 1.1538005 1.1538026 1.1538035 1.1538038 1.1538040 1.1538040 1.1538040
## [20] 1.1538040
plot(1:120, y, type = 'none')
lines(1:120, y)
lines(101:120, y.predicted$pred, col = 'red')
lines(101:120, y.predicted$pred + 2*y.predicted$se, col = 'red', lty = 'dotted')
lines(101:120, y.predicted$pred - 2*y.predicted$se, col = 'red', lty = 'dotted')解析結果をグラフにしたものが次のようになる。黒線は乱数により作られたシミュレーションデータである。赤の実線は予測した平均値で、赤の点線は予測された平均値 ± 2×標準偏差である。予測結果を見ると、近い未来であれば、ある程度予測がうまく行っていることが確認できる。
