時系列解析の基礎(04/XX) - 「大人の教養・知識・気付き」を伸ばすブログ

　以下の書籍

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

作者:竜義, 沖本
朝倉書店

Amazon

を中心に時系列解析を勉強していきます。

前回

power-of-awareness.com

前回
3.　ARMA過程
4.　予測
次回
参考文献

3.　ARMA過程

3.7　AR過程における母数推定

　 $\mathrm{AR}$ 過程における母数推定方法は、(1)最小二乗法、(2)最尤法が一般的である。 $(p,q)$ 次 $\mathrm{AR}$ 過程

$\begin{aligned} y_t=c+\displaystyle{\sum_{i=1}^{p}\phi_{i}y_{t-i}}+\varepsilon_{t},\varepsilon_{t}\sim W.N.(\sigma^2) \end{aligned}$

において推定対象となるのは係数 $\phi_{i},i=1,\cdots, p$ および定数項 $c$ 、さらには分散 $\sigma^2$ である。ここでは係数 $\phi_{i},i=1,\cdots, p$ および定数項 $c$ を知りたいとしよう。

3.7.1　最小二乗法

　最小二乗法では通常の回帰モデルと同様に残差平方和 $SSR$ を最小化するような $\phi_{i},i=1,\cdots, p$ および $c$ として推定量を得る。残差平方和 $SSR$ は以下の通りである：

$\begin{aligned} SSR=\displaystyle{\sum_{t=1}^{T}{e_{t}^2}}=\displaystyle{\sum_{t=1}^{T}\left\{y_t-\left(\tilde{c}+\sum_{i=1}^{p}\tilde{\phi}_{i}y_{t-i}\right)\right\}^2} \end{aligned}$

二次関数( $y=x^2$ )は任意の区間で下に凸であるから $SSR$ の極小値を与えるようなものを求めればよいということとなる。この計算には $p$ 個の初期値が必要となる点に注意したい。また上式からも明らかなように、最小二乗法で $\phi_{i},i=1,\cdots, p$ および $c$ を推定するには $\sigma^2$ を考慮しなくてもよい点が利点である。

　最小二乗推定量には3つの特徴がある*1：

	一致性	最小二乗推定量は一致推定量である。
	漸近正規性	最小二乗推定量を基準化したものは漸近的に正規分布に従う。
	正規性の下での漸近有効性	$\varepsilon_{t}\sim i.i.d. N(0,\sigma^2)$ のとき、最小二乗推定量は一致推定量の中で漸近的に最小の分散共分散行列をもつ。

なお $\mathrm{AR}$ 過程は過去の誤差項と相関を有するため、最小二乗推定量は不偏性を持たない点に留意されたい。

3.7.2　最尤法

　 $\mathrm{AR}$ 過程はシンプルなモデルであるために最小二乗法が可能であった一方で、 $\mathrm{ARMA}$ 過程などより複雑なモデルであれば最尤法を用いる。簡単のため $\mathrm{AR}(1)$ モデルについて、確率変数 $X,Y$ （およびその確率密度関数 $f_{X}(x), f_{Y}(y)$ ）に関する $\mathrm{Bayes}$ の定理

$\begin{aligned} f_{X,Y}(x,y) = f_{X|Y}(x|y)\cdot f_{Y}(y) \end{aligned}$

を適用すると、確率変数 $Y_0, Y_1,\cdots,Y_{T}$ および母数ベクトル $\boldsymbol{\theta}$ について

$\begin{aligned} f_{Y_T,\cdots,Y_1|Y_0}(y_T,\cdots,y_1|y_0;\boldsymbol{\theta}) = \prod_{t=1}^{T}f_{Y_t|Y_{t-1}}(y_t|y_{t-1};\boldsymbol{\theta}) \end{aligned}$

となる。したがって対数尤度 $\mathcal{l}(\boldsymbol{\theta})$ は

$\begin{aligned} \mathcal{l}(\boldsymbol{\theta}) = \displaystyle{\sum_{t=1}^{T}\log{f_{Y_t|Y_{t-1}}(y_t|y_{t-1};\boldsymbol{\theta})}} \end{aligned}$

であり、求めたかった推定量 $\hat{\boldsymbol{\theta}}$ は

$\begin{aligned} \hat{\boldsymbol{\theta}} = \mathrm{arg}\max_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\{\mathcal{l}(\boldsymbol{\theta})\} \end{aligned}$

