Bayesian analysis of seasonally cointegrated VAR model

Bayesian analysis of seasonally cointegrated VAR model

\thispagestyle{fancy} \fancyhead

Keywords: seasonal cointegration, reduced rank regression, error correction model, Bayesian analysis, Bayesian model comparison;

JEL Classification: C11, C32, C53;

Introduction

Many macroeconomic quarterly or monthly time series display strong both trend and seasonal behavior. The idea of cointegration at zero frequency (assuming common stochastic trend controlling the long-run behavior of the series) is well known and oft employed in the empirical analyses. Methods enabling estimation of parameters of vector error correction models are well known within both classical (see e.g. Johansen1995) and Bayesian paradigm (see e.g. Koop_al2006). However, although the idea of cointegration at seasonal frequencies was introduced in 1990 (Hylleberg_al1990) it got much fewer interests. In most researches, the seasonally adjusted data are analyzed or the seasonality is modeled via seasonal dummies. However, there are papers presenting the unwilling results of seasonal adjustment. Such procedures may for example change both the short- and long-run behavior of the series (e.g. Nerlove1964, Cubadda1999, Hecq1998, GrangerSiklos1995).\footnote{Note that the very same problems appear in the case of filtering non-seasonal unit roots (see e.g. MeyerWinker2005, Hamilton2018).} Abeysinghe1994 discusses the problems which may occur in researchers where seasonality is modeled by seasonal dummies in cases when such behavior is generated by seasonal unit roots. Moreover, LofFranses2001 shows that taking into account seasonal cointegration improves forecast abilities of models.\\ The likelihood method enabling estimation of seasonal VEC models was proposed by Johansen_Schaumburg1999 (see also Cubadda_Omtzigt2005 for further discussion), but to the best of our knowledge, there are no Bayesian seasonal VEC models. This paper is aimed to fulfill this gap by building a Bayesian model for the quarterly seasonally cointegrated time series. Priors for the parameters are imposed. The set of full conditional posterior distributions is obtained. Additionally, the point estimation of the cointegrating spaces is discussed. The proposed methods are illustrated by a small simulation experiment and the empirical analysis of money and prices in the Polish economy conducted in the four-dimensional seasonal system.

Bayesian seasonal VEC model

We start the analysis with the assumption that the $n-$dimensional quarterly time series $\{y_t\}$ has the following VAR($k$) representation:

where $D_t$ contains deterministic components such as a constant, a trend or seasonal dummies. The initial conditions - $y_0, y_{-1}, \dots, y_{-k+1}$ - are fixed, $A(L)=I_n-A_1L-\dots-A_kL^k$ is the polynomial matrix of the process ((ref)).\\ According to Lagrange expansion the polynomial $A(z)$ around the points $z_1,z_2,\dots,z_S$ may be written as (Johansen_Schaumburg1999, Kotlowski2005):

where $p(z)=\prod_{s=1}^S(1-\bar{z}_sz)$, $p_j(z)=\prod_{s\neq j}^S(1-\bar{z}_sz)=\frac{p(z)}{1-\bar{z}_jz},\ z\neq z_j$ and $A_0(z)$ is a polynomial matrix. If additionally $z_s$ is a root of the characteristic polynomial of the process ((ref)), i.e. $|A(z_s)|=0$ then $A(z_s)$ is of the reduced rank and can be decomposed as the product of full column rank matrices $A(z_s)=a_sb_s'$.\\ The expansion ((ref)) leads for processes cointegrated at zero and quarterly frequencies (i.e. $|A(z)|=0$ for $z_1=1$, $z_2=-1$, $z_3=i$, $z_4=\bar{z}_3=-i$) to the following seasonal error correction representation (see e.g. Hylleberg_al1990, Johansen_Schaumburg1999, Cubadda_Omtzigt2005, Kotlowski2005):

where $\Delta_4y_t=(1-L^4)y_t=y_t-y_{t-4}$ with $L$ denoting the lag operator, $\Gamma_j=-\sum_{l=1}^{[(k-j)/4]}A_{j+4l},\ j=1,2,\dots,k-4$. Vectors $\tilde{y}_t^{(\cdot)}$ may contain deterministic components, so $\tilde{y}_t^{(\cdot)}=\left(

\right)$. In particular, let consider the deterministic component of the following form $\Phi D_t=\mu+a\cos\left(\frac{\pi}{2}t\right)+b\sin\left(\frac{\pi}{2}t\right)+c\cos\left(\pi t\right)+\gamma t$ (see \citealp{Franses_Kunst1999} and \citealp{Kotlowski2005}).\\ Employing the polynomials $p_{\cdot}(L)$ to the stochastic and deterministic parts of the process (\ref{eqn:VAR}) leads to the vectors $\tilde{y}_t^{(\cdot)}$ of the following forms:

