Identification and Estimation of SVARMA models with Independent and Non-Gaussian Inputs

Identification and Estimation of SVARMA models with Independent and Non-Gaussian Inputs

\thispagestyle{empty}

Proposed Running Head

Non-Gaussian SVARMA Identification

Affiliation

and

E-mail

[email removed]

\thispagestyle{empty}

Abstract

This paper analyzes identifiability properties of structural vector autoregressive moving average (SVARMA) models driven by independent and non-Gaussian shocks. It is well known, that SVARMA models driven by Gaussian errors are not identified without imposing further identifying restrictions on the parameters. Even in reduced form and assuming stability and invertibility, vector autoregressive moving average models are in general not identified without requiring certain parameter matrices to be non-singular. Independence and non-Gaussianity of the shocks is used to show that they are identified up to permutations and scalings. In this way, typically imposed identifying restrictions are made testable. Furthermore, we introduce a maximum-likelihood estimator of the non-Gaussian SVARMA model which is consistent and asymptotically normally distributed.

Keywords: Structural vector autoregressive moving-average models, non-Gaussianity, Identifiability

JEL classification: C32, C51, E52

Introduction

\setcounter{page}{1}

Recently, LMS_svarIdent16 and GourierouxZakoianRenne17 have shown that structural vector autoregressive (SVAR) models driven by independent non-Gaussian components are identified up to scaling and permutations which makes the typically imposed identifying restrictions testable. If the error terms driving the economy are Gaussian or (cross-sectionally) merely uncorrelated (as opposed to independent), one has to resort to identifying restrictions in order to conclude on the fundamental shocks driving the economy. From analysis in terms of second moments, the true shocks can be identified only up to multiplication with orthogonal matrices (which all lead to the same second moments of the observed process). Non-Gaussianity combined with cross-sectional independence, however, allows to identify the shocks up to permutations and scalings. In particular, infinitely many linear combinations of shocks generating the same second moments are reduced to a finite set of linear combinations generating the same distributional outcome. It is thus possible to employ a data-driven approach instead of a story-telling approach. Most importantly, the identifying (story-imposed) restrictions are made testable when using the (data-driven) non-Gaussian SVARMA approach.

Structural econometric analysis is usually conducted with SVAR models. The situation for structural VARMA models driven by independent non-Gaussian shocks is more complicated because one has to take additional identifiability restrictions on the parameter space into account. In this paper, spectral factorization techniques are employed to generalize the SVAR results by LMS_svarIdent16 to the SVARMA case. While the literature on SVAR models is abundant, see KilianLut17 and references therein, the contributions regarding SVARMA models are easier to keep track of, see, e.g., BoubacarFrancq11 and GourierouxMR_svarma19. In structural econometric analysis, the impulse response function (IRF) and variance decompositions are the primary objects of interest luet05,KilianLut17. Especially in macroeconometrics, where data is sometimes available only at quarterly instances, it is of paramount importance to use a parsimoniously parameterized models (like e.g. SVARMA models) for which the IRF and other can be obtained straight-forwardly. It is widely known that SVARMA models are superior to SVAR models in this respect, see, e.g., HannanDeistler12. Moreover, the articles Poskitt16, PoskittYao17, RaghavanAthSilvapulle16, AthVahid08, and AthVahid_JTSA_08 provide ample evidence and make a strong point for using VARMA models instead of VAR models for econometric analysis.

In a recent contribution, GourierouxMR_svarma19 consider the dynamic identification problem in SVARMA models. While their focus is a general treatment of whether it is possible to identify the root location of determinantal roots of the associated MA polynomial matrix in the structural VARMA case, we focus here on the precise derivation of the properties of the maximum likelihood (ML) estimator of the fundamental representation, including the first and second partial derivatives with respect to all system and noise parameters.

One (perceived) disadvantage of VARMA models is increased complexity of the estimation procedure compared to VAR models. Two rebuttals are in order. First, there are many sophisticated (e.g. non-linear threshold) VAR models whose estimation is arguably more involved than the one of VARMA models. Second, there are many stable and openly available software implementations which should put the complexities of estimation of VAR and VARMA models on the same level. Examples for implementations in the R software environment Rcore are Scherrer_rldm, Tsay13,Tsay_R_MTS and Gilbert_R_dse, see also ScherrerDeistler2019_handbook for a comparison and further comments on these packages, and in MATLAB Gomez_matlab_15,Gomez16. The estimation procedure described in this article is implemented in R and can be installed with the command devtools::install_github(“bfunovits/svarma_id”) in the R console.

The rest of the paper is structured as follows. In section 2, the SVARMA model is introduced. In section 3, the identification result is stated and proved. In section 4, the maximum likelihood (ML) estimator is derived and shown to be consistent and asymptotically normal. In section 5, we illustrate the method. Proofs and technical details are available in the Online Appendix.