で得られる。
　 $\mathrm{ARMA}$ 過程に限らない一般のモデルについては時点 $t-1$ の情報集合を $\Omega_{t-1}$ とすると

$\begin{aligned} \mathcal{l}(\theta) = \sum_{t=1}^{T}\log{f_{Y_t|\Omega_{t-1}}(y_t|\Omega_{t-1};\Theta)} \end{aligned}$

で表現できる。
　 $\mathrm{ARMA}$ 過程の最尤推定量は3つの性質をもつ*2：

	一致性	最尤推定量は一致推定量である。
	漸近正規性	基準化された最尤推定量は漸近的に正規分布に従う。
	漸近有効性	最尤推定量はすべての一致推定量の中でも漸近的に最小の分散共分散行列をもつ。

3.8　モデル選択

　定常かつ反転可能な $\mathrm{ARMA}(p,q)$ 過程が真のモデルであるとき、観測値を基に適当なモデルを構築する方法を考えたい。

　最初、効率的な探索を行うべく標本自己相関ないし標本偏自己相関を用いてモデル候補を絞り込む。 $\mathrm{AR}$ 過程と $\mathrm{ARMA}$ 過程の自己相関関数の絶対値は指数的に減衰していった。これに対して $\mathrm{MA}(q)$ 過程の自己相関は $q+1$ 次以降で $0$ になるために標本自己相関が $q+1$ 次以降で急遽 $0$ 近くになるかどうかが1つの判断材料になるからである。 $\mathrm{AR}$ 過程と $\mathrm{ARMA}$ 過程との判別は偏自己相関を用いて判断することができる。

　確認のために再提示しておくと、任意の $t,k$ について自己相関係数 $\rho_{k}$ は

$\begin{aligned} \rho_{k} =\frac{\gamma_{k,t}}{\sqrt{\gamma_{0,t}\gamma_{0,t-k}}}=\frac{\gamma_{k}}{\gamma_{0}} \end{aligned}$

である。
　これに対して偏自己相関 $\alpha_{k}^{(m)}$ は直近 $m$ 時点への $y_t$ の線型射影における係数であり、自己共分散に関する以下の連立方程式の解として得られる*3：

$\begin{pmatrix} \alpha_{1}^{(m)}\\ \alpha_{2}^{(m)}\\ \vdots\\ \alpha_{m}^{(m)} \end{pmatrix}=\begin{pmatrix} \gamma_{0} & \gamma_{1} &\cdots& \gamma_{m-1}\\ \gamma_{1} &\gamma_{0}&\cdots& \gamma_{m-2}\\ \vdots&\vdots&\cdots& \vdots\\ \gamma_{m-1} & \gamma_{m-2} & \cdots & \gamma_{0} \end{pmatrix}^{-1}\begin{pmatrix} \gamma_{1}\\ \gamma_{2}\\ \vdots\\ \gamma_{m}\\ \end{pmatrix}$

　 $\mathrm{AR}(p)$ モデルであれば $p+1$ 次以降の偏自己相関は $0$ となる一方で $\mathrm{MA}$ 過程であれば無限の偏自己相関を有するものの、その絶対値は指数的に減衰していく。
　整理すると、

モデル	自己相関	偏自己相関
$\mathrm{AR}(p)$ モデル	減衰していく	$p+1$ 次以降で $0$
$\mathrm{MA}(q)$ モデル	$q+1$ 次以降で $0$	減衰していく
$\mathrm{ARMA}$ モデル	減衰していく	減衰していく

3.8.1　情報量基準

　 $T$ をモデルの推定に用いた標本数としたときに $p(T)$ を $T$ に関する関数、 $\mathcal{l}(\hat{\theta})$ を最大対数尤度*4、 $k$ を推定した母数の数として、情報量基準 $IC$ を一般に

$\begin{aligned} IC =-2\mathcal{l}(\hat{\theta})+p(T)k \end{aligned}$

と定義する。
　たいていは $AIC$ ないし $SIC$ 、すなわち

$\begin{aligned} AIC&=-2\mathcal{l}(\hat{\theta})+2k,\\ SIC&=-2\mathcal{l}(\hat{\theta})+\log(T)k \end{aligned}$

を用いる*5。これを通じて適当な母数の選択を行う*6。

3.9　モデルの診断

