The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
66,987 characters
Identification and Estimation of SVARMA models with Independent and Non-Gaussian Inputs
\title{Identification and Estimation of SVARMA models with Independent and
Non-Gaussian Inputs}
\author{Bernd Funovits}
\maketitle
\thispagestyle{empty}
\section*{Proposed Running Head}
Non-Gaussian SVARMA Identification
\section*{Affiliation}
\begin{singlespace}
\textbf{University of Helsinki}
Faculty of Social Sciences
Discipline of Economics
P. O. Box 17 (Arkadiankatu7)
FIN-00014 University of Helsinki
\end{singlespace}
and
\begin{singlespace}
\textbf{TU Wien}
Institute of Statistics and Mathematical Methods in Economics
Econometrics and System Theory
Wiedner Hauptstr. 8
A-1040 Vienna
\end{singlespace}
\section*{E-mail}
[email removed]
\pagebreak{}
\thispagestyle{empty}
\section*{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
\pagebreak{}
\section{Introduction}
\setcounter{page}{1}
Recently, \citet{LMS_svarIdent16} and \citet{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 \citet{LMS_svarIdent16} to the SVARMA case. While the literature
on SVAR models is abundant, see \citet{KilianLut17} and references
therein, the contributions regarding SVARMA models are easier to keep
track of, see, e.g., \citet{BoubacarFrancq11} and \citet{GourierouxMR_svarma19}.
In structural econometric analysis, the impulse response function
(IRF) and variance decompositions are the primary objects of interest
\citep{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., \citet{HannanDeistler12}.
Moreover, the articles \citet{Poskitt16}, \citet{PoskittYao17},
\citet{RaghavanAthSilvapulle16}, \citet{AthVahid08}, and \citet{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, \citet{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 \citet{Rcore} are
\citet{Scherrer_rldm}, \citet{Tsay13,Tsay_R_MTS} and \citet{Gilbert_R_dse},
see also \citet{ScherrerDeistler2019_handbook} for a comparison and
further comments on these packages, and in MATLAB \citet{Gomez_matlab_15,Gomez16}.
The estimation procedure described in this article is implemented
in R and can be installed with the command \texttt{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.
\section{\label{sec:Model}Model}
We start from an $n$-dimensional VARMA system
\begin{equation}
\underbrace{\left(I_{n}-a_{1}z-\cdots a_{p}z^{p}\right)}_{=a(z)}y_{t}=\underbrace{\left(I_{n}+b_{1}z+\cdots+b_{q}z^{q}\right)}_{=b(z)}B\varepsilon_{t},\quad a_{i},b_{i}\in\mathbb{R}^{n\times n}.\label{eq:system}
\end{equation}
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 \eqref{sec:Model} 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
\begin{equation}
\det\left(a(z)\right)\neq0,\ \left|z\right|\leq1,\label{eq:stability}
\end{equation}
and the strict invertibility condition
\begin{equation}
\det\left(b(z)\right)\neq0,\ \left|z\right|\leq1\label{eq:invertibility}
\end{equation}
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 \citet{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 \citet{Hannan71}
or \citet{HannanDeistler12}, Chapter 2.7 on page 77.
The stationary solution $\left(y_{t}\right)_{t\in\mathbb{Z}}$ of
the system \eqref{eq:system} is called an ARMA process.
We follow \citet{Rothenberg71} to define identifiability of parametric
models. The external characteristic of the stationary solution $\left(y_{t}\right)_{t\in\mathbb{Z}}$
of \eqref{eq:system} is the probability distribution function (or
a subset of corresponding moments). A particular system \eqref{eq:system}
is described by the parameters of \eqref{eq:system} which satisfy
assumptions \eqref{eq:stability} and \eqref{eq:invertibility} 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 \eqref{eq:system} 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 \citet{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 \eqref{eq:system}. 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 \citep[page 66]{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 \eqref{eq:stability} and \eqref{eq:invertibility}
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 \eqref{eq:system} 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 \eqref{eq:system} 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.
\section{\label{sec:identification_scheme}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 \citet{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.
\begin{assumption}
\label{assu:non_gaussianIID}We assume the following.
\begin{enumerate}
\item The error process $\varepsilon_{t}=\left(\varepsilon_{1,t},\ldots,\varepsilon_{n,t}\right)$
is a sequence of i.i.d. random vectors. Each component $\varepsilon_{i,t},\ i\in\left\{ 1,\ldots,n\right\} $
has zero mean and positive variance.
\item For any (fixed) point in time, the components of $\varepsilon_{t}$
are mutually independent and at most one of the components has a Gaussian
marginal distribution.
\end{enumerate}
\end{assumption}
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(\begin{smallmatrix}1 & 1\\
1 & -1
\end{smallmatrix}\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.
\subsection{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
\begin{lem}[\citet{Kagan73}, Theorem 3.1.1]
\label{lem:Kagan} Let $X_{1},\ldots X_{n}$ be independent (not
necessarily identically distributed) random variables, and define
$Y_{1}=\sum_{i=1}^{n}a_{i}X_{i}$ and $Y_{2}=\sum_{i=1}^{n}b_{i}X_{i}$
where $a_{i}$ and $b_{i}$ are constants. If $Y_{1}$ and $Y_{2}$
are independent, then the random variables $X_{j}$ for which $a_{j}b_{j}\neq0$
are all normally distributed.
\end{lem}
In the following, Lemma \ref{lem:Kagan} 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}\begin{pmatrix}\varepsilon_{1,t}^{*} & \cdots & \varepsilon_{n,t}^{*}\end{pmatrix}^{'}$
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
\begin{thm}
The set of observationally equivalent ARMA systems of the form in
section \ref{sec:Model} is described by the set of matrices $PD$
where $P$ is a permutation matrix and $D$ a diagonal matrix with
non-zero diagonal entries.
\end{thm}
\begin{proof}
Consider two systems \eqref{eq:system}, say $\left(a(z),b(z);B,\Sigma\right)$
and $\left(a^{*}(z),b^{*}(z);B^{*},\Sigma^{*}\right)$ whose stationary
solutions have the same spectral density (or equivalently the same
second moments), in particular $B\Sigma^{2}B'=B^{*}\Sigma^{*2}B^{*'}$.
Written differently, we consider
\[
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}
\]
and
\[
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}^{*},
\]
post-multiply $\left(x_{t-1}',s_{t-1}'\right)$ and $\left(x_{t-1}',s_{t-1}^{*'}\right)$
respectively, where $s_{t-1}^{*'}=\left(\varepsilon_{t-1}^{*'}B^{*'},\ldots,\varepsilon_{t-q}^{*'}B^{*'}\right)$,
and take expectations such that
\begin{gather}
\left(\gamma_{1},\ldots,\gamma_{p},k_{1}BD^{2}B',\ldots,k_{q}BD^{2}B'\right)=\nonumber \\
\begin{pmatrix}a_{1} & \cdots & a_{p} & b_{1} & \cdots & b_{q}\end{pmatrix}\mathbb{E}\left(\begin{pmatrix}x_{t-1}\\
s_{t-1}
\end{pmatrix}\begin{pmatrix}x_{t-1}' & s_{t-1}'\end{pmatrix}\right)\label{eq:YW_sys1}
\end{gather}
and
\begin{gather}
\left(\gamma_{1},\ldots,\gamma_{p},k_{1}B^{*}D^{*2}B^{*'},\ldots,k_{q}B^{*}D^{*2}B^{*'}\right)-\cdots\nonumber \\
\cdots-\begin{pmatrix}a_{1}^{*} & \cdots & a_{p}^{*} & b_{1}^{*} & \cdots & b_{q}^{*}\end{pmatrix}\mathbb{E}\left(\begin{pmatrix}x_{t-1}\\
s_{t-1}^{*}
\end{pmatrix}\begin{pmatrix}x_{t-1}' & s_{t-1}^{*'}\end{pmatrix}\right)=0.\label{eq:YW_sys2}
\end{gather}
Since the stationary solution $\left(y_{t}\right)_{t\in\mathbb{Z}}$
of \eqref{eq:system} depends only on past inputs, we obtain that
the right-hand-side of the equation is zero. The square matrices in
\eqref{eq:YW_sys1} and \eqref{eq:YW_sys2} are non-singular due to
the coprimeness assumption on $\left(a(z),b(z)\right)$ and the full-rank
assumption on $\left(a_{p},b_{q}\right)$, compare \citet{deistler83}.
The elements in this matrix correspond either to autocovariances or
can be obtained as, e.g., $\mathbb{E}\left(y_{t-1}\varepsilon_{t-1}^{'}B'\right)=\mathbb{E}\left[\left(\sum_{j=0}^{\infty}k_{j}B\varepsilon_{t-1-j}\right)\varepsilon_{t-1}^{'}B'\right]=k_{0}B\Sigma^{2}B'$.
Now, it follows that $\begin{pmatrix}a_{1} & \cdots & a_{p} & b_{1} & \cdots & b_{q}\end{pmatrix}=\begin{pmatrix}a_{1}^{*} & \cdots & a_{p}^{*} & b_{1}^{*} & \cdots & b_{q}^{*}\end{pmatrix}$
and $B\varepsilon_{t}=B^{*}\varepsilon_{t}^{*}$ because both equation
system involve the same second moments (in particular $B\Sigma^{2}B'=B^{*}\Sigma^{*2}B^{*'}$).
The remainder of the proof follows from what was discussed below
Lemma \ref{lem:Kagan}.
\end{proof}
While the proof above is easily understandable for readers who know
the paper \citet{LMS_svarIdent16}, the following proof uses less
matrix algebra.
\begin{proof}
A different way to prove this theorem uses spectral factorization
arguments (in the guise of linear projections and the Wold representation
theorem). Starting from the stationary solution $\left(y_{t}\right)_{t\in\mathbb{Z}}$
of \eqref{eq:system}, we project $y_{t}$ on its infinite past in
order to obtain the linear innovation $v_{t}$, i.e. $y_{t}-Proj\left(y_{t}|y_{t-1},y_{t-2},\ldots\right)=v_{t}$.
Note that $Proj\left(y_{t}|y_{t-1},y_{t-2},\ldots\right)=Proj\left(y_{t}|B\varepsilon_{t-1},B\varepsilon_{t-2},\ldots\right)=Proj\left(y_{t}|v_{t-1},v_{t-2},\ldots\right)$
because the linear space spanned by the components of $\left\{ y_{t-1},y_{t-2},\ldots\right\} $
coincides with the linear space spanned by the components of $\left\{ v_{t-1},v_{t-2},\ldots\right\} $
and the linear space spanned by the components of $\left\{ B\varepsilon_{t-1},B\varepsilon_{t-2},\ldots\right\} $\footnote{Note that not only the projection is unique (as follows from the projection
theorem) but also the representation in the given basis $\left(y_{t-1},\ldots,y_{t-p},B\varepsilon_{t-1},\ldots,B\varepsilon_{t-q}\right)$
because of the assumptions that $\left(a(z),b(z)\right)$ be left-coprime
and that $\left(a_{p},b_{q}\right)$ be of full rank.}. Now, knowing that the inputs $\varepsilon_{t}$ to \eqref{eq:system}
are not only uncorrelated but also independent and non-Gaussian, we
factorize the covariance matrix of $v_{t}$ as $\mathbb{E}\left(v_{t}v_{t}'\right)=B\mathbb{E}\left(\varepsilon_{t}\varepsilon_{t}'\right)B'$
where $\mathbb{E}\left(\varepsilon_{t}\varepsilon_{t}'\right)=\Sigma^{2}$
is diagonal. It is obvious that it is impossible to distinguish between
$B^{*}\varepsilon_{t}^{*}=B\Sigma^{*}Q\Sigma^{-1}\varepsilon_{t}$
and $B\varepsilon_{t}$ for any orthogonal matrix $Q$ by second moments
only. The rest of the proof is the same as above.
\end{proof}
The difference in these two proofs is as follows. In the first proof,
we use model \eqref{eq:system} 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.
\subsection{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 \citet{LMS_svarIdent16}
(which are in turn based on \citet{IlmonenPaindaveine11} and \citet{HallinMehta15}).
The \textit{first identification} scheme, which is convenient for
deriving asymptotic properties and which we refer to as \textbf{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 \textit{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 \citet{LMS_svarIdent16} including an example
of a matrix or which the above identification schemes are not defined.
The \textit{third identification} \textit{scheme}, similar to the
one in \citet{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.
\section{Parameter Estimation}
In this section, we treat local ML estimation of \eqref{eq:system}.
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 \citet{LMS_svarIdent16} are deferred to the Online Appendix.
\subsection{\label{subsec:asy_lowlevel}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{assu:densities}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
\begin{align}
L_{T}\left(\theta\right) & =\frac{1}{T}\sum_{t=1}^{T}l_{t}\left(\varepsilon_{t}(\theta).\theta\right)\label{eq:likelihood}
\end{align}
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 \citet{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
\begin{thm}
\label{thm:generic_mle}For \eqref{eq:likelihood}, the following
conditions are assumed to be true:
\begin{enumerate}
\item \label{enu:generic_mle_solution}There exists a sequence of estimators
$\left(\hat{\theta}_{T}\right)$ converging almost surely to an interior
point $\theta_{0}\in\Theta$ for which $L_{\theta,T}\left(\hat{\theta}_{n}\right)=o_{P}\left(\frac{1}{\sqrt{T}}\right)$.
\item \label{enu:generic_mle_stat_ergod}$\left(\varepsilon_{t}\right)$
is stationary and ergodic with density $f\left(x,\theta_{0}\right)$
\item \label{enu:generic_mle_mds}$l_{\theta,t}\left(\varepsilon_{t}(\theta).\theta\right)$
is an MDS.
\item \label{enu:generic_mle_density}For the parametric family $\left\{ f\left(x.\theta\right)\ |\ \theta\in\Theta\subseteq\mathbb{R}^{k},\ x\in\mathbb{R}^{n}\right\} $
of densities it holds that $f\left(x,\theta\right)>0$ for all $\left(x,\theta\right)$
and that $f\left(x,\theta\right)$ is twice continuously differentiable
with respect to $\theta$ in an open neighborhood centered at $\theta_{0}$
for all $x$.
\item \label{enu:generic_mle_clt_opg}The individual contributions to the
standardized log-likelihood function satisfy $\mathbb{E}\left(\left\Vert l_{\theta,t}\left(\varepsilon_{t}(\theta_{0}),\theta_{0}\right)\right\Vert ^{2}\right)<\infty$.
\item \label{enu:generic_mle_ulln}There exists a (non-singleton) compact
set $\Theta_{0}$ such that $\mathbb{E}\left(\sup_{\theta\in\Theta_{0}}\left\Vert l_{\theta\theta,t}\left(\varepsilon_{t}(\theta),\theta\right)\right\Vert \right)<\infty$.
\item \label{enu:generic_mle_non_singular}The Hessian matrix $\mathbb{E}\left(l_{\theta\theta,t}\left(\varepsilon_{t}(\theta_{0}),\theta_{0}\right)\right)$
is non-singular.
\item \label{enu:generic_mle_hess_equal2_opg}At the true parameter value,
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, i.e. $\mathbb{E}\left(l_{\theta\theta,t}\left(\varepsilon_{t}(\theta_{0}),\theta_{0}\right)\right)=-\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)$
holds.
\end{enumerate}
Under 1) to 8), we obtain that
\[
\sqrt{T}\left(\hat{\theta}_{T}-\theta_{0}\right)\xrightarrow{d}\mathcal{N}\left(0,\left[\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]^{-1}\right).
\]
\end{thm}
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{enu:generic_mle_solution})
is necessary for the existence of a consistent sequence $\left(\hat{\theta}_{T}\right)$
and is, together with point \ref{enu:generic_mle_stat_ergod}), \ref{enu:generic_mle_mds}),
\ref{enu:generic_mle_density}), and \ref{enu:generic_mle_clt_opg}),
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{enu:generic_mle_ulln}), \ref{enu:generic_mle_non_singular})
and \ref{enu:generic_mle_hess_equal2_opg})\footnote{Point \ref{enu:generic_mle_hess_equal2_opg}) 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:
\begin{thm}
If in addition to the assumptions of the above Theorem \ref{thm:generic_mle},
we assume that
\[
\mathbb{E}\left(\sup_{\theta\in\Theta_{0}}\left\Vert l_{\theta,t}\left(\varepsilon_{t}(\theta),\theta\right)\right\Vert ^{2}\right)<\infty
\]
then we obtain that
\[
B_{T}=\frac{1}{T}\sum_{t=1}^{T}\left[l_{\theta,t}\left(\varepsilon_{t}\left(\hat{\theta}_{T}\right),\hat{\theta}_{T}\right)l_{\theta,t}\left(\varepsilon_{t}\left(\hat{\theta}_{T}\right),\hat{\theta}_{T}\right)'\right]\xrightarrow{p}\mathbb{E}\left[l_{\theta,t}\left(\varepsilon_{t}\left(\hat{\theta}_{T}\right),\hat{\theta}_{T}\right)l_{\theta,t}\left(\varepsilon_{t}\left(\hat{\theta}_{T}\right),\hat{\theta}_{T}\right)'\right].
\]
\end{thm}
\subsection{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{assu:CLT}
and \ref{assu:ULLN} below (together with Assumptions \ref{assu:non_gaussianIID},
\ref{assu:densities}, and \ref{assu:paramSpace}) imply the ones
in the generic Theorem \ref{thm:generic_mle} 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 \citet{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.
\begin{assumption}
\label{assu:paramSpace}The true parameter value $\theta_{0}$ belongs
to the permissible parameter space $\Theta=\Theta_{\pi}\times\Theta_{\beta}\times\Theta_{\sigma}\times\Theta_{\lambda},$
where
\begin{enumerate}
\item $\Theta_{\pi}=\Theta_{\pi_{2}}\times\Theta_{\pi_{3}}$ with $\Theta_{\pi_{2}}\subseteq\mathbb{R}^{n^{2}p}$
and $\Theta_{\pi_{3}}\subseteq\mathbb{R}^{n^{2}q}$ are such that
condition \eqref{eq:stability}, \eqref{eq:invertibility}, the coprimeness
assumption and the full rank assumption on $\left(a_{p},b_{q}\right)$
are satisfied, and
\item $\Theta_{\beta}=vecd\text{�}\left(\mathcal{B}\right)=\left\{ \beta\in\mathbb{R}^{n(n-1)}\,|\,\beta=vecd\text{�}\left(B\right)\text{ for some }B\in\mathcal{B}\right\} $.
The vector $\beta$ collects the off-diagonal elements of $B$.
\item For the scalings, $\Theta_{\sigma}=\mathbb{R}_{+}^{n}$ holds, and
\item for the additional parameters appearing in the component densities,
we have $\Theta_{\lambda}=\Theta_{\lambda_{1}}\times\cdots\times\Theta_{\lambda_{n}}\subseteq\mathbb{R}^{d}$
with $\Theta_{\lambda_{i}}\subseteq\mathbb{R}^{d_{i}}$ open for every
$i\in\left\{ 1,\ldots,n\right\} $ and $d=d_{1}+\cdots+d_{n}$.
\end{enumerate}
\end{assumption}
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
\begin{assumption}
\label{assu:densities}For each $i\in\left\{ 1,\ldots,n\right\} $
the distribution of the error term $\varepsilon_{i,t}$ has a (Lebesgue)
density $f_{i,\sigma_{i}}\left(x;\lambda_{i}\right)=\sigma_{i}^{-1}f_{i}\left(\sigma_{i}^{-1}x;\lambda_{i}\right)$
which may also depend on a parameter vector $\lambda_{i}\in\mathbb{R}^{d_{i}}$.
\end{assumption}
Thus, the individual contributions in the (standardized) log-likelihood
function \eqref{eq:likelihood} are
\begin{equation}
l_{t}\left(\theta\right)=\sum_{i=1}^{n}\log\left[f_{i}\left(\sigma_{i}^{-1}\iota_{i}^{'}B\left(\beta\right)^{-1}u_{t}\left(\theta\right);\lambda_{i}\right)\right]-\log\left\{ \left|\det\left[B\left(\beta\right)\right]\right|\right\} -\sum_{i=1}^{n}\log\left(\sigma_{i}\right),\label{eq:likelihood_individual}
\end{equation}
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
\begin{align*}
\frac{\partial l_{t}\left(\theta\right)}{\partial\pi_{2}} & =-x_{b,t-1}\left(\theta\right)B'\left(\beta\right)^{-1}\Sigma^{-1}e_{x,t}\left(\theta\right)\\
\frac{\partial l_{t}\left(\theta\right)}{\partial\pi_{3}} & =-w_{b,t-1}\left(\theta\right)^{'}\Sigma^{-1}e_{x,t}\left(\theta\right).\\
\frac{\partial l_{t}\left(\theta\right)}{\partial\beta} & =-H'\sum_{i=1}^{q}\left(B\left(\beta\right)^{-1}u_{t-i}\left(\theta\right)\otimes b_{i}'B'\left(\beta\right)^{-1}\Sigma^{-1}e_{x,t}\left(\theta\right)\right)\\
& \qquad-H'\left(B\left(\beta\right)^{-1}u_{t}\left(\theta\right)\otimes B'\left(\beta\right)^{-1}\Sigma^{-1}e_{x,t}\left(\theta\right)\right)\\
& \qquad-H'vec\left(B'\left(\beta\right)^{-1}\right)\\
\frac{\partial}{\partial\sigma}l_{t}\left(\theta\right) & =-\Sigma^{-2}\left[e_{x,t}\left(\theta\right)\odot\varepsilon_{t}\left(\theta\right)+\sigma\right]\\
\frac{\partial}{\partial\lambda}l_{t}\left(\theta\right) & =e_{\lambda,t}\left(\theta\right)
\end{align*}
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.
\begin{assumption}
\label{assu:CLT}The following conditions hold for $i\in\left\{ 1,\ldots,n\right\} $.
\begin{enumerate}
\item For all $x\in\mathbb{R}$ and all $\lambda_{i}\in\Theta_{0,\lambda_{i}},\ f_{i}\left(x;\lambda_{i}\right)>0$
and $f_{i}\left(x;\lambda_{i}\right)$ is twice continuously differentiable
with respect to $\left(x;\lambda_{i}\right)$.
\item The function $f_{i,x}\left(x;\lambda_{i,0}\right)$ is integrable
with respect to x, i.e., $\int\left|f_{i,x}\left(x;\lambda_{i,0}\right)\right|dx<\infty$
.
\item For all $x\in\mathbb{R}$
\[
x^{2}\frac{f_{i,x}^{2}\left(x;\lambda_{i}\right)}{f_{i}^{2}\left(x;\lambda_{i}\right)}\ and\ \frac{\left\Vert f_{i,\lambda_{i}}\left(x;\lambda_{i}\right)\right\Vert ^{2}}{f_{i}^{2}\left(x;\lambda_{i}\right)}
\]
are dominated by $c_{1}\left(1+\left|x\right|^{c_{2}}\right)$ with
$c_{1},c_{2}\geq0$ and $\int\left|x\right|^{c_{2}}f_{i}\left(x;\lambda_{i,0}\right)dx<\infty$
\item $\int\sup_{\lambda_{i}\in\Theta_{0,\lambda_{i}}}\left\Vert f_{i,\lambda_{i}}\left(x;\lambda_{i,0}\right)\right\Vert dx<\infty$.
\end{enumerate}
\end{assumption}
and
\begin{assumption}
\label{assu:ULLN}The following conditions hold for $i\in\left\{ 1,\ldots,n\right\} $.
\begin{enumerate}
\item The functions $f_{i,xx}\left(x;\lambda_{i,0}\right)$ and $f_{i,x\lambda_{i}}\left(x;\lambda_{i,0}\right)$
are integrable with respect to $x$, i.e.,
\[
\int\left|f_{i,xx}\left(x;\lambda_{i,0}\right)\right|dx<\infty\ and\ \int\left\Vert f_{i,x\lambda_{i}}\left(x;\lambda_{i,0}\right)\right\Vert dx<\infty.
\]
\item $\int\sup_{\lambda_{i}\in\Theta_{0,\lambda_{i}}}\left\Vert f_{i,\lambda_{i}\lambda_{i}}\left(x;\lambda_{i,0}\right)\right\Vert dx<\infty$
\item For all $x\in\mathbb{R}$ and all $\lambda_{i}\in\Theta_{0,\lambda_{i}}$,
\[
\frac{f_{i,x}^{2}\left(x;\lambda_{i}\right)}{f_{i}^{2}\left(x;\lambda_{i}\right)}\text{ and }\left|\frac{f_{i,xx}\left(x;\lambda_{i}\right)}{f_{i}\left(x;\lambda_{i}\right)}\right|
\]
are dominated by $a_{0}\left(1+\left|x\right|^{a_{1}}\right)$,
\[
\left\Vert \frac{f_{i,x\lambda_{i}}\left(x;\lambda_{i}\right)}{f_{i}\left(x;\lambda_{i}\right)}\right\Vert \text{and }\left\Vert \frac{f_{i,x}\left(x;\lambda_{i}\right)}{f_{i}\left(x;\lambda_{i}\right)}\frac{f_{i,\lambda_{i}}\left(x;\lambda_{i}\right)}{f_{i}\left(x;\lambda_{i}\right)}\right\Vert
\]
are dominated by $a_{0}\left(1+\left|x\right|^{a_{2}}\right)$,
\[
\left\Vert \frac{f_{i,\lambda_{i}}\left(x;\lambda_{i}\right)}{f_{i}\left(x;\lambda_{i}\right)}\right\Vert ^{2}\text{and }\left\Vert \frac{f_{i,\lambda_{i}\lambda_{i}}\left(x;\lambda_{i}\right)}{f_{i}\left(x;\lambda_{i}\right)}\right\Vert
\]
are dominated by $a_{0}\left(1+\left|x\right|^{a_{3}}\right),$ with
$a_{0},a_{1},a_{2},a_{3}\geq0$ such that $\int\left(\left|x\right|^{2+a_{1}}+\left|x\right|^{1+a_{2}}+\left|x\right|^{a_{3}}\right)f_{i}\left(x;\lambda_{i,0}\right)dx<\infty$.
\end{enumerate}
\end{assumption}
In combination, these assumptions allow to prove (in the Online
Appendix)
\begin{thm}
\label{thm:MLE}Under assumptions 2-5, there exists a sequence of
maximizers $\hat{\theta}_{T}$ of (\ref{eq:likelihood}) such that
$\sqrt{T}\left(\hat{\theta}_{T}-\theta_{0}\right)$ converges in distribution
to $\mathcal{N}\left\{ 0,\mathbb{E}\left[l_{\theta,t}\left(\theta_{0}\right)l_{\theta,t}'\left(\theta_{0}\right)\right]^{-1}\right\} $.
\end{thm}
\section{Empirical Application}
\subsection{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. \citet[Chapter 9.4 and 11.7 ]{luet05}
and \citet[Chapter 4]{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 \citep[page 15]{HannanDeistler12}. Then,
the impulse responses and variance decompositions are obtained in
the same way as in the SVAR case, see e.g. \citep[page 108]{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 \citet{Ravenna07} and
\citet{PoskittYao17}. \citet{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. \citet{PoskittYao17} decompose the truncation
error introduced in \citet{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.
\subsection{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.
\subsubsection{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.
\begin{center}
\begin{figure}[h]
\includegraphics[width=17cm]{cf_plot_data_wrap}\caption{Raw data}
\end{figure}
\par\end{center}
\subsubsection{Estimation Procedure}
In order to select appropriate integer-valued parameters $p$ and
$q$, we estimate a number of VARMA models with the \texttt{MTS}-package
and select the one with the smallest AIC value. It is worth noting
that neither the \texttt{dse}-package nor the \texttt{MTS}-package
enforces the stability condition \eqref{eq:stability} or the invertibility
condition \eqref{eq:invertibility}. 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 \citet{MahdiMcLeod_R_portes}, we choose $p=2$
and $q=2$. Moreover, the distribution of the shocks of the initial
model in Figure \ref{fig:shocks_initial_model} and Figure \ref{fig:qqplot_initial_model}
suggest, and the Jarque-Bera test indicates that the individual series
are not normally distributed.
\textcolor{red}{}
\begin{figure}
\textcolor{red}{\includegraphics[height=6cm]{cfa_hist_shocks}}
\caption{\textcolor{red}{\label{fig:shocks_initial_model}}Histogram of shocks
$\hat{\varepsilon}_{t}=\Sigma^{-1}B^{-1}\hat{a}(z)^{-1}\hat{b}(z)y_{t}$
based on initial model}
\end{figure}
\begin{figure}
\includegraphics[height=6cm]{cfa_qq_shocks}
\caption{\label{fig:qqplot_initial_model}Quantile-quantile plot of shocks
$\hat{\varepsilon}_{t}=\Sigma^{-1}B^{-1}\hat{a}(z)^{-1}\hat{b}(z)y_{t}$
based on initial model}
\end{figure}
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 \citet{LMS_svarIdent16} for fixing
a particular permutation and scaling. The estimates for $B$ and $\sigma$
are
\[
\hat{B}=\begin{pmatrix}1 & 0.1224 & -0.1282\\
-0.0168 & 1 & 0.0107\\
0.0280 & 0.175 & 1
\end{pmatrix}
\]
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)=\begin{pmatrix}0 & 0.1204 & 0.01496\\
0.02816 & 0 & 0.0156\\
0.06903 & 0.0928 & 0
\end{pmatrix}
\]
and $\hat{\mathbb{V}}_{bs}\left(\hat{\sigma}\right)=\left(0.00427,\,0.00199,\,0.01419\right)$.
\subsubsection{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}{}
\begin{figure}[H]
\includegraphics[width=16cm]{cea_irf_alternative}
\caption{Impulse responses}
\end{figure}
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.
\section{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 \textit{Finnish Grid and Cloud
Infrastructure} with persistent identifier \textit{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.
\section{Conclusion}
In this article, we showed that stable and invertible SVARMA models
\eqref{eq:system} 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 \citet{LMS_svarIdent16}.
SVARMA models capture (macroeconomic) dynamics more parsimoniously
than SVAR models and are therefore advantageous in situations with
relatively small sample sizes.
\begin{thebibliography}{38}
\expandafter\ifx\csname urlstyle\endcsname\relax
\else
\fi
\bibitem[Athanasopoulos and Vahid(2008{a})]{AthVahid08}
George Athanasopoulos and Farshid Vahid.
\newblock Varma versus var for macroeconomic forecasting.
\newblock \emph{Journal of Business \& Economic Statistics}, 26\penalty0
(2):\penalty0 237--252, 2008{a}.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1198/073500107000000313}.
\bibitem[Athanasopoulos and Vahid(2008{b})]{AthVahid_JTSA_08}
George Athanasopoulos and Farshid Vahid.
\newblock A complete varma modelling methodology based on scalar components.
\newblock \emph{Journal of Time Series Analysis}, 29\penalty0 (3):\penalty0
533--554, 2008{b}.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1111/j.1467-9892.2007.00568.x}.
\bibitem[Boubacar~Mainassara and Francq(2011)]{BoubacarFrancq11}
Yacouba Boubacar~Mainassara and Christian Francq.
\newblock Estimating structural varma models with uncorrelated but
non-independent error terms.
\newblock \emph{Journal of Multivariate Analysis}, 102\penalty0 (3):\penalty0
496 -- 505, 2011.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1016/j.jmva.2010.10.009}.
\bibitem[Chen and Bickel(2005)]{ChenBickel05}
Aiyou Chen and Peter~J. Bickel.
\newblock Consistent independent component analysis and prewhitening.
\newblock \emph{IEEE Transactions on Signal Processing}, 53:\penalty0
3625--3632, 2005.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1109/TSP.2005.855098}.
\bibitem[Deistler(1983)]{deistler83}
Manfred Deistler.
\newblock {The Properties of the Parameterization of ARMAX Systems and Their
Relevance for Structural Estimation and Dynamic Specification}.
\newblock \emph{Econometrica}, 51\penalty0 (4):\penalty0 1187--1207, July 1983.
\newblock URL \texttt{http://www.jstor.org/stable/1912058}.
\bibitem[Deistler and Seifert(1978)]{DeistlerSeifert78}
Manfred Deistler and Hans-G\"unther Seifert.
\newblock {Identifiability and Consistent Estimability in Econometric Models}.
\newblock \emph{Econometrica}, 46\penalty0 (6):\penalty0 969--980, July 1978.
\newblock URL \texttt{http://www.jstor.org/stable/1909759}.
\bibitem[Gilbert(2015)]{Gilbert_R_dse}
Paul~D. Gilbert.
\newblock \emph{dse: Dynamic Systems Estimation (Time Series Package)}, 2015.
\newblock URL \texttt{https://cran.r-project.org/package=dse}.
\bibitem[Golub and Van~Loan(2013)]{GvL13}
Gene~H. Golub and Charles~F. Van~Loan.
\newblock \emph{Matrix Computations}.
\newblock The Johns Hopkins University Press, 4th edition, 2013.
\bibitem[Gomez(2015)]{Gomez_matlab_15}
Victor Gomez.
\newblock {SSMMATLAB}: A set of matlab programs for the statistical analysis of
state space models.
\newblock \emph{Journal of Statistical Software}, 66\penalty0 (9):\penalty0
1--37, 2015.
\newblock ISSN 1548-7660.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.18637/jss.v066.i09}.
\bibitem[Gomez(2016)]{Gomez16}
Victor Gomez.
\newblock \emph{{Multivariate Time Series With Linear State Space Structure}}.
\newblock Springer, 2016.
\bibitem[Gourieroux et~al.(2017)Gourieroux, Monfort, and
Renne]{GourierouxZakoianRenne17}
Christian Gourieroux, Alain Monfort, and Jean-Paul Renne.
\newblock {Statistical inference for independent component analysis:
Application to structural VAR models}.
\newblock \emph{Journal of Econometrics}, 196:\penalty0 111--126, 2017.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1016/j.jeconom.2016.09.007}.
\bibitem[Gourieroux et~al.(2019)Gourieroux, Monfort, and
Renne]{GourierouxMR_svarma19}
Christian Gourieroux, Alain Monfort, and Jean-Paul Renne.
\newblock {Identification and Estimation in Non-Fundamental Structural VARMA
Models}.
\newblock \emph{Review of Economic Studies}, pages 1--39, 2019.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.193/restud/rdz028}.
\bibitem[Hallin and Mehta(2015)]{HallinMehta15}
Marc Hallin and Chintan Mehta.
\newblock R-estimation for asymmetric independent component analysis.
\newblock \emph{Journal of the American Statistical Association}, 110\penalty0
(509):\penalty0 218--232, 2015.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1080/01621459.2014.909316}.
\bibitem[Hannan(1970)]{Hannan70}
Edward~J. Hannan.
\newblock \emph{Multiple Time Series}.
\newblock Wiley, 1970.
\bibitem[Hannan(1971)]{Hannan71}
Edward~J. Hannan.
\newblock The identification problem for multiple equation systems with moving
average errors.
\newblock \emph{Econometrica}, 39\penalty0 (5):\penalty0 751--765, September
1971.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.2307/1909577}.
\bibitem[Hannan and Deistler(2012)]{HannanDeistler12}
Edward~J. Hannan and Manfred Deistler.
\newblock \emph{The Statistical Theory of Linear Systems}.
\newblock SIAM Classics in Applied Mathematics, Philadelphia, 2012.
\bibitem[Harville(1997)]{Harville97}
David~A. Harville.
\newblock \emph{{Matrix Algebra From a Statistician's Perspective}}.
\newblock Springer, 1997.
\bibitem[Ilmonen and Paindaveine(2011)]{IlmonenPaindaveine11}
Pauliina Ilmonen and Davy Paindaveine.
\newblock Semiparametrically efficient inference based on signed ranks in
symmetric independent component models.
\newblock \emph{Annals of Statistics}, 39\penalty0 (5):\penalty0 2448--2476,
2011.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1214/11-AOS906}.
\bibitem[Kagan et~al.(1973)Kagan, Linnik, and Rao]{Kagan73}
Abram~M. Kagan, Yuri~V. Linnik, and Calyampudi~R. Rao.
\newblock \emph{{Characterization Problems in Mathematical Statistics}}.
\newblock John Wiley \& Sons, 1973.
\bibitem[Kilian and L\"utkepohl(2017)]{KilianLut17}
Lutz Kilian and Helmut L\"utkepohl.
\newblock \emph{Structural Vector Autoregressive Analysis}.
\newblock Cambridge University Press, 2017.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1017/9781108164818}.
\bibitem[Lanne et~al.(2017)Lanne, Meitz, and Saikkonen]{LMS_svarIdent16}
Markku Lanne, Mika Meitz, and Pentti Saikkonen.
\newblock Identification and estimation of non-gaussian structural vector
autoregressions.
\newblock \emph{Journal of Econometrics}, 196\penalty0 (2):\penalty0 288 --
304, 2017.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1016/j.jeconom.2016.06.002}.
\bibitem[L\"utkepohl(1996)]{luet_mat96}
Helmut L\"utkepohl.
\newblock \emph{Handbook of Matrices}.
\newblock John Wiley \& Sons Ltd., 1996.
\bibitem[L\"utkepohl(2005)]{luet05}
Helmut L\"utkepohl.
\newblock \emph{New Introduction to Multiple Time Series Analysis}.
\newblock Springer Berlin, 2005.
\bibitem[Magnus and Neudecker(2007)]{MagnusNeudecker07}
Jan~R. Magnus and Heinz Neudecker.
\newblock \emph{{Matrix Differential Calculus with Applications in Statistics
and Econometrics}}.
\newblock John Wiley \& Sons, 2007.
\bibitem[Mahdi and McLeod(2018)]{MahdiMcLeod_R_portes}
Esam Mahdi and A.~Ian McLeod.
\newblock \emph{portes: Portmanteau Tests for Univariate and Multivariate Time
Series Models}, 2018.
\newblock R package version 3.0.
\bibitem[Meitz and Saikkonen(2013)]{meitzSaikkonen13}
Mika Meitz and Pentti Saikkonen.
\newblock Maximum likelihood estimation of a noninvertible arma model with
autoregressive conditional heteroskedasticity.
\newblock \emph{Journal of Multivariate Analysis}, 114:\penalty0 227 -- 255,
2013.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1016/j.jmva.2012.07.015}.
\bibitem[Poskitt(2016)]{Poskitt16}
Donald~S. Poskitt.
\newblock Vector autoregressive moving average identification for macroeconomic
modeling: A new methodology.
\newblock \emph{Journal of Econometrics}, 192:\penalty0 468--484, June 2016.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1016/j.jeconom.2016.02.011}.
\bibitem[Poskitt and Yao(2017)]{PoskittYao17}
Donald~S. Poskitt and Wenying Yao.
\newblock Vector autoregressions and macroeconomic modeling: An error taxonomy.
\newblock \emph{Journal of Business \& Economic Statistics}, 35\penalty0
(3):\penalty0 407--419, 2017.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1080/07350015.2015.1077139}.
\bibitem[P\"otscher and Prucha(1997)]{poetpruch97}
Benedikt P\"otscher and Ingmar Prucha.
\newblock \emph{Dynamic Nonlinear Econometric Models, Asymptotic Theory}.
\newblock Springer Berlin, 1997.
\bibitem[{R Core Team}(2019)]{Rcore}
{R Core Team}.
\newblock \emph{R: A Language and Environment for Statistical Computing}.
\newblock R Foundation for Statistical Computing, Vienna, Austria, 2019.
\newblock URL \texttt{https://www.R-project.org/}.
\bibitem[Raghavan et~al.(2016)Raghavan, Athanasopoulos, and
Silvapulle]{RaghavanAthSilvapulle16}
Mala Raghavan, George Athanasopoulos, and Param Silvapulle.
\newblock Canadian monetary policy analysis using a structural varma model.
\newblock \emph{Canadian Journal of Economics}, 49\penalty0 (1):\penalty0
347--373, 2016.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1111/caje.12200}.
\bibitem[Ravenna(2007)]{Ravenna07}
Federico Ravenna.
\newblock Vector autoregressions and reduced form representations of \{DSGE\}
models.
\newblock \emph{Journal of Monetary Economics}, 54\penalty0 (7):\penalty0 2048
-- 2064, 2007.
\newblock ISSN 0304-3932.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.1016/j.jmoneco.2006.09.002}.
\bibitem[Rothenberg(1971)]{Rothenberg71}
Thomas~J. Rothenberg.
\newblock {Identification in Parametric Models}.
\newblock \emph{Econometric Theory}, 39\penalty0 (3):\penalty0 577--591, May
1971.
\newblock doi: \begingroup \urlstyle{rm}\Url{10.2307/1913267}.
\bibitem[Scherrer and Deistler(2019)]{ScherrerDeistler2019_handbook}
Wolfgang Scherrer and Manfred Deistler.
\newblock Vector autoregressive moving average models.
\newblock In Hrishikesh~D. Vinod and C.~R. Rao, editors, \emph{Handbook of
statistics 41}, volume~41. North-Holland, 2019.
\bibitem[Scherrer and Funovits(2019)]{Scherrer_rldm}
Wolfgang Scherrer and Bernd Funovits.
\newblock \emph{rldm: A package for modeling of time series with a rational
spectral density}, 2019.
\bibitem[Seber(2008)]{seber08}
George A.~F. Seber.
\newblock \emph{{A Matrix Handbook for Statisticians}}.
\newblock John Wiley \& Sons, 2008.
\bibitem[Tsay(2013)]{Tsay13}
Ruey~S. Tsay.
\newblock \emph{Multivariate Time Series Analysis With R and Financial
Applications}.
\newblock John Wiley \& Sons, Inc., Hoboken, New Jersey, 2013.
\bibitem[Tsay and Wood(2018)]{Tsay_R_MTS}
Ruey~S. Tsay and David Wood.
\newblock \emph{MTS: All-Purpose Toolkit for Analyzing Multivariate Time Series
(MTS) and Estimating Multivariate Volatility Models}, 2018.
\newblock R package version 1.0.
\end{thebibliography}
\pagebreak{}