2.1. A Gaussian linear state-space model

For the purpose of explanation, we start with a Gaussian linear state-space model, which is extended to the heteroskedastic $t$ -QVAR model in this paper. The observable variables are collected into the vector $Y_{t}$ $(\mathcal{N}\times 1)$ and the state variables are collected into the vector $X_{t}$ $(\mathcal{N}\times 1)$ . The first representation of the Gaussian linear state-space model is the following system of equations:

(1) \begin{equation} X_{t}=A X_{t-1}+B \epsilon _{t} \end{equation}

(2) \begin{equation} Y_{t}=C X_{t-1}+D \epsilon _{t} \end{equation}

where $\epsilon _{t}$ $(\mathcal{N}\times 1)$ has an i.i.d. multivariate normal distribution $\epsilon _{t}\sim N(0,\Omega \Omega^{\prime})$ and a diagonal covariance matrix $\Omega \Omega^{\prime}$ . The parameter matrices are $A$ , $B$ , $C$ , and $D$ , which are square matrices with dimensions $(\mathcal{N} \times \mathcal{N})$ . If matrix $D$ is non-singular, then, by using the reduced-form error $v_{t}\equiv D \epsilon _{t}$ , an equivalent representation of the system of Equations (1) and (2) is

(3) \begin{equation} X_{t}=A X_{t-1}+B D^{-1} v_{t} \end{equation}

(4) \begin{equation} Y_{t}=C X_{t-1}+v_{t} \end{equation}

where $v_{t}\sim N(0,\Sigma )\equiv N[0,D\Omega \Omega^{\prime}D^{\prime}]$ is the multivariate i.i.d. reduced-form error term. Moreover, if $C$ is non-singular, then we get the following Gaussian-VARMA(1,1) representation:

(5) \begin{equation} Y_{t}=C A C^{-1} Y_{t-1} + \left(C B D^{-1} -C A C^{-1}\right)v_{t-1}+v_{t} \end{equation}

The reduced-form representation of Equations (3) and (4) is extended to the heteroskedastic $t$ -QVAR model according to the following points. In the heteroskedastic $t$ -QVAR model, (i) $v_{t}$ has a multivariate $t$ -distribution with degrees of freedom $\nu$ , (ii) $v_{t}$ is conditionally heteroskedastic, and (iii) the updating term $v_{t}$ in Equation (3) is replaced by the score function $u_{t}$ , for which we have that $u_{t} \rightarrow _{p} v_{t}$ if $\nu \rightarrow \infty$ . The resulting heteroskedastic $t$ -QVAR model is presented in the following section.

2.2. Heteroskedastic score-driven $t$ -QVAR model

The first representation of the heteroskedastic $t$ -QVAR model is

(6) \begin{equation} X_{t}=A X_{t-1}+B D^{-1} u_{t} \end{equation}

(7) \begin{equation} Y_{t}=C X_{t-1}+D \epsilon _{t} \end{equation}

where $D$ is non-singular and $\epsilon _{t}|\mathcal{F}_{t-1}\sim t(0,\Omega _{t}\Omega^{\prime}_{t},\nu )$ is a multivariate $t$ -distribution with scale matrix $\Omega _{t}\Omega^{\prime}_{t}$ and degrees of freedom $\nu \gt 2$ , where the sigma algebra is $\mathcal{F}_{t-1}=\sigma (Y_{1},\ldots,Y_{t-1},X_{1},\Omega _{1})$ . We use parameter vector $X_{0}$ , which is jointly estimated with the rest of the time-invariant parameters, to initialize the state variables $X_{t}$ . Technical details for the updating term $u_{t}$ (i.e., scaled score function) are presented in Section 2.3. By using the reduced-form error $v_{t}\equiv D \epsilon _{t}$ , the second representation is

(8) \begin{equation} X_{t}=A X_{t-1}+B D^{-1} u_{t} \end{equation}

(9) \begin{equation} Y_{t}=C X_{t-1}+v_{t} \end{equation}

where $v_{t}|\mathcal{F}_{t-1}\sim t(0,\Sigma _{t},\nu )\equiv t(0,D \Omega _{t}\Omega^{\prime}_{t}D^{\prime},\nu )$ . By comparing the models of Equations (8) and (9) and Equations (3) and (4), the mathematical formulas indicate that $v_{t}$ of Equation (3) is replaced by $u_{t}$ , and the error term with normal distribution of Equation (4) is replaced by an error term with $t$ -distribution. If $C$ is non-singular, then we get a nonlinear $t$ -VARMA(1,1) with heteroskedastic errors representation:

(10) \begin{equation} Y_{t}=C A C^{-1} Y_{t-1} + C B D^{-1} u_{t-1} - C A C^{-1} v_{t-1}+v_{t} \end{equation}

2.3. Score functions

The observation-driven updates of filter $X_{t}$ and the scale matrix $\Omega _{t}\Omega^{\prime}_{t}$ of the $t$ -distribution are defined as follows, respectively: First, for the multivariate $t$ -distribution, the log conditional density of $Y_{t}$ is