　モデル構築を終えた後、その正当性を判断しなければその妥当性を担保できない。
　 $\mathrm{ARMA}$ モデルであれば、回帰分析でいう残差に相当する誤差項 $\varepsilon_{t}$ についてホワイトノイズであるとの仮定を置いている。そこでこれが自己相関を有するか否かを検定すればよい。このとき母数の数だけ自由度が減るため、 $\mathrm{ARMA}(p,q)$ モデルであればかばん統計量の値と $\chi^2(m-p-q)$ 分布の $100\alpha$ %点の値とを比較する。

4.　予測

4.1　予測の基本的な考え方

　これまでで時系列モデルを構築してきたが、それを用いて予測を立てることが今回の課題である。

　推定方法が様々あるように予測の考え方も様々なものがある。そのために如何なる予測が適切なものかを判断するのは難しく、その判断のための尺度が必要となる。
　一つには最小二乗法との考え方を敷衍し、予測誤差 $\hat{e}_{t+1|t}=y_{t+1}-y_{t+1|t}$ を最小化するような予測を良い予測であると見なす考え方がある。とはいえ予測誤差は確率変数の差であるからそれ自身も確率変数であり、どのような意味で最小化するものが適当なのかを定めなければ十分とは言えない。これにも様々な考え方があるが、もっとも頻用されるのは平均二乗誤差

$\begin{aligned} MSE =E[{\hat{e}_{t+1|t}}^2]=E[(y_{t+1}-y_{t+1|t})^2] \end{aligned}$

を最小化していくことである。このような平均二乗誤差を最小化することで得られる予測を最適予測（optimal forecast）と呼ぶ。

4.2　より厳密な考え方

　時点 $t$ において利用可能な情報(観測値)からなる情報集合 $\Omega_{t}=\{y_t,y_{t-1},\cdots,y_{1}\}$ を用いて、 $h\gt0$ 期先の値 $y_{t+h}$ を予測したいとする。 $\Omega_{t}$ を用いて得られる予測量を $\hat{y}_{t+h|t}$ とする。
　このとき前述の考えに従うと

$\begin{aligned} MSE\left(\hat{e}_{t+h|t}\right)=E\left[\hat{e}_{t+h|t}^2|\Omega_{t}\right] =E\left[(y_{t+h}-\hat{y}_{t+h|t})^2|\Omega_{t}\right] \end{aligned}$

を考えることとなる。
　これを展開すると、

$\begin{aligned} MSE(\hat{e}_{t+h|t})=E\left[\hat{e}_{t+h|t}^2|\Omega_{t}\right]=E\left[(y_{t+h}-\hat{y}_{t+h|t})^2|\Omega_{t}\right] \end{aligned}$

である。ここで $\mu_{t+h|t}=E\left[\hat{y}_{t+h}|\Omega_{t}\right]$ とおくと

$\begin{aligned} MSE(\hat{e}_{t+h|t})=&E\left[(y_{t+h}-\hat{y}_{t+h|t})^2|\Omega_{t}\right]\\ 　　　　　　　　 =&E\left[(y_{t+h}-\mu_{t+h|t}+\mu_{t+h|t}+\hat{y}_{t+h|t})^2|\Omega_{t}\right]\\ 　　　　　　　　 =&E\left[\left(y_{t+h}-\mu_{t+h|t}\right)^2|\Omega_{t}\right]\\ 　　　　　　　　 &-2E\left[(y_{t+h}-\mu_{t+h|t})(\hat{y}_{t+h|t}-\mu_{t+h|t})|\Omega_{t}\right]\\ 　　　　　　　　 &+E\left[\left(\hat{y}_{t+h|t}-\mu_{t+h|t}\right)^2|\Omega_{t}\right]\\ 　　　　　　　　 =&E\left[\left(y_{t+h}-\mu_{t+h|t}\right)^2|\Omega_{t}\right]+E\left[\left(\mu_{t+h|t}-\hat{y}_{t+h|t}\right)^2|\Omega_{t}\right] \end{aligned}$

が成り立つ。右辺第2項は非負であり、 $MSE(\hat{e}_{t+h|t})$ は $\hat{y}_{t+h|t}=\mu_{t+h|t}$ で最小化できる。以上の考察から、分布に特定の仮定を与えることで計算を可能（ないし容易）にした条件付き期待値を利用するのが妥当であると結論付けることができた。

4.3　AR過程の予測