We use $z$ as a complex variable as well as the backward shift operator on a stochastic process, i.e. $z\left(y_{t}\right)_{t\in\mathbb{Z}}=\left(y_{t-1}\right)_{t\in\mathbb{Z}}$ and define $i=\sqrt{-1}$. The transpose of an $\left(m\times n\right)$-dimensional matrix $A$ is denoted by $A'$. The column-wise vectorization of $A\in\mathbb{R}^{m\times n}$ is denoted by $vec\left(A\right)\in\mathbb{R}^{mn\times1}$ and for a square matrix $B\in\mathbb{R}^{n\times n}$ we denote with $vecd{^\circ}\left(B\right)\in\mathbb{R}^{n(n-1)}$ the vectorization where the diagonal elements of $B$ are left out. The $n$-dimensional identity matrix is denoted by $I_{n}$, an $n$-dimensional diagonal matrix with diagonal elements $\left(a_{1},\ldots,a_{n}\right)$ is denoted by $\text{diag}\left(a_{1},\ldots,a_{n}\right)$, and the inequality $">0"$ means positive definiteness in the context of matrices. The column vector $\iota_{i}$ has a one at positions $i$ and zeros everywhere else. The expectation of a random variable with respect to a given probability space is denoted by $\mathbb{E}\left(\cdot\right)$. Convergence in probability and in distribution are denoted by $\xrightarrow{p}$ and $\xrightarrow{d}$, respectively. Partial derivatives $\left.\frac{\partial f(x)}{\partial x}\right|_{x=x_{0}}$ of a real-valued function $f(x)$ evaluated at a point $x_{0}\in\mathbb{R}^{k}$ are denoted by $f_{x}\left(x_{0}\right)$ and considered columns.

Model

We start from an $n$-dimensional VARMA system

The shocks $\left(\varepsilon_{t}\right)_{t\in\mathbb{Z}}$ driving the system are identically and independently distributed (i.i.d.) in cross-section and time, have zero mean, and diagonal covariance matrix $\Sigma^{2}$ with positive diagonal elements $\sigma_{i}^{2}$, whose positive square root is in turn denoted by $\sigma_{i}$ . To simplify presentation, we also introduce the column vector $\sigma=\left(\sigma_{1},\ldots,\sigma_{n}\right)'$ and $\Sigma=\text{diag}\left(\sigma_{1},\ldots,\sigma_{n}\right)$, as well as $x_{t-1}'=\left(y_{t-1}',\ldots,y_{t-p}'\right)$ and $s_{t-1}'=\left(\varepsilon_{t-1}'B',\ldots,\varepsilon_{t-q}'B'\right)$ such that equation (ref) can be written as \[ y_{t}=\left(a_{1},\ldots,a_{p}\right)x_{t-1}+\left(b_{1},\ldots,b_{q}\right)s_{t-1}+B\varepsilon_{t}. \] We assume that the stability condition

and the strict invertibility condition

hold, and that $B$ is invertible and has ones on its diagonal. Furthermore, we assume that the polynomial matrices $a(z)$ and $b(z)$ are left-coprime\footnote{Two matrix polynomials are called left-coprime if $\left(a(z),b(z)\right)$ is of full row rank for all $z\in\mathbb{C}$. For equivalent definitions see HannanDeistler12 Lemma 2.2.1 on page 40.} and that $\left(a_{p},b_{q}\right)$ is of full rank\footnote{The stability, invertibility, coprimeness, and full-rank assumptions on the parameters in $a(z)$ and $b(z)$ can be relaxed. Imposing them, allows us to focus on the essential part of this contribution: To reduce the class of observational equivalence in terms of second moments from the orthogonal matrices to permutation matrices in the context of SVARMA models.}. Note that this full rank assumption is over-identifying in the sense that some rational transfer function cannot be parameterized by any VARMA(p,q) system which satisfies this assumption, see Hannan71 or HannanDeistler12, Chapter 2.7 on page 77.

The stationary solution $\left(y_{t}\right)_{t\in\mathbb{Z}}$ of the system (ref) is called an ARMA process.

We follow Rothenberg71 to define identifiability of parametric models. The external characteristic of the stationary solution $\left(y_{t}\right)_{t\in\mathbb{Z}}$ of (ref) is the probability distribution function (or a subset of corresponding moments). A particular system (ref) is described by the parameters of (ref) which satisfy assumptions (ref) and (ref) as well as the coprimeness assumption, the full rank assumption and the assumptions on $B$ and $D^{2}$. The model is then characterized by the set of all a priori possible systems which we will call internal characteristics. Two systems of the form (ref) are called observationally equivalent if they imply the same external characteristics of $\left(y_{t}\right)_{t\in\mathbb{Z}}$. A system is identifiable if there is no other observationally equivalent system. The identifiability problem is concerned with the existence of an injective function from the internal characteristics to the external characteristics\footnote{The inverse of this function, i.e. from the external to the internal characteristics, is called the identifying function.}, see DeistlerSeifert78 for a more detailed discussion.