(11) \begin{align} \ln f(Y_{t}|\mathcal{F}_{t-1},\Theta ) & = \ln \Gamma \left (\frac{\nu +\mathcal{N}}{2}\right ) -\ln \Gamma \left (\frac{\nu }{2}\right ) -\frac{\mathcal{N}}{2}\ln (\pi \nu )\nonumber\\& -\frac {1}{2}\ln |\Sigma _{t}|-\frac {\nu +\mathcal {N}}{2}\ln \left (1+\frac {v^{\prime}_{t}\Sigma _{t}^{-1}v_{t}}{\nu }\right ) \end{align}

where $v_{t}=Y_{t}-C X_{t-1}$ and $\Theta$ is the vector of the time-invariant parameters with dimension $S$ . The derivative of the log conditional density, with respect to the conditional mean $C X_{t-1}$ , is [Harvey (Reference Harvey2013)]:

(12) \begin{equation} \frac{\partial \ln f(Y_{t}|\mathcal{F}_{t-1},\Theta )}{\partial [C X_{t-1}]}= \frac{\nu +\mathcal{N}}{\nu }\Sigma _{t}^{-1} \times \left (1+\frac{v^{\prime}_{t}\Sigma _{t}^{-1}v_{t}}{\nu }\right )^{-1} v_{t}\equiv \frac{\nu +\mathcal{N}}{\nu }\Sigma _{t}^{-1} \times u_{t} \end{equation}

where the last equality defines the scaled score function $u_{t}$ with respect to location, that is,

(13) \begin{equation} u_{t}=\left (1+\frac{v^{\prime}_{t}\Sigma _{t}^{-1}v_{t}}{\nu }\right )^{-1} v_{t} =\left [1+\frac{\epsilon^{\prime}_{t} D^{\prime}\Sigma _{t}^{-1}D\epsilon _{t}}{\nu }\right ]^{-1} D\epsilon _{t} \end{equation}

which, asymptotically and at the true values of parameters $\Theta _{0}$ , is an uncorrelated white noise vector; the conditions of which are presented in the Supplementary Material.

Second, the diagonal of the scale matrix $\Omega _{t}\Omega^{\prime}_{t}$ is driven by score function $e_{t}$ . For the heteroskedastic $t$ -QVAR model, $\Omega _{t}\Omega^{\prime}_{t}$ is formulated as follows:

(14) \begin{equation} \Omega _{t}\Omega^{\prime}_{t}= \left [ \begin{array}{c@{\quad}c@{\quad}c} \exp (2\lambda _{1,t})&0&0\\[4pt] 0&\ddots &0\\[4pt] 0&0&\exp (2\lambda _{\mathcal{N},t}) \end{array}\right ] \end{equation}

where the filters within the diagonal of $\Omega _{t}\Omega^{\prime}_{t}$ are score-driven, and specified as follows:

(15) \begin{equation} \lambda _{i,t}= \omega _{i}+\beta _{i}\lambda _{i,t-1}+\alpha _{i}e_{i,t-1}+\alpha _{i}^{*}\text{sgn}({-}\epsilon _{i,t-1})(e_{i,t-1}+1) \equiv \omega _{i}+\beta _{i}\lambda _{i,t-1}+g(e_{i,t-1}) \end{equation}

for $i=1,\ldots,\mathcal{N}$ , where $\text{sgn}(x)$ is the signum function, and $\alpha _{i}^{*}$ measures asymmetric effects in log-scale. The specification of each element in the diagonal of $\Omega _{t}$ is the Beta- $t$ -EGARCH(1,1) with leverage effects model [Harvey and Chakravarty (Reference Harvey and Chakravarty2008)]. The score function $e_{i,t}$ is defined in what follows.

The marginal distribution of the multivariate $t$ -distribution is the univariate Student’s $t$ -distribution, where the degrees of freedom parameters coincide. The partial derivatives of the log marginal densities of the univariate Student’s $t$ -distributions $\ln f_{i}(Y_{i,t}|\mathcal{F}_{t-1},\Theta )$ , with respect to $\lambda _{i,t}$ , for $i=1,\ldots,\mathcal{N}$ , are the following score functions with respect to log-scale [Harvey and Chakravarty (Reference Harvey and Chakravarty2008)]:

(16) \begin{equation} e_{i,t}= \frac{\partial \ln f_{i}(Y_{i,t}|\mathcal{F}_{t-1},\Theta )}{\partial \lambda _{i,t}} =\frac{(\nu +1)v_{i,t}^{2}}{\nu \exp (2\lambda _{i,t})+v_{i,t}^{2}}-1 \end{equation}

for $i=1,\ldots,\mathcal{N}$ . These score functions update $\lambda _{i,t}$ within the scale matrix $\Omega _{t}\Omega^{\prime}_{t}$ . The updating terms $e_{i,t}$ and $g(e_{i,t})$ for $i=1,\ldots,\mathcal{N}$ in Equation (15), asymptotically and at the true values of parameters $\Theta _{0}$ , are uncorrelated white noise; the conditions of which are presented in the Supplementary Material.

If $\nu \rightarrow \infty$ , then $e_{i,t}$ of Equation (16) converges in probability to