　 $\mathrm{MSE}$ を基準とすれば、最適予測は条件付き期待値で与えられる。 $\varepsilon_t$ が正規ホワイトノイズだと仮定すると、

$\begin{aligned} E\left[y_t|\Omega_{\tau}\right]&=y_{\tau},&\tau\leq t\\ E\left[\varepsilon_{t+k}|\Omega_{t}\right]&=0,&k\gt0 \end{aligned}$

を用いればよい。

例：AR(1)モデル
　 $\mathrm{AR}(1)$ モデル
$\begin{aligned} y_t=c+\phi_1 y_{t-1}+\varepsilon_t,\ \varepsilon_t\sim i.i.d.N(0,\sigma^2) \end{aligned}$
の最適予測とその性質を考える。
　まず1期先 $t=t+1$ の予測 $y_{t+1}$ を考えると、
$\begin{aligned} E\left[y_t|\Omega_{\tau}\right]&=y_{\tau},&\tau\leq t\\ E\left[\varepsilon_{t+k}|\Omega_{t}\right]&=0,&k\gt0 \end{aligned}$
を用いることで
$\begin{aligned} y_{t+1|t}=E[c+\phi_1 y_t+\varepsilon_{t+1}]=c+\phi_1 y_t \end{aligned}$
で与えられる。またこのときの予測誤差は
$\begin{aligned} \hat{e}_{t+1|t}=y_{t+1}-\hat{y}_{t+1|t}=\varepsilon_{t+1} \end{aligned}$
であり、その $\mathrm{MSE}$ は $\varepsilon_{t}$ が独立であることに注意すれば、
$\begin{aligned} MSE(\hat{y}_{t+1|t})=E[\varepsilon_{t+1}^2]=\sigma^2 \end{aligned}$
を得る。
　次に $2$ 期先予測を考えれば、
$\begin{aligned} y_{t+2}&=c+\phi_1 y_{t+1}+\varepsilon_{t+2}=c+\phi_1(c+\phi_1 y_t)+\varepsilon_{t+2}\\ &=(1+\phi_1)c+\phi_1^2 y_t+\varepsilon_{t+2}+\phi\varepsilon_{t+1} \end{aligned}$
であるから、
$\begin{aligned} \hat{y}_{t+2|t}&=(1+\phi_1)c+\phi_1^2 y_t,\\ MSE(\hat{y}_{t+2|t})&=E\left[\left(y_{t+2}-\hat{y}_{t+2|t}\right)^2\right]=E\left[\left(\varepsilon_{t+2}+\phi_1\varepsilon_{t+1}\right)^2\right]\\ &=\sigma^2+\phi_1^2\sigma^2=(1+\phi_1^2)\sigma^2 \end{aligned}$
を得る。同様にして $h\gt0$ 期先予測を考えると、
$\begin{aligned} y_{t+h}&=(1+\phi_1+\phi_1^2+\cdots+\phi_1^{h-1})+\phi_1^2 y_t+\varepsilon_{t+h}+\phi_1\varepsilon_{t+h-1}+\cdots+\phi_1^{h-1}\varepsilon_{t+1}\\ &=c\displaystyle{\sum_{k=1}^{h}\phi_1^{k-1}}+\phi_1^h y_t+\displaystyle{\sum_{k=0}^{h-1}\phi_1^k\varepsilon_{t+h-k}} \end{aligned}$
と書き直すことができる。したがって
$\begin{aligned} \hat{y}_{t+h|t}&=c\displaystyle{\sum_{k=1}^{h}\phi_1^{k-1}}+\phi_1^h y_t=\displaystyle{\frac{(1-\phi_1^h)c}{1-\phi_1}}+\phi_1^h y_t,\\ MSE(\hat{y}_{t+h|t})&=E\left[\left(\displaystyle{\sum_{k=0}^{h-1}\phi_1^k\varepsilon_{t+h-k}}\right)^2\right]\\ &=\sigma^2\displaystyle{\sum_{k=0}^{h-1}\phi_1^{2k}}\\ &=\displaystyle{\frac{(1-\phi_1^{2h})\sigma^2}{1-\phi_1^2}} \end{aligned}$
である。
　もし $h\rightarrow\infty$ ならば定常条件 $|\phi_1|\lt1$ より
$\begin{aligned} \hat{y}_{t+h|t}\rightarrow\displaystyle{\frac{c}{1-\phi_1}}(h\rightarrow\infty) \end{aligned}$
が得られる。これは予測期間が長くなるにつれて、現在ある情報が徐々に有用でなくなり、その結果、長期的には予測値が過程の期待値とほぼ一致することを意味する。またこの結果は定常 $\mathrm{AR}(1)$ 過程は長期的に過程の平均の方向に戻ることが期待できることを意味する。これを平均回帰的(mean reverting)という。
　 $\mathrm{MSE}$ を考えると、予測期間が長くなればなるほど単調増加するため、予測が困難になる。また $\phi_1^h\rightarrow0(n\rightarrow\infty)$ であるから、 $\mathrm{MSE}$ は予測期間が長くなるにつれて過程の分散
$\begin{aligned} \displaystyle{\frac{\sigma^2}{1-\phi_1^2}} \end{aligned}$
に近づいていく。すなわち定常過程ならば平均回帰的であるために、 $\mathrm{MSE}$ は上限を持ち、それは過程の分散に等しい。