The classical (non-)identifiability issues where the external characteristics are described by the second moments of $\left(y_{t}\right)_{t\in\mathbb{Z}}$ are best understood in terms of the spectral density of the stationary solution of (ref). The spectral density, i.e. the Fourier transform of the autocovariance function $\gamma(s)=\mathbb{E}\left(y_{t}y_{t-s}'\right),\ s\in\mathbb{Z},$ of $\left(y_{t}\right)_{t\in\mathbb{Z}}$ , is \[ f(z)=a(z)^{-1}b(z)B\Sigma^{2}B'b'\left(\frac{1}{z}\right)a'\left(\frac{1}{z}\right)^{-1}=k(z)\left(B\Sigma^{2}B'\right)k'\left(\frac{1}{z}\right), \] evaluated at $z=e^{-i\lambda}$, $\lambda\in\left[-\pi,\pi\right]$, where $k(z)=a(z)^{-1}b(z)=\sum_{j=0}^{\infty}k_{j}z^{j},\ k(0)=I_{n}$, and $k(z)B$ corresponds to the transfer function relating the output $y_{t}$ to the $\left(\varepsilon_{t}\right)_{t\in\mathbb{Z}}$ .

On the one hand, transforming the pair $\left(B,\Sigma\right)$ with an orthogonal matrix\footnote{A square matrix is orthogonal if $QQ'=Q'Q=I_{n}$.} $Q$ to $\left(B\Sigma Q\Sigma_{1}^{-1},\Sigma_{1}\right)$, where $\Sigma_{1}$ is a diagonal matrix such that the diagonal elements of $B$ are equal to one, generates the same spectral density because $B_{1}\Sigma_{1}^{2}B_{1}'=B\Sigma^{2}B'$ where $B_{1}=B\Sigma Q\Sigma_{1}^{-1}$. Hence, the class of observational equivalence is at least $\frac{n(n-1)}{2}$-dimensional. On the other hand, it is easy to see Hannan70 that two spectral factors\footnote{A spectral factor $l(z)$ is a rational matrix function for which $l(z)l'\left(\frac{1}{z}\right)$, evaluated at the unit circle, is equal to the spectral density.} of the form $k(z)B\Sigma=a(z)^{-1}b(z)B\Sigma$, where $a(z)$ and $b(z)$ satisfy (ref) and (ref) as well as the coprimeness assumption, the full rank assumption and where $B$ and $\Sigma$ satisfy the assumptions outlined above, obtained from the spectral density corresponding to the stationary solution of (ref) are related through orthogonal matrices. This means that any other spectral factor is of the form $a(z)^{-1}b(z)B\Sigma Q$ where $Q$ is an orthogonal matrix. By normalizing the diagonal elements of $B\Sigma Q$, we obtain a new pair $\left(B_{1},\Sigma_{1}\right)$ of the required form. Hence, the class of observational equivalence is $\frac{n(n-1)}{2}$-dimensional. This result, however, only uses second moment information and not the full distribution of the stochastic process $\left(\varepsilon_{t}\right)_{t\in\mathbb{Z}}$.

We will show in the next section that if the inputs $\left(\varepsilon_{t}\right)$ to (ref) are non-Gaussian and independent, the spectral factors are related by permutation matrices (modulo sign). Thus, we reduce the class of observational equivalence from the group of orthogonal matrices to the group of (signed) permutations.

Identification of the Instantaneous Shock Transmission

In this section, we first use the cross-sectional independence and non-Gaussianity of the components of the shocks $\varepsilon_{t}$ for identifying the matrix $B$ up to permutation and scaling of its columns.Finally, we discuss advantages and disadvantages of various rules for choosing a particular permutation and scaling.

The assumptions on the error term $\varepsilon_{t}=\left(\varepsilon_{1,t},\ldots,\varepsilon_{n,t}\right)$ are the same as in LMS_svarIdent16, the essential one being that the components (at one point in time) are mutually independent and that at most one of them has a Gaussian marginal distribution.

In order to strengthen intuition as to how non-Gaussianity and independence help reducing the size of the class of observational equivalence, consider the following example featuring two identically and independently uniformly distributed random variables. Rotating these two variables 45 degrees (with rotation matrix $\frac{1}{\sqrt{2}}\left(

\right)$) leads to marginal distributions which are “more Gaussian” (e.g. measured by the absolute value of the excess kurtosis) than the original variables. This suggests that searching for linear combinations that lead to “maximally non-Gaussian” variables might pin down a rotation. In the following, we present a formal approach.

Fixing a Rotation

The theoretical background for reducing the class of observational equivalence from orthogonal matrices to (signed) permutations is provided by the following lemma. It allows to conclude from the independence of the sums of independent variables on the distribution of the underlying summands. In particular, it is useful to conclude on the coefficients pertaining to the summands if one makes additional assumptions on the distribution of the summands.

We use

In the following, Lemma (ref) is used to conclude on the columns of $M$ in $\varepsilon_{t}=M\varepsilon_{t}^{*}$, where $M=B^{-1}B^{*}$, where both $\varepsilon_{t}$ and $\varepsilon_{t}^{*}$ are assumed to be (cross-sectionally) independent and non-Gaussian. The components of $\varepsilon_{t}$ correspond to $Y_{1},\ Y_{2}$, the components of $\varepsilon_{t}^{*}$ correspond to $X_{1},\ldots,X_{n}$. E.g., for component 1 and 2 of $\varepsilon_{t}$ we have $\varepsilon_{1,t}=\left(m_{11},\ldots,m_{1n}\right)\varepsilon_{t}^{*}$ and $\varepsilon_{2,t}=\left(m_{21},\ldots,m_{2n}\right)\varepsilon_{t}^{*}$. If any pair of coefficients $\left(m_{1k},m_{2k}\right)$ satisfies $m_{1k}m_{2k}\neq0$, then the corresponding component $\varepsilon_{k,t}^{*}$ is Gaussian according to the Lemma. By Assumption 1, at most one component of $\varepsilon_{t}^{*}$ is allowed to have a Gaussian marginal distribution. It follows that there cannot be another pair $\left(m_{1l},m_{2l}\right),\ l\neq k,$ that satisfies $m_{1l}m_{2l}\neq0$. In particular, there is (at most) one non-zero coefficient in the scalar product $\left\langle m_{1,\bullet},m_{2,\bullet}\right\rangle =m_{1k}m_{2k}\neq0$, where $m_{i,\bullet}$ denotes the $i$-th row of $M$. If $\left\langle m_{1,\bullet},m_{2,\bullet}\right\rangle =m_{1k}m_{2k}\neq0$, we obtain a contradiction to the assumption that $\mathbb{E}\left(\varepsilon_{1,t}\varepsilon_{2,t}\right)=0$ because from the fact that one (exactly one) component $\varepsilon_{k,t}^{*}$ is Gaussian and $\varepsilon_{i,t}=m_{i,\bullet}

^{'}$ we obtain that $\mathbb{E}\left(\varepsilon_{1,t}\varepsilon_{2,t}\right)=m_{1,\bullet}D^{*}m_{2,\bullet}'=d_{k}^{*}m_{1k}m_{2k}\neq0$. It thus follows that all pairs $\left(m_{1k},m_{2k}\right)$ satisfy $m_{1k}m_{2k}=0$. Since this argument holds for all pairs in $\varepsilon_{1,t},\ldots,\varepsilon_{n,t}$, it follows that every column contains at most one non-zero element. Finally, non-singularity implies that every column contains exactly one non-zero element.

Now we are ready to prove

While the proof above is easily understandable for readers who know the paper LMS_svarIdent16, the following proof uses less matrix algebra.

The difference in these two proofs is as follows. In the first proof, we use model (ref) together with its assumptions earlier, i.e. we write down the ARMA equation, take expectations, and obtain that any two independent error terms satisfy $\varepsilon_{t}^{*}=\left(B^{*}\right)^{-1}B\varepsilon_{t}$. In the second proof, we focus firstly on the linear innovations $v_{t}$ and only use the fact that the error terms are independent when it comes to parameterizing the covariance matrix of the innovations. Note that the second proof suggests that as soon as one can identify the true inputs (irrespective of the model), one may use cross-sectional independence and non-Gaussianity to reduce the equivalence class of orthogonal matrices to the one of (signed) permutation matrices.

Identification Scheme: Choosing a Unique Permutation and Scaling

In this section, we describe how to pick one particular permutation and scaling from the class of observational equivalence described in the previous section. In order to do this, we describe different identification schemes, i.e. rules for choosing a particular permutation and scaling of the matrix $B$.

We start by repeating two identification schemes presented in LMS_svarIdent16 (which are in turn based on IlmonenPaindaveine11 and HallinMehta15). The first identification scheme, which is convenient for deriving asymptotic properties and which we refer to as identification scheme A, consists in firstly scaling all columns of $B$ such that their norm is equal to one, secondly, permutating the columns such that the absolute value of each diagonal element is larger than the absolute value of all elements in the same row with a higher column index, and finally scaling all columns of $B$ such that the diagonal elements are equal to one\footnote{Note that in the derivation of the ML estimator, we impose only that the diagonal elements of $B$ be equal to one. Thus, the restrictions, in general, do not suffice to pin down the particular permutation and scaling for $B$. However, the fact that the observationally equivalent points in the parameter space are discrete ensures the existence of a consistent root, i.e. the solution of the first order conditions obtained from taking derivatives of the standardized log-likelihood function. Should the gradient descent algorithm return a $B$ matrix which does not satisfy the identification scheme, it can be easily transformed such that the identification scheme is satisfied. The companion R-package to this article transforms the $B$ matrix such that all restrictions described here are satisfied.}. The second identification scheme consists of the same first two steps but instead of scaling the columns in the last step such that their diagonal elements are equal to one, it is required that the diagonal elements are positive. Sometimes, the second identification scheme turns out to be more flexible, for example when testing hypotheses involving diagonal elements. Regarding the derivation of asymptotic properties, however, one would need to maximize the constrained (log-) likelihood function where the restrictions that the columns of $B$ have length one are taken into account.

It is important to realize that the transformations used in the identification schemes described above, exist not on the whole parameter space but only on a topologically large set in the parameter set. For details, see Proposition 2 in LMS_svarIdent16 including an example of a matrix or which the above identification schemes are not defined. The third identification scheme, similar to the one in ChenBickel05 on page 3626, does not exclude any non-singular matrix $B$ and is defined by the following transformations. Firstly, the columns of $B$ are scaled to have norm equal to one. Secondly, in each column, the element with largest absolute value is made positive. Finally, the columns are ordered according to $\prec$ such that $c\prec d$ for two columns $c,d$ of $B$ if and only if there exists a $k\in\left\{ 1,\ldots,n\right\} $ such that $c_{k}<d_{k}$ and $c_{j}=d_{j}$ for all $j\in\left\{ 1,\ldots,k-1\right\} $.

Now that we have firstly obtained a discrete set of observationally equivalent SVARMA systems and secondly provided different rules to select a unique representative, we may proceed to local ML estimation of the true underlying parameter.

Parameter Estimation

In this section, we treat local ML estimation of (ref). In particular, we prove local consistency and asymptotic normality of the ML estimator (MLE).

In order to separate the essential ideas from technicalities, we start by stating a theorem for local asymptotic normality of the MLE in terms of (easily understandable and intuitive) high-level assumptions on the densities of i.i.d. shocks. Next, we discuss (component-wise) the densities of the i.i.d. shocks $\left(\varepsilon_{t}\right)$, the admissible parameter space, and the (standardized) log-likelihood function of our problem at hand. Last, we state a theorem for local asymptotic normality of the MLE in terms of low-level integrability and differentiability assumptions on the densities and verify the high-level assumptions. The proofs and many technicalities (e.g. partial derivatives of the likelihood function) which are similar to the ones in LMS_svarIdent16 are deferred to the Online Appendix.

Local Asymptotic Normality in terms of High-Level Assumptions

For the sake of clarity, and in order to understand where the low-level assumptions on the densities that we will introduce in Assumption (ref)below come into play, we state a theorem proving local asymptotic normality in terms of high-level assumptions. In the Online Appendix, we show how the low-level assumptions imply the high-level assumptions. The (standardized) log-likelihood function to be maximized is

where $l_{t}\left(\varepsilon_{t}(\theta).\theta\right)=\log\left(f\left(\varepsilon_{t}(\theta).\theta\right)\right)$ are the individual contributions to the log-likelihood function and $f(\cdot)$ is the (joint) density of a residuals $\varepsilon_{t}\left(\theta\right)$ which are obtained from a parametric model with parameter $\theta\in\Theta\subseteq\mathbb{R}^{k}$.

The following discussion builds on poetpruch97. Firstly, the existence of a sequence of solutions $\left(\hat{\theta}_{T}\right)$ of the first order condition $L_{\theta,T}\left(\theta\right)=0$ of the standardized log-likelihood function which converges almost surely towards $\theta_{0}$ is required. This is essentially guaranteed by the identification result in the previous section (and some technical conditions), showing that the observationally equivalent points in the parameter space are discrete (in the sense that there exist disjoint open sets around each point of this kind). Furthermore, the score of the individual contributions $l_{t}\left(\theta\right)$ to $L_{T}\left(\theta\right)$ has to satisfy a Central Limit Theorem (CLT) for martingale difference sequences (MDS) and the Hessian of the individual contributions has to satisfy a Uniform Law of Large Numbers (ULLN). If these conditions are satisfied, the sequence $\sqrt{T}\left(\hat{\theta}_{T}-\theta_{0}\right)$ is asymptotically normal. To make this discussion more precise, we state

The basic idea consists in applying (component-wise) the mean value theorem to $\left(L_{\theta,T}\left(\hat{\theta}_{T}\right),L_{\theta,T}\left(\theta_{0}\right)\right)$ such that one obtains asymptotically $\sqrt{T}\left(L_{\theta,T}\left(\hat{\theta}_{T}\right)-L_{\theta,T}\left(\theta_{0}\right)\right)=\bar{A}_{T}\sqrt{T}\left(\hat{\theta}_{T}-\theta_{0}\right)$ where the matrix $\bar{A}_{T}$ corresponds to the Hessian whose rows are evaluated at the respective mean values. Point (ref)) is necessary for the existence of a consistent sequence $\left(\hat{\theta}_{T}\right)$ and is, together with point (ref)), (ref)), (ref)), and (ref)), required for the CLT for MDS. It follows that $\sqrt{T}L_{\theta,T}\left(\theta_{0}\right)$ and $-\bar{A}_{T}\sqrt{T}\left(\hat{\theta}_{T}-\theta_{0}\right)$ are asymptotically normal with the same asymptotic distribution, i.e. $\mathcal{N}\left(0,\mathbb{E}\left(l_{\theta,t}\left(\varepsilon_{t}(\theta_{0}),\theta_{0}\right)l_{\theta,t}\left(\varepsilon_{t}(\theta_{0}),\theta_{0}\right)'\right)\right)$. Moreover, one needs to ensure that the Hessian satisfies a ULLN\footnote{This means that $\sup_{\theta\in\Theta_{0}}\left\Vert \frac{1}{T}\sum_{t=1}^{T}l_{\theta\theta,t}(\theta)-\mathbb{E}\left(l_{\theta\theta,t}(\theta)\right)\right\Vert \rightarrow0,\ a.s.,$ and as a byproduct $\mathbb{E}\left(l_{\theta\theta,t}\left(\theta\right)\right)$ is continuous at $\theta_{0}$.} such that $\bar{A}_{T}$ converges towards the non-singular expectation of the Hessian evaluated at the true parameter value which is moreover equal to the negative of the expectation of the outer product of the score. This is ensured by points (ref)), (ref)) and (ref))\footnote{Point (ref)) is, e.g., implied by requiring that $\int\sup_{\theta\in\Theta_{0}}\left\Vert l_{\theta,t}\left(\varepsilon_{t}(\theta),\theta\right)\right\Vert dx<\infty$ and $\int\sup_{\theta\in\Theta_{0}}\left\Vert l_{\theta\theta,t}\left(\varepsilon_{t}(\theta),\theta\right)\right\Vert dx<\infty$ but can also be obtained by less stringent assumptions.}, respectively.

Note that the covariance matrix can be consistently estimated by $-A_{T}^{-1}$, where $A_{T}=\frac{1}{T}\sum_{t=1}^{T}\left(l_{\theta\theta,t}\left(\varepsilon_{t}\left(\hat{\theta}_{T}\right),\hat{\theta}_{T}\right)\right)$. Under an additional condition, the outer product of the score can be used as well:

Parameter Space, Log-Likelihood Function, and Low-Level Assumptions

In this section, we specialize the generic theorem stated in the previous section for the problem at hand. First, we describe the parameter space on which we optimize the log-likelihood function. Second, we make assumptions on the densities of the components of $\varepsilon_{t}$. This allows us to provide explicit expressions for the individual contributions to the standardized log-likelihood function and its first partial derivatives\footnote{The expression for the second partial derivatives as well as the tedious but straightforward derivations are deferred to the Online Appendix.}. Third, we state integrability and dominance conditions on the first and second partial derivatives of the densities of the components of $\varepsilon_{t}$. Last, we verify that Assumptions (ref) and (ref) below (together with Assumptions (ref), (ref), and (ref)) imply the ones in the generic Theorem (ref) and are thus sufficient for consistency and local asymptotic normality of the MLE.

In order to introduce the parameter space for the SVARMA parameters, we define $\pi=\left(\pi_{2},\pi_{3}\right)$ where $\pi_{2}=vec\left(a_{1},\ldots,a_{p}\right)$, and $\pi_{3}=vec\left(b_{1},\ldots,b_{q}\right)$. Compared with LMS_svarIdent16 there is an additional sub-vector $\pi_{3}$ for the MA parameters and we abstract in our model from the mean by setting it equal to zero.

We also introduce the non-singleton compact and convex subset $\Theta_{0}=\Theta_{0,\pi}\times\Theta_{0,\beta}\times\Theta_{0,\sigma}\times\Theta_{0,\lambda}$ of the interior of $\Theta$ which contains the true parameter value $\theta_{0}$.

Regarding the component densities of the i.i.d. shock process $\left(\varepsilon_{t}\right)$, we have

Thus, the individual contributions in the (standardized) log-likelihood function (ref) are

where $u_{t}\left(\theta\right)=y_{t}-a_{1}y_{t-1}-\cdots-a_{p}y_{t-p}-b_{1}B\left(\beta\right)\varepsilon_{t-1}\left(\theta\right)-\cdots-b_{q}B\left(\beta\right)\varepsilon_{t-q}\left(\theta\right)$.

The expressions for the partial derivatives of the individual contributions to the standardized log-likelihood function are given as

where $x_{b,t-1}\left(\theta\right)=\left(x_{t-1}\otimes b'(z)^{-1}\right)$, $w_{b,t-1}\left(\theta\right)=\left(w_{t-1}\left(\theta\right)\otimes b'(z)^{-1}\right)$, $w_{t-1}\left(\theta\right)=\left(u'_{t-1}\left(\theta\right),\ldots,u'_{t-q}\left(\theta\right)\right)^{'}$, the matrix $H\in\mathbb{R}^{n^{2}\times n(n-1)}$ consisting of zeros and ones is implicitly defined by $vec\left(B(\beta)\right)=H\beta+vec\left(I_{n}\right)$ for $B$ in $\mathcal{B}$, and \[ e_{i,x,t}(\theta)=\frac{\partial}{\partial x}\log\left[f_{i}\left(\sigma_{i}^{-1}\iota_{i}^{'}B\left(\beta\right)^{-1}u_{t}\left(\theta\right);\lambda_{i}\right)\right]=\frac{f_{i,x}\left(\sigma_{i}^{-1}\varepsilon_{i,t}\left(\theta\right);\lambda_{i}\right)}{f_{i}\left(\sigma_{i}^{-1}\varepsilon_{i,t}\left(\theta\right);\lambda_{i}\right)} \] and \[ e_{i,\lambda_{i},t}(\theta)=\frac{\partial}{\partial\lambda_{i}}\log\left[f_{i}\left(\sigma_{i}^{-1}\iota_{i}^{'}B\left(\beta\right)^{-1}u_{t}\left(\theta\right);\lambda_{i}\right)\right]=\frac{f_{i,\lambda}\left(\sigma_{i}^{-1}\varepsilon_{i,t}\left(\theta\right);\lambda_{i}\right)}{f_{i}\left(\sigma_{i}^{-1}\varepsilon_{i,t}\left(\theta\right);\lambda_{i}\right)}, \] with $f_{i,x}\left(x;\lambda_{i}\right)=\frac{\partial}{\partial x}f_{i}\left(x;\lambda_{i}\right)$ and $f_{i,\lambda_{i}}\left(x;\lambda_{i}\right)=\frac{\partial}{\partial\lambda_{i}}f_{i}\left(x;\lambda_{i}\right)$. Evaluated at the truth, i.e. $\theta=\theta_{0}$, we have that $\varepsilon_{i,t}\left(\theta_{0}\right)=\varepsilon_{i,t}$ and \[ e_{i,x,t}=e_{i,x,t}(\theta_{0})=\left.\frac{\partial}{\partial x}\log\left[f_{i}\left(\sigma_{i}^{-1}\iota_{i}^{'}B\left(\beta\right)^{-1}u_{t}\left(\pi\right);\lambda_{i}\right)\right]\right|_{\theta=\theta_{0}}=\frac{f_{i,x}\left(\sigma_{i}^{-1}\varepsilon_{i,t};\lambda_{i,0}\right)}{f_{i}\left(\sigma_{i,0}^{-1}\varepsilon_{i,t};\lambda_{i,0}\right)}. \]