(17) \begin{equation} e_{i,t}= \frac{\partial \ln f_{i}(Y_{i,t}|\mathcal{F}_{t-1},\Theta )}{\partial \lambda _{i,t}} =\frac{v_{i,t}^{2}}{\exp (2\lambda _{i,t})}-1 \end{equation}

for $i=1,\ldots,\mathcal{N}$ , which are quadratic transformations of $v_{i,t}$ as in the GARCH model [Bollerslev (Reference Bollerslev1986)].

If $\nu \rightarrow \infty$ , then we also have that $u_{t} \rightarrow _{p} v_{t}$ , $\epsilon _{t} \sim t(0,\Omega _{t}\Omega^{\prime}_{t},\nu ) \rightarrow _{d} N(0,\Omega _{t}\Omega^{\prime}_{t})$ , and $v_{t}\sim t(0,\Sigma _{t},\nu ) \rightarrow _{d} N(0,\Sigma _{t})$ . Moreover, if we also assume homoskedastic innovations, that is $\Sigma _{t}=\Sigma$ , then we obtain a classical Gaussian linear state-space model (Section 2.1).

2.4. Impulse responses

For the reduced-form error term $\text{Var}(v_{t}|\mathcal{F}_{t-1})=\Sigma _{t} \times \nu/(\nu -2)$ , which is factorized as

(18) \begin{equation} \text{Var}(v_{t}|\mathcal{F}_{t-1})= \Sigma _{t} \times \frac{\nu }{\nu -2}= \left (\frac{\nu }{\nu -2}\right )^{1/2} \times D\Omega _{t}\Omega^{\prime}_{t}D^{\prime} \times \left (\frac{\nu }{\nu -2}\right )^{1/2} \end{equation}

For the IRF analysis, the multivariate i.i.d. structural-form error term $\tilde{\epsilon }_{t}$ is introduced as

(19) \begin{equation} \tilde{\epsilon }_{t}=\left (\frac{\nu }{\nu -2}\right )^{1/2}D^{-1}\Omega _{t}^{-1}v_{t} =\left (\frac{\nu }{\nu -2}\right )^{1/2}D^{-1}\Omega _{t}^{-1}D\epsilon _{t} \end{equation}

where $\tilde{\epsilon }_{t} \sim t[0,I_{\mathcal{N}}\times (\nu -2)/\nu,\nu ]$ is an i.i.d. multivariate $t$ -distribution with zero mean and an identity covariance matrix. From Equations (12) to (19), the $\tilde{\epsilon }_{t}$ error term-based representation of $u_{t}$ is

(20) \begin{equation} u_{t}=[(\nu -2)\nu ]^{1/2} D\Omega _{t} \times \frac{\tilde{\epsilon }_{t}}{\nu -2+\tilde{\epsilon}^{\prime}_{t}\tilde{\epsilon }_{t}} \end{equation}

The nonlinear MA( $\infty$ ) representation of $Y_{t}$ is

(21) \begin{equation} Y_{t}= \left \{\sum _{j=1}^{\infty }C A^{j-1} B [(\nu -2)\nu ]^{1/2} \Omega _{t-j}\frac{\tilde{\epsilon }_{t-j}}{\nu -2+\tilde{\epsilon }_{t-j}^{\prime}\tilde{\epsilon }_{t-j}}\right \} +\left (\frac{\nu }{\nu -2}\right )^{1/2} D \Omega _{t}\tilde{\epsilon }_{t} \end{equation}

The IRFs for the heteroskedastic $t$ -QVAR model are

(22) \begin{equation} \frac{\partial Y_{t+j}}{\partial \tilde{\epsilon }_{t}}=\left (\frac{\nu }{\nu -2}\right )^{1/2} D\Omega _{t} \quad \text{for} \quad j=0 \end{equation}

(23) \begin{equation} \frac{\partial Y_{t+j}}{\partial \tilde{\epsilon }_{t}}=C A^{j-1} B [(\nu -2)\nu ]^{1/2}\Omega _{t} \tilde{D}_{t} \quad \text{for} \quad j=1,\ldots,\infty \end{equation}

where $\tilde{D}_{t}$ is

(24) \begin{align} \tilde{D}_{t} & = \frac{\partial \frac{\tilde{\epsilon }_{t}}{\nu -2+\tilde{\epsilon}^{\prime}_{t}\tilde{\epsilon }_{t}}}{\partial \tilde{\epsilon }_{t}}= \left [ \begin{array}{c@{\quad}c@{\quad}c} d_{11,t}&\cdots &d_{1\mathcal{N},t}\\ \cdots &\cdots &\cdots \\ d_{\mathcal{N}1,t}&\cdots &d_{\mathcal{N}\mathcal{N},t}\\ \end{array} \right ]=\nonumber\\[7pt] & =\left [ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \frac {\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t}-2\tilde {\epsilon }_{1t}^{2}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}}& \frac {-2\tilde {\epsilon }_{1t}\tilde {\epsilon }_{2t}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}}&\cdots & \frac {-2\tilde {\epsilon }_{1t}\tilde {\epsilon }_{\mathcal {N}t}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}}\\[5pt] \frac {-2\tilde {\epsilon }_{2t}\tilde {\epsilon }_{1t}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}}& \frac {\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t}-2\tilde {\epsilon }_{2t}^{2}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}}&\cdots &\cdots \\ \cdots &\cdots &\cdots &\cdots \\ \frac {-2\tilde {\epsilon }_{\mathcal {N}t}\tilde {\epsilon }_{1t}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}} &\cdots &\cdots &\frac {\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t}-2\tilde {\epsilon }_{\mathcal {N}t}^{2}}{(\nu -2+\tilde {\epsilon}^{\prime}_{t}\tilde {\epsilon }_{t})^{2}} \end {array} \right ] \end{align}