Note that the unrestricted constant occurring in each of the above considered frequencies results form the linear trend assumed for the level of the analyzed process ((ref)). If there is no linear trend, but only constant, there is neither trend restricted to the cointegration space at the zero frequency nor an unrestricted constant in the representation ((ref)), but there is a constant restricted to the cointegration spaces of the zero frequency ($d_t^{(1)}=1$), see e.g. Juselius2006 for the discussion of the meaning of dummies gathered in the vectors $D_t$ and $\tilde{D}_t$, and also the relations between them.\\ Note also that $\alpha_{\star}\bar{\beta}_{\star}'\tilde{y}_t^{(3)}$ and $\bar{\alpha}_{\star}\beta_{\star}'\bar{\tilde{y}}_t^{(3)}$ are complex conjugate matrices, so their sum gives their real part multiplied by 2:

where $\alpha_R$, $\beta_R$ denote the real parts of $\alpha_{\star}$ and $\beta_{\star}$ respectively, whereas $\alpha_I$, $\beta_I$ - their imaginary parts ($\alpha_{\star}=\alpha_R+i\alpha_I$, $\beta_{\star}=\beta_R+i\beta_I$).\\ These leads to the more commonly used representation of seasonally cointegrated quarterly VAR process (see e.g. Hylleberg_al1990, Johansen_Schaumburg1999, Cubadda_Omtzigt2005, Kotlowski2005):

where $\Pi_1=\alpha_1\beta_1'$, $\Pi_2=\alpha_2\beta_2'$, $\Pi_3=-2(\alpha_R\beta_R'+\alpha_I\beta_I')$, $\Pi_4=2(\alpha_I\beta_R'-\alpha_R\beta_I')$, $\tilde{y}_t^{(31)}=\left(

\right)$, where $y_t^{(31)}=(1-L^2)Ly_t=y_{t-1}-y_{t-3}$, $d_t^{(31)}=sin(\frac{\pi t}{2})$, $\tilde{y}_t^{(32)}=\left(

\right)$, where $y_t^{(32)}=(1-L^2)L^2y_t=y_{t-2}-y_{t-3}$, $d_t^{(32)}=cos(\frac{\pi t}{2})$.

To save on notation we introduce the matrix form of the model ((ref)).

where $Z_0=\left(

\right)'$, $Z_1=\left(

\right)'$, $Z_2=\left(

\right)'$, $Z_3=\left(

\right)'=-Z_{32}-iZ_{31}$, $Z_{31}=\left(

\right)'$, $Z_{32}=\left(

\right)'$, $Z_4=\left(

\right)'$, $z_t'=\left(

\right)$, $\Gamma=\left(

\right)'$, $E=\left(

\right)'$.\\ As the analyzed data inform only about the cointegration space not the cointegration vectors, during the estimation we have to deal with the non-identification occurring in the products: $\alpha_1\beta_1'$, $\alpha_2\beta_2'$, $\bar{\beta}_{\star}\alpha_{\star}'$, $\beta_{\star}\bar{\alpha}_{\star}'$. Therefore, we employ the methods proposed by \citet{Koop_al2009}. Following their ideas we will consider two observationally equivalent representation for each of the considered products. In the $A-B$ representations it is assumed that the matrices belong to the $\mathbb{R}^{\cdot}$ or $\mathbb{C}^{\cdot}$ spaces of appropriate dimensions, whereas in the $\alpha-\beta$ representations $\beta$s have orthonormal columns and $\alpha$s still belong tho the $\mathbb{R}^{\cdot}$ or $\mathbb{C}^{\cdot}$ spaces:

