knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(warn = -1) 

\newcommand{\bs}{\boldsymbol} \newcommand{\mb}{\mathbb} \newcommand{\bm}{\mathbf} \newcommand{\mc}{\mathcal}

1. 들어가며 (Introduction)

bayesVAR 패키지는 베이지안 벡터자기회귀모형에 대한 추론을 위한 패키지이다. 다음과 같이 $k$개의 내생변수와 $p$개의 시차, $m$​개의 외생변수를 갖는 벡터자귀회귀(Vector Autoregressive: VAR) 모형을 고려해 보자. $$ \scriptsize \begin{aligned} \begin{pmatrix} y_{1,t}\ y_{2,t}\ \vdots\ y_{k,t} \end{pmatrix} = & \begin{pmatrix} a_{11}^1&a_{12}^1&\cdots&a_{1k}^1\ a_{21}^1&a_{22}^1&\cdots&a_{2k}^1\ \vdots&\vdots&\ddots&\vdots\ a_{k1}^1&a_{k2}^1&\cdots&a_{kk}^1 \end{pmatrix} \begin{pmatrix} y_{1,t-1}\ y_{2,t-1}\ \vdots\ y_{k,t-1} \end{pmatrix}+\cdots+ \begin{pmatrix} a_{11}^p&a_{12}^p&\cdots&a_{1k}^p\ a_{21}^p&a_{22}^p&\cdots&a_{2k}^p\ \vdots&\vdots&\ddots&\vdots\ a_{k1}^p&a_{k2}^p&\cdots&a_{kk}^p \end{pmatrix} \begin{pmatrix} y_{1,t-p}\ y_{2,t-p}\ \vdots\ y_{k,t-p} \end{pmatrix} +\begin{pmatrix} c_{11}&a_{12}&\cdots&c_{1m}\ c_{21}&a_{22}&\cdots&c_{2m}\ \vdots&\vdots&\ddots&\vdots\ c_{k1}&a_{n2}&\cdots&c_{km} \end{pmatrix} \begin{pmatrix} x_{1,t}\ x_{2,t}\ \vdots\ x_{m,t} \end{pmatrix}+ \begin{pmatrix} \varepsilon_{1,t}\ \varepsilon_{2,t}\ \vdots\ \varepsilon_{k,t} \end{pmatrix} \end{aligned} $$ 보다 간략히 나타내면, $$ \bm y_t = \bm A_1 \bm y_{t-1}+\bm A_2\bm y_{t-2}+\cdots+\bm A_p\bm y_{t-p}+\bm C\bm x_t+\bs\varepsilon_t $$ 여기서 $\bm y_t$ 는 내생변수들로 이루어진 $k\times 1$ 벡터, $\bm A_1,\cdots,\bm A_p$ 는 $k\times k$ 행렬, $\bm x_t$ 는 외생변수들로 이루어진 $m\times 1$ 벡터이고, 오차항 $\bs\varepsilon_t$ 는 각 시점에 대하여 독립이고, 평균이 $\bm 0$ 벡터이고, 분산이 $\bs\Sigma$인 다변량 정규분포를 따른다. 즉, $$ \bs\varepsilon_t\sim \text{MVN}(\bm 0,\bs\Sigma),\quad E(\bs\varepsilon_t\bs\varepsilon_s')=\bm0,\quad t\neq s. $$ 여기서 우리가 추정해야 하는 계수의 수는 총 $k(kp+m)$ 개이다.

(예시 1.1) 내생변수 2개, 시차 2, 외생변수 1개인 VAR 모형

실질GDP(RGDP), 원달러환율(EX)가 내생변수이고, 콜금리(CALL)가 외생변수인 시차 2인 VAR모형을 고려해 보자. $$\small \begin{aligned} RGDP_t &= a^1_{11}RGDP_{t-1}+a^1_{12}EX_{t-1}+a^2_{11}RGDP_{t-2}+a^2_{12}EX_{t-2}+c_{11}+c_{12}CALL_{t}+\varepsilon_{1t}\ EX_t&=a^1_{21}RGDP_{t-1}+a^1_{22}EX_{t-1}+a^2_{21}RGDP_{t-2}+a^2_{22}EX_{t-2}+c_{21}+c_{22}CALL_{t}+\varepsilon_{2t} \end{aligned} $$ 행렬형태로 표현하면, $$\scriptsize \begin{pmatrix}RGDP_t\EX_t\end{pmatrix}=\begin{pmatrix}a_{11}^1&a_{12}^1\a_{21}^1&a_{22}^1\end{pmatrix}\begin{pmatrix}RGDP_{t-1}\EX_{t-1}\end{pmatrix}+\begin{pmatrix}a_{11}^2&a_{12}^2\a_{21}^2&a_{22}^2\end{pmatrix}\begin{pmatrix}RGDP_{t-2}\EX_{t-2}\end{pmatrix}+\begin{pmatrix}c_{11}&c_{12}\c_{22}&c_{22}\end{pmatrix}\begin{pmatrix}1\CALL_{t}\end{pmatrix}+\begin{pmatrix}\varepsilon_{1t}\\varepsilon_{2t}\end{pmatrix} $$

벡터자귀회귀모형의 각 변수에 대하여 $T$ 개의 관측치를 가지고 있다고 가정해 보자. 즉, $\bm y_t,\; t=1,\dots,T$ 를 $t$시점의 관측치 벡터라고 하고, 위 VAR 모형의 기본형태에 전치행렬을 취하여 다음과 같이 표현해 보자. $$ \bm y_t' = \bm y_{t-1}'\bm A_1'+\bm y_{t-2}'\bm A_2'+\cdots+\bm y_{t-p}'\bm A_p'+\bm x_t'\bm C'+\bs\varepsilon_t' $$ 모든 관측치와 계수를 단일 행렬형태로 나타내면 다음과 같다. $$\scriptsize \begin{pmatrix} \bm y_1'\ \bm y_2'\ \vdots\ \bm y_T' \end{pmatrix}= \begin{pmatrix} \bm y_0'& \bm y_{-1}'&\cdots &\bm y_{1-p}' &\bm x_1'\ \bm y_1'&\bm y_0'&\cdots &\bm y_{2-p}'&\bm x_2'\ \vdots&\vdots&&\vdots&\vdots\ \bm y_{T-1}'&\bm y_{T-2}'&\cdots&\bm y_{T-p}'&\bm x_T'\end{pmatrix} \begin{pmatrix} \bm A_1'\\bm A_2'\\vdots\\bm A_p'\\bm C'\end{pmatrix}+\begin{pmatrix}\bs\varepsilon_1'\\bs\varepsilon_2'\\vdots\\bs\varepsilon_T'\end{pmatrix} $$ 보다 단순하게 표현하면 $$ \mathbf{Y}=\mathbf{X}\mathbf{B}+\mathcal{E}. $$ 자귀회귀모형의 계수 행렬 $\mathbf B$ 와 공분산 행렬 $\Sigma$ 의 최소제곱추정량은 다음과 같이 주어진다. $$ \hat{\mathbf B}=(\mathbf X^T \mathbf X)^{-1}\mathbf X^T \mathbf Y,\quad \hat{\mathcal E} = \bm Y -\bm X\hat{\bm B},\quad \hat{\bs\Sigma}= \frac{1}{T-k(kp+m)-1}\hat{\mathcal E}^T \hat{\mathcal E} $$

(예시 1.2) 실질GDP와 환율 VAR 모형의 OLS 추정

2001년 1분기 ~ 2020년 4분기 기간중 분기별 실질GDP와 원달러환율, 콜금리 실제 데이터를 이용하여 총 10개(2$\times$(2$\times$2+1)) 개의 회귀계수와 공분산을 추정해 보자. 데이터를 행렬식 $\mathbf{Y}=\mathbf{X}\mathbf{B}+\mathcal{E}$ 로 나타내고 최소제곱추정법을 이용하여 추정량 $\hat {\bm A}_1$, $\hat {\bm A}_2$, $\hat {\bm C}$ 및 $\hat{\bs\Sigma}$를 각각 구하면 $$\scriptsize\begin{aligned} \hat{\bm A}_1=\begin{pmatrix}0.911 & -0.102\-0.004&1.168\end{pmatrix},&\quad \hat{\bm A}_2=\begin{pmatrix}-0.254 & 2.743\0.0054&-0.323\end{pmatrix},\ \hat{\bm C}=\begin{pmatrix}-1.540&162.480 \0.2939& 0.213\end{pmatrix},&\quad \hat{\bs\Sigma}=\begin{pmatrix}1.785 & -25.025\-25.025& 2423.148\end{pmatrix} \end{aligned}$$ 과 같다. 다음 그림은 실제 시계열 데이터(흑색 실선)와 추정한 회귀계수를 이용하여 실질 GDP와 환율을 적합한 시계열(적색 점선)을 보여주고 있다.

VAR모형의 또 다른 표현 방법은 겹층 형태(stacked form)의 벡터와 행렬을 이용하는 방법이다. 즉, $$\scriptsize \begin{pmatrix}y_{1,1}\\vdots\y_{1,T}\\vdots\y_{k,1}\\vdots\y_{k,T}\end{pmatrix}= \begin{pmatrix} y_0'&y_{-1}'&\cdots &y_{1-p}' &x_1'&&0&&\cdots&&0\ y_1'&y_0'&\cdots &y_{2-p}'&x_2'&&&&&&\ \vdots&\vdots&&\vdots&\vdots&&\vdots&&\cdots&&\vdots\ y_{T-1}'&y_{T-2}'&\cdots&y_{T-p}'&x_T'&&0&&\cdots&&0&\ &&&&&\ddots&&&\ 0&&\cdots&&0&&y_0'&y_{-1}'&\cdots &y_{1-p}' &x_1'\ &&&&&&y_1'&y_0'&\cdots &y_{2-p}'&x_2'\ \vdots&&\cdots&&\vdots&&\vdots&\vdots&&\vdots&\vdots\ 0&&\cdots&&0&&y_{T-1}'&y_{T-2}'&\cdots&y_{T-p}'&x_T' \end{pmatrix} \begin{pmatrix}A_{1}^{(1)}\\vdots\A_{p}^{(1)}\C^{(1)}\\vdots\A_{1}^{(k)}\\vdots\A_{p}^{(k)}\C^{(k)}\end{pmatrix}+ \begin{pmatrix}\varepsilon_{1,1}\\vdots\\varepsilon_{1,T}\\vdots\\varepsilon_{k,1}\\vdots\\varepsilon_{k,T}\end{pmatrix} $$ 여기서, $A_i^{(j)}$ 와 $C^{(j)}$는 행렬 $\bm A_i$와 $\bm C$의 $j$번째 행의 전치벡터이다. 이를 간단히 표현하면 $$ \mathbf y_s= \mathbf X_s\boldsymbol \beta + \boldsymbol\varepsilon $$ 여기서, $\bm y_s=vec(\bm Y)$, $\bm X_s= \bm I_k\otimes\bm X$, $\bs\beta=vec(\bm B)$이고 $\otimes$는 크로넥커곱을 나타낸다. 한편, 오차항은 $\bs\varepsilon=vec(\mathcal E)$이고,
$$ \bs\varepsilon\sim \text{N}(\bm 0, \bs\Sigma_s),\quad \bs\Sigma_s=\bs\Sigma\otimes \bm I_T $$ 최소제곱추정량 $\hat{\bs\beta}$는 다음과 같이 구할 수 있다. $$ \hat{\bs\beta}=(\bm X_s^T\bm X_s)^{-1}\bm X_s^T\bm y_s $$

2. 베이지안 접근법 (Bayesian Inference)

VAR모형의 관측치 $\bm y_s$는 $\bm y_s\sim \text N(\bs\beta,\bs\Sigma_s)$ 를 따르고, 이에 대한 우도함수는 다음과 같다. $$ L(\bm y\mid \bs\beta,\bs\Sigma_s)=(2\pi)^{-nT/2}|\bs\Sigma_s|^{-1/2}\exp\left[-\frac{1}{2}(\bm y -\bm X_s\bs\beta)^T\bs\Sigma_s^{-1}(\bm y -\bm X_s \bs\beta)\right] $$ 여기서 정규분포에 대한 파라미터는 $\bs\beta$ 와 $\bs\Sigma_s$ 이다. 우선 $\bs\Sigma_s$는 주어진 값으로 보고 파라미터 $\bs\beta$의 사후분포를 구하는 방법에 대하여 알아보고, 그 다음 두개의 파라미터를 모두 확률변수로 간주하여 사후분포를 구하는 법에 대하여 알아보자.

2.1 경험 베이즈 접근법과 Minnesota 사전분포

먼저 $\bs\Sigma_s$ 를 알고 있다고 가정하는 경우 사용할 $\bs\Sigma_s$의 추정치 $\hat{\bs\Sigma}_s$ 을 구해보자. 가장 일반적인 접근법은 위에서 고려한 VAR모형의 최소제곱추정법(OLS) 잔차항으로 $\bs\Sigma_s$을 추정하는 것이다. 즉, $$ \hat{\bs\Sigma}_s=\frac{(\bm y - \bm X_s\hat{\bs\beta})^T(\bm y - \bm X_s\hat{\bs\beta})}{T-k(kp+m)-1}=\frac{\hat{\bs\varepsilon}^T\hat{\bs\varepsilon}}{T-k(kp+m)-1}. $$ 이제 확률변수 $\bs\beta$ 의 사후분포 $\pi(\bs\beta\mid \bm y)$ 를 구하기 위해서는 $\bs\beta$의 사전분포를 부여하여야 한다. $\bs\beta$의 사전분포로 다음과 같은 다변량 정규분포로 가정해 보자. 즉, $$ \begin{aligned} \bs\beta&\sim\text{MVN}(\bs\beta_0,\bs\Omega_0)\ p(\bs\beta) &\propto \exp \left[-\frac{1}{2}(\bs\beta-\bs\beta_0)^T\bs\Omega_0^{-1}(\bs\beta-\bs\beta_0)\right] \end{aligned} $$ 정규분포인 $\bs\beta$의 사전분포에서 우리는 평균과 공분산 두개의 하이퍼 파라미터 $\bs\beta_0,\bs\Omega_0$를 갖게 되는데, Littermann(1986)은 두개의 하이퍼 파라미터에 대하여 MInnesota 파라미터 가정을 제시하였다.

먼저 $\bs\beta_0$는 대부분의 거시변수들은 단위근을 가지고 있기 때문에 각각의 내생변수들은 자신의 첫번째 시차에는 단위근을 갖고, 이후 시차나 다른 변수의 계수는 모두 0으로 만든다. 또한, 외생변수에 대해서도 변수의 영향을 중립적인 것으로 가정하여 마찬가지고 모든 계수를 0으로 놓는다. 예를 들어 내생변수 2개, 시차 2, 외생변수 1개인 VAR 모형에서 $\bs\beta_0$는 $$\small \bs\beta_0=\begin{pmatrix}1\0\0\0\0\-\0\1\0\0\0\end{pmatrix} $$ $k$개의 내생변수와 $p$개 시차, $m$개의 외생변수를 가진 모형의 경우 $\bs\beta_0$는 $k(kp+m)\times 1$ 벡터의 형태가 된다.

공분산 행렬 $\bs\Omega_0$에 대해서는 $\bs\beta$ 의 서로 다른 변수 간에는 공분산이 0이고, 같은 변수라도 시차가 증가할 수록 분산은 작아진다고 가정 하였다. 한편, 외생변수에 대해서는 정보가 없기 때문에 분산에 매우 큰 값을 부여하였다. 이를 정리하면

  1. 내생변수의 자기자신의 시차에 대해서 분산은 $$ \sigma_{a_{ii}}^2=\left(\frac{\lambda_1}{l^{\lambda_3}}\right) $$ 여기서 $l$은 시차를 의미하며, $\lambda_1$은 모형의 전반적인 분산을, $\lambda_3$는 시차가 1보다 커질 때 분산이 0으로 수렴하는 속도를 조절하는 파라미터이다.

  2. 다른 변수의 시차에 대하여 분산은 $$ \sigma_{a_{ij}}^2=\left(\frac{\sigma_i^2}{\sigma_j^2}\right)\left(\frac{\lambda_1\lambda_2}{l^{\lambda_3}}\right)^2 $$ 여기서 $\sigma_i^2$ 와 $\sigma_j^2$ 는 $i$와 $j$ 변수의 OLS 자귀회귀모형(AR)모형 잔차의 분산을 나타내고, $\lambda_2$는 다른 변수의 분산 크기를 조절하는 파라미터이다.

  3. 외생변수의 분산은 $$ \sigma_{c_i}=\sigma_i^2(\lambda_1\lambda_4)^2 $$ 여기서 $\lambda_4$는 무한대 또는 매우 큰 값을 갖는 파라미터이다.

예를 들어 내생변수 2개, 시차 2, 외생변수 1개인 VAR 모형에서 공분산 행렬 $\bs\Omega_0$ 는 $$\scriptsize \bs\Omega_0=\begin{pmatrix} \lambda_1^2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\ 0 & \left(\frac{\sigma_1^2}{\sigma_2^2}\right)(\lambda_1\lambda_2)^2&0 &0 & 0 & 0 & 0 & 0 & 0 & 0\ 0& 0& \left(\frac{\lambda_1}{2^{\lambda_3}}\right)^2& 0 & 0 & 0 & 0 & 0 & 0 & 0\ 0& 0& 0 &\left(\frac{\sigma_1^2}{\sigma_2^2}\right)\left(\frac{\lambda_1\lambda_2}{2^{\lambda_3}}\right)^2 & 0 & 0& 0 & 0 & 0 & 0\ 0& 0 & 0 & 0 & \sigma_1^2(\lambda_1\lambda_4)^2& 0 & 0 & 0 & 0 & 0\ 0 & 0 & 0 & 0 & 0& \left(\frac{\sigma_2^2}{\sigma_1^2}\right)(\lambda_1\lambda_2)^2&0 &0 & 0 & 0 \ 0 & 0 &0 & 0 & 0 & 0 &\lambda_1^2 & 0 & 0 & 0 \ 0& 0 & 0 & 0 & 0&0& 0& \left(\frac{\sigma_2^2}{\sigma_1^2}\right)\left(\frac{\lambda_1\lambda_2}{2^{\lambda_3}}\right)^2 & 0 &0\ 0 & 0 & 0 & 0 & 0 & 0 & 0& 0& \left(\frac{\lambda_1}{2^{\lambda_3}}\right)^2& 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma_2^2(\lambda_1\lambda_4)^2\ \end{pmatrix} $$

파라미터 $\lambda_1,\;\lambda_2,\;\lambda_3,\;\lambda_4$ 에 대하여 보통 다음과 같은 값을 이용한다. $$ \lambda_1=0.1,\quad \lambda_2=1/2,\quad \lambda_3= 1\;\text{or}\; 2,\quad \lambda_4=10^2 $$ 이제 $\bs\beta$의 사후분포 $\pi(\bs\beta\mid \bm y)$ 는 다음과 우도함수와 사전분포의 곱으로 $$\small \begin{aligned} \pi(\bs\beta\mid \bm y) &\propto L(\bm y\mid \bs\beta)p(\bs\beta)\ &\propto\exp\left[-\frac{1}{2}\left{(\bm y-\bm X_s\bs\beta)^T\hat{\bs\Sigma}_s^{-1}(\bm y-\bm X_s \bs\beta)+(\bs\beta-\bs\beta_0)^{-1}\bs\Omega_0^{-1}(\bs\beta-\bs\beta_0)\right}\right]\ &\propto \exp\left[-\frac{1}{2}\left{(\bs\beta-\tilde{\bs\beta_0})^{-1}\tilde{\bs\Omega}_0^{-1}(\bs\beta-\tilde{\bs\beta_0})\right}\right] \end{aligned} $$ 여기서, $\tilde{\bs\Omega}_0=[\bs\Omega_0^{-1}+\bs\Sigma^{-1}\otimes\bm X^T\bm X]^{-1}$이고, $\tilde{\bs\beta}_0=\tilde{\bs\Omega}_0[\;\bs\Omega_0^{-1}\bs\beta_0+(\bs\Sigma\otimes\bm X^{-1})\bm y\;]$ 이다. 최종적으로 $\bs\beta$ 의 사후분포는 정규분포로 다음과 같이 정리할 수 있다. $$ \bs\beta\mid \bm y \sim \text{N}(\tilde{\bs\beta}_0,\tilde{\bs\Omega}_0) $$

2.2. 완전베이즈 접근법: Normal-Invese Wishart 사전분포

이번 절에서는 파라미터 $\bs\beta$ 와 $\bs\Sigma_s$ 모두 확률변수로 간주하고, $\bs\beta$와 $\bs\Sigma_s$의 사후분포를 구해보자. 여기서 문제는 공분산 행렬 $\Sigma_s$의 차원이 너무 크기 때문에 $\Sigma_s$를 직접 추출하는 하는 것보다 차원이 작은 $\bs\Sigma$를 이용하는 방법을 고려한다. 따라서, 우도함수 $L(\bm y\mid \bs\beta,\bs\Sigma)$를 파라미터 $\bs\beta$와 $\bs\Sigma$로 다시 나타내면 $$\small \begin{aligned} L(\bm y\mid \bs\beta,\bs\Sigma)&\propto |\bs\Sigma|^{-(kp+m)/2}\exp\left[-\frac{1}{2}(\bs\beta -\hat{\bs\beta})^T(\bs\Sigma\otimes(\bm X ^T \bm X)^{-1})^{-1}(\bs\beta -\hat{\bs\beta})\right]\ &\times|\bs\Sigma|^{-[(T-(kp+m)-k-1)+k+1]/2}\exp\left[-\frac{1}{2}\text{tr}\left{\bs\Sigma^{-1}(\bm Y-\bm X \hat{\bm B})^T(\bm Y-\bm X \hat{\bm B})\right}\right]\ &\propto |\bs\Sigma|^{-T/2}\exp\left[-\frac{1}{2}(\bs\beta -\hat{\bs\beta})^T(\bs\Sigma\otimes(\bm X ^T \bm X)^{-1})^{-1}(\bs\beta -\hat{\bs\beta})\right]\times\exp\left[-\frac{1}{2}\text{tr}\left{\bs\Sigma^{-1}(\bm Y-\bm X \hat{\bm B})^T(\bm Y-\bm X \hat{\bm B})\right}\right]\ \end{aligned} $$ 파라미터 $\bs\beta$에 대하여 Minnesota 정규 사전분포 $\text{N}(\bs\beta_0, \bs\Omega_0^{-1})$ 를 부여하고, 파라미터 $\bs\Sigma$ 대하여 사전분포로 역위샤트분포 $\text{IW}(S_0,\alpha_0)$를 부여하기로 하자. 역위샤트분포 $\text{IW}(S_0,\alpha_0)$의 밀도함수는 다음과 같다. $$ p(\bs\Sigma)\propto |\bs\Sigma|^{-(\alpha_0+n+1)/2}\exp\left{-\frac{1}{2}\text{tr}(\bs\Sigma^{-1}S_0)\right} $$ 결합사후분포 $\pi(\bs\beta,\bs\Sigma\mid \bm y)$ 는 우도함수에 두개의 사전분포의 밀도함수를 곱한 형태에 비례하게 된다. 즉, $$\small \begin{aligned} \pi(\bs\beta,\bs\Sigma\mid \bm y)\propto & L(\bm y \mid \bs\beta,\bs\Sigma)\cdot p(\bs\beta)\cdot p(\bs\Sigma)\ \propto&|\bs\Sigma|^{-T/2}\exp\left[-\frac{1}{2}(\bs\beta -\hat{\bs\beta})^T(\bs\Sigma\otimes(\bm X ^T \bm X)^{-1})^{-1}(\bs\beta -\hat{\bs\beta})\right]\times\exp\left[-\frac{1}{2}\text{tr}\left{\bs\Sigma^{-1}(\bm Y-\bm X \hat{\bm B})^T(\bm Y-\bm X \hat{\bm B})\right}\right]\ &\times \exp\left[-\frac{1}{2}(\bs\beta-\bs\beta_0)^T\bs\Omega_0^{-1}(\bs\beta-\bs\beta_0)\right]\times|\bs\Sigma|^{-(\alpha_0+n+1)/2}\exp\left{-\frac{1}{2}\text{tr}(\bs\Sigma^{-1}S_0)\right} \end{aligned} $$

결합사후분포의 형태로부터 $\bs\beta$의 완전조건부 사후분포 $\pi(\bs\beta\mid \bs\Sigma,\bm y)$는 경험베이즈 추정과 같이 $$ \bs\beta\mid \bs\Sigma, \bm y \sim \text{N}(\tilde{\bs\beta}_0,\tilde{\bs\Omega}_0) $$ 여기서 $\tilde{\bs\Omega}_0=[\bs\Omega_0^{-1}+\bs\Sigma^{-1}\otimes\bm X^T\bm X]^{-1}$이고, $\tilde{\bs\beta}_0=\tilde{\bs\Omega}_0[\;\bs\Omega_0^{-1}\bs\beta_0+(\bs\Sigma\otimes\bm X^{-1})\bm y\;]$ 이다.

$\bs\Sigma$의 완전조건부 사후분포 $\pi(\bs\Sigma\mid \bs\beta,\bm y)$는 $$\small \pi(\bs\Sigma\mid \bs\beta,\bm y)\propto|\bs\Sigma|^{-(\alpha_0+n+1)/2}\exp\left[-\frac{1}{2}\text{tr}\left[ \bs\Sigma^{-1}\left{(\bm Y-\bm X \hat{\bm B})^T(\bm Y-\bm X \hat{\bm B})+S_0\right}\right]\right] $$ 로 주어지고, 역위샤트 분포의 커널형태가 되어 파라미터 $\hat{\bm S}$ 와 $\hat\alpha$ 를 갖는 역쉬샤트 분포를 샘플링 할 수 있다. $$ \bs\Sigma\mid \bs\beta, \bm y \sim \text{IW}(\hat{\bm S}, \hat\alpha) $$ 여기서, $\hat{\bm S}=(\bm Y-\bm X\hat{\bm B})^T(\bm Y-\bm X\hat{\bm B})+\bm S_0$ 이고, $\hat\alpha=T+\alpha_0$이다. 위 두 완전조건부 사후분포를 이용한 Gibbs 샘플링 추출방법은 다음과 같다.

(Algorithm 1) 회귀계수 $\bs\beta$ 와 공분산 $\bs\Sigma$ 의 Gibbs 샘플링

  1. 최소제곱추정량 $\hat{\bs\beta}$ 와 $\hat{\bs\Sigma}$를 초기값으로 놓는다.
  2. 현재 $\bs\Sigma^{(c)}$ 를 이용하여 완전 조건부분포 $\bs\beta\mid \bs\Sigma,\bm y\sim\text N(\tilde{\bs\beta}_0,\tilde{\bs\Omega}_0)$ 에서 $\bs\beta$값을 샘플링한다.
  3. 현재 $\bs\beta^{(c)}$ 를 이용하여 완전 조건부분포 $\bs\Sigma\mid \bs\beta, \bm y \sim \text{IW}(\hat{\bm S}, \hat\alpha)$ 에서 $\bs\Sigma$값을 샘플링한다.
  4. 반복 추출을 통하여 분포가 목표 분포에 충분히 수렴하면 Monte Carlo 기법을 이용하여 파라미터의 분포를 구한다.

3. 예시 (Example)

bayesVAR 패키지는 거시변수 데이터 $macrodata$를 제공하는데, 실질GDP성장율(RGDP_Q), 원달러환율변화율(USDKRW_G) 등 2000년 3분기부터 2021년 2분기까지 주요 분기별 자료를 $\texttt{xts}$ 시계열 형태로 제공한다.

suppressMessages(library(bayesVAR))
data(macrodata)
class(macrodata)
tail(macrodata)

이제 실질GDP증가율(RGDP_G), 원달러환율변화율(USDKRW_G)이 내생변수이고, 콜금리(INT_CALL)가 외생변수인 시차 $p=4$인 VAR모형을 bayesVAR 패키지를 이용하여 추론해 보자. 세 거시변수의 움직임은 다음과 같다.

data <- macrodata[,c("RGDP_Q_G","USDKRW_G","INT_CALL")]
plot(data,multi.panel=TRUE,yaxis.same = FALSE)

$VAR_bayes$함수는 bayesVAR 패키지에서 베이지안 추론을 위한 기본 함수이며, Minnesota 사전분포와 Normal-Wishart 사전분포를 결합하여 사후분포를 Gibbs 샘플링 기법을 이용하여 추정한다. 추정 결과로 $VAR_bayes$함수는 $bayesVAR$ 객체를 제공하는데, 이를 개체를 이용하여 기본적인 파라미터에 대한 추론, (조건부) 예측 및 충격-반응(impulse-response) 함수 등을 구할 수 있다. 여기서 $p$는 시차, $N$은 Gibbs sampling의 표본수, $waramup$은 분포의 수렴을 위해 최종 표본에서 제외할 표본수, $exos$는 외생변수, $minnesota_par$는 Minnesota 사전분포의 각 파라미터 값을 의미한다. $bayesVAR$ 객체에 $plot$ 함수를 적용하면 추론 결과 얻어진 각 내생변수의 사후분포 평균과 95% 확률구간을 보여준다.

# Bayesian VAR
object <- VAR_bayes(data,p=4, N=150, warmup = 50, exos="INT_CALL",
                   minnesota_par=list(rho=0.8, lambdas = c(0.1, 0.5, 2, 100)))
class(object)
plot(object, ncol.fig = 2)

VAR모형은 주로 관측되지 않은 시계열에 대한 예측을 위하여 이용된다. 미레 시점에 대한 예측은 크게 외생변수가 주어졌을 때 내생변수 전체를 예측하는 문제(unconditional prediction)와 내생변수 중 일부에 대하여 특정 값을 전제하고 다른 내생변수들의 움직임을 예측하는 문제(conditional prediction)로 나눌 수 있다. 예측에 사용되는 함수는 일반적인 $predict$ 함수를 이용한다. 이 때 bayesVAR 패키지에서는 외생변수의 미래값과 내셩변수의 미래 조건값 모두 $\texttt{xts}$ 형태로 입력해야 한다. 예측의 결과는 각 내생변수의 미래 시점의 사후분포를 제공하고, $plot$함수를 적용하면 예측시점의 사후분포 평균 및 95% 확률구간을 제공한다.

# unconditional forecasting
# newdata: exogenous variable
newdata <- xts(data.frame(INT_CALL=c(0.81,0.91,1.0,1.2)),
               order.by=as.yearqtr(seq(2021.5,2022.25,by=0.25)))
pred2 <- bayesVAR:::predict.bayesVAR(object, newdata=newdata)
plot(pred2,ncol.fig=2)
pred2$fitted

조건부 예측은 무조건부 예측과 유사하는 아래와 같이 특정 시계열이 이용자가 주어진 경로를 따라 움직일 때, 다른 변수들의 움직을 추정한다는 점에서 차이가 있다.

# conditional forecasting
condition <- xts(data.frame(RGDP_Q_G = c(4.3,0.5,-4.5,-2.3)),
                 order.by=as.yearqtr(seq(2021.5,2022.25,by=0.25)))

pred3 <- predict(object, newdata=newdata, condition=condition)
plot(pred3,ncol.fig=2)
pred3$fitted


morner75/bayesVAR documentation built on April 6, 2024, 7:25 a.m.