Motivated by the work of Herwartz and Lütkepohl (Reference Herwartz and Lütkepohl2000), we replace the matrices $\Omega _{t}$ and $\Omega _{t} \tilde{D}_{t}$ by $E(\Omega _{t})$ and $E(\Omega _{t} \tilde{D}_{t})$ in Equations (22) and (23), respectively, which are estimated by using sample averages of the full sample period. In the literature, there are alternative approaches for the estimation of nonlinear IRFs [Lütkepohl (Reference Lütkepohl2005)]. Future empirical applications could consider other nonlinear IRF estimation approaches, which may be more appropriate than our approach for specific policy analyses.

2.5. Statistical inference

We use the maximum likelihood (ML) method [Creal et al. (Reference Creal, Koopman and Lucas2013), Harvey (Reference Harvey2013), Blasques et al. (Reference Blasques, Gorgi, Koopman and Wintenberger2018, Reference Blasques, van Brummelen, Koopman and Lucas2022)]:

(25) \begin{equation} \hat{\Theta }=\arg \max _{\Theta } \text{LL}(Y_{1},\ldots,Y_{T},\Theta )= \arg \max _{\Theta } \sum _{t=1}^{T} \ln f(Y_{t}|\mathcal{F}_{t-1},\Theta ) \end{equation}

where the log conditional density $\ln f(Y_{t}|\mathcal{F}_{t-1},\Theta )$ is Equation (11). The LL function is maximized numerically, by using alternative start values of parameters. The standard errors are estimated using the inverse information matrix: $\{(1/T)\sum _{t=1}^{T}[G_{t}(\hat{\Theta })^{\prime}G_{t}(\hat{\Theta })]\}^{-1}$ , where $G_{t}(\hat{\Theta })$ is the gradient.

With respect to the sufficient conditions of the asymptotic properties of the ML estimator for score-driven location plus score-driven scale models, such as the heteroskedastic $t$ -QVAR model, we refer to the work of Blazsek et al. (Reference Blazsek, Escribano and Licht2022), which is directly related to the present work from the perspective of statistical inference. The heteroskedastic $t$ -QVAR model is an extension of the bivariate $t$ -QVAR model of Blazsek et al. (Reference Blazsek, Escribano and Licht2022), because in the present paper we have multivariate location plus multivariate scale filters, while in the work of those authors the model includes univariate location plus univariate scale filters. Although most of the ML assumptions coincide for the two models, due to the different levels of technical complexity at some points of model formulation, we report the asymptotic theory for the heteroskedastic $t$ -QVAR model in the Supplementary Material.

In the remainder of this section, we present some of those conditions for covariance stationarity and invertibility. For the covariance stationarity of $X_{t}$ , asymptotically and at the true values of parameters $\Theta _{0}$ , it is required that the maximum modulus of the eigenvalues of $A$ , that is $\text{Stat}_{\mu }$ , is less than one. For the invertibility of $X_{t}$ , that is $X_{t}$ converges almost surely (a.s.) to a unique strictly stationary vector sequence for all $\Theta \in \tilde{\Theta }$ where $\tilde{\Theta }$ is the parameter set, one of the conditions of invertibility of $X_{t}$ is

(26) \begin{equation} \text{Inv}_{\mu }= \inf _{n \geq 1}\left \{n^{-1} E\left (\ln \left |\left |\prod _{t=1}^{n} A+B D^{-1} \frac{\partial u_{t-1}}{\partial (X_{t-1})^{\prime}} \right |\right |_{2}\right )\right \}\lt 0 \end{equation}

where the matrix norm $||\mathcal{W}||_{2}$ is the spectral norm. We present empirical estimates of $\text{Stat}_{\mu }$ and $\text{Inv}_{\mu }$ in the application.

Covariance stationarity of $\lambda _{i,t}$ , for $i=1,\ldots,\mathcal{N}$ , asymptotically and at the true values of parameters $\Theta _{0}$ , requires $|\beta _{i}|\lt 1$ for $i=1,\ldots,\mathcal{N}$ , respectively, which are denoted by $\text{Stat}_{\lambda,i}$ for $i=1,\ldots,\mathcal{N}$ . Hence, $E(\lambda _{i,t})=\omega _{i}/(1-\beta _{i})$ , for $i=1,\ldots,\mathcal{N}$ , which is applied to initialize $\lambda _{i,t}$ for $i=1,\ldots,\mathcal{N}$ , respectively. For the invertibility of $\lambda _{i,t}$ for $i=1,\ldots,\mathcal{N}$ , that is $\lambda _{i,t}$ for $i=1,\ldots,\mathcal{N}$ converge a.s. to a unique strictly stationary sequences for all $\Theta \in \tilde{\Theta }$ , one of the conditions of the invertibility of $\lambda _{i,t}$ is