where $\mathbb{V}_{r_j,m_j},\ j=1,2$ denotes the Stiefel manifold, i.e. the set of $m_j\times r_j$ matrices with orthonormal columns ($\mathbb{V}_{r_j,m_j}=\left\{X\in\mathbb{R}^{m_j\times r_j}: X'X=I_{r_j}\right\}$ ), $\mathbb{V^C}_{r_3,m_3}$ stands for the complex Stiefel manifold, i.e. the set of $m_3\times r_3$ semi-unitary matrices ($\mathbb{V^C}_{r_3,m_3}=\left\{X\in\mathbb{C}^{m_3\times r_3}: \bar{X}'X=I_{r_3}\right\}$).\\ Note that by this approach the non-identification issue is only partially solved as there is many-to-one relationship between the Stiefel manifolds and the Grassmann manifolds \footnote{$\mathbb{G}_{r_j,m_j-r_j}, j=1,2,\ \mathbb{G^C}_{r_3,m_3-r_3}$, collecting $r_j$-dimensional planes, passing through the origin, in the real ($j=1,2$)/complex ($j=3$) vector $m_j$-dimensional space, see e.g. James1954, Chern_Wolfson1987} to which belong the cointegration spaces: if $X$ is the element of the (complex) Stiefel manifold and the $r_j\times r_j$ ($j=1,2,3$) matrix $O$ is the element of the (complex) orthonormal group ($O'O=OO'=I_{r_j},\ j=1,2$, $\bar{O}'O=O\bar{O}'=I_{r_3}$) than $XO$ is the element of the same (complex) Stiefel manifold and they span the same spaces (the projection matrices are equal, i.e. $XX'=XOO'X'$ in the real case and $X\bar{X}'=XO\bar{O}'\bar{X}'$ in the complex case).\\ The first two products (i.e. $\alpha_1\beta_1$ and $\alpha_2\beta_2$) involve only matrices with real numbers, so can be treated exactly as proposed by Koop_al2009. While in the third case we have to adjust their procedures to the complex matrices and spaces.\\ We start the analysis with the $A-B$ parameterization and impose the following prior distributions over the model parameters:

The joint prior distribution is truncated by the non-explosive condition taking into account the appropriate numbers of unit roots at each frequency.

The above proposed prior distributions together with the likelihood function lead to the joint posterior distribution with the following kernel:

where $\theta=(\Sigma,\Gamma,A_1,A_2,A_{\star},B_1,B_2,B_{\star})$ collects all model's parameters, $l$ denotes the number of deterministic components gathered in $\tilde{D}_t$, $I_{[a,b]}(\cdot)$ is the indicator function of the interval $[a,b]$ and $\lambda$ stands for the vector of the eigenvalues of the companion matrix, that is the matrix of the form:

where $I_n$ is an $n$-dimensional identity matrix, $A_1=\Pi_1+\Pi_2+\Pi_3+\Gamma_1$, $A_2=\Pi_1-\Pi_2+\Pi_4+\Gamma_2$, $A_3=\Pi_1+\Pi_2-\Pi_3+\Gamma_3$, $A_4=I_n+\Pi_1-\Pi_2-\Pi_4+\Gamma_4$, $A_i=\Gamma_i-\Gamma_{i-4}$ for $i=5,6,\dots$, $\Pi_1=\alpha_1\beta_1'\equiv A_1B_1'$, $\Pi_2=\alpha_2\beta_2'\equiv A_1B_1'$, $\Pi_3=2(\alpha_I\beta_R'-\alpha_R\beta_I')\equiv2(A_IB_R'-A_RB_I')$, $\Pi_4=-2(\alpha_R\beta_R'+\alpha_I\beta_I')\equiv-2(A_RB_R'+A_IB_I')$ and $\Gamma_i=0$ for $i>k-4$. The matrix A makes it possible to write the analyzed process in the VAR(1) form (see e.g. Lutkepohl2005).\\ From equation ((ref)) we obtain the set of full conditional posterior distributions for model's parameters:

Having the set of full conditional posterior distributions, the pseudo-random sample from the joint posterior distribution may be obtained with the help of the Gibbs sampler, similarly as Koop_al2009 in CI(1,1) case.\\ In the first step the initial values are proposed - $\Sigma^{(0)},\ \nu^{(0)},\ \Gamma^{(0)},\ A_1^{(0)},\ B_1^{(0)},\ A_2^{(0)},\ B_2^{(0)},\ A_{\star}^{(0)},\ B_{\star}^{(0)}$, then the following steps are reiterated:

