計量時系列分析
いわゆる龍義本。リンクはこちらから。
pdfも無料で落とせるらしい。
1. 時系列分析の基礎概念
1.1 時系列分析の基礎
1.1.1 時系列分析の目的
時系列データ解析の主要な目的は、時系列データに関しての何らかの予測である。また、変数間の動的な関係を明らかにすることも重要である。
1.1.2 時系列データの種類
データそのものを原系列という。
1.1.3 基本統計量と時系列モデル
特にファイナンスの世界では標準偏差のことをボラティリティvolatilityということがある。時系列分析で特有の統計量に自己共分散autocovarianceがある。これは同一の時系列データにおける異時点間の共分散であり、以下のように定義される。$${\gamma}_{kt} = \text{Cov}(y_t, y_{t-k}) = E[(y_t-\mu_t)(y_{t-k}-\mu_{t-k})] \tag{1.2}\label{1.2}$$ここで、\(\mu_{t} = E(y_t)\)である。自己共分散を\(k\)の関数とみたものを自己共分散関数という。これは正定値になることが知られている。
この定義の問題としては、単位依存性である。そこで、以下の自己相関係数autocorrelation coefficientを考える。$$\rho_{kt} = \text{Corr}(y_t, y_{t-k}) = \frac{\text{Cov}(y_t, y_{t-k})}{\sqrt{\text{Var}(y_t)\text{Var}(y_{t-k})}} = \frac{\gamma_{kt}}{\sqrt{\gamma_{0t}\gamma_{0, t-k}}} \tag{1.3}\label{1.3}$$この定義から\(\rho_{0t} = 1\)であることは明らかである。また、\(k\geq 1\)において\(|\rho_{kt}|\leq 1\)も成立する。これを\(k\)の関数とみたものは自己相関関数と呼ばれる。\(x\)軸に\(k\)、\(y\)軸に自己相関関数を描いたものはコレログラムcorrelogramと言われる。
時系列分析の大きな問題としては、データの観測は一度しかできず、また単一のデータから分散や自己共分散を計算することができないことが挙げられる。このため、時系列データは確率変数\(\{y_t\}_{t=-\infty}^{\infty}\)からの一つの実現値\(\{y_t\}_{t=1}^{T}\)とみなし、確率変数列の生成過程において何らかの性質や構造を仮定する。この確率変数列は確率過程stochastic processあるいはデータ生成過程DGP: data generating processと呼ばれる。
1.2 定常性
定常性stationarityについて。定常性は弱定常性weak stationarityと強定常性strict stationarityに分類される。
定義1.1（弱定常性）任意の\(t, k\)に対して$$\begin{eqnarray}E(y_t) & =& \mu \tag{1.4}\label{1.4}\\ \text{Cov}(y_t, t_{t-k}) & = & E[(y_t-\mu)(y_{t-k}-\mu)] = {\gamma}_k \tag{1.5}\label{1.5}\end{eqnarray}$$が成立する場合、過程は弱定常(weak stationary)といわれる。
\eqref{1.5}から弱定常過程においては自己共分散は時点には依存せず、時間差\(k\)のみに依存する。したがって、任意の\(k\)に対して\(\gamma_{k} = \gamma_{-k}\)が成立する。また、このとき$$\begin{eqnarray}\text{Corr}(y_t, y_{t-k}) & = & \frac{\gamma_{kt}}{\sqrt{\gamma_{0t}\gamma_{0,t-k}}}\\ & = & \frac{\gamma_{k}}{\gamma_0}\\ & = & \rho_{k}\end{eqnarray}$$であり、自己相関も時点に依存しない。さらに\(\rho_{k} = \rho_{-k}\)も成立する。弱定常性は共分散定常性(covariance stationarity)と呼ばれることもある。これに対して、強定常性は同時分布が不変であることを要求する。
定義1.2（強定常性） 任意の\(t, k\)に対して、\((y_t, y_{t+1}, \cdots, y_{t+k})^{\prime}\)の同時分布が同一となる場合、過程は強定常(strict stationary)といわれる。
1.3 ホワイトノイズ
定義1.3 （iid系列）各時点のデータが互いに独立でかつ同一の分布に従う系列はiid系列(iid sequence)と呼ばれる。
iidはindependently and identically distributedの略である。\(y_t\)が期待値\(\mu\)、分散\({\sigma}^2\)（標準偏差\(\sigma\)）のiid系列であることを\(y_t\sim iid(\mu, {\sigma}^2)\)と表す。このiid系列自体を時系列モデルとするより、期待値\(0\)のiid系列を時系列モデルの撹乱項(innovation, disturbance term)として用いる。より弱い仮定をおいたものが、以下のホワイトノイズになる。
定義1.4（ホワイトノイズ） すべての時点\(t\)において$$\begin{eqnarray}E({\epsilon}_t) & = & 0 \tag{1.6}\label{1.6}\\ \gamma_k & = & E({\epsilon}_t{\epsilon}_{t-k}) = \begin{cases}{\sigma}^2, k = 0\\ 0, k\ne 0\end{cases} \tag{1.7}\label{1.7}\end{eqnarray}$$が成立するとき、\({\epsilon}_t\)はホワイトノイズ(white noise)と呼ばれる。
以下\({\epsilon}_t\)が分散\({\sigma}^2\)のホワイトノイズであることを\({\epsilon}_t \sim \text{W.N}({\sigma}^2)\)と表記する。この定義から、ホワイトノイズはすべての時点において期待値が\(0\)で、かつ分散が一定であり、自己相関を持たない。これからホワイトノイズが弱定常過程であることは明らかである。最も基本的な弱定常過程はホワイトノイズを用いて$$y_t = \mu + {\epsilon}_t, \ \ {\epsilon}_t \sim \text{W.N}({\sigma}^2) \tag{1.8}\label{1.8}$$と表される。
様々な\(\mu, \sigma\)に対して系列をplotしたものが以下になる。Python3のコードと合わせて記載しておく。
import numpy as np
import matplotlib.pyplot as plt
mu = [0, 2, -2]
std = [1, 2, 3]
dlen = 124
fig, ax = plt.subplots(3, 3, sharex = True, sharey = False, figsize = (10, 10))
for i in range(3):
for j in range(3):
mu_temp = mu[i]
std_temp = std[j]
y = np.random.normal(mu_temp, std_temp, dlen)
ax[i, j].plot(y)
ax[i, j].plot(0, 0, label = "$\mu$ = {:3.0f}\n$\ \sigma$ = {:3.0f}".format(mu_temp, std_temp))
ax[i, j].legend(fontsize = 8)
plt.tight_layout()1.4 自己相関の検定
時系列データに対して期待値、自己共分散、自己相関の推定量は以下で与えられる。$$\begin{eqnarray}\bar{y} & = & \frac{1}{T}\sum_{t=1}^{T}{y_t}\\ {\hat{{\gamma}_k}} & = & \frac{1}{T}\sum_{t = k+1}^{T}{(y_t-\bar{y})(y_{t-k}-\bar{y})}, k = 0, 1, 2, \cdots \\ \bar{{\rho}_k} & = & \frac{\hat{\gamma}_k}{\hat{\gamma}_0}, k = 1, 2, 3, \cdots \tag{1.10}\label{1.10} \end{eqnarray}$$この自己相関係数\(\hat{\rho}_k\)を用いて\(H_{0}:\hat{\rho}_k = 0\)という帰無仮説を\(H_1: \hat{\rho}_k\ne 0\)という対立仮説に対して検定する。\(y_t\)がiid系列の場合には\(\hat{\rho}_k\)は漸近的に平均\(0\)、分散\(\displaystyle \frac{1}{T}\)の正規分布に従うことが知られている。複数の自己相関係数がすべて\(0\)であるという帰無仮説を検定するときは、かばん検定portmanteau testといわれ、Ljung and Box (1978)の統計量が知られている。$$Q(m) = T(T+2)\sum_{k=1}^{m}{\frac{{\hat{\rho}_k}^2}{T-k}} \sim {\chi}^2(m)\ \tag{1.11}\label{1.11}$$つまり、\(H_0: {\rho}_1 = {\rho}_2 = \cdots = {\rho}_m = 0\)という帰無仮説を、\(H_1: \)少なくとも\(1\)つの\(k\in [1, m]\)において\({\rho}_k\ne 0\)という対立仮説に対して検定する。ここで\({\chi}^2(m)\)は自由度\(m\)のカイ二乗分布を表す。この\(m\)の選択について、小さい\(m\)では高次の自己相関を見逃す可能性があり、 大きい\(m\)だと検出力が小さくなるという問題がある。適切な\(m\)について明確な基準はないが、\(m \approx \log{T}\)が一つの目安となる。
以下はドル円為替終値の過去\(100\)日のコレログラムである。lagは50日目までとしている。
これを見ると、前日からのデータとは高い相関があり、減少傾向にあることなどが言えよう。参考までに、Python3のコードは以下になる。
#Correlogram
import statsmodels.graphics.api as smg
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(usdjpy_d["closing"][-100:], lags = 50);偏自己相関を求める場合は、上の”plot_acf”を”plot_pacf”にすれば良い。
https://www.statsmodels.org/dev/generated/statsmodels.graphics.tsaplots.plot_acf.html
https://www.statsmodels.org/dev/generated/statsmodels.graphics.tsaplots.plot_pacf.html
2. ARMA過程
ARMA過程 (AutoRegressive Moving Average)過程はより複雑なモデルの基礎となる。
2.1 ARMA過程の性質
1章から、時系列データは自己相関を持つものが多く、モデルに自己相関を組み込む必要がある。方法は\(2\)つあり、ここでは\(1\)次の自己相関を考える。まずは、\(y_t\)と\(y_{t-1}\)に共通の成分を含める方法で、$$\begin{cases}y_t = a+b\\ y_{t-1} = b+c\end{cases}$$のようにモデル化できる。もう一つの方法が$$y_{t} = ay_{t-1} + b$$のようなモデル化で、前者が移動平均 (Moving Average)モデルで、後者が自己回帰 (Auto Regressive)モデルである。これら\(2\)つを組み合わせたものがARMA過程になる。
2.1.1 MA過程
移動平均(MA)過程はホワイトノイズを拡張したものである。具体的にはホワイトノイズの線形和で表される。\(1\)次のMA過程は$$y_t = \mu + {\epsilon}_t + {\theta}_1{\epsilon}_{t-1} \ \ {\epsilon}_t\sim \text{W.N.}({\sigma}^2) \tag{2.1}\label{2.1}$$と定義され、これをMA(\(1\))などとする。\(y_t\)がMA\(1\)に従うことは\(y_t\sim MA(1)\)と表記される。このとき、\(y_{t-1}\)を考えると、$$y_{t-1} = \mu + {\epsilon}_{t-1} + {\theta}_1{\epsilon}_{t-2}$$となり、\(y_t\)と\(y_{t-1}\)は\({\epsilon}_{t-1}\)という共通項を持つ。\eqref{1.8}と比べると\({\theta}_1\)というパラメータが増えており、このパラメータが\(1\)次自己相関の強さを決定する。せっかく上でホワイトノイズのグラフを描いているので、ここではいくつかのパラメータの組み合わせでMA\(1\)過程を描いてみよう。
コードはこちら。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
mu = [0, 2, -2]
theta = [[0.8, 0.5, 0.3], [-0.3, -0.5, -0.8]]
std = [1, 0.5, 2]
dlen = 124
fig, ax = plt.subplots(3, 2, sharex = True, sharey = False, figsize = (10, 10))
for j in range(2):
for i in range(3):
mu_temp = mu[i]
theta_temp = theta[j][i]
std_temp = std[i]
wn = np.random.normal(0, std_temp, dlen)
y = mu_temp + wn[1:] + theta_temp*wn[:1]
ax[i, j].plot(y)
ax[i, j].plot(0, 0, label = "$\ \mu$ = {0} \n$\ \\theta$ = {1} \n$\ \sigma$ = {2}".format(mu_temp, theta_temp, std_temp))
ax[i, j].legend(fontsize = 8)
plt.tight_layout()細かい点だが、labelの\(\theta\)部分で
\tはタブと捉えられてしまうので、
\\thetaとバックスラッシュを\(2\)個重ねる必要がある。
このグラフをみると、MA\(1\)過程が\(\mu\)の周りを変動していることがわかる。つまり、期待値が\(\mu\)であることが予想される。このことを示すには、$$\begin{eqnarray}E(y_i) & = & E(\mu + {\epsilon}_t + \theta_1{\epsilon}_{t-1})\\ & = & E(\mu) + E({\epsilon}_t) + E(\theta_1{\epsilon}_{t-1})\\ & = & \mu\end{eqnarray}$$を見ると良い。最後の等号はホワイトノイズの性質\eqref{1.6}から明らかである。
次に、\(MA(1)\)過程の分散を考慮する。$$\begin{eqnarray}{\rho}_0 & = & \text{Var}(y_t)\\ & = & \text{Var}(\mu + {\epsilon}_t + \theta_1{\epsilon}_{t-1})\\ & = & \text{Var}({\epsilon}_t)+{\theta_1}^2\text{Var}({\epsilon}_{t-1})+2\theta_1\text{Cov}({\epsilon}_t, {\epsilon}_{t-1})\\ & = & (1+{\theta_1}^2){\sigma}^2\end{eqnarray}$$である。最後の等号はやはりホワイトノイズの性質\eqref{1.7}から。つまり、\eqref{1.8}のモデルと比べて、\(MA(1)\)過程の分散は\({\theta_1}^2{\sigma}^2\)の分だけ大きくなる。この事実をもって\eqref{1.8}のグラフと\(MA(1)\)のグラフを眺めて見るとよい。また、\(MA(1)\)のグラフ単体を見ると、\(\theta_1\)が大きくなるにつれてグラフがなめらかに見えることに気がつく。これは\(MA(1)\)が\(0\)でない自己相関を持つためである。つまり、\(\theta_1\)が正の場合、\(y_t\)に含まれる\(\theta_1{\epsilon}_{t-1}\)と\(y_{t-1}\)に含まれる\({\epsilon}_{t-1}\)が同じ符号を持つため、\(y_t\)と\(y_{t-1}\)が同じ方向に動く傾向がある。その結果、\(1\)次の正の自己相関が生じ、それは\(\theta_1\)が\(1\)に近づくほど強くなるので、グラフがより滑らかになる。逆に\(\theta_1\)が負の場合、\(y_t\)と\(y_{t-1}\)は傾向として逆に動き、\(1\)次の負の自己相関が生まれる。これは\(\theta_1\)が\(-1\)に近づくほど強くなるので、グラフがギザギザとした形になる。
そこで\(MA(1)\)過程の自己相関の値を求める。$$\begin{eqnarray}{\gamma}_1 & = & \text{Cov}(y_t, y_{t-1})\\ & = & \text{Cov}(\mu + {\epsilon}_t + {\theta}_1{\epsilon}_{t-1}, \mu + {\epsilon}_{t-1}+\theta_1{\epsilon}_{t-2})\\ & = & \text{Cov}({\epsilon}_t, {\epsilon}_{t-1}) + \text{Cov}(\epsilon_t, \theta_1{\epsilon}_{t-2}) + \text{Cov}(\theta_1\epsilon_{t-1}, \epsilon_{t-1}) + \text{Cov}(\theta_1\epsilon_{t-1}, \theta_1\epsilon_{t-2})\\ & = & \theta_1\text{Cov}(\epsilon_{t-1}, \epsilon_{t-1})\\ & = & \theta_1{\sigma}^2\end{eqnarray}$$になる。これから\(MA(1)\)過程の\(1\)次自己相関が$${\rho}_1 = \frac{{\gamma}_1}{{\gamma}_0} = \frac{\theta_1}{1+{\theta_1}^2} \tag{2.2}\label{2.2}$$で与えられることがわかる。この絶対値は\(\theta_1 = \pm 1\)のときに最大値\(\displaystyle \frac{1}{2}\)を取る。この事実は、絶対値が\(\displaystyle \frac{1}{2}\)よりも大きい過程はどのようにしても\(MA(1)\)過程としてモデル化することが難しいということを示している。
最後に\(MA(1)\)過程の\(2\)次以降の自己共分散と自己相関について考える。\(1\)次共分散のときと同様に、\(2\)次以降の自己共分散は\(k\geq 2\)として$$\begin{eqnarray}{\gamma}_k & = & \text{Cov}(y_t, y_{t-k})\\ & = & \text{Cov}(\mu + \epsilon_{t} + \theta_1 \epsilon_{t-1}, \mu + \epsilon_{t-k}+\theta_1\epsilon_{t-k})\\ & = & 0\end{eqnarray}$$となる。これは\(MA(1)\)過程の\(2\)次以降の自己相関が\(0\)となることを意味する。これも上と同様に、\(MA(1)\)過程は\(1\)次自己相関をモデル化することはできるが、\(2\)次以降の自己相関を記述できないということになる。
また、\(MA(1)\)過程の期待値も自己共分散も時刻\(t\)に依存しないという定常性\eqref{1.4}, \eqref{1.5}の性質を満たすので、\(MA(1)\)過程はパラメータの値に関わらず定常過程である。
以下に\(MA(1)\)過程および\(MA(2)\)過程のコレログラムを示す。\(\theta_1, \theta_2\)の正負とコレログラムの形状をよく確認されたい。
以下にこのコレログラムを描写したPython3のコードを示す。
import numpy as np
from matplotlib import pylab as plt
import seaborn as sns
sns.set()
import statsmodels.graphics.api as smg
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
theta1 = [[0.8, 0.8, 0.8], [-0.8, 0.3, -0.8]]
theta2 = [[0, 0.3, -0.3], [0, 0.8, -0.3]]
mu = 0
dlen = 10000
fig, ax = plt.subplots(3, 2, sharex = True, sharey = False, figsize = (10, 10))
for j in range(2):
for i in range(3):
theta1_temp = theta1[j][i]
theta2_temp = theta2[j][i]
wn = np.random.normal(0, 1, dlen)
y = wn[:-2] + theta1_temp*wn[1:-1] + theta2_temp*wn[2:]
plot_acf(y, ax = ax[i, j])
ax[i, j].plot(0, 0, label = "$\\theta_1$ = {0} \n$\ \\theta_2$ = {1}".format(theta1_temp, theta2_temp))
ax[i, j].legend(fontsize = 8)
plt.tight_layout()一般に\(MA(q)\)過程は\(q\)次の自己相関しかモデル化できないが、\(q\)を増やせば良いという発想もある。これが\(q\)次の移動平均過程で、$$y_t = \mu + \epsilon_t + \theta_1\epsilon_{t-1}+\theta_2\epsilon_{t-2}+\cdots + \theta_q\epsilon_{t-q}, \epsilon_t \sim \text{W.N.}(\sigma^2)\ \ \tag{2.3}\label{2.3}$$これは\(MA(q)\)過程と表記される。性質は、
定理2.1 （\(MA(q)\)過程の性質）\(MA(q)\)過程\eqref{2.3}は以下の性質を持つ。
\((1)\) \(E(y_t) = \mu\)
\((2)\) \({\rho}_0 = \text{Var}(y_t) = (1+{\theta_1}^2 + {\theta_2}^2+\cdots + {\theta_q}^2){\sigma}^2\)
\((3)\) $${\rho}_k = \begin{cases}(\theta_k + \theta_1\theta_{k+1}+\cdots + \theta_{q-k}\theta_q){\sigma}^2, \ \ & & 1\leq k\leq q\\ 0, & & k\geq q+1\end{cases}$$
\((4)\) \(MA)\)過程は常に定常である。
\((5)\) $${\rho}_k = \begin{cases}\displaystyle \frac{\theta_k + {\theta}_1{\theta}_{k+1}+\cdots + {\theta}_{q-k}\theta_{q}}{1+{\theta_1}^2+{\theta_2}^2+\cdots + {\theta_q}^2},\ \ & &1\leq k\leq q\\ 0, & & k\geq q+1\end{cases}$$\(MA\)過程は常に定常であること、および\(MA(q)\)過程の\(q+1\)次以降の自己相関は\(0\)になることを示す\((4), (5)\)は特に重要である。
\(MA(q)\)過程は\(q\)次の自己相関をモデル化するのに\(q\)個のパラメータを要する。これから、あまり長期間の自己相関をモデル化するのには向かないモデルである。また、ホワイトノイズは観測できないが、それが式に含まれるのでモデルの解釈がやや難しくなる。このことはモデル推定や予測がより複雑になるということであり、大きな弱点である。そこで以下に示す自己回帰過程（ARモデル）が登場する。
2.1.2 AR過程
自己回帰(AR)過程 (autoregressive processe)は過程が自身の過去に回帰された形で表現される過程である。\(1\)次AR過程（AR(1)）は$$y_t = c+\phi_1y_{t-1} + \epsilon_t, \epsilon_t \sim \text{W.N}({\sigma}^2)\tag{2.4}\label{2.4}$$で定義される。\eqref{2.4}は\eqref{1.8}に\(\phi_1y_{t-1}\)が追加された形であり、\(y_t\)と\(y_{t-1}\)が相関を持つことがすぐに分かる。具体的な例として、以下があげられている。$$y_t = 1+0.5y_{t-1}+\epsilon_t, \epsilon_t \sim \text{W.N}(1)$$まず、ホワイトノイズ\(\epsilon_t\)が決まり、それに続いて\(y\)の値が順次決まっていく。
以下の図は、いくつかのパラメータの組み合わせに対して、\(\text{AR}(1)\)過程をグラフ化したものである。
このグラフを描写するPythonコードは以下となる。ただし、上の4つのグラフに関しては初期値が\(\displaystyle \text{N}\left(\frac{c}{1-\phi}, \frac{{\sigma}^2}{1-{\phi}^2}\right)\)に従うことを用いた。また、下の2つのグラフの初期値は\(0\)としている。
import numpy as np
from matplotlib import pylab as plt
import seaborn as sns
sns.set()
import statsmodels.graphics.api as smg
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
const = np.array([[[2, 0.8, 1], [-2, 0.3, 0.5]], [[0, -0.3, 2], [-2, -0.8, 1]], [[2, 1, 1], [2, 1.1, 1]]])
dlen = 10000
fig, ax = plt.subplots(3, 2, sharex = True, sharey = False, figsize = (10, 10))
for i in range(3):
for j in range(2):
c_temp = const[i][j][0]
phi_temp = const[i][j][1]
sigma_temp = const[i][j][2]
wn = np.random.normal(0, sigma_temp, dlen)
y = np.zeros(100)
if 0 <= i <= 1: y[0] = np.random.normal(c_temp/(1-phi_temp), sigma_temp**2/(1-phi_temp**2))
else: y[0] = 0
for t in range(2, 100):
y[t] = c_temp + phi_temp*y[t-1] + np.random.normal(0, sigma_temp)
ax[i, j].plot(y)
ax[i, j].plot(0, 0, label = "$c$ = {0} \n$\ \\phi$ = {1} \n$\ \\sigma$ = {2}".format(c_temp, phi_temp, sigma_temp))
ax[i, j].legend(fontsize = 8)
plt.tight_layout()グラフを診ると、下の2つは上4つと大きく異なる。これは\(AR\)過程の定常性に起因するものである。\(AR(1)\)の場合、\(|\phi| < 1\)のときに過程は定常となり、それ以外では過程は指数的に上昇しており、爆発的(explosive)と言われる。
以下では\(AR(1)\)は定常的、すなわち\(|\phi|< 1\)であるとして、上記のグラフからその性質を考えていく。その際、\(MA\)過程のグラフも合わせて見ると理解が深まる。
まず、過程の期待値についてであるが、\(MA\)過程と異なり定数ではない。例えば左上のグラフでは、\(c = 2\)であるのに、グラフは\(10\)の周りを変動している。実際、\(AR\)過程の期待値が\(\displaystyle \mu = \frac{c}{1-\phi}\)で与えられることを確認しよう。\eqref{2.4}の両辺の期待値を考え、$$E(y_t) = E(c+\phi y_{t-1} + \epsilon_t) = c + \phi E(y_{t-1})$$となり、\(y\)は定常なので、\(\mu = E(y_t) = E(y_{t-1})\)であるから、$$\mu = c + \phi \mu$$であり、これより期待値\(\displaystyle \mu = \frac{c}{1-\phi}\)を得る。
また、\(MA\)過程と同様に、過程の分散が撹乱項の分散\({\sigma}^2\)と異なるのも図からわかる。\(MA(1)\)過程と比べると、\(AR(1)\)過程の分散は撹乱項の分散よりも大きいことが予想される。これも定常性の仮定を用いることで容易に証明される。再び\eqref{2.4}から、両辺の分散を考えると、$$\begin{eqnarray}\text{Var}(y_t) & = & \text{Var}(c + \phi y_{t-1} + \epsilon_t) \\ & = & {\phi}^2 \text{Var}(y_{t-1}) + \text{Var}(\epsilon_t) + 2\text{Cov}(y_{t-1}, \epsilon_t) \\ & = & {\phi}^2\text{Var}(y_{t-1}) + {\sigma}^2\end{eqnarray}$$となる。ただし、最後の統合は撹乱項\(\epsilon_t\)が過去の\(y\)とは無相関であるので、\(\text{Cov}(y_{t-1}, \epsilon_t) = 0\)であることから成立する。\(y_t\)が定常のとき、\(\gamma_0 = \text{Var}(y_t) = \text{Var}(y_{t-1})\)であるから、結局\(AR(1)\)過程の分散は\(\displaystyle \gamma_0 = \frac{{\sigma}^2}{1-{\phi}^2}\)で与えられることがわかる。
最後に、\(AR(1)\)過程の自己相関を考えよう。\(MA\)過程のグラフと見比べると、同じ\(\theta = \phi = 0.8\)であっても、\(AR(1)\)過程の方が滑らかに見える。これは、\(\theta_1 = \phi_1\)のとき、一次自己相関の絶対値が\(AR(1)\)の方が\(MA(1)\)過程よりも大きくなること、および\(\phi_1 > 0\)のとき、\(AR(1)\)過程は二次以降の自己相関もすべて正になることに起因している。これを確認するため、\(k\)次の自己共分散を考えると、$$\begin{eqnarray}\gamma_k & = & \text{Cov}(y_t, y_{t-k})\\ & = & \text{Cov}(\phi_1 y_{t-1} + \epsilon_t, y_{t-k}) \\ & = & \text{Cov}(\phi_1y_{t-1}, y_{t-k}) + \text{Cov}(\epsilon_t, y_{t-k}) \\ & = & \phi_1 \gamma_{k-1}\end{eqnarray}$$が得られる。両辺を\(\gamma_0\)で割ると、$$\rho_k = \phi_1 \rho_{k-1} \tag{2.8}\label{2.8}$$が得られる。これをユール・ウォーカー方程式(Yule-Walker equation)という。これは、\(AR\)過程の自己相関が、\(y_t\)が従う\(AR\)過程と同一の係数をもつ差分方程式に従うことを示すものとなる。この方程式と、\(\rho_0 = 1\)を用いて、順次\(AR\)過程の自己相関を求める事ができる。\(AR(1)\)の場合は容易で、$$\rho_k = {\phi_1}^2{\rho_{k-2}} = {\phi_1}^3{\rho_{k-3}} = \cdots = {\phi_1}^k{\rho_0} = {\phi_1}^k$$となる。\(|\phi_1| < 1\)であるから、\(AR(1)\)過程の自己相関の絶対値は指数的に減衰していく。どういうことかと言うと、\(AR(1)\)過程は二次以降の自己相関も\(0\)にはならないが、コレログラムは制約的であり、より複雑な自己相関構造を記述するためには、一般的なモデルを要するということである。
そこで用意されるのが、\(AR(p)\)過程であり、以下のように記述される。$$y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p + y_{t-p} + \epsilon_t, \epsilon_t \sim \text{W.N.}({\sigma}^2) \tag{2.9}\label{2.9}$$\(AR(p)\)過程は、\(AR(1)\)過程と同様、常に定常になるとは限らないが、以下は定常であるとして、性質をまとめる。
定理2.2（定常\(AR\)）過程の性質）定常\(AR(p)\)過程\eqref{2.9}は以下の性質をもつ。
\((1)\) $$\mu = E(y_t) = \frac{c}{1-\phi_1-\phi_2-\cdots \phi_p} \tag{2.10}\label{2.10}$$
\((2)\) $$\gamma_0 = \text{Var}(y_t) = \frac{{\sigma}^2}{1-\phi_1 \rho_1 -\phi_2 \rho_2 -\cdots -\phi_p \rho_p}$$
\((3)\) 自己共分散と自己相関は\(y_t\)が従う\(AR\)過程と同一の係数をもつ以下の\(p\)次差分方程式に従う。$$\gamma_k = \phi_1\gamma_{k-1} + \phi_2 \gamma_{k-2} + \cdots + \phi_p \gamma_{k-p}, \ \ \ k\geq 1 \tag{2.11}\label{2.11}$$ $$\rho_k = \phi_1 \rho_{k-1} + \phi_2 \rho_{k-2} + \cdots + \phi_p + \rho_{k-p}, \ \ k\geq 1 \tag{2.12}\label{2.12}$$ここで、\eqref{2.12}はユール・ウォーカー方程式と呼ばれる。
\((4)\) \(AR\)過程の自己相関は指数関数的に減衰する。
重要なのは自己相関に関するもので、\(AR\)過程の自己相関はユール・ウォーカー方程式\eqref{2.12}を落ちいて求めることができ、また自己相関は指数関数的に減衰する。
一般の\(AR(p)\)過程の場合、\(\rho_1, \cdots, \rho_p\)はユール・ウォーカー方程式と\(\rho_0 = 1, \rho_k = \rho_{-k}\)という性質を用いて得られる\(p-1\\)次の連立方程式を解いて求めることができる。\(\rho_p\)以降は、やはりユール・ウォーカー方程式を利用して逐次的に求めることができる。
\(AR(1), AR(2)\)過程のコレログラムを見ておく。
最上段が\(AR(1)\)、残りの4つは\(AR(2)\)のコレログラムとなる。注目すべき点としては、\(AR(2)\)が循環する自己相関を表現可能な点である。理論的には、方程式$$1-\phi_1z-\phi_2z^2 = 0$$が共役な複素数\(a \pm bi\)を解に持つとき、自己相関は循環的になり、またその周期は\(\displaystyle \frac{2\pi}{\cos^{-1}{\left(\frac{a}{\sqrt{a^2+b^2}}\right)}}\)で与えられる。この結果は、\(AR(p)\)にも容易に拡張でき、$$1-\phi_1 z – \cdots -\phi_p z^p = 0 \tag{2.15}\label{2.15}$$という方程式が共役な複素数\(a \pm bi\)を解にもつとき、\(AR(p)\)過程の自己相関は、\(\displaystyle \frac{2\pi}{\cos^{-1}{\left(\frac{a}{\sqrt{a^2+b^2}}\right)}}\)に等しい循環成分を持つ。したがって、\(AR(p)\)過程の自己相関は、最大\(p/2\)個の循環成分を有することができる。景気循環というように、経済データは循環的な変動を示す。このことは、\(AR\)過程をより魅力的なものとする理由になるであろう。
上記のPythonコードも掲載しておく。
import numpy as np
from matplotlib import pylab as plt
import seaborn as sns
sns.set()
import statsmodels.graphics.api as smg
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
phi1 = [[0.8, 0.5, 0.5], [-0.8, 0.1, 0.9]]
phi2 = [[0, 0.35, -0.8], [0, 0.5, -0.8]]
mu = 0
dlen = 10000
fig, ax = plt.subplots(3, 2, sharex = True, sharey = False, figsize = (10, 10))
for j in range(2):
for i in range(3):
phi1_temp = phi1[j][i]
phi2_temp = phi2[j][i]
sigma_temp = 1.0
y = np.zeros(100)
if i == 0:
y[0] = np.random.normal(2/(1-phi1_temp), sigma_temp**2/(1-phi1_temp**2))
for t in range(2, 100):
y[t] = phi1_temp*y[t-1] + np.random.normal(0, sigma_temp)
else:
y[0] = np.random.normal(2/(1-phi1_temp), sigma_temp**2/(1-phi1_temp**2))
y[1] = np.random.normal(2/(1-phi1_temp), sigma_temp**2/(1-phi1_temp**2))
for t in range(2, 100):
y[t] = phi1_temp*y[t-1] + phi2_temp*y[t-2]
plot_acf(y, ax = ax[i, j])
ax[i, j].plot(0, 0, label = "$\\phi_1$ = {0} \n$\ \\phi_2$ = {1}".format(phi1_temp, phi2_temp))
ax[i, j].legend(fontsize = 8)
plt.tight_layout()2.1.3 ARMA過程
\(ARMA\)過程は、自己回帰移動平均過程である(autoregressive moving average process)。今までの自己回帰項と、移動平均項の両方を含んでいる。\((p, q)\)次の\(ARMA\)過程は、$$y_t = c + \phi_1y_{t-1} + \cdots + \phi_py_{t-p} + \epsilon_t + \theta_1\epsilon_1 + \cdots + \theta_q\epsilon_{t-q}, \ \ \epsilon_t \sim \text{W.N.}(\sigma^2) \tag{2.16}\label{2.16}$$と定義される。この過程は\(AR\)過程と\(MA\)過程の性質を併せ持っており、両家庭の性質の強い方が\(ARMA\)過程の性質になる。例えば、\(MA\)過程は常に定常であるが、\(AR\)過程は定常であるとは限らない。この場合は、\(ARMA\)過程も\(AR\)過程と同様に、常に定常であるかは不明ということになる。
定理2.3（定常\(ARMA(p, q)\)過程の性質）定常\(ARMA(p, q)\)過程\eqref{2.16}は、以下の性質を有する。
\((1)\) \(\displaystyle \mu = E(y_t) = \frac{c}{1-\phi_1phi_2-\cdots -\phi_p}\)
\((2)\)\(q+1\)次以降の自己共分散と自己相関は\(y_t\)が従う\(ARMA\)過程の\(AR\)部分と同一の係数を持つ以下の\(p\)次差分方程式（ユール・ウォーカー方程式）に従う。$$\begin{eqnarray}\gamma_k & = & \phi_1\gamma_{k-1} + \phi_2\gamma_{k-2} + \cdots + \phi_p\gamma_{k-p}, \ \ k\geq q+1 \\ \sigma_k & = & \phi_1\sigma_{k-1} + \phi_2 \sigma_{k-2} + \cdots + \phi_p \sigma_{k-p}, \ \ k\geq q+1\end{eqnarray}$$
\((3)\) \(ARMA\)過程の自己相関は指数的に減衰する。
\((2)\)から分かるように、\(ARMA(p, q)\)過程の場合、\(q+1\)次以降の自己共分散と自己相関は、ユール・ウォーカー方程式を用いて逐次求めていくことが可能である。\(q\)次までの自己共分散と自己相関は、\(MA\)項の影響があるため、一般的な表現は難しい。
2.2 \(ARMA\)過程の定常性と反転可能性
2.2.1 \(AR\)過程の定常性
\(AR(p)\)過程\eqref{2.9}の定常条件は\eqref{2.15}の方程式を用いて述べることができる。\eqref{2.15}のすべての解の絶対値が\(1\)より大きいとき、\(AR\)過程は定常である。\eqref{2.15}は\(AR\)特性方程式と呼ばれ、特性方程式の左辺の多項式は\(AR\)多項式と呼ばれる。以下で簡単な例を考えてみる。
\(AR(1)\)過程$$y_t = c + \phi_1y_{t-1} + \epsilon, \ \ \epsilon_t \sim \text{W.N.}(\sigma^2) \tag{2.17}\label{2.17}$$の定常条件を考える。特性方程式は\(1-\phi_1z = 0\)となるので、解は\(z = \phi_1^{-1}\)である。したがって、定常条件は\(|\phi_1| < 1\)である。
\(AR(2)\)過程では具体的な数値を入れてみる。$$y_t = 0.5y_{t-1} + 0.5y_{t-2} + \epsilon_{t}, \ \ \epsilon_t \sim \text{W.N.}(\sigma^2)$$の特性方程式は$$0 = 1-0.5z-0.5z^2 = (1-z)(1+0.5z)$$となるので、解は\(z = -2, 1\)である。したがって、この\(AR(2)\)過程は非定常である。より一般的には、$$y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + \epsilon_{t}, \ \ \epsilon_t \sim \text{W.N.}(\sigma^2) \tag{2.18}\label{2.18}$$の定常条件は、$$\begin{cases}\phi_2 > -1 \\ \phi_2 < 1 + \phi_1\\ \phi_2 < 1-\phi_1\end{cases} \tag{2.19}\label{2.19}$$で与えられることが知られている。
定常\(AR\)過程の性質について、\(MA(\infty)\)で書き直すことができる。正確に言うと、\(AR\)過程が定常であることと、\(AR\)過程を\(MA\)過程に書き直すことができるというのは同値になる。$$y_t = \phi_1y_{t-1} + \epsilon_t, \ \ \epsilon_t \sim \text{W.N.}(\sigma^2) \tag{2.20}\label{2.20}$$という\(AR(1)\)過程についてこれを確認しておく。\eqref{2.20}を繰り返し利用して、$$\begin{eqnarray}y_t & = & \phi_1y_{t-1} + \epsilon_t \\ & = & \phi_1 (\phi_1y_{t-2} + \epsilon_{t-1}) + \epsilon_t \\ & = & {\phi_1}^2 y_{t-2} + \epsilon_t + \phi_1 \epsilon_{t-1} \\ & = & {\phi_1}^2 (\phi_1y_{t-3} + \epsilon_{t-2}) + \epsilon_t + \phi_1\epsilon_{t-1} \\ & = & {\phi_1}^3 y_{t-3} + \epsilon_t + \phi_1\epsilon_{t-1} + {\phi_1}^2\epsilon_{t-2} \\ & \cdots & \\ & = & {\phi_1}^m y_{t-m} + \sum_{k = 0}^{m-1}{{\phi_1}^k\epsilon_{t-k}} \end{eqnarray}$$となる。\(|\phi_1| < 1\)ならば、\(m\to \infty\)のとき\({\phi_1}^my_{t-m}\to 0\)が成り立つので、$$y_{t} = \sum_{k = 0}^{\infty}{{\phi_1}^k \epsilon_{t-k}}$$と書き直すことができる。逆に\(|\phi_1| \geq 1\)のとき、この\(AR\)過程は\(MA\)過程では表現できない。したがって、\(AR\)過程の定常条件と\(AR\)過程が\(MA\)表現できる条件が同一になることが分かる。また、定常\(AR\)過程を\(MA(\infty)\)過程に書き直したとき、\(MA\)係数が指数的に減衰していくことが、上記の表現から分かる。
2.2.2 \(MA\)過程の反転可能性
\(MA\)過程は常に定常である（2.2.1）。したがって、\(MA\)過程の定常性については何も心配する必要がないのであるが、実は以下の問題がある。すなわち、任意の\(MA\)過程について、同一の期待値と自己相関構造を持つ、異なる\(MA\)過程が複数存在するのである。以下の\(MA(1)\)過程、$$y_t = \epsilon_t + \theta \epsilon_{t-1}, \ \ \epsilon_t \sim \text{W.N.}({\sigma}^2) \tag{2.21}\label{2.21}$$を用いてこのことを確認する。2.1.1から\eqref{2.21}の期待値と自己共分散は$$\begin{cases} \text{E}(y_t) & = & 0 \\ \gamma_0 & = & (1+{\theta}^2){\sigma}^2 \\ \gamma_1 & = & \theta{\sigma}^2 \\ \gamma_k & = & 0, k\geq 2\end{cases} \tag{2.22}\label{2.22}$$となる。次に、$$y_t = \bar{\epsilon_t} + \frac{1}{\theta}\bar{\epsilon_t} \tag{2.23}\label{2.23}$$という別の\(MA(1)\)過程を考える。\eqref{2.22}の結果を基に、\(MA(1)\)過程\eqref{2.23}の期待値と自己共分散を求めると、$$\begin{cases}\text{E}(y_t) & = & 0 \\ \gamma_0 & = & \left(1 + \frac{1}{\theta^2}\right)\theta^2{\sigma}^2 \\ & = & (1+{\theta}^2){\sigma}^2 \\ \gamma_1 & = & \frac{1}{\theta}\cdot \theta^2{\sigma}^2\\ & = & \theta{\sigma}^2 \\ \gamma_{k} & = & 0, k\geq 2\end{cases}$$となる。つまり、\eqref{2.21}と\eqref{2.23}という\(2\)つの\(MA(1)\)過程は同一の期待値と自己相関係数を有する。時系列モデルの導入の目的の一つは、データの平均的な挙動と自己相関構造をモデル化することになる。しかし、この観点からすると、同一の期待値と同一の自己相関構造を有する\(MA\)過程が複数存在することから、どの\(MA\)過程を採用すべきか明らかでない。一つの基準として、\(MA\)過程の反転可能性(invertibility)がある。
定義2.1 （\(MA\)過程の反転可能性） \(MA\)過程が\(AR(\infty)\)に書き直せるとき、\(MA\)過程は反転可能(invertible)といわれる。
\(MA\)過程が反転可能のとき、撹乱項\(\epsilon_t\)は過去の\(y_t\)の関数として表現され、さらに過去の\(y\)を用いて\(y_t\)を予測したときの予測誤差として解釈可能である。このため、反転可能表現に伴う\(\epsilon_t\)を\(y_t\)の本源的な撹乱項 (fundamental innovation)と呼ぶことがある。この本源的な撹乱項を使った\(MA\)過程を用いると、パラメータの推定や予測を自然な形で行うことができる。すなわち、同一の期待値と自己相関構造を有する複数の\(MA\)過程のうち、反転可能な表現を選択することが望ましい。同一の期待値と自己相関構造を有する\(MA(q)\)過程は一般に\(2^q\)個存在する。そのうち、反転可能なものは一つだけとなる。
\(MA(q)\)過程\eqref{2.3}の反転可能条件は、\(AR\)過程の定常条件と同様に、$$1+\theta_1 z + \theta_2 z^2 + \cdots + \theta_p z^p = 0 \tag{2.24}\label{2.24}$$という\(MA\)特性方程式を用いて述べることが可能である。具体的には、\eqref{2.24}のすべての解の絶対値が\(1\)よりも大きいとき、\(MA\)過程は反転可能となる。
以下でいくつか例を見てみる。\(MA(1)\)過程\eqref{2.21}の特性方程式は\(1 + \theta z = 0\)となるので、反転可能条件は\(|z| = |\theta^{-1}| > 1\)、すなわち\(|\theta| < 1\)で与えられる。
以下では\(MA(1)\)過程\eqref{2.21}が反転可能なとき、\(AR\)過程で書き直せることを確認する。定常\(AR(1)\)過程を\(MA(\infty)\)に書き直したときと同様に、\(\epsilon_t\)を\(AR\)過程で書き直すと、$$\begin{eqnarray}\epsilon_t & = & -\theta \epsilon_{t-1} + y_t \\ & = & (-\theta)^m \epsilon_{t-m} + \sum_{k = 0}^{m-1}{(-\theta)^{k}y_{t-k}} \\ & \to & \sum_{k = 0}^{\infty}{(-\theta)^{k} y_{t-k}}\end{eqnarray}$$となる。反転可能条件\(|\theta| < 1\)から、\(m \to \infty\)のとき\((-\theta)^{m}\epsilon_{t-m}\to 0\)が成立することに注意。したがって、\(MA(1)\)過程\eqref{2.21}は、反転可能なとき$$y_t = -\sum_{k = 1}^{\infty}{(-\theta)^{k}y_{t-k}} + \epsilon_t$$という\(AR(\infty)\)過程で書き直すことができる。
2.2.3 \(ARMA\)過程の定常・反転可能性
\(MA\)過程は常に定常なので、\(ARMA\)過程の定常性を考えるとき、\(MA\)過程の部分を無視できる。したがって、\(AR\)過程の特性方程式\eqref{2.15}を考え、そのすべての解の絶対値が\(1\)よりも大きければ、\(ARMA\)過程は定常となる。
同様に、\(ARMA\)過程の反転可能性の条件は、\(AR\)過程はすでに\(AR\)過程で表現できているので、\(MA\)過程の部分が反転可能であればいいということになる。したがって、\(MA\)過程の特性方程式\eqref{2.24}を考え、そのすべての解の絶対値が\(1\)よりも大きければ、\(ARMA\)過程は反転可能となる。
2.3 \(ARMA\)モデルの推定
2.3.1 最小二乗法
最小二乗法(OLS: ordinary least squares)について述べる。考え方はシンプルで、モデルに含まれる説明できない部分（誤差）を最小にしようというものである。以下では\(AR(1)\)モデル$$y_t = c + \phi y_{t-1} + \epsilon_t, \ \ \epsilon_t \sim \text{iid}(0, {\sigma}^2) \tag{2.25}\label{2.25}$$を用いてOLSについて理解する。OLSでは推定するモデル\eqref{2.25}のことを回帰モデル（regression model）と呼び、左辺の変数\(y_t\)のことを被説明変数(dependent variable)と言う。また、\(\epsilon_t\)のことを誤差項(error term)という。誤差項以外の右辺の変数を説明変数(independent variable)という。この場合、\(c, y_{t-1}\)が説明変数である。\eqref{2.25}は説明変数が一つしか無いので単回帰モデル(simple regression model)であり、\(p > 1\)の\(AR(p)\)モデルのように定数以外の説明変数が複数ある場合は、重回帰モデル(multiple regression model)である。\eqref{2.25}における\(c, \phi\)の任意の推定量を\(\bar{c}, \bar{\phi}\)とすると、OLSでモデルが説明できない誤差は$$e_t = y_t-\tilde{c}-\tilde{\phi}y_{t-1}$$であり、これを残差(residual)という。今は、残差平方和(SSR: sum of squared residuals)$$\text{SSR} = \sum_{t = 1}^{T}{{e_t}^2} = \sum_{t = 1}^{T}{(y_t-\tilde{c}-\tilde{\phi}y_{t-1})^2} \tag{2.26}\label{2.26}$$が最小となるような\(\bar{c}, \bar{\phi}\)を選択する。\eqref{2.26}の両辺を\(\bar{c}, \bar{\phi}\)で偏微分して、$$\begin{cases}\displaystyle \frac{\partial \text{SSR}}{\partial \tilde{c}} & = & \displaystyle -2\sum_{t = 1}^{T}{(y-\tilde{c}-\tilde{\phi}y_{t-1})} \\ \displaystyle \frac{\partial \text{SSR}}{\partial \tilde{\phi}} & = & \displaystyle -2\sum_{t = 1}^{T}{y_{t-1}(y_t-\tilde{c}-\tilde{\phi}y_{t-1})}\end{cases} \tag{2.27}\label{2.27}$$が得られる。\eqref{2.27}の2式が\(0\)になるような\(\hat{c}, \hat{\phi}\)がOLS推定量であるから、$$\begin{cases}\displaystyle \sum_{t = 1}^{T}{(y_t-\hat{c}-\hat{\phi}y_{t-1})} & = & 0 \\ \displaystyle \sum_{t = 1}^{T}{y_{t-1}(y_t-\hat{c}-\hat{\phi}y_{t-1})} & = & 0\end{cases}$$を満たす。この方程式を正規方程式(normal equations)という。正規方程式を解くと、$$\begin{eqnarray}\hat{c} & = & \bar{y_t}-\hat{\phi}\bar{y_{t-1}} \tag{2.28}\label{2.28} \\ \hat{\phi} & = & \frac{\sum_{t=1}^{T}{(y_t-\bar{y_t})(y_{t-1}-\overline{y_{t-1}})}}{\sum_{t = 1}^{T}{(y_{t-1}-\overline{y_{t-1}})^2}} \tag{2.29}\label{2.29}\end{eqnarray}$$で与えられることが分かる。ただし、\(\displaystyle \bar{y_t} = \frac{1}{T}\sum_{t = 1}^{T}{y_t}\)である。この結果から分かるように、OLS推定量の計算には初期値\(y_0\)が必要である。また、撹乱項の分散\({\sigma}^2\)のOLS推定量\(s^2\)は以下のOLS残差に基づいた不偏標本分散で与えられる。$$s^2 = \frac{1}{T-2}\sum_{t = 1}^{T}{(y_t-\hat{c}-\hat{\phi}y_{t-1})^2} \tag{2.30}\label{2.30}$$以上の議論は容易に\(AR(p)\)過程に拡張できる。定理として以下にまとめる。
定理2.4（OLS推定量の性質） ARモデルにおけるOLS推定量は以下の性質を有する。
\((1)\) OLS推定量は一致推定量(consistent estimator)である。
\((2)\) OLS推定量を基準化したものは漸近的に正規分布に従う。
\((3)\) \(\epsilon_t \sim \text{iid N}(0, {\sigma}^2)\)のとき、OLS推定量は一致推定量の中で漸近的に最小の分散共分散行列を持つ。
性質\((1)\)は標本数が十分大きいとき、OLS推定量が真の値に限りなく近づくことを意味する（\(\hat{\phi}\)が\(\phi\)に収束する）。性質\((2)\)は仮説検定に用いることができる（第5章で導出方が述べられる）。また、性質\((3)\)は正規性の下でのOLSの漸近有効性(asymptotic efficiency)を表す。
回帰分析の標準的な仮定では、以上に加えてOLS推定量が不偏推定量(unbiased estimator)であり、すべての線形不偏推定量の中で、最小の分散共分散行列を持つ最良線形不偏推定量(BLUE: best linear unbiased estimator)となることが知られている（ガウス・マルコフの定理）。しかし、\(AR\)過程は説明編数が過去の誤差項と相関を有するため、OLS推定量は不偏推定量とはならない。
このように、\(AR\)モデルに対してはOLS推定量は（一部望ましくない性質を持つものの）統計的推測を行うのに十分である。しかも、\(\epsilon_t\)の分布を仮定しなくてもよいという利点もある。しかし、\(ARMA\)モデルや\(GARCH\)モデル、マルコフ転換モデルなどはOLSを単純に適用することはできない。そのような場合、OLSではなく、以下に述べる最尤法を利用する。
2.3.2 最尤法
OLSではモデルが説明できない誤差を最小にするようなパラメータ選択が行われたが、最尤法の場合得られた観測値をモデルが最も実現しやすくするようなパラメータ選択が行われる。
以下では\(AR(1)\)モデルに対して最尤法を適用してみる。ここで注意点として、最尤法では尤度の評価に\(\epsilon_t\)の分布が必要である。簡便さから正規分布が使われることが多いが、その他の分布を用いることも可能である。ここでは\(\epsilon_t \sim \text{iid N} (0, {\sigma}^2)\)と仮定しておく。
\(\theta\)を推定したいパラメータとすると、\(AR(1)\)モデル\eqref{2.25}の場合は\(\theta = (c, \phi, {\sigma}^2)^{\prime}\)となる。\(\epsilon_t\)に連続な分布を仮定した場合は、ある\(1\)点をとる確率は\(0\)となるので、確率の評価はできない。したがって、この場合は確率密度関数を\(\theta\)の関数としてみたものを尤度と呼び、尤度を最大にするパラメータが最尤推定量となる。一般的に、時系列データの尤度の計算には、条件付き分布を用いた同時密度の分解を用いることが多い。ベイズの法則(Bayes rule)により、同時密度に関して$$f_{X, Y}(x, y) = f_{X|Y}(x|y)\cdot f_{Y}(y)$$が成立するので、\(y_1\)と\(y_2\)の同時密度に関して、$$\begin{eqnarray}f_{Y_2, Y_1| Y_0}(y_2, y_2 | y_0; \theta) & = & f_{Y_2|Y_1, Y_0}(y_2|y_1, y_0; \theta)\cdot f_{Y_1|Y_0}(y_1|y_0; \theta)\\ & = & f_{Y_2|Y_1}(y_2|y_1; \theta)\cdot f_{Y_1|Y_0}(y_1|y_0; \theta)\end{eqnarray}$$が成立する。最後の等号が成立するのは\(AR(1)\)過程の場合、直近\(1\)期の値にしか条件付き分布が依存しないことによる。次に、\(y_1, y_2, y_3\)の同時密度を考えると、$$\begin{eqnarray}f_{Y_3, Y_2, Y_1 | Y_0}(y_3, y_2, y_1 | y_0; \theta) & = & f_{Y_3 | Y_2, Y_1, Y_0}(y_3 | y_2, y_1, y_0; \theta)\cdot f_{Y_2, Y_1 | Y_0}(y_2, y_1 | y_0; \theta) \\ & = & f_{Y_3 | Y_2}(y_3 | y_2; \theta) \cdot f_{Y_2 | Y_1}(y_2 | y_1; \theta)\cdot f_{Y_1 | Y_0}(y_1 | y_0; \theta)\end{eqnarray}$$が得られる。同様にすると、$$f_{Y_T, Y_{T-1}, \cdots, Y_1 | Y_0}(y_T, y_{T-1}, \cdots, y_1 | y_0; \theta) = \prod_{t = 1}^{T}{f_{Y_t} | Y_{t-1}(y_t | y_{t-1}; \theta)}$$が成立することが分かる。したがって、対数尤度は$$L(\theta) = \sum_{t = 1}^{T}{\log{f_{Y_t | Y_{t-1}}(y_t | y_{t-1}; \theta)}} \tag{2.31}\label{2.31}$$と書ける。詳しい説明は後ほど出てくるが、\(\epsilon_t \sim \text{iid N}(0, {\sigma}^2)\)のとき\(y_t | y_{t-1} \sim \text{N}(c + \phi y_{t-1}, {\sigma}^2)\)であるので、$$f_{Y_t | Y_{t-1}}(y_t | y_{t-1}; \theta) = \frac{1}{\sqrt{2\pi {\sigma}^2}}\exp{\left(\frac{-(y_t-c-\phi y_{t-1})^2}{2{\sigma}^2}\right)} \tag{2.32}\label{2.32}$$となることに注意する。\eqref{2.32}を\eqref{2.31}に代入し、$$L(\theta) = -\frac{T}{2}\log{(2\pi)} – \frac{T}{2}\log{({\sigma}^2)} – \sum_{t = 1}^{T}{\left(\frac{(y_t-c-\phi y_{t-1})^2}{2{\sigma}^2}\right)}$$となることが分かる。これを最大にするような\(c, \phi\)は\((y_t – c -\phi y_{t-1})^2\)を最小とする\(c, \phi\)であるが、これはOLS推定量そのものである。したがって、\(AR(1)\)モデルにおける係数の最尤推定量とOLS推定量は一致する。また、撹乱項の分散\({\sigma}^2\)の最尤推定量は$$\hat{{\sigma}^2} = \frac{1}{T}\sum_{t = 1}^{T}{(y_t – \hat{c}-\hat{\phi}y_{t-1})^2} \tag{2.33}\label{2.33}$$で与えられる。最尤法においても\(y_1\)の条件付き分布の計算に必要な\(y_0\)は初期値、つまり所与の値として扱うことがほとんどである。実際には\(y_0\)も確率変数として扱い、\(y_0\)も含めた尤度として計算したほうが正確ではあるが、計算が複雑になる上に、標本数\(T\)が大きければ初期値の影響は無視しても構わない。
以上の議論は\(ARMA, GARCH\)モデルに容易に一般化でき、基本的には\eqref{2.31}における\(g_{Y_t | Y_{t-1}} (y_t | y_{t-1}; \theta)\)を適切な条件付き分布に置き換えるだけで良い。具体的には、\(t\)期までに利用可能な情報の集合、時点\(t-1\)の情報集合(information set)を\(\Omega_{t-1}\)で表すと、対数尤度は一般的に$$L(\theta) = \sum_{t = 1}^{T}{\log{f_{Y_t | \Omega_{t-1}}(y_t | \Omega_{t-1}; \theta)}} \tag{2.34}\label{2.34}$$という形で書くことができるので、\eqref{2.34}を最大化するような\(\theta\)が最尤推定量となる。ただし、\(MA\)過程や\(ARMA\)過程の場合は、\(y\)だけではなく\(\epsilon\)の初期値も必要となり、\(\epsilon\)の初期値としては、期待値の\(0\)を用いる。最尤推定量の性質は以下のようになる。
定理 2.5（最尤推定量の性質） \(ARMA\)過程における最尤推定量は以下の性質を有する。
\((1)\) 最尤推定量は一致推定量である。
\((2)\) 最尤推定量を基準化したものは漸近的に正規分布に従う。
\((3)\) 最尤推定量は一致推定量の中で漸近的に最小の分散共分散行列を持つ。
性質\((1), (2)\)はOLSと同様の望ましい性質である。性質\((3)\)は最尤法の漸近有効性を示しており、分布の仮定が正しければ、最尤法が最も有効な方法であることを意味する。
2.4 \(ARMA\)モデルの選択
本節の過程は、「真のモデルが定常かつ反転可能な\(ARMA(p, q)\)過程である」というものになる。この仮定のもとで、与えられたデータに対して最適な\(ARMA\)モデルを決定する方法を議論する。これらの基礎にはBox and Jenkinsの一連の仕事がある。
2.4.1 モデル候補の選択
\(ARMA\)モデルの選択は、標本自己相関や標本偏自己相関を用いてモデルの候補を選択することから始めることが多い。最終的には、情報量基準を用いて行うのであるが、\(ARMA\)モデルの場合は選択すべき変数が\(2\)つあるので、最初に絞った方が効率的であるからである。
ここで重要なのが、\eqref{1.10}の標本自己相関関数である。まず、\(AR\)過程と\(ARMA\)過程の自己相関関数の絶対値は指数的に減衰する。対して、\(MA(q)\)過程の自己相関は\(q+1\)次以降\(0\)になる。したがって、標本自己相関が\(q+1\)次で切断され、それ以降\(0\)に近い値をとっていることが確認されれば、\(MA(q)\)過程の可能性が高いということである。しかし、標本自己相関だけからは、\(AR\)過程と\(ARMA\)過程の区別は難しい。その区別に役立つのが、偏自己相関である。
偏自己相関(partial autocorrelation)は、例えば\(y_t\)と\(y_{t-k}\)の自己相関を計算するときに、\(y_{t-1}, \cdots, y_{t-k+1}\)の影響を取り除いたものである。すなわち、\(y_t\)と\(y_{t-k}\)は\(y_{t-1}\)や\(y_{t-2}\)を通じて相関している可能性があるので、間にあるデータの影響を取り除いたうえで、\(y_t\)と\(y_{t-k}\)の相関を評価したものが\(k\)次偏自己相関である。より正確には、\(k\)次偏自己相関\(\alpha_k\)は\(y_t\)を定数と\(y_{t-1}, \cdots, y_{t-k}\)に線形射影したときの\(y_{t-k}\)の射影係数である。また偏自己相関を\(k\)の関数と見たものは偏自己相関関数と言われる。
\(AR\)モデルの偏自己相関について考えてみる。まず、\(AR(p)\)であるが、\(y_t\)は\(y_{t-1}, \cdots, y_{t-p}\)の線形関数で書ける。言い換えると、\(y_{t}\)から\(y_{t-1}, \cdots, y_{t-p}\)の影響を取り除くと、残るのは自己相関を持たないホワイトノイズ\(\epsilon_t\)の部分だけとなる。したがって、\(AR(p)\)モデルの場合、\(p+1\)次以降の偏自己相関は\(0\)となる。それに対して、\(MA\)過程は\(AR(\infty)\)過程で書き直すことができる。したがって、\(y_t\)から\(y_{t-1}, \cdots, y_{t-k+1}\)の影響を取り除いたとしても、\(y_{t}\)と\(y_{t-k}\)が相関を持つような\(k\)は限りなく存在する。しかし、2.2.2でも述べた通り、\(MA\)過程を\(AR(\infty)\)で書き直したときの\(AR\)係数の絶対値は指数的に減衰していくことが知られている。故に、\(MA\)過程は無限の偏自己相関を持つが、その絶対値は指数的に減衰していく。同様に、\(ARMA\)過程の偏自己相関も指数的に減衰していく。
標本偏自己相関関数は与えられたデータから偏自己相関を計算したものであり、\(k\)次標本自己相関\(\hat{\alpha_k}\)は\(y_t\)を定数と\(y_{t-1}, y_{t-2}, \cdots, y_{t-k}\)に回帰したときの\(y_{t-k}\)のOLS係数推定量で与えられる。標本自己相関と同様に、\(y_t\)が\(\text{iid}\)系列の場合、\(k\)次標本偏自己相関\(\hat{\alpha_k}\)は漸近的に平均\(0\)、分散\(1/T\)の正規分布に従うことが知られており、この結果を基に、\(\alpha_k = 0\)を検定することができる。また、標本偏自己相関関数は偏自己相関関数と同様の性質を持つことが要求されるので、標本偏自己相関関数の値が\(p + 1\)次で切断され、それ以降で\(0\)に近い値をとっていることが確認されれば、\(AR(p)\)過程の可能性が高いということになる。しかしながら、\(MA\)過程と\(ARMA\)過程の区別は偏自己相関からは難しい。
以上をまとめたものが、表2.2になる。
| モデル | 自己相関 | 偏自己相関 |
| \(AR(p)\)モデル | 減衰していく。 | \(p+1\)次以降は\(0\) |
| \(MA(q)\)モデル | \(q+1\)次以降\(0\) | 減衰していく。 |
| \(ARMA\)モデル | 減衰していく。 | 減衰していく。 |
これをみると、標本自己相関と標本偏自己相関を用いてある程度最適なモデルを選択できる。しかし、いずれも主観的なものになりがちであり、以下で述べる情報量基準の登場になる。
2.4.2 情報量基準
情報量基準(information criterion)は最尤法の推定結果を基に最適なモデル選択を行う客観的基準である。一般に$$\text{IC} = -2L(\hat{\theta}) + p(T)k \tag{2.35}\label{2.35}$$で定義される。ここで、\(L(\hat{\theta})\)は対数尤度を最尤推定値で評価した最大対数尤度、\(T\)はモデルの推定に用いた標本数、\(p(T)\)は\(T\)の何らかの関数、\(k\)は推定に用いたパラメータの数である。一般的に、情報量基準は\(2\)つの部分からなり、第\(1\)項はモデルの当てはまりの良さであり、第\(2\)項はモデルが複雑になることに対するペナルティとみなすことができる。モデルを複雑にすればするほど、説明力は必ず上昇するので、第\(1\)項は必ず小さくなるが、その分パラメータの推定精度が下がるため、ペナルティを課すというのが第\(2\)項である。そうすることで、パラメータの追加に対して第\(1\)項と第\(2\)項との間でトレードオフが生じるので、モデルのoverfittingを防ぐことができる。つまり、パラメータを増やすペナルティ以上にモデルが改善しなくなったとき、モデルは最適になるので、情報量基準を最小にするようなモデルが最適である。
ペナルティ関数の選択によって、複数の情報量基準が存在する。最も頻繁に用いられるものが、赤池情報量基準(\(AIC\))と呼ばれるものであり、$$\text{AIC} = -2L(\hat{\theta}) + 2k$$で定義される。これは\eqref{2.35}において\(p(T) = 2\)としたものである。もう一つは、Schwarz情報量基準(\(SIC\))もしくはベイズ情報量規準(\(BIC\))と呼ばれるものであり、\(SIC(BIC)\)は$$\text{SIC} = -2L(\hat{\theta}) + \log{(T)}k$$で定義される。つまり、\eqref{2.35}において\(p(T) = \log{T}\)としたものである。\(\log{T} > 2\)すなわち\(T\geq 8\)であれば、\(SIC\)は\(AIC\)よりもペナルティが大きいことになる。実際の解析ではほぼ\(T\geq 8\)となるので、\(SIC\)は\(AIC\)よりも小さいモデルを選ぶ傾向があり、必ずしも\(2\)つの情報量基準が同一のモデルを選択することにはならない。
真のモデルを\(AR(p)\)とし、\(p\)が有限であるとすると、\(AIC\)が選択する\(p\)に関して、$$\begin{eqnarray}\lim_{T\to\infty}{P(\hat{p} < p)} & = & 0\\ \lim_{T\to\infty}{P(\hat{p} > p)} & > & 0\end{eqnarray}$$であることが知られている。つまり、標本数\(T\)が大きくなるとき、\(AIC\)はモデル次数を過小に評価することはないが、過大なモデルを選択する確率が\(0\)とはならない。\(AIC\)はこの意味で一貫性を持たない。それに対して、\(SIC\)が選択する\(\hat{p}\)に関しては、$$\lim_{T\to\infty}{P(\hat{p} = p)} = 1$$が成立し、\(SIC\)は一貫性を有する。したがって、\(SIC\)は標本数が大きければ、ほぼ必ず正しいモデルを選択することができる。この事実からすると、\(AIC\)よりも\(SIC\)が優れているように見えるが、必ずしもそうではない。例えば、真のモデルが\(AR(\infty)\)の場合、\(AIC\)はある種の漸近最適正を持つことが知られている。また、\(AIC\)が過大なモデルを選択したとしても、標本数が大きくなるにつれて、過大な部分のパラメータの最尤推定量は真の値である\(0\)に収束するので、過大なモデルの選択はそれほど問題にならない。
注意点としては、\(AIC, SIC\)は\(ARMA\)モデルだけではなく、より複雑なモデルに対しても適用可能であるが、その場合どのような性質があるかは一般的にはよくわかっておらず、慣例として用いられることが多いことである。
2.4.3 モデルの診断
最適なモデルを選択した後、そのモデルが適切かどうかを診断する必要がある。真のモデルにおいて、\(\epsilon_t\)はホワイトノイズであるので、自己相関を持たない。したがって、モデルが適切であれば、モデルから推定された\(\hat{\epsilon_t}\)も自己相関を持たないはずである。1.4節で述べた方法を用いて、\(\hat{\epsilon_t}\)の自己相関を検定することによって、モデルの診断を行うことができる。\(ARMA(p, q)\)モデルを推定した後、\(m\)次までの自己相関を用いたかばん検定の漸近分布は\({\chi}^2(m)\)ではなく、\({\chi}^2(m-p-q)\)となる。しまり、推定した\(ARMA\)パラメータの分だけ自由度が下がる。したがって、\eqref{1.11}のかばん統計量\(Q(m)\)の値と\({\chi}^2(m-p-q)\)の\(95%\)分位点を比較し、\(Q(m)\)の方が大きければ、帰無仮説を棄却する。帰無仮説が棄却されたとき、モデルは最良のものではないので、何らかの修正を考える必要がある。
以下では、本書とは異なる例を用いてみよう。上記は、直近\(100\)日のドル円の終値の標本自己相関および標本偏自己相関である。自己相関が減衰していき、\(2\)次以降の偏自己相関が\(0\)に漸近していることが見て取れる。
PythonのLibraryを用いると、簡単に\(AIC\)や\(BIC\)を表示することができる。例えば、このデータに対して、\(AR(3)\)モデルを選択した場合、\(AIC/BIC\)は以下のようになる。
AutoReg Model Results
==============================================================================
Dep. Variable: Close No. Observations: 100
Model: AutoReg(3) Log Likelihood -122.699
Method: Conditional MLE S.D. of innovations 0.857
Date: Thu, 18 Jul 2024 AIC 255.398
Time: 16:17:27 BIC 268.272
Sample: 03-05-2024 HQIC 260.603
- 07-17-2024
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 5.7629 3.636 1.585 0.113 -1.364 12.890
Close.L1 0.8410 0.101 8.286 0.000 0.642 1.040
Close.L2 0.3596 0.130 2.773 0.006 0.105 0.614
Close.L3 -0.2374 0.101 -2.351 0.019 -0.435 -0.039
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 -1.7870 +0.0000j 1.7870 0.5000
AR.2 1.0443 +0.0000j 1.0443 0.0000
AR.3 2.2575 +0.0000j 2.2575 0.0000
-----------------------------------------------------------------------------
\(ARMA(1, 1)\)モデルは以下である。
SARIMAX Results
==============================================================================
Dep. Variable: Close No. Observations: 100
Model: ARIMA(1, 0, 1) Log Likelihood -130.958
Date: Thu, 18 Jul 2024 AIC 269.915
Time: 16:18:19 BIC 280.336
Sample: 02-29-2024 HQIC 274.133
- 07-17-2024
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 154.1993 2.555 60.342 0.000 149.191 159.208
ar.L1 0.9772 0.024 40.006 0.000 0.929 1.025
ma.L1 -0.1077 0.111 -0.968 0.333 -0.326 0.110
sigma2 0.7806 0.073 10.741 0.000 0.638 0.923
===================================================================================
Ljung-Box (L1) (Q): 0.17 Jarque-Bera (JB): 249.40
Prob(Q): 0.68 Prob(JB): 0.00
Heteroskedasticity (H): 1.85 Skew: -1.66
Prob(H) (two-sided): 0.08 Kurtosis: 9.99
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
\(MA(3)\)モデルも以下のように簡単に計算できる。
SARIMAX Results
==============================================================================
Dep. Variable: Close No. Observations: 100
Model: ARIMA(0, 0, 3) Log Likelihood -172.084
Date: Thu, 18 Jul 2024 AIC 354.168
Time: 16:22:58 BIC 367.194
Sample: 02-29-2024 HQIC 359.440
- 07-17-2024
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 154.9232 0.506 306.434 0.000 153.932 155.914
ma.L1 1.2371 0.094 13.216 0.000 1.054 1.421
ma.L2 1.1307 0.133 8.530 0.000 0.871 1.391
ma.L3 0.4677 0.089 5.280 0.000 0.294 0.641
sigma2 1.7847 0.318 5.619 0.000 1.162 2.407
===================================================================================
Ljung-Box (L1) (Q): 6.25 Jarque-Bera (JB): 1.26
Prob(Q): 0.01 Prob(JB): 0.53
Heteroskedasticity (H): 1.13 Skew: -0.14
Prob(H) (two-sided): 0.73 Kurtosis: 2.53
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
様々なモデルに対して、\(AIC/BIC\)を計算し、選択を行うといいのであるが、便利なコマンドもあり、
from statsmodels.tsa.stattools import arma_order_select_ic
arma_order_select_ic(usdjpy_d["Close"], max_ar = 4, max_ma = 4, ic = "bic", trend = "c", model_kw = None, fit_kw = None)
{'bic': 0 1 2 3 4
0 557.517439 467.641242 395.812279 367.193637 339.229120
1 277.246709 280.335990 280.110010 282.985528 286.959784
2 279.588840 278.594397 281.557818 286.148812 290.611416
3 278.417202 281.596497 286.152789 290.766598 295.097879
4 281.583746 286.131964 290.718944 295.312598 295.611528,
'bic_min_order': (1, 0)}のコマンドで\(AIC, BIC\)を基に次数を決定してくれる。上記の計算に用いたPythonコードは以下。
#Correlogram
import statsmodels.graphics.api as smg
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.ar_model import AutoReg, ar_select_order
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import arma_order_select_ic
import matplotlib.pyplot as plt
import pandas as pd
usdjpy_d = pd.read_csv("/Users/Documents/FX/foreign_exchange_historical_data/usdjpy_d.csv")[-100:]
usdjpy_d = usdjpy_d.set_index("Date")
# Plot
fig, axs = plt.subplots(1, 2, figsize = (15, 5))
plot_acf(usdjpy_d["Close"], lags = 50, ax = axs[0])
plot_pacf(usdjpy_d["Close"], lags = 40, ax = axs[1])
plt.show()
# AR
ar_model = AutoReg(usdjpy_d["Close"], 3, old_names = False)
print(ar_model.fit().summary())
# ARIMA
arima_model = ARIMA(usdjpy_d["Close"], order = (1, 0, 1))
print(arima_model.fit().summary())
# order select by IC
arma_order_select_ic(usdjpy_d["Close"], max_ar = 4, max_ma = 4, ic = "bic", trend = "c", model_kw = None, fit_kw = None)
# MA
ma_model = ARIMA(usdjpy_d["Close"], order = (0, 0, 3))
print(ma_model.fit().summary())
\(ARIMA\)モデルのorder(a, b, c)の部分は、順に\(AR\)次数、\(I\)次数、および\(MA\)次数となる。\(I\)次数については後で説明する。予測を行いたいときは、例えば\(AR\)モデルなら、
ar_model = AutoReg(usdjpy_d["Close"], 3, old_names = False)
res = ar_model.fit()
pred = res.predict()とすると、predに予測値が収納される。以下を参照。
コメント