(27) \begin{equation} \text{Inv}_{\lambda,i} =\inf _{n \geq 1}\left \{n^{-1}E\left [\ln \left |\prod _{t=1}^{n} \beta _{i}+[\alpha _{i}+\alpha _{i}^{*}\text{sgn}({-}v_{i,t})] \frac{\partial e_{i,t}}{\partial \lambda _{i,t}}\right |\right ]\right \} \end{equation}

for $i=1,\ldots,\mathcal{N}$ . For the invertibility of filters $\lambda _{i,t}$ , it is required that $\text{Inv}_{\lambda,i}\lt 0$ for $i=1,\ldots,\mathcal{N}$ . We present the estimates of $\text{Stat}_{\lambda,i}$ and $\text{Inv}_{\lambda,i}$ , for $i=1,\ldots,\mathcal{N}$ , in the empirical application.

3.2. Heteroskedastic $t$ -ABCD representation

We present the heteroskedastic $t$ -ABCD representation which is adapted to the DSGE model of An and Schorfheide (Reference An and Schorfheide2007). We specify the ABCD representation according to Morris (Reference Morris2014), as follows:

(34) \begin{equation} \underbrace{ \left [ \begin{array}{c} z_{t}\\ g_{t}\\ r_{t} \end{array} \right ]}_{X_{t}} = \underbrace{ \left [ \begin{array}{c@{\quad}c@{\quad}c} \rho _{z}&0&0\\ 0&\rho _{g}&0\\ c_{r,z}&0&c_{r,r} \end{array} \right ]}_{A} \underbrace{ \left [ \begin{array}{c} z_{t-1}\\ g_{t-1}\\ r_{t-1} \end{array} \right ]}_{X_{t-1}} + \underbrace{ \left [ \begin{array}{c@{\quad}c@{\quad}c} 1&0&0\\ 0&1&0\\ c_{r,z}/\rho _{z}&0&c_{r,r}/\rho _{r} \end{array} \right ]}_{B} \underbrace{ \left [ \begin{array}{c@{\quad}c@{\quad}c} c_{r,z}/\rho _{z}&0&c_{r,r}/\rho _{r}\\ c_{y,z}/\rho _{z}&1&c_{y,r}/\rho _{r}\\ c_{\pi,z}/\rho _{z}&0&c_{\pi,r}/\rho _{r} \end{array} \right ]^{-1}}_{D^{-1}} \underbrace{ \left [ \begin{array}{c} u_{z,t}\\ u_{g,t}\\ u_{r,t} \end{array} \right ]}_{u_{t}} \end{equation}

(35) \begin{equation} \underbrace{ \left [ \begin{array}{c} r_{t}\\ y_{t}\\ \pi _{t} \end{array} \right ]}_{Y_{t}} = \underbrace{ \left [ \begin{array}{c@{\quad}c@{\quad}c} c_{r,z}&0&c_{r,r}\\ c_{y,z}&\rho _{g}&c_{y,r}\\ c_{\pi,z}&0&c_{\pi,r} \end{array} \right ]}_{C} \underbrace{ \left [ \begin{array}{c} z_{t-1}\\ g_{t-1}\\ r_{t-1} \end{array} \right ]}_{X_{t-1}} + \underbrace{ \left [ \begin{array}{c} v_{z,t}\\ v_{g,t}\\ v_{r,t} \end{array} \right ]}_{v_{t}} \qquad\qquad\qquad\qquad\qquad\quad\qquad\end{equation}

These equations correspond to Equations (8) and (9) of the heteroskedastic $t$ -QVAR model. The parameters of the ABCD-matrices are identified, and the IRFs can be estimated according to the nonlinear MA( $\infty$ ) representation of the heteroskedastic $t$ -QVAR model.

For the heteroskedastic $t$ -ABCD representation, the following parameters are estimated:

(36) \begin{align} \Theta & = (\rho _{z},\rho _{g},\rho _{r},c_{r,z},c_{r,r},c_{y,z},c_{y,r},c_{\pi,z},c_{\pi,r},X_{0,z},X_{0,g},X_{0,r}, \omega _{z},\omega _{g},\omega _{r},\beta _{z},\beta _{g},\beta _{r},\nonumber\\& \qquad \alpha _{z},\alpha _{g},\alpha _{r},\alpha _{z}^{*}, \alpha _{g}^{*},\alpha _{r}^{*},\nu )^{\prime} \end{align}