The square root of the complex Hermitian matrix $(\bar{B}_{\star}^{(s)'}B_{\star}^{(s)})^{\frac{1}{2}}$, may be obtained with the Newton's method proposed by Highman1986.

Point estimation of the cointegration space

Information on the cointegration spaces obtained from the data may be summarized in the point estimate of these spaces and the measure of their posterior distributions' dispersion. Villani2006 proposed to employ the Frobenius (Hilbert-Schmidt) matrix norm to build the loss function needed to point estimation of the real cointegration space. The same approach can be used to estimate complex spaces (see e.g. Srivastava2000).\\ Employing the Frobenius matrix norm $\|A\|_F=(tr(\bar{A}'A))^{\frac{1}{2}}$ to the projection matrices we can built the loss function of the following form:

where $r$ denotes the number of cointegrating vectors.\\ The loss function ((ref)) reaches its minimum in:

where $\nu_i$ ($i=1,2,\dots,r$) is the eigenvector of the matrix $E(\beta\bar{\beta}')$ corresponding to its $i$th largest eigenvalue (see Chikuse2003, Villani2006).\\ The numerical realization of $\hat{\beta}$ is obtained with the use of the pseudo-random sample from the posterior distribution of $\beta$, $\{\beta^{(s)},\ s=1,2,\dots,S\}$, by approximating $E(\beta\bar{\beta}')$ as $\frac{1}{S}\sum_{s=1}^S\beta^{(s)}\bar{\beta}^{(s)'}$.\\ Following Villani2006 we use the projective Frobenius span variation:

where $\lambda_i$ is the $i$th largest eigenvalue of $E(\beta\bar{\beta}')$. The measure $\tau_{sp(\beta)}^2$ reaches its minimum value when the distribution is degenerated, whereas it hits the maximum value for the uniform distribution over the complex Grassmann manifold ($\mathbb{G}_{r,m-r}$).\\

Simulation experiment

As a first illustration of the proposed methods, we perform a small simulation study. We use one of the data generating processes proposed by Cubadda_Omtzigt2005. We simulate 250 data points from:

where $B_1=B_2=\left(

\right)$, $B_{\star}=\left(

\right)+i\left(

\right)$, $A_1=\left(

\right)$, $A_2=\left(

\right)$, $A_{\star}=i\left(

\right)$, $\Gamma=\left(

\right)$ and $\Sigma=\left(

\right)$.\\ Initial values are set to zeros, then the first 50 points are discarded, so we are left with 200 modeled data points.\\ We impose the following priors:

Table (ref) gathered the results of the Bayesian model comparison (see appendix for more information about the employed methods).

There are 3 models with the posterior probability higher than 0.001 and they gathered 0.997 of probability mass. The true model is on the second place with the posterior probability equal 0.383. All the models displayed in the Table (ref) have proper numbers of cointegrating vectors, have no seasonal dummies and differ only in the type of deterministic components.\\ In Table (ref) we present the marginal posterior probabilities of models' features.

The results of models comparison generally correctly point to proper model's characteristics with the exception of the type of deterministic components because the models with the constant restricted to the cointegrating space at zero frequency gathered 0.606 of the posterior probability which is approximately 1.5 more than the whole posterior probability of models without a constant, i.e. the true specification. In the next step of this small simulation experiment we estimate the true model and present the distances between obtained cointegrating spaces and the assumed ones (Table (ref)).

At the first glance one may notice significant differences between point estimate of the cointegration vectors at the annual frequency and the assumed one, but it should be remembered that the data contain information only about the cointegrating spaces, not the vectors, and the distance between the true and estimated space is very low (see the last column of Table (ref)). Generally, the results of the performed simulation experiment proofed that the proposed procedures works well, as the differences between true and estimated spaces at all frequencies are negligible. Moreover, values of the measure $\tau_{sp(\beta)}^2$ are close to zero, so the posterior distributions of the cointegrating spaces are almost degenerated, which can be expected during the analysis of artificial data. Now we proceed to the employment of the proposed model in the real data analysis.

Empirical illustration

In the empirical analysis, we will consider the four-dimensional time series consisted of GDP in constant prices from 2010, consumer price index ($2015=100$), the broad monetary aggregate M3, and the spread between long- and short-term interest rates approximated as the difference between 10-Year Bond Yield and 3-months WIBOR. The quarterly data cover the period 2002Q1 - 2019Q4.\\ A similar model was analyzed by Kotlowski2005.\\ Figure (ref) depicts the analyzed time series. The seasonality may be observed in the paths of all the considered time series, whereas the strongest seasonal pattern is visible in GDP. The trending behavior of the series is also noticeable.

Following Franses1994 and Granger_al1993 we present also the unit transformation of the analyzed time series (Figure (ref)), where the trending behavior and seasonal patterns are more visible. The seasonal variations of GDP differ from the changes observed at the rest of the series.

To fully define the Bayesian seasonally cointegrated VAR model we impose the following priors:

In order to check the nature of the analyzed series we start the analysis by the Bayesian model comparison. The models may differ in the number of cointegrating relations ($r_j\in\{0,1,2,3,4\}$ for $j=1,2,3$) at zero ($j=1$), $\pi$ ($j=2$) and $\frac{\pi}{2}$, $\frac{3\pi}{2}$ ($j=3$) frequencies. We consider models with a linear trend restricted to the cointegration space at zero frequency and with an unrestricted constant ($d=1$), models with an unrestricted constant ($d=2$), specifications with a constant restricted to the cointegration space at zero frequency ($d=3$), and models without constant ($d=4$). The models without ($s=0$) and with ($s=1$) seasonal dummies were examined. Each of the considered specifications has five lags in VAR representation. After excluding non-possible feature combinations and leaving in the set one representation of the observationally equivalent models we are left with 784 pairwise different models. We assume equal prior probability of each specification, i.e. $p(M_{d,s,r_1,r_2,r_3})=\frac{1}{784}\approx 0.0013$. Note that imposing uniform probability on the models' space does not lead to uniform prior distribution for the model features (see the numbers in parentheses in Table (ref)).\\ The models with posterior probability higher than 0.01 are displayed in Table (ref). Almost all of the listed models assume that the analyzed times series may be treated as a realization of the process with 3 or 4 bi-annual relations (note that 4 means stability at $\pi$ frequency). The models ranked at the first and second place assume one long-run relationship and an unrestricted constant. They differ only in the number of cointegrating vectors at the bi-annual frequency (4 or 3). Note that according to the results displayed in Table (ref), models with two relations at 0 frequency are the most probable and gathered 0.409 of the posterior probability, whereas models with one relation 0.363. Models assuming stability at bi-annual frequency together obtained 0.432 of the probability mass, so there is still evidence of cointegration at this frequency as the whole posterior probability of models assuming it is higher and equals 0.562. The posterior probability of the number of cointegrating relations at annual frequency is also diffused, but the whole probability of cointegration equals 0.949, so there is a strong confirmation of the existence of cointegration of this type.

The discussion of nature of the analyzed time series will be complemented by point estimation of the cointegration spaces in the most probable model - $M_{2,0,1,4,3}$.\\

The posterior distribution of the cointegration space at zero frequency is quite diffuse, as the measure $\tau^2_{sp(\beta_1)}$ is roughly in the middle of the interval $[0, 1]$. The visual inspection of deviations from the obtained cointegrating relation (graphed at Figure (ref)) confirms its stationarity.\\ The estimation's results of the cointegration space at the annual frequency are a bit surprising. The measure $\tau^2_{sp(\beta_{\star})}$ is just over 0.9, so the posterior distribution of the cointegration space is almost uniform.

Deviations from the obtained relations seem to be stationary, but by a closer look at the first relation we can notice that its real and imaginary parts are almost the same \footnote{Similar remark applies to the $3^{rd}$ relation}. Moreover, these paths resemble the dynamics of transformations for SPREAD, i.e. $SPREAD_t-SPREAD_{t-2}$ (see Figure (ref)). Such results may indicate for stationarity of SPREAD at the annual frequency. This hypothesis may be checked by the Bayesian model comparison between models with such stationarity restriction imposed and without it, but this is left for further research.

Conclusions

In this paper, the Bayesian seasonally cointegrated vector error correction model for quarterly data was introduced. The empirical usefulness of the discussed methods was illustrated by the analysis of four-dimensional time series. Results of the model comparison indicate that data support the hypothesis of cointegration at zero and annual frequency. There is also evidence of cointegration at bi-annual frequency.\\ This paper focuses only on Bayesian model comparison and point estimation in the most probable model, but as the posterior distribution of the model's specification is diffused it would be useful to take advantage of the Bayesian knowledge pooling in the set of the most probable models. Moreover, as there are papers demonstrating risk in omitting seasonal cointegration (see Introduction), the presented research can also be extended by comparison of forecasts performed in the set of models taking into account relationships at seasonal frequencies and those allowing for only long-run relations. A similar comparison may be performed for structural analyses. Such examination is left for further research.