In order to show that the scores are MDS, that the resulting covariance matrix is finite at the true parameter point, that the expectation of the supremum on $\Theta_{0}$ of the Hessian is finite, and that the expectation of the outer product of the score is equal to the negative expectation of the Hessian of the individual contribution to the likelihood, we need the following two assumptions on the first and second derivatives of the component densities.

and

In combination, these assumptions allow to prove (in the Online Appendix)

Empirical Application

Impulse Response Functions

Often, the goal of macroeconometric analyses is gaining an understanding of the impact of structural economic shocks on the observable variables. This is usually done through analysis of the impulse response function (IRF) or the analysis of variance decompositions (see e.g. luet05 and KilianLut17).

We comment on the differences between obtaining them from SVARMA or SVAR representations. First, note that calculation of the IRF is as straightforward as in the SVAR case after one has obtained the estimates of the structural parameters. One possibility is representing the system in state space form HannanDeistler12. Then, the impulse responses and variance decompositions are obtained in the same way as in the SVAR case, see e.g. KilianLut17.

If, furthermore, the object of interest is the impulse response function it is hard to come up with reasons favoring SVAR models over SVARMA models. While in theory one may approximate SVARMA models (or even “infinite VAR models”) by SVAR models, it is well known that the approximation is bad in many practically relevant cases. This was emphasized in a macroeconometric context by Ravenna07 and PoskittYao17. Ravenna07 decomposes the error when SVARMA models are approximated by SVAR models into a truncation error and an identification error, pertaining to the parameters describing the economic shocks. PoskittYao17 decompose the truncation error introduced in Ravenna07 further into an estimation and approximation error and argue that both are large for commonly used lag lengths and sample sizes. They conclude that “using VAR($n$) may not be justified unless $n$ and {[}the sample size{]} $T$ are enormous”. Obviously, these errors carry over to the IRF which is a non-linear transformation of the structural parameters.