The heteroskedastic $t$ -ABCD representation is a generalization of the homoskedastic Gaussian-ABCD representation of Fernández-Villaverde et al. (Reference Fernández-Villaverde, Rubio-Ramírez, Sargent and Watson2007) and the heteroskedastic score-driven Gaussian-ABCD representation of Angelini and Gorgi (Reference Angelini and Gorgi2018). By assuming that $\Omega _{t}=\Omega$ and $\nu \rightarrow \infty$ , we get the homoskedastic Gaussian-ABCD representation [Fernández-Villaverde et al. (Reference Fernández-Villaverde, Rubio-Ramírez, Sargent and Watson2007)]. Moreover, by assuming that $\Omega _{t}=\Omega$ , we get the new homoskedastic $t$ -ABCD representation. In addition, we note that we provide the following contributions to the work of Angelini and Gorgi (Reference Angelini and Gorgi2018): (i) Angelini and Gorgi (Reference Angelini and Gorgi2018) assume that only the scale filter is score-driven. We use score-driven transition equation and scale filters. (ii) Angelini and Gorgi (Reference Angelini and Gorgi2018) use a score-driven multivariate Gaussian distribution, which we extend to the score-driven multivariate $t$ -distribution. (iii) We extend the score-driven scale filter of Angelini and Gorgi (Reference Angelini and Gorgi2018) by adding asymmetric effects to the updates.

3.3. Data

All quarterly time series data of this paper are from Federal Reserve Economic Data (FRED). We use the following variables for the period of 1954 Q3 to 2022 Q1: (a) not seasonally adjusted Effective Federal Funds Rate (in % points), (b) seasonally adjusted US GDP level $\tilde{y}_{t}$ , and (c) seasonally adjusted US CPI (consumer price index) for all urban consumers. We use the maximum period for which data are available from FRED for these variables. The sample period includes the recent and on-going crisis periods of the coronavirus disease of 2019 (COVID-19) pandemic and the Russian invasion of Ukraine.

The US inflation rate series is computed as quarterly log percentage change (in % points) of the US CPI. The seasonally adjusted US GDP level is detrended by using the following model:

(37) \begin{equation} \ln \tilde{y}_{t}=\ln A_{t}+\ln y_{t}^{\dagger } \end{equation}

(38) \begin{equation} \ln A_{t}=\gamma +\ln A_{t-1}+ z_{t} \end{equation}

where log aggregate productivity $\ln A_{t}$ is the unobservable stochastic trend for US GDP level [Giacomini (Reference Giacomini2013)], and the model is estimated as a Gaussian linear state-space model by using the Kalman filter.

In this paper, the steady state of each observable dependent variable is estimated by using the sample average [Morris (Reference Morris2014)]. There are some negative observations for the US inflation time series. Hence, we do not use the logarithmic percentage difference between the value of the variables and their steady states [Giacomini (Reference Giacomini2013)], but the sample average is subtracted from each variable [Morris (Reference Morris2014)]. Hence, the sample averages are subtracted from the US inflation rate and the Effective Federal Funds Rate in the computations of $\pi _{t}$ and $r_{t}$ , respectively (both are measured in % points). For percentage deviation of US GDP from its trend, we use $y_{t}\equiv 100 \times \ln y_{t}^{\dagger }$ (hence $y_{t}$ is in % points), and $E(y_{t})=0$ due to Equations (37) and (38).

Table 1. Descriptive statistics

Notes: Augmented Dickey–Fuller generalized least squares (ADF–GLS); autoregressive conditional heteroskedasticity (ARCH). For the ADF–GLS test with constant, the modified Bayesian information criterion (BIC) is used for lag-order selection. For the ARCH test, lag-order $4$ is used. $\text{H}_{0}$ for the Shapiro–Wilk test [Shapiro and Wilk (Reference Shapiro and Wilk1965)] is the normal distribution (rejected for all variables at all levels of significance). $\text{H}_{0}$ for the ADF–GFS test is integration of order one (not able to reject for any variable at the 10% level of significance). $\text{H}_{0}$ for the ARCH test is no ARCH effect is present (rejected for all variables at all levels of significance).

Descriptive statistics are presented in Table 1. The evolution of the observable dependent variables $r_{t}$ , $y_{t}$ , and $\pi _{t}$ , for the period of 1954 Q3 to 2022 Q1, is presented in Figure 1.

3.4. Parameter estimates and model diagnostics

In the empirical application, we estimate the following models: (i) a homoskedastic Gaussian-QVAR model which is called the homoskedastic Gaussian-ABCD representation; (ii) a heteroskedastic Gaussian-QVAR model which is called the heteroskedastic score-driven Gaussian-ABCD representation; (iii) a homoskedastic $t$ -QVAR model which is called the homoskedastic score-driven $t$ -ABCD representation; and (iv) a heteroskedastic $t$ -QVAR model which is called the heteroskedastic score-driven $t$ -ABCD representation. As aforementioned, models (i), (ii), and (iii) are special cases of model (iv). In Table 2, the ML parameter estimates, likelihood-based model performance metrics, and stationarity and invertibility statistics are presented for all ABCD representations.

Figure 1. Observable dependent variables for the period of 1954 Q3 to 2022 Q1.

Table 2. Parameter estimates and model diagnostics