　 $p\gt0$ として一般の定常 $\mathrm{AR}(p)$ 過程

$\begin{aligned} y_t=c+\phi_1 y_{t-1}+\cdots+\phi_p y_{t-p}+\varepsilon_t,\varepsilon_t\sim i.i.d. N(0,\sigma^2) \end{aligned}$

の予測を考える。

$\begin{aligned} \hat{y}_{t+1|t}=c+\phi_1 y_t+\phi_2 y_{t-1}+\cdots+\phi_p y_{t-p+1} \end{aligned}$

で与えられる。また

$\begin{aligned} MSE(\hat{y}_{t+1|t})=E[\varepsilon_{t+1}^2]=\sigma^2 \end{aligned}$

を得る。
　2期先予測は

$\begin{aligned} \hat{y}_{t+2|t}=(1+\phi_1)c+(\phi_1^2+\phi_2)y_t+\cdots+\phi_1\phi_p y_{t-p+1} \end{aligned}$

である。同様にすればより長期の予測もできるが非常に煩雑になる。
　そこで逐次予測を用いるのがよい。このとき

$\begin{aligned} \hat{y}_{t+h|t}=c+\phi_1\hat{y}_{t+h-1|t}+\phi_2\hat{y}_{t+h-2}+\cdots+\phi_p \hat{y}_{t+h-p|t} \end{aligned}$

で得られる。ここで $\tau\leq t$ ならば $\hat{y}_{\tau|t}=y_{\tau}$ である。他方でこの方法では $\mathrm{MSE}$ の計算は困難である。
　証明無しに一般の定常 $\mathrm{AR}(p)$ 過程の性質を述べる。

定常 $\mathrm{AR}(p)$ 過程の性質　定常 $\mathrm{AR}(p)$ 過程( $p\gt0$ )は以下の性質をもつ：

最適予測は直近 $p$ 期の $y,$ すなわち $y_t,y_{t-1},\cdots,y_{t-p+1}$ のみに依存する。
予測期間が長くなるにつれて最適予測における過去の観測値から受ける影響は指数的に減衰し、過程の期待値に収束する。
予測期間が長くなるにつれて $\mathrm{MSE}$ は単調増加し、過程の分散に収束する。

次回

power-of-awareness.com

参考文献

沖本竜義(2010)「経済・ファイナンスデータの計量時系列分析」(朝倉書店)
北川源四郎(2020)「Rによる時系列モデリング入門」（岩波書店）
柴田里程(2017)「時系列解析」(共立出版)
白石博(2022)「時系列データ解析」(森北出版)
萩原淳一郎，瓜生真也，牧山幸史[著]，石田基広[監修](2018)「基礎からわかる時系列分析　Rで実践するカルマンフィルタ・MCMC・粒子フィルタ」(技術評論社)

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

*2:同前掲書　P46

*3:Hamilton, James D.(1994), "Time Series Analysis", Princeton University Press, PP.111-112参照

*4:対数尤度を最尤推定値で評価した値

*5:いずれを用いるかは難しい問題である。

*6:とはいえ実際の変数選択はしらみつぶしになりかねず、手間がかかる。

前回

3. ARMA過程

3.7 AR過程における母数推定

3.7.1 最小二乗法

3.7.2 最尤法

3.8 モデル選択

3.8.1 情報量基準

3.9 モデルの診断

4. 予測

4.1 予測の基本的な考え方

4.2 より厳密な考え方

4.3 AR過程の予測

次回