Empirical Application

To illustrate the developed methods, we estimate a three equation macroeconomic model and analyze its impulse response function. A more detailed analysis can be found in the vignette of the R-package associated with this article.

Data

We use the FRED database of the Federal Reserve Bank of St. Louis and retrieve series for the unemployment gap $n_{t}$, i.e. we subtract the unemployment rate (UNRATE) from the natural unemployment rate (NROU), inflation $\pi_{t}$ (lagged differences of GDPDEF), and the effective federal funds rate $R_{t}$ (FEDFUNDS). The observation period starts with Q3 1954 and ends with Q1 2019, thus there are 259 observations.

Estimation Procedure

In order to select appropriate integer-valued parameters $p$ and $q$, we estimate a number of VARMA models with the MTS-package and select the one with the smallest AIC value. It is worth noting that neither the dse-package nor the MTS-package enforces the stability condition (ref) or the invertibility condition (ref). In case of unstable or non-invertible determinantal roots of the $a(z)$ or $b(z)$ matrix polynomials, we mirror these roots outside the unit circle and adjust the error covariance accordingly. Based on the AIC value, the Ljung-Box test, and the McLeod-Li test MahdiMcLeod_R_portes, we choose $p=2$ and $q=2$. Moreover, the distribution of the shocks of the initial model in Figure (ref) and Figure (ref) suggest, and the Jarque-Bera test indicates that the individual series are not normally distributed.