Notes: AIC, Akaike information criterion; BIC, Bayesian information criterion; HQC, Hannan–Quinn criterion; LL, log-likelihood. Standard errors are reported in parentheses. The bold likelihood-based model performance metrics indicate superior statistical performance. The heteroskedastic score-driven Gaussian-ABCD representation is the DSGE model of Angelini and Gorgi (Reference Angelini and Gorgi2018). The Gaussian-ABCD representations have VAR(1) representations, hence $\text{Inv}_{\mu }=0$ . $^{***}$ , $^{**}$ , $^{*}$ , and $^{+}$ indicate parameter significance at the 1%, 5%, 10%, and 15% levels, respectively.

First, almost all parameters of the ABCD-matrices for all representations are significantly different from zero in Table 2. By comparing the ML results for alternative ABCD representations, the estimate of the dynamic parameter $\rho _{r}$ is lower (though statistically significant) for the heteroskedastic score-driven Gaussian-ABCD and heteroskedastic score-driven $t$ -ABCD representations ( $\hat{\rho }_{r}=0.0077^{***}$ and $\hat{\rho }_{r}=0.0373^{***}$ , respectively), than for the homoskedastic Gaussian-ABCD and score-driven homoskedastic $t$ -ABCD representations ( $\hat{\rho }_{r}=0.7388^{***}$ and $\hat{\rho }_{r}=0.3384^{***}$ , respectively). As these results may suggest that the dynamic effects for $r_{t}$ are different for the homoskedastic and heteroskedastic ABCD representations, we return to this point in Section 3.5 where the IRF estimates are presented. Moreover, the parameter estimates for the score-driven heteroskedastic $t$ -ABCD representation show that scale dynamics, that is $\beta _{z}$ , $\beta _{g}$ , $\beta _{r}$ and $\alpha _{z}$ , $\alpha _{g}$ , $\alpha _{r}$ , $\alpha _{z}^{*}$ , $\alpha _{g}^{*}$ , $\alpha _{r}^{*}$ , are significantly different from zero. In addition, the degrees of freedom parameter $\nu$ estimates for both $t$ -ABCD representations justify the use of the Student’s $t$ -distribution, as their estimates are significantly lower than 30.

Second, the statistical performances of the ABCD representations are compared by using the LL, Akaike information criterion, Bayesian information criterion, and Hannan–Quinn criterion metrics in Table 2. We find that the performances of the Gaussian-ABCD representations are improved by using the score-driven $t$ -ABCD representations. We also find that the best statistical performance is provided by the heteroskedastic score-driven $t$ -ABCD representation.

Third, the empirical estimates of the covariance stationarity and Lyapunov exponent statistics are also presented in Table 2. For the Gaussian-ABCD representations, covariance stationarity is supported by the $\text{Stat}_{\mu }$ statistic. For the $t$ -ABCD representations, covariance stationarity and invertibility are supported by the $\text{Stat}_{\mu }$ and $\text{Inv}_{\mu }$ statistics, respectively. For both heteroskedastic ABCD representations, covariance stationarity and invertibility are supported by the $\text{Stat}_{\lambda,z}$ , $\text{Stat}_{\lambda,g}$ , $\text{Stat}_{\lambda,r}$ , and $\text{Inv}_{\lambda,z}$ , $\text{Inv}_{\lambda,g}$ , $\text{Inv}_{\lambda,r}$ , respectively. We note that both Gaussian-ABCD representations of the DSGE model of An and Schorfheide (Reference An and Schorfheide2007) have VAR(1) forms [see Morris (Reference Morris2014)], hence $\text{Inv}_{\mu }=0$ for those models.

Fourth, in Figure 2, as a graphical illustration of heteroskedasticity for the DSGE model, the evolution of $\epsilon _{i,t} \pm \sigma _{i,t}$ (i.e., error $\pm$ one standard deviation) is presented for the period of 1954 Q3 to 2022 Q1, by using the estimates for the best-performing heteroskedastic score-driven $t$ -ABCD representation. The figure shows that the standard deviations of the error terms are time-varying for the DSGE model.

3.5. Impulse response analysis

IRFs are identified using sign restrictions, based on 10,000 Monte Carlo simulations of the matrix $\Omega$ , in accordance with the procedure of Rubio-Ramírez et al. (Reference Rubio-Ramírez, Waggoner and Zha2010). First, the ML estimates of $\Omega$ are used. Second, an $\mathcal{N} \times \mathcal{N}$ matrix $\tilde{K}$ of i.i.d. $N(0,1)$ numbers is simulated. Third, the QR decomposition of $\tilde{K}$ is performed, and the resulting matrices are denoted $\tilde{Q}$ and $\tilde{R}$ . Fourth, we define $\tilde{\Omega }\equiv \Omega \times \tilde{Q}^{\prime}$ for each simulation. The parameter matrices $\Omega$ and $\Omega _{t}$ are replaced by $\tilde{\Omega }$ in the IRFs.

For each simulation of $\tilde{\Omega }$ , sign restrictions are used in accordance with Table 3, for which we refer to the works of Leeper et al. (Reference Leeper, Walker and Yang2009) and Cogan et al. (Reference Cogan, Cwick, Taylor and Wieland2010). For the simulations of $\tilde{\Omega }$ that satisfy the sign restrictions of Table 3, we report the 5%, 50%, and 95% percentiles of the IRFs up to 20 leads in Figures 3 and 4 for the classical homoskedastic Gaussian-ABCD and the new heteroskedastic score-driven $t$ -ABCD representations, respectively.

Figure 2. Evolution of $\epsilon _{t}$ and $\epsilon _{t} \pm \sigma _{t}$ for score-driven heteroskedastic $t$ -ABCD (1954 Q3 to 2022 Q1).

Notes: $\sigma _{i,t}=\exp (\lambda _{i,t})[\nu/(\nu -2)]^{1/2}$ for $i=z,g,r$ .

Figure 3. IRFs for the homoskedastic Gaussian-ABCD representation (5%, 50%, and 95% percentiles).

Figure 4. IRFs for the heteroskedastic score-driven $t$ -ABCD representation (5%, 50%, and 95% percentiles).

Figure 5. Comparison of IRFs estimates for the periods of 1954 Q3 to 2019 Q4 and 1954 Q3 to 2022 Q1 (5%, 50%, and 95% percentiles).

We highlight the following IRF results: First, the relatively low estimate of $\rho _{r}$ for the heteroskedastic ABCD representations, compared to the homoskedastic ABCD representations (Section 3.4), might be interpreted as a lower effect of an interest rate shock on subsequent interest rates, but, as is shown in Panel (c) of Figures 3 and 4, that is not the case. The persistence of interest rate shocks is similar for all ABCD representations [Panel (c) of Figures 3 and 4].

Second, a difference between the IRFs of the two ABCD representations, which may have policy implications, is in relation to the effects of an aggregate productivity shock on the US inflation rate $\tilde{\epsilon }_{z,t}\rightarrow \pi _{t+j}$ . For the homoskedastic Gaussian-ABCD representation, $\tilde{\epsilon }_{z,t}\rightarrow \pi _{t+j}$ is negative, and its contemporaneous value is $-0.41[-0.11,-0.62]$ [Panel (g) of Figure 3]. For the heteroskedastic $t$ -ABCD representation, $\tilde{\epsilon }_{z,t}\rightarrow \pi _{t+j}$ is also negative, but its effect is lower in absolute value, that is, its contemporaneous value is $-0.10[-0.02,-0.14]$ [Panel (g) of Figure 4].

Third, another difference between the IRFs of different ABCD representations is in relation to the effects of government spending shocks on the US inflation rate $\tilde{\epsilon }_{g,t}\rightarrow \pi _{t+j}$ [Panel (h) of Figures 3 and 4]. For the homoskedastic Gaussian-ABCD representation, $\tilde{\epsilon }_{g,t}\rightarrow \pi _{t+j}$ is positive contemporaneously, and its value is $0.28[0.02,0.61]$ [Panel (h) of Figure 3]. For the heteroskedastic $t$ -ABCD representation, $\tilde{\epsilon }_{z,t}\rightarrow \pi _{t+j}$ is also positive, but its effect is lower in absolute value, that is, its contemporaneous value is $0.07[0.01,0.14]$ [Panel (h) of Figure 4].

Fourth, another difference between the IRFs of different ABCD representations is in relation to monetary policy shocks, that is, effects of interest rate shocks on the US inflation rate $\tilde{\epsilon }_{r,t}\rightarrow \pi _{t+j}$ [Panel (i) of Figures 3 and 4]. For the homoskedastic Gaussian-ABCD representation, $\tilde{\epsilon }_{r,t}\rightarrow \pi _{t+j}$ is negative, and its contemporaneous value is $-0.25[-0.02,-0.61]$ [Panel (i) of Figure 3]. For the heteroskedastic $t$ -ABCD representation, $\tilde{\epsilon }_{r,t}\rightarrow \pi _{t+j}$ is also negative, but its effect is lower in absolute value, that is, its contemporaneous value is $-0.06[-0.01,-0.13]$ [Panel (i) of Figure 4].

Fifth, as measured by the IRF confidence intervals, the estimation precision of the IRF for US GDP growth [Panels (d)−(f) of Figures 3 and 4] and US inflation rate [Panels (g)−(i) of Figures 3 and 4] is higher for the heteroskedastic score-driven $t$ -ABCD representation than for the homoskedastic Gaussian-ABCD representation. This indicates that the score-driven updates of the dependent variable may improve the statistical inferences of the impulse responses.

Table 3. Sign restrictions on impact responses

Finally, in Figure 5, we compare some of the IRF estimates for the homoskedastic Gaussian-ABCD representation and the heteroskedastic score-driven $t$ -ABCD representation, by using data for the periods of 1954 Q3 to 2019 Q4 and 1954 Q3 to 2022 Q1 (i.e., data without and with the COVID-19 pandemic, respectively). The IRF estimates show that the heteroskedastic score-driven $t$ -ABCD representation is robust to the COVID-19 pandemic, while for the homoskedastic Gaussian-ABCD representation some of the IRF estimates are significantly different for the two samples. These results indicate the robustness of the score-driven $t$ -ABCD representation to extreme observations.