\textcolor{red}

We proceed thus under the assumption that the errors of the VARMA(2,2) model are independent and not normally distributed. The individual series seem leptokurtic and we assume that they follow a Laplace distribution. The full conditional likelihood is subsequently estimated using the same identification scheme as in LMS_svarIdent16 for fixing a particular permutation and scaling. The estimates for $B$ and $\sigma$ are \[ \hat{B}=

\] and $\hat{\sigma}=\left(0.0685,\,0.0315,\,0.14\right)$, the respective (bootstrapped) standard deviations are \[ \hat{\mathbb{V}}_{bs}\left(\hat{B}\right)=

\] and $\hat{\mathbb{V}}_{bs}\left(\hat{\sigma}\right)=\left(0.00427,\,0.00199,\,0.01419\right)$.

Impulse Response Function

In order to interpret the results, we calculate the impulse response function together with bootstrapped confidence intervals (1000 bootstrap replications).\textcolor{red}

We interpret the third shock as a (in the longer term) contractionary monetary policy shock since it is the only one to which the interest rate reacts significantly (judging from the confidence bands). In the medium term, the unemployment gap widens (the economy shrinks) and initially there is a positive response of inflation to Shock 3. The first shock is interpreted as negative demand shock because the economy shrinks (widening of the unemployment gap) and prices fall slightly. We interpret the remaining Shock 2 as negative supply shock since the economy shrinks (widening unemployment gap) with increasing prices.

Acknowledgements

Financial support by the Research Funds of the University of Helsinki as well as by funds of the Oesterreichische Nationalbank (Austrian Central Bank, Anniversary Fund, project number: 17646) is gratefully acknowledged. For bootstrapping, the Finnish Grid and Cloud Infrastructure with persistent identifier urn:nbn:fi:research-infras-2016072533 was used. Juho Koistinen, Mika Meitz, Markku Lanne, and Wolfgang Scherrer provided helpful comments on various versions of this article.

Conclusion

In this article, we showed that stable and invertible SVARMA models (ref) driven by independent and non-Gaussian shocks are identifiable up to permutation and scaling. This result extends the identifiability results regarding structural VAR models in LMS_svarIdent16. SVARMA models capture (macroeconomic) dynamics more parsimoniously than SVAR models and are therefore advantageous in situations with relatively small sample sizes.