Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,272 characters · 15 sections · 82 citation commands
Robust Estimation of the non-Gaussian Dimension in Structural Linear Models
{\it Keywords:} SVARMA models; Normality Test; Matrix Rank; Higher Order Spectrum.
Central to assessing and quantifying the effects of economic shocks over real activity variables and other macroeconomic outcomes is the identification of such disturbances. For achieving this, it is usual to employ structural vector auto-regressive models and to apply a particular identification strategy. Most of these strategies resort on imposing external restrictions, which may come from a theoretical model, a particular data frequency or by exploiting granular institutional information (see sims1980macroeconomics,blanchard1989dynamic,bernanke1998measuring,blanchard2002empirical,gali1992well as examples of imposing zero, linear and non-linear identification restrictions, economically motivated, over VAR parameters; uhlig2005effects,canova2007price,arias2018inference for agnostic economically motivated restrictions; romer1989does,romer2010macroeconomic for narrative identification strategy; and mertens2013dynamic,mertens2014reconciliation,stock2018identification for instrumental variables identification approach).
Alternatively, the statistical identification strategy (hereafter SIS) of structural VAR models has appeared on the scene as part of a data-driven identification approach and has gained relevance, because this strategy permits to identify a SVAR model without employing external information and makes feasible to evaluate any economically motivated restriction. The seminal work of comon1994independent states that a linear transformation of a vector whose components are mutually independent and with at most one Gaussian distributed element is identified up to signed permutation (see lanne2017identification,maxand2020identification,guay2021identification for a direct application within the context of a causal SVAR model). For general linear, stationary processes, the non-Gaussian requirement needs to be satisfied by each of the structural shocks. Particularly, chan2006note show that identification of a linear, possibly non-causal, stationary process needs that structural shocks vector must be fully independent \oldFootnote{Being fully independent implies that a shocks is serial and contemporaneously independent.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi and each shock must have non-zero and finite third and fourth order cumulants, respectively. This result allows is the cornerstone to statistically identify possibly non-fundamental SVARMA models \oldFootnote{Non-fundamental behavior is equivalent to AR and MA lag polynomials have some roots inside the unit circle}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi (see lanne2010structural,gourieroux2020identification for applications with likelihood methods and velasco2022identification for an equivalent identification result through higher-order cumulant conditions).
In contrast to other identification strategies, the SIS permits to evaluate the validity of its identification assumptions. For instance, hong1999hypothesis develops a procedure for testing the serial independence assumption. Also, amengual2022moment provide a framework for evaluating the contemporaneous independence and the non-Gaussianity of structural shocks. In this paper, we focus on developing a procedure for evaluating the non-Gaussian identification condition, since even though serial and mutual independence hold, non-Gaussian behavior of structural shocks is decisive for making feasible application of the SIS. Besides, as it is detailed below, our proposal seeks to be implemented without requiring estimation of structural model. Thus, the procedure can be performed before applying the SIS.
There is a voluminous literature coping with the problem of determining whether a random variable (or vector) is Gaussian distributed or not. Some procedures are based on the empirical distribution of data (see kolmogorov1933sulla,smirnov1948table,massey1951kolmogorov,shapiro1965analysis); other approaches employ the characteristic function (see hall1983test,epps1983test); while others exploit third and fourth order centered moments (see d1971omnibus,bera1982model,lobato2004simple). These procedures only assess the hypothesis of joint Gaussianity, which in case of no rejection, would make infeasible to apply the SIS to identify structural parameters in a SVARMA model. Nonetheless, if joint Gaussianity were rejected, this would support only the existence of at least one non-Gaussian component in the system, and this conclusion would not be sufficient for justifying the application of the SIS.
The literature of estimating the non-Gaussian dimension or the number of non-Gaussian distributed shocks inside a random vector is less prolific. nordhausen2017asymptotic propose both asymptotic and bootstrap tests to estimate the non-Gaussian dimension of an unobservable random vector by analyzing the number of non-zero eigenvalues of a convenient transformation of the scatter matrix \oldFootnote{The scatter matrix contains fourth-order moments of a vector.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi constructed using information from an observable random vector, which is an affine transformation of the unobservable one. \oldFootnote{Let $\mathbf{x}$ and $\boldsymbol{\varepsilon}$ be an observable and unobservable $d$-dimensional random vectors, respectively. An affine transformation is define as $\mathbf{x} = \boldsymbol{b}_{0} + \boldsymbol{B}\boldsymbol{\varepsilon}$, with $\boldsymbol{b}_{0}$ possibly non-zero vector and $\boldsymbol{B}$ time-invariant, full rank, square matrix.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi They find the null eigenvalues of the scatter matrix is equal to the number of Gaussian distributed components in the vector of unobservable shocks. \oldFootnote{To be completely precise, there exists $r$ eigenvalues equal to $d+2$, but nordhausen2017asymptotic employ a normalized version of the scatter matrix.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi Within the context of a causal SVAR models, maxand2020identification applies nordhausen2017asymptotic approach directly to estimated structural shocks, because according to her, the non-Gaussian block remains identified even though there is more than one Gaussian shock in the system. Following a similar approach, amengual2022moment estimates structural VAR model by pseudo MLE and applies Normality test to each estimated structural error. Alternatively, guay2021identification approach does not require to estimate the structural model, but only the reduced-form VAR model. He relates the non-Gaussian dimension to the rank of the matrix version of fourth order cumulants of the reduced-form errors. \oldFootnote{When third and fourth order moments exist, the third order cumulant is the same as the asymmetry coefficient, while the fourth order cumulant represents the excess of kurtosis. Gaussian distribution is the only one with all cumulants of third or higher order equal to zero. See marcinkiewicz1939propriete for technical details.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi Guay's approach reduces the original task to a matrix rank estimation problem, thus he directly applies robin2000tests sequential procedure.
These approaches are helpful but restrictive for some reasons. First, they consider only the particular case of having a causal structural VAR model as the true, underlying data generating process. However, many structural macroeconomic models can be represented and fitted more accurately by an SVARMA model rather than an SVAR one (see ravenna2007vector,fernandez2007abcs). Second, as discussed in velasco2022identification,gourieroux2020identification and the references therein, the fundamentalness assumption is necessary for avoiding the dynamic identification problem -i.e., for any stationary linear process, its fundamental and non-fundamental representations are observationally equivalent when only second-order information is exploited. And, as surveyed by alessi2011non, non-fundamentalness may be consistent with many macroeconomic models and data features.
Consequently, this article attempts to fill this gap in the literature and proposes a procedure for determining the non-Gaussian dimension in a SVARMA model, without requiring prior knowledge about the roots location of the lag polynomials. Nonetheless, this is a quite challenging task, since the reduced-form (hereafter RF) errors from a fundamental VARMA approximation -or RF-VARMA model-, when the roots location of lag polynomials is unknown, might not be a simple static rotation of structural shocks, but a dynamic filter of them. This characteristic makes invalid guay2021identification, maxand2020identification and amengual2022moment approaches, because the rank or the number of non-zero eigenvalues of a matrix constructed from contemporaneous fourth-order cumulants of the RF errors are not related one-to-one to the non-Gaussian dimension, unless extreme assumptions are made. In contrast, we exploit third or fourth-order cumulant spectrums of RF errors for constructing a matrix whose rank unveils the non-Gaussian dimension in the vector of structural shocks. Our work can be seen as the extension of lobato2004simple approach to a multivariate context and the generalization of guay2021identification and amengual2022moment works to situations in which roots location is unknown.
\textcolor{red}{CHANGE THIS PART!!} For estimating the rank of a matrix, there exist several strategies in the literature. We follow kleibergen2006generalized (hereafter KP) approach. The KP statistic is built from the singular value decomposition of the matrix of interest. The asymptotic distribution of the statistic is a standard chi-square whose degrees of freedom change depending on the null hypothesis. Unlike KP work, our context changes the asymptotic distribution of the test statistic for some particular rank values, specifically under joint Gaussianity. For other null hypotheses, the asymptotic distribution is chi-square, but the degrees of freedom are generally unknown. Thus, we propose a bootstrap strategy. This path implies another challenge: to impose the null hypothesis in the resampled data. Montecarlo exercises were performed to analyze the size and power of the bootstrap test; these show that our strategy estimates satisfactorily the non-Gaussian dimension. We apply our procedure to two well-known macroeconomic datasets. Our proposal detects a skewed structural shock in the system described by blanchard1989dynamic (hereafter BQ). Using blanchard2002empirical (hereafter BP), our approach detects at least two skewed and non-mesokurtic structural errors, unlike guay2021identification whose procedure only could detect one non-mesokurtic structural shock. In the case of imposing roots outside the unit circle, this latter result implies that the SIS can be applied to BQ or BP datasets.
The remainder of this paper is structured as follows: Section 2 describes our time series model and states the main assumptions. Section 3 shows the connection between the number of non-Gaussian structural shocks and the rank of a matrix constructed from third and fourth order spectrum of reduced-form errors. In Section 4, the test procedure and estimation are detailed. Section 5 shows simulation results and the empirical application. We conclude in Section 6.
Let $\boldsymbol{y}_{t}$ be a $d$-dimensional stationary, zero-mean process generated by a structural vector autoregression moving-average (SVARMA) model \oldFootnote{See gourieroux2020identification,velasco2022identification,mainassara2011estimating for alternative representations}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi
where $\boldsymbol{\Phi}(L)=\operatorname{\boldsymbol{I}}_{d}-\sum_{j=1}^{p}{\boldsymbol{\Phi}_{j}L^{j}}$, $\boldsymbol{\Theta}(L)=\operatorname{\boldsymbol{I}}_{d}+\sum_{j=1}^{q}{\boldsymbol{\Theta}_{j}L^{j}}$. $L$ is the lag (back-shift) operator, i.e. $L^{j}\boldsymbol{y}_{t}=\boldsymbol{y}_{t-j}$ for any $j\in\mathbb{Z}$. $\boldsymbol{B}$ is a time-invariant, full-rank, squared matrix, also known as the contemporaneous effects matrix. $\boldsymbol{\varepsilon}_{t}$ is a $d$-dimensional vector, representing the structural shocks. \oldFootnote{These shocks have economic interpretation.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi
The $\boldsymbol{\Phi}(z)$ and $\boldsymbol{\Theta}(z)$ represents the autoregressive (AR) and moving-average (MA) polynomials, respectively. These hold
Condition ((ref)) rules out only the presence of unit roots in both AR and MA polynomials, otherwise a stationary solution to equation ((ref)) does not exist. Unlike current approaches in the literature, this requirement does not impose an exact location of AR and MA polynomials roots \oldFootnote{Compares this to alternative works such as maxand2020identification,guay2021identification,amengual2022moment, where they assume that condition in ((ref)) holds for any $z\in\mathbb{T}_{+}=\left\{ z\in\mathbb{C}\;|\;|z|\leq1 \right\}$.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi, implying that the model in ((ref)) allows for possibly non-fundamental behavior.
The structural shocks, $\boldsymbol{\varepsilon}_{t}$, have a specific behavior summarized by the following assumption.
Assumption ((ref).(i)) is standard in the literature \oldFootnote{See for instance lanne2010structural,lanne2017identification,gourieroux2020identification,velasco2022identification,guay2021identification.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi and makes the analysis simpler. However, it is a restrictive requisite since it implies structural shocks do not exhibit any linear or nonlinear serial dependence \oldFootnote{This condition rules out some characteristics that may be relevant in the empirical analysis of macroeconomic outcomes such as conditional heteroskedasticity.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. Assumption ((ref).(ii)) is similar to the one stated in comon1994independent. It is the basis of Independent Component Analysis (hereafter ICA). This requirement is not as restricting as it seems, e.g., it is standard in macroeconomics to assume that productivity shock is independent of monetary or fiscal policy shocks.
Assumption ((ref).(iii)) does not impose a particular number of non-Gaussian shocks, because this will be determined empirically. This implies the structural model cannot be estimated, since it is unidentified. Besides, if $d_{ng}$ is the number of non-Gaussian components in $\boldsymbol{\varepsilon}_{t}$, the rest $d-d_{ng}$ structural errors are Gaussian distributed. Finally, Assumption ((ref).(iv)) imposes structural errors are centered and standardized with finite moments up to eighth order. This last requirement is necessary, because we will characterize non-Gaussian behavior using third and fourth order information.
The vector of structural parameters of the model in equation ((ref)) can be represented by $\boldsymbol{\vartheta}_{0}$, a $K$-dimensional column vector. We split $\boldsymbol{\vartheta}_{0}$ into two blocks: one related to dynamic behavior ($\boldsymbol{\vartheta}_{0,1}$) and another that governs the static component ($\boldsymbol{\vartheta}_{0,2}$). Hence, the structural model can be written as $\boldsymbol{\Phi}(L,\boldsymbol{\vartheta}_{0,1})\boldsymbol{y}_{t} = \boldsymbol{\Theta}(L,\boldsymbol{\vartheta}_{0,1})\boldsymbol{B}(\boldsymbol{\vartheta}_{0,2})\boldsymbol{\varepsilon}_{t}$. The moving-average representation of the observable variables is $\boldsymbol{y}_{t}=\boldsymbol{\Psi}(L,\boldsymbol{\vartheta}_{0,1})\boldsymbol{B}(\boldsymbol{\vartheta}_{0,2})\boldsymbol{\varepsilon}_{t}$, where $\boldsymbol{\Psi}(L,\boldsymbol{\vartheta}_{0,1})\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\boldsymbol{\Phi}^{-1}(L,\boldsymbol{\vartheta}_{0,1})\boldsymbol{\Theta}(L,\boldsymbol{\vartheta}_{0,1})$ is a possibly non-causal filter. \oldFootnote{A non-causal filter is of the form $\boldsymbol{\Psi}(L)=\sum_{j=-\infty}^{\infty}{\boldsymbol{\Psi}_{j}L^{j}}$ with non-zero $\boldsymbol{\Psi}_{j}$ for some $j<0$.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi The vector of structural parameters, $\boldsymbol{\vartheta}_{0}$, is not identified, because Assumption ((ref)) does not assure the non-Gaussian behavior of the structural shocks. Identification assuming that all roots of AR and MA polynomials are outside unit circle requires $d_{ng}\geq{d-1}$ (comon1994independent); while, when roots location is unknown and only unit roots are discarded $d_{ng}=d$ is needed (chan2006note,velasco2020identification).
Instead of working with structural estimates, which requires identification of $\boldsymbol{\vartheta}_{0}$, in this research we propose a procedure that exploits only reduced-form estimates, for which structural parameters identification is unnecessary. According to Wold's Decomposition Theorem (hereafter WDT), any stationary process $\boldsymbol{y}_{t}$ can be represented as a square summable, infinite, causal moving-average of serially uncorrelated errors (see anderson2011statistical for more details). This means that $\boldsymbol{y}_{t}=\tilde{\boldsymbol{\Psi}}(L)\boldsymbol{u}_{t}$ with $\tilde{\boldsymbol{\Psi}}(L)=\sum_{j=0}^{\infty}{\tilde{\boldsymbol{\Psi}}_{j}L^{j}}$ such that $\sum_{j=0}^{\infty}{\left\|{\tilde{\boldsymbol{\Psi}}_{j}}\right\|^{2} }<\infty$, $\tilde{\boldsymbol{\Psi}}_{0}=\boldsymbol{I}_{d}$ and $\operatorname{\mathbb{C}ov}\left( \boldsymbol{u}_{t}, \boldsymbol{u}_{t-k} \right)=\boldsymbol{0}$ for all $\left|{k}\right|\geq1$. Without loss of generality, we can assume that $\tilde{\boldsymbol{\Psi}}(L)=\tilde{\boldsymbol{\Phi}}^{-1}(L)\tilde{\boldsymbol{\Theta}}(L)$ with $\det\left(\tilde{\boldsymbol{\Phi}}(z)\right)\det\left(\tilde{\boldsymbol{\Theta}}(z)\right)\neq{0}$ for all $\left|{z}\right|\leq1$. Furthermore, let $\boldsymbol{\vartheta}_{f}$ denotes the parameters vector of the fundamental approximation of $\boldsymbol{y}_{t}$, thus $\boldsymbol{y}_{t}=\boldsymbol{\Psi}(L,\boldsymbol{\vartheta}_{f})\boldsymbol{u}_{t}=\boldsymbol{\Phi}^{-1}(L,\boldsymbol{\vartheta}_{f})\boldsymbol{\Theta}(L,\boldsymbol{\vartheta}_{f})\boldsymbol{u}_{t}$. From this and equation ((ref)), the serially uncorrelated RF errors, $\boldsymbol{u}_{t}$, are equal to
Notice that if the structural model in equation ((ref)) were fundamental, then $\boldsymbol{\delta}(L,\boldsymbol{\vartheta}_{f},\boldsymbol{\vartheta}_{0})=\boldsymbol{B}(\boldsymbol{\vartheta}_{0,2})$, which means that RF errors are a simply rotation of contemporaneous, unobserved structural shocks. Nonetheless, when only unit roots are discarded from AR or MA polynomials and the structural model is possibly non-fundamental, then $\boldsymbol{\delta}(L,\boldsymbol{\vartheta}_{f},\boldsymbol{\vartheta}_{0})$ is a non-causal filter. This implies that RF errors are not a static rotation of structural shocks, but a linear combination of leads and lags of the structural errors.
For a simpler exposition, we write $\boldsymbol{\delta}(L,\boldsymbol{\vartheta}_{f},\boldsymbol{\vartheta}_{0})=\boldsymbol{\delta}(L,\boldsymbol{\vartheta}_{f})$. The no serial correlation of RF errors implies that
where $\boldsymbol{\Omega}$ is positive definite, symmetric, time-invariant matrix; $\boldsymbol{\delta}^{*}(e^{i\lambda},\boldsymbol{\vartheta}_{f})$ denotes the transposed, conjugate of matrix polynomial $\boldsymbol{\delta}(e^{i\lambda},\boldsymbol{\vartheta}_{f})$. velasco2022identification,baggio2018parametrization and the references therein denominate filters satisfying condition ((ref)) as all-pass filters. \oldFootnote{To be completely precise, an all-pass filter requires that $\boldsymbol{\Omega}=\boldsymbol{I}_{d}$. However, this is not a problem in our analysis because $\boldsymbol{\Omega}$ is a constant matrix.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi
In order to determine the non-Gaussian dimension in the vectors of structural shocks $\boldsymbol{\varepsilon}_{t}$ and RF errors vector $\boldsymbol{u}_{t}$ it is necessary to characterize non-Gaussian behavior. Like jarque1987test,lobato2004simple,guay2021identification, we typify non-Gaussianity as deviations of the higher order information respect to the their values under Gaussianity. And for assessing higher order information we exploit third and fourth order cumulants.
The cumulant generating function of vector $\boldsymbol{\varepsilon}_{t}$ is $\boldsymbol{\kappa}_{\boldsymbol{\varepsilon}}(\boldsymbol{\tau})\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\log\left(\phi_{\boldsymbol{\varepsilon}}(\boldsymbol{\tau})\right)$, where $\phi_{\boldsymbol{\varepsilon}}(\boldsymbol{\tau})\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\mathbb{E}(\exp{(i\boldsymbol{\tau}^{\prime}\boldsymbol{\varepsilon})})$ represents the characteristic function of $\boldsymbol{\varepsilon}_{t}$ and $\boldsymbol{\tau}=(\tau_{1},\dots,\tau_{d})^{\prime}$ a $d$-dimensional real-valued vector. A cumulant of order $k$ (or the cumulant of a $k$-tuple $\{\varepsilon_{j_{m},t}\}_{m=1}^{k}$ ) is defined as
where $\sum_{j=1}^{d}{k_{j}}=k\geq1$ and $k_{j}\geq{0}$. $j_{m}\in\{1,\dots,d\}$ for all $m=1,\dots,k$. When moments (up to order $k$) exist, jammalamadaka2006higher shows that the expression in equation ((ref)) is equal to
where $\mathfrak{p}$ is a partition of set $\{j_{1},\dots,j_{k}\}$, $\left|{\mathfrak{p}}\right|$ represents the number of parts in partition $\mathfrak{p}$. $P$ is an element of $\mathfrak{p}$. Although the expression in ((ref)) seems difficult to handle it, let see some simple examples. When $k=2$ and we choose the duplet $(\varepsilon_{1,t},\varepsilon_{2,t})$, then there exists only two partitions of set $\{1,2\}$, which are $\mathfrak{p}_{1}=\{1,2\}$ and $\mathfrak{p}_{2}=\{\{1\},\{2\}\}$. Thus $\operatorname{\mathbb{C}um}(\varepsilon_{1,t},\varepsilon_{2,t})=\operatorname{\mathbb{E}}[\varepsilon_{1,t}\varepsilon_{2,t}]-\operatorname{\mathbb{E}}[\varepsilon_{1,t}]\operatorname{\mathbb{E}}[\varepsilon_{2,t}]=\operatorname{\mathbb{C}ov}(\varepsilon_{1,t},\varepsilon_{2,t})$.
When $k=3$ and the chosen triplet is $(\varepsilon_{1,t},\varepsilon_{2,t},\varepsilon_{3,t})$, there are five possible partitions $\mathfrak{p}_{1}=\{1,2,3\}$, $\mathfrak{p}_{2}=\{\{1\},\{2,3\}\}$, $\mathfrak{p}_{3}=\{\{2\},\{1,3\}\}$, $\mathfrak{p}_{4}=\{\{3\},\{1,2\}\}$ and $\mathfrak{p}_{5}=\{\{1\},\{2\},\{3\}\}$. Thus, $\operatorname{\mathbb{C}um}(\varepsilon_{1,t},\varepsilon_{2,t},\varepsilon_{3,t})=\operatorname{\mathbb{E}}[\varepsilon_{1,t}\varepsilon_{2,t}\varepsilon_{3,t}]-\operatorname{\mathbb{E}}[\varepsilon_{1,t}]\operatorname{\mathbb{E}}[\varepsilon_{2,t}\varepsilon_{3,t}]-\operatorname{\mathbb{E}}[\varepsilon_{2,t}]\operatorname{\mathbb{E}}[\varepsilon_{1,t}\varepsilon_{3,t}]-\operatorname{\mathbb{E}}[\varepsilon_{3,t}]\operatorname{\mathbb{E}}[\varepsilon_{1,t}\varepsilon_{2,t}]+2\operatorname{\mathbb{E}}[\varepsilon_{1,t}]\operatorname{\mathbb{E}}[\varepsilon_{2,t}]\operatorname{\mathbb{E}}[\varepsilon_{3,t}]=\operatorname{\mathbb{E}}[\varepsilon_{1,t}\varepsilon_{2,t}\varepsilon_{3,t}]$.
When the elements in the $k$-tuple satisfies $j_{1}=j_{2}=\cdots=j_{k}=m$, then we call $\operatorname{\mathbb{C}um}\left( \varepsilon_{m,t},\dots,\varepsilon_{m,t} \right)$ the $m$-th marginal cumulant of order $k$, and it is represented by $\kappa_{k,m}^{\varepsilon}$. Collecting all the cumulants of order $k$ in a single-column vector, we obtain
where $\boldsymbol{\upsigma}_{k}(\{1,\dots,d\})$ is the set of all permutations of length $k$ formed with numbers in the set $\{1,\dots,d\}$. For instance, if $d=2$ and $k=2$, $\boldsymbol{\upsigma}_{2}(\{1,\dots,2\})=\left\{ \{1,1\}; \{1,2\}; \{2,1\}; \{2,2\} \right\}$. From expression in ((ref)), it clear that $\boldsymbol{\kappa}_{k}(\boldsymbol{\varepsilon})$ is $d^{k}$-dimensional real-valued vector. The number of non-repeated elements in $\boldsymbol{\kappa}_{k}(\boldsymbol{\varepsilon})$ is $\binom{d+k-1}{k}=\frac{(d+k-1)!}{k!(d-1)!}$.
Let $\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}$ be a $d^{k-1}\times{d}$ matrix such that $\boldsymbol{\kappa}_{k}(\boldsymbol{\varepsilon})=\operatorname{\text{vec}}(\mathbf{v}_{k}^{\boldsymbol{\varepsilon}})$, i.e, $\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}$ is a matrix version of vector $\boldsymbol{\kappa}_{k}(\boldsymbol{\varepsilon})$. By assumptions ((ref).(ii)-(iv)), all non-marginal cumulants of $\boldsymbol{\varepsilon}_{t}$ are zero, then $\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}$ has the following structure
where $\boldsymbol{e}_{j}$ is the $j$-th canonical vector in $\mathbb{R}^{d}$.
Since we are characterizing non-Gaussian behavior by its deviation from Gaussianity using third and fourth order information, then the number of non-Gaussian distributed components in the vector of structural shocks $\boldsymbol{\varepsilon}$ is equal to number of non-zero third or fourth order marginal cumulants $\kappa_{k,m}^{\varepsilon}$. And, from expression ((ref)), the number of non-zero third order cumulants is equal to the rank of matrix $\mathbf{v}_{3}^{\boldsymbol{\varepsilon}}$; while the number of non-zero fourth order cumulants is $\operatorname{\text{rank}}(\mathbf{v}_{4}^{\boldsymbol{\varepsilon}})$. Besides, the rank of the block matrix $\mathbf{v}_{34}^{\boldsymbol{\varepsilon}}=
$ gives the joint number of non-zero third and fourth order cumulants. Therefore, it can be said that $\operatorname{rank}(\mathbf{v}_{3}^{\boldsymbol{\varepsilon}})$ reveals the amount of asymmetric shocks, $\operatorname{rank}(\mathbf{v}_{4}^{\boldsymbol{\varepsilon}})$ provides the number of non-mesokurtic shocks.
This result is behind the approach followed by maxand2020identification work. Given that her work restricts the focus to a causal (fundamental) SVAR model, RF errors are a simple rotation of contemporaneous value of structural shocks. And, even though the whole vector of structural errors is not identified, she demonstrates that the non-Gaussian block it is. Therefore, she employs nordhausen2017asymptotic strategy for estimating the rank or the non-Gaussian dimension.
guay2021identification works with a similar context as maxand2020identification, i.e., assuming a causal (fundamental) SVAR as the generating process for observables. But, unlike Maxand's work, he does not exploit $\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}$ directly. As we explained above, the RF errors, under fundamentalness, are a linear rotation of structural shocks, that is $\boldsymbol{u}_{t}=\boldsymbol{B}(\boldsymbol{\vartheta}_{0,2})\boldsymbol{\varepsilon}_{t}$ where $\boldsymbol{\vartheta}_{0,2}$ governs the contemporaneous or static part of the fundamental SVAR model and is not identified. According to jammalamadaka2006higher, the vector of $k$-th order cumulant of the vector of RF errors is $\boldsymbol{\kappa}_{k}^{\boldsymbol{u}} = \boldsymbol{B}^{\otimes{k}}\boldsymbol{\kappa}_{k}({\boldsymbol{\varepsilon}})$ where $\boldsymbol{B}^{\otimes{k}}$ denotes the $k$-th Kronecker power. \oldFootnote{Notice that we are employing $\boldsymbol{B}\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\boldsymbol{B}(\boldsymbol{\vartheta}_{0,2})$.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi Then, similarly to $\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}$ we define
Since $\boldsymbol{B}$ is a square, full rank matrix, then $\operatorname{\text{rank}}\left( \mathbf{v}_{k}^{\boldsymbol{u}} \right)=\operatorname{\text{rank}}\left( \mathbf{v}_{k}^{\boldsymbol{\varepsilon}} \right)$. Therefore, it can be concluded that under fundamentalness, the non-Gaussian dimension in the structural model can be identified by calculating the rank of a matrix constructed with third or fourth cumulants of RF errors.
Nonetheless, under a more general structural model like the one considered in this paper, maxand2020identification and guay2021identification approaches are invalid. First, a possibly non-fundamental structural VARMA is only identified if all shocks are non-Gaussian, which we ignore. Hence, it would be no reliable any approach that employs estimated structural shocks. Second, as it showed in the expression ((ref)), RF errors from a fundamental VARMA model are a filtered version of structural shocks. In this case, $\boldsymbol{\kappa}_{k}^{\boldsymbol{u}}$ and $\mathbf{v}_{k}^{\boldsymbol{u}}$ are
From this latter expression, it is visible that $\operatorname{\text{rank}}\left( \mathbf{v}_{k}^{\boldsymbol{u}} \right)\neq \operatorname{\text{rank}}\left( \mathbf{v}_{k}^{\boldsymbol{\varepsilon}} \right)$ in general, because even though the rank of $\boldsymbol{\delta}_{j}^{\otimes{k-1}}\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}\boldsymbol{\delta}_{j}^{\prime}$ may coincide with $\operatorname{\text{rank}}\left(\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}\right)$, this term may be positive or negative definite. Thus, the rank of the summation could be higher or lower than $\operatorname{\text{rank}}\left(\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}\right)$. \oldFootnote{Notice a simple example of having two matrices $A=
$ and $B=
$. Both are full rank matrices, but the sum $C=A+B=
$ is rank-deficient.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi Consequently, the rank of $\mathbf{v}_{k}^{\boldsymbol{u}}$ does not identify anymore our parameter of interest, $\operatorname{rank}(\mathbf{v}_{k}^{\boldsymbol{\varepsilon}})$.
To overcome this issue and instead of using only third or fourth information of contemporaneous RF errors, we exploit the third and fourth information at all leads and lags of RF errors. Hence, we employ the $k$-th order cumulant spectral density of the RF errors. According to brillinger2001time, the cumulant spectrum of order $k$ is
where $\bigotimes_{j=1}^{n}{\boldsymbol{A}_{j}}\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\boldsymbol{A}_{1}\otimes\boldsymbol{A}_{2}\otimes\cdots\otimes{\boldsymbol{A}_{n}}$.
Similarly to $\mathbf{v}_{k}^{\boldsymbol{u}}$, we define the $d^{k-1}\times{d}$ matrix $G_{k}^{u}(\boldsymbol{\lambda})$, such that $g_{k}^{u}(\boldsymbol{\lambda})=\operatorname{\text{vec}}\left( G_{k}^{u}(\boldsymbol{\lambda}) \right)$. Thus
Based on $G_{k}^{u}(\boldsymbol{\lambda})$, let define the $d\times{d}$ matrix $G_{k}^{u,2}(\lambda_{k})\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\left[G_{k}^{u}(\boldsymbol{\lambda})\right]^{*}G_{k}^{u}(\boldsymbol{\lambda})$ where $\boldsymbol{A}^{*}$ denotes the conjugate transpose of the complex-valued matrix $A$. Since $\boldsymbol{\delta}(z,\boldsymbol{\vartheta}_{f})$ is an all-pass filter and using $(A\otimes{B})(C\otimes{D})=\left(AC\otimes{BD}\right)$, provided that $AC$ and $BD$ are conformable for multiplication, we obtain
with $\lambda_{k}=-\sum_{m=1}^{k-1}{{\lambda}_{m}}$. Given that $G_{k}^{u,2}(\lambda_{k})$ only depends on the scalar frequency $\lambda_{k}\in[-\pi,\pi]$, thus, from now on, we can drop the sub-index $k$.
Since $\boldsymbol{\delta}(z,\boldsymbol{\vartheta}_{f})$ is full rank for any $z\in\mathbb{C}$ such that $\left|{z}\right|=1$, then it is clear that $\operatorname{\text{rank}}\left( G_{k}^{u,2}(\lambda) \right) = \operatorname{\text{rank}}\left( \mathbf{v}_{k}^{\boldsymbol{\varepsilon}\prime}\boldsymbol{\Omega}^{\otimes{k-1}}\mathbf{v}_{k}^{\boldsymbol{\varepsilon}} \right)=\operatorname{\text{rank}}\left(\mathbf{v}_{k}^{\boldsymbol{\varepsilon}\prime}\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}\right)$. In consequence, the non-Gaussian dimension in $\boldsymbol{\varepsilon}_{t}$ can be identified by $\operatorname{\text{rank}}\left( G_{k}^{u,2}(\lambda) \right)$. The following proposition summarizes our main finding.
Proposition ((ref)) can be interpreted as follows: the non-Gaussian dimension in the vector of structural shocks is equal to the rank of a matrix constructed from higher order cumulant spectrum of RF errors at a given frequency. In particular, if only the third-order cumulant spectrum is employed, we obtain the non-Gaussian dimension delivered by asymmetric non-Gaussian shocks. If only the fourth-order spectrum is utilized, we capture the non-Gaussian dimension spanned by non-mesokurtic shocks.
Besides, notice that in case the condition in ((ref)) holds for all $z\in\mathbb{T}_{+}=\{x\in\mathbb{C}\;|\;\left|{x}\right|\leq{1}\}$, higher order cumulant spectrums are constant. Hence, array $G^{u,2}_{k}(\lambda)$ is constant, i.e., the same for any frequency. In particular, $G^{u,2}_{k}(\lambda)=\boldsymbol{B}(\boldsymbol{\vartheta}_{0,2})\left(\mathbf{v}_{k}^{\boldsymbol{\varepsilon}\prime}\boldsymbol{\Omega}^{\otimes{k-1}}\mathbf{v}_{k}^{\boldsymbol{\varepsilon}}\right)\boldsymbol{B}^{\prime}(\boldsymbol{\vartheta}_{0,2})=\mathbf{v}_{k}^{\boldsymbol{u}\prime}\mathbf{v}_{k}^{\boldsymbol{u}}$ for any $\lambda\in[-\pi,\pi]$. This setup is the case analyzed in guay2021identification. Under this situation, proposition ((ref)) holds as well. In other words, if a fundamental structural VARMA model generates $\boldsymbol{y}_{t}$, a rectangular array constructed from cumulant spectrums of order $3$ or $4$ based on RF errors identifies the non-Gaussian dimension.
From the discussion above, to determine the non-Gaussian dimension in the structural shocks vector, we need to find the rank of a rectangular array constructed from the cumulant spectrum of order $k=3,4$. Thus, our empirical problem becomes the estimation of the rank of a matrix. In the literature, this problem has been dealt with in two approaches: sequential testing or estimation by information criteria (see camba2009statistical for a detailed review). We opt for the former approach because it may be problematic to construct a pseudo-likelihood function for a non-causal filter. Additionally, the parametric estimation of our matrix of interest is not feasible because the possibly non-fundamental filter $\boldsymbol{\delta}(z,\boldsymbol{\vartheta}_{f},\boldsymbol{\vartheta}_{0,1})=\boldsymbol{\Psi}^{-1}(z,\boldsymbol{\vartheta}_{f})\boldsymbol{\Psi}(z,\boldsymbol{\vartheta}_{0,1})$ is not identified unless we impose fundamentalness. Consequently, we decide to estimate $g_{k}(\boldsymbol{\lambda})$ non-parametrically.
Dealing with $G_{k}^{u,2}(\lambda)$ may be problematic since it implies working with complex-valued terms. For simplifying the analysis, we use $\text{Re}(G_{k}^{2}(\lambda))$ instead. At frequency zero, there is no problem since $G_{k}^{u,2}(0)=\text{Re}(G_{k}^{u,2}(0))$. At other frequencies, $\text{Re}(G_{k}^{u,2}(\lambda))$ is the sum of two positive semidefinite quadratic forms; thus, we cannot lose rank, and we gain rank only if both quadratic forms has linearly independent columns.
We now describe our estimation of the non-Gaussian dimension in structural shocks through a sequential hypothesis test. At step $s=1,2,\dots,d$ in the sequential procedure, the null and alternative hypotheses are
where $r_{s}$ denotes the rank under the null hypothesis at step $s$.
We start the sequential procedure (at step $s=1$) by imposing $r_{1}=0$, i.e., joint Gaussianity of the structural shocks. If this is rejected, we continue with the next step, $s=2$, and the null hypothesis is updated to $r_{2}=1$, i.e., only one asymmetric or non-mesokurtic structural shock. We continue this way until we cannot reject a null hypothesis or reach the final step $s=d$, where $r_{d}=d-1$. If this is not rejected, the non-Gaussian dimension in the structural shocks vector is $d-1$; on the contrary, the non-Gaussian dimension is $d$.
The literature on rank estimation via hypothesis testing is extensive (see al2017unifying for a complete survey). We center our approach on kleibergen2006generalized (hereafter KP) proposal. They employ the singular value decomposition (hereafter SVD) of a matrix of interest, $\boldsymbol{\Pi}$. The SVD of matrix $\boldsymbol{\Pi}_{m\times{n}}$ consists on finding squared orthonormal matrices $\boldsymbol{R}_{1}$ and $\boldsymbol{R}_{2}$ of dimension $m$ and $n$, respectively; and a quasi-diagonal rectangular matrix $\boldsymbol{\mathcal{L}}_{m\times{n}}$ such that \oldFootnote{See golub2012matrix for more details on SVD of a matrix.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi
where $\boldsymbol{\mathcal{L}}$ is a rectangular array with a squared block $\boldsymbol{\mathcal{L}}_{\overline{m}}=\text{diag}(l_{1},\dots,l_{\overline{m}})$ with $\overline{m}=\min\{m,n\}$ and $l_{1}\geq{l_{2}}\geq\dots\geq{l_{\overline{m}}}\geq{0}$ are the singular values of $\boldsymbol{\Pi}$. If $\boldsymbol{\Pi}$ is squared, then $\boldsymbol{\mathcal{L}}=\boldsymbol{\mathcal{L}}_{\overline{m}}$. Assuming that $m>n$, $\boldsymbol{\mathcal{L}}=
^{\prime}$.
From decomposition in ((ref)) and selecting an integer $r\in\{0,\dots,\overline{m}-1\}$, kleibergen2006generalized obtain that
where $\boldsymbol{C}_{r}\boldsymbol{D}_{r}=\underbrace{
}_{m\times{r}}\underbrace{\boldsymbol{\mathcal{L}}_{1}}_{r\times{r}}\underbrace{
}_{r\times{n}}$; $\boldsymbol{C}_{r,\perp}\boldsymbol{\mathcal{L}}_{r}\boldsymbol{D}_{r,\perp}=\underbrace{
}_{m\times{(m-r)}}\underbrace{\boldsymbol{\mathcal{L}}_{2}}_{(m-r)\times(n-r)}\underbrace{
}_{(n-r)\times{n}}$. $\boldsymbol{C}_{r,\perp}$ and $\boldsymbol{D}_{r,\perp}$ are the orthogonal complements of $\boldsymbol{C}_{r}$ and $\boldsymbol{D}_{r}$, respectively. Thus, \citet{kleibergen2006generalized} approach consists in decomposing the matrix $\boldsymbol{\Pi}$ into two linear independent parts: one that is made of the multiplication of two matrices that are full column rank, $\boldsymbol{C}_{r}\boldsymbol{D}_{r}$, and another with null rank $\boldsymbol{C}_{r,\perp}\boldsymbol{\mathcal{L}}_{r}\boldsymbol{D}_{r,\perp}$. Therefore, $\operatorname{rank}(\boldsymbol{\Pi})=\operatorname{rank}(\boldsymbol{C}_{r}\boldsymbol{D}_{r})=r$. This result is only achieved if we place all the non-zero singular values into $\boldsymbol{\mathcal{L}}_{1}$. Accordingly, the null hypothesis of $\operatorname{rank}(\boldsymbol{\Pi})=r$ is equivalent to $\boldsymbol{\mathcal{L}}_{r}=\boldsymbol{0}$.
Since $\boldsymbol{\Pi}$ is not observable, the SVD is applied to its estimator $\boldsymbol{\widehat{\Pi}}$. Thus, we can write $\boldsymbol{\widehat{\Pi}}=\boldsymbol{\widehat{C}}_{r}\boldsymbol{\widehat{D}}_{r} + \boldsymbol{\widehat{C}}_{r,\perp}\boldsymbol{\widehat{\mathcal{L}}}_{r}\boldsymbol{\widehat{D}}_{r,\perp}$. Under the assumption of $\sqrt{T}(\operatorname{\text{vec}}{[\boldsymbol{\widehat{\Pi}}]}-\operatorname{\text{vec}}[\boldsymbol{\Pi}]) \xrightarrow{d} \mathcal{N}_{mn}(\boldsymbol{0},\boldsymbol{\Xi})$ with $\boldsymbol{\Xi}$ positive definite, then under the null hypothesis of $\operatorname{\text{rank}}(\boldsymbol{\Pi})=r$ the asymptotic distribution of $\boldsymbol{\widehat{\ell}}_{r}=\operatorname{\text{vec}}\left(\boldsymbol{\widehat{\mathcal{L}}}_{r}\right)$ is $\sqrt{T}\boldsymbol{\widehat{\ell}}_{r} \xrightarrow{d} \mathcal{N}_{(m-r)(n-r)}\left(\boldsymbol{0}, \widetilde{\boldsymbol{\Xi}}_{r}\right)$ where $\boldsymbol{\widetilde{\Xi}}_{r}=(\boldsymbol{D}_{r,\perp}\otimes{\boldsymbol{C}_{r,\perp}^{\prime}})\boldsymbol{\Xi}(\boldsymbol{D}_{r,\perp}\otimes{\boldsymbol{C}_{r,\perp}^{\prime}})$. The KP-statistic is
where $\widehat{\boldsymbol{\widetilde{\Xi}}}_{r}$ is the consistent estimator of the asymptotic variance, ${\boldsymbol{\widetilde{\Xi}}}_{r}$.
In some dimensions, our problem departs from kleibergen2006generalized context. First, our matrix of interest $\boldsymbol{\Pi}=\text{Re}({G}_{k}^{u,2}(\lambda))$ is symmetric, i.e. it contains repeated elements. Hence, the asymptotic variance of unrestricted estimator $\boldsymbol{\widehat{\Pi}}$, $\boldsymbol{\Xi}$, is only positive semi-definite. Second, as we detail in the section where it is discussed the estimation of higher order spectrum, the asymptotic distribution of the unrestricted estimator of our matrix of interest changes under the null hypothesis of joint Gaussianity, i.e., under the null of $\operatorname{\text{rank}}(\boldsymbol{\Pi})=0$. Finally, our convergence rates for the asymptotic distribution of the statistic are lower in comparison to the standard speed of convergence, $\sqrt{T}$.
\paragraph{Test Statistic at first step}\newline For the first step in the sequential procedure, the null hypothesis is $\operatorname{\text{rank}}(\boldsymbol{\Pi})=0$, i.e., the vector of structural shocks is an uncorrelated Gaussian process. The KP statistic under this null hypothesis is
where $\boldsymbol{A}^{\dagger}$ denotes the Moore-Penrose inverse of matrix $\boldsymbol{A}$, $\boldsymbol{Q}$ is the vectorized form of a linear combination of an inverse Wishart distribution, which parameters depends on the asymptotic variance of spectrum estimates (see Section $4.3$ below), and $ \overline{\boldsymbol{D}}_{0,\perp}$ and $\overline{\boldsymbol{C}}^{\prime}_{0,\perp}$ represent the limiting values of $\widehat{\boldsymbol{D}}_{0,\perp}$ and $\widehat{\boldsymbol{C}}_{0,\perp}$. $a_{T}$ is the convergence rate of spectrum estimates.
\paragraph{Test Statistic under mixed cases}\newline Now, when the null hypothesis is $\operatorname{\text{rank}}(\boldsymbol{\Pi})=r$ with $r>{0}$, i.e., the vector of structural shocks is an i.i.d process with mutually independent components and there are $r$ non-Gaussian distributed shocks. In this case, the asymptotic distribution of $\boldsymbol{\widehat{\Pi}}$ is Normal, and the KP statistic adopts a similar form as in kleibergen2006generalized,
The statistic in ((ref)) is asymptotically distributed as a chi-square with degrees of freedom $\nu=\operatorname{\text{rank}}{\left[(\widehat{\boldsymbol{D}}_{r,\perp}\otimes{\widehat{\boldsymbol{C}}_{r,\perp}^{\prime}})\widehat{\boldsymbol{\Xi}}(\widehat{\boldsymbol{D}}_{r,\perp}\otimes{\widehat{\boldsymbol{C}}_{r,\perp}^{\prime}})^{\prime}\right]}$. Unfortunately, the exact value for $\nu$ is not easy to determine, because -although $\widehat{\boldsymbol{C}}_{r,\perp}$ and $\widehat{\boldsymbol{D}}^{\prime}_{r,\perp}$ are full column rank matrices of dimension $d\times{(d-r)}$- the matrix $(\widehat{\boldsymbol{D}}_{r,\perp}\otimes{\widehat{\boldsymbol{C}}_{r,\perp}^{\prime}})$ is full row rank matrix of dimension $(d-r)^{2}\times{d^{2}}$ and the matrix $\widehat{\boldsymbol{\Xi}}$ is only positive semidefinite, i.e. it is not full rank. From this discussion, we can bound the degrees of freedom, $1\leq \nu\leq \min\{ \operatorname{\text{rank}}(\boldsymbol{\Xi}), (d-r_{s})^{2} \}$.
Because of the difficulties characterizing the asymptotic distributions of the KP statistic for each step in the sequential procedure, we opt for a bootstrap strategy. Nonetheless, this path entails other challenges. The most problematic is to accomplish that the bootstrap sample appropriately reflects the null hypothesis of each step. Otherwise, it may severely compromise the size and power of the test (see hall1991two,portier2014bootstrap).
It is a well-known result from the SVD of a matrix $\boldsymbol{\Pi}_{m\times{n}}$ with $m>n$ and $\operatorname{\text{rank}}(\boldsymbol{\Pi})=r$, that $\text{null}(\boldsymbol{\Pi})$ is spanned by the last $n-r$ columns of $\boldsymbol{R}_{2}$ and the range $\text{ran}(\boldsymbol{\Pi})$ is spanned by the first $r$ columns of $\boldsymbol{R}_{1}$. Thus, when $\operatorname{\text{rank}}(\boldsymbol{\Pi})=r$, the first $r$ columns of $\boldsymbol{R}_{1}$ span the non-Gaussian dimension, while last $n-r$ columns of $\boldsymbol{R}_{2}$ span the Gaussian one. Thus, we follow nordhausen2017asymptotic approach and use the estimates of $\boldsymbol{R}_{2}$ for a projection of residuals into Gaussian and non-Gaussian dimensions.
Given the sample $\{\boldsymbol{y}_{t}\}_{t=1}^{T}$, we estimate a RF-VARMA model and compute its RF residuals $\{\hat{\boldsymbol{u}}_{t}\}_{t=1}^{T}$. Based on these estimates, we construct our matrix of interest $\widehat{\boldsymbol{\Pi}}=\text{Re}\left(\widehat{G}_{k}^{\Hat{u},2,(T)}(\lambda)\right)$ and the matrices $\widehat{\boldsymbol{\mathcal{L}}}$, $\widehat{\boldsymbol{R}}_{1}$ and $\widehat{\boldsymbol{R}}_{2}$. Under $H_{0,r}:\,\operatorname{\text{rank}}(\boldsymbol{\Pi})=r$ and according to nordhausen2017asymptotic, we construct $\widehat{\boldsymbol{R}}_{2,r}\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\left[ \hat{\mathbf{r}}_{2,r+1}\cdots{\hat{\mathbf{r}}_{2,d} } \right]$ and the projection matrix $\mathbf{M}_{r}=\boldsymbol{I}-\hat{\Sigma}_{\hat{u}}^{1/2}\widehat{\boldsymbol{R}}_{2,r}\widehat{\boldsymbol{R}}_{2,r}^{\prime}\hat{\Sigma}_{\hat{u}}^{-1/2}$, where $\hat{\Sigma}_{\hat{u}}$ is variance of RF residuals vector. $\mathbf{M}_{r}$ is the projection matrix into the orthogonal space to the one spaned by $\widehat{\boldsymbol{R}}_{2,r}$, i.e., $\mathbf{M}_{r}$ is a projection into non-Gaussian space. Based on these inputs, at any step $s$ in the sequential testing procedure, each bootstrap sample is created following Algorithm 1.
It is worth mentioning that the unrestricted bootstrap sample of RF residuals is not performed by the usual independent bootstrap procedure stated by efron1979bootstrap, but by the proposed method in politis1994stationary, called the stationary bootstrap. According to the authors, this procedure is suitable for stationary weakly dependent time series because, unlike other proposals such as kunsch1989jackknife,politis1992general, it exhibits the desirable property that the resampled time series obtained are stationary conditional on the original data. This type of bootstrap sampling is necessary because RF errors are only serially uncorrelated, but they are not independent unless the fundamentalness of data is imposed or Gaussian is assumed. Besides, for obtaining the restricted bootstrap sample for observable data, we employ a parametric approach; otherwise is not possible to impose the null hypothesis in the data.
Once the restricted bootstrap sample has been obtained and we set a nominal size $\alpha\in(0,1)$ for all the steps in the sequence, the test proceeds as follows:
This part briefly discusses some details of estimating higher order cumulant spectrum. The discussion follows closely brillinger1967asymptotic. The sample periodogram for a $k$-tuple of RF errors is
where $\boldsymbol{c}=\left(c_{1},\dots,c_{k}\right)\in\upsigma_{k}\left(\{1,\dots,d\}\right)$. $\lambda_{t_{j}}=2\pi \frac{t_{j}}{T}$ for $t_{j}=1,\dots,T-1$, $j=1,\dots,k-1$ and $\sum_{j=1}^{k}{\lambda_{t_{j}}}=0[\text{mod}(2\pi)]$. Besides, $z_{c_{j}}^{(T)}(\lambda_{t_{j}})=\sum_{t=0}^{T-1}{u_{c_{j},t}e^{-i\lambda_{t_{j}}{t}}}$ is the discrete Fourier transform (DFT) of $c_{j}$-th reduced-form error, $u_{c_{j},t}$.
A consistent estimator of the cumulant spectrum of order $k$ for a $\boldsymbol{c}$-tuple of RF errors is
with $W_{T}(\boldsymbol{a})=H_{T}^{-(k-1)}\sum_{\boldsymbol{j}=(j_{1},\dots,j_{k}):\sum_{m=1}^{k}{(j_{m}+a_{m})}=0}{ W\left( \frac{1}{H_{T}}(\boldsymbol{a}+2\pi\boldsymbol{j}) \right) }$ and $H_{T}$ satisfies $\lim\limits_{T\to{\infty}}{H_{T}}=0$ and $\lim\limits_{T\to\infty}{TH_{T}^{k-1}}=\infty$. The weighting function, $W(\boldsymbol{a})$, is symmetric around $\boldsymbol{0}$ and satisfies conditions stated in brillinger1967asymptotic.
The consistent estimator of $g_{k}^{u}(\boldsymbol{\lambda})$ is the collection of all $\hat{g}_{\boldsymbol{c},k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$ in a vector, which we call $\hat{g}_{k}^{\hat{u},(T)}(\boldsymbol{\lambda})$. Thus, we can obtain the $d^{k-1}\times{d}$ matrix $\hat{G}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$ such that $\hat{g}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda}) = \operatorname{\text{vec}}\left( \hat{G}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda}) \right)$. Besides, $\hat{G}_{k}^{\Hat{u},2,(T)}(\lambda)= \left[\hat{G}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda})\right]^{*}\hat{G}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$. Finally, our matrix of interest is $\text{Re}\left( \hat{G}_{k}^{\Hat{u},2,(T)}(\lambda) \right)$.
brillinger1967asymptotic show that
where $\mathcal{N}_{2d^{k}}$ represents a multivariate Gaussian random $2d^{k}$-dimensional column-vector with mean $\boldsymbol{0}$ and variance $\boldsymbol{\mathcal{V}}(\boldsymbol{\lambda})$; $\boldsymbol{\mathcal{V}}(\boldsymbol{\lambda})=
$ and $\boldsymbol{\Lambda}^{k}(\boldsymbol{\lambda},\boldsymbol{\mu}) = \left[ \Lambda_{(\boldsymbol{c},\boldsymbol{b})}^{k}(\boldsymbol{\lambda},\boldsymbol{\mu}) \right]_{(\boldsymbol{c},\boldsymbol{b})\in[\upsigma_{k}\left(\{1,\dots,d\}\right)]^{2} }$ with
where $\upsigma\boldsymbol{b}=(b_{\upsigma_{1}},\dots,b_{\upsigma_{k}})$, $\upsigma=(\upsigma_{1},\dots,\upsigma_{k})\in\upsigma\left(\{1,\dots,k\}\right)$, $\sum_{m=1}^{k}{\lambda_{m}}=0[\text{mod}(2\pi)]$, $\eta\left({x}\right)=\sum_{j=-\infty}^{\infty}{\delta_{D}(x+2\pi{j})}$ and $\delta_{D}(x)$ is the delta-Dirac function.
The estimator of our matrix of interest is $\widehat{\boldsymbol{\Pi}}=\text{Re}\left\{\hat{G}_{k}^{\Hat{u},2,(T)}(\boldsymbol{\lambda})\right\}$. Besides, given that $\hat{G}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda}) = \operatorname{\text{vec}}^{-1}\left( \hat{g}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda}) \right)$, $\widehat{\boldsymbol{\Pi}}$ is a function of real and imaginary parts of $\hat{g}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$ as the following equation shows
where $\hat{g}_{\text{Re},k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$ and $\hat{g}_{\text{Im},k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$ denote the real and imaginary parts of vector $\hat{g}_{k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$, respectively; and $\mathbf{I}(d)=\operatorname{\text{vec}}(\boldsymbol{I}_{d})\operatorname{\text{vec}}(\boldsymbol{I}_{d})^{\prime}\otimes\boldsymbol{I}_{d^{k-1}}$.
Consequently, the asymptotic distribution of $\operatorname{\text{vec}}\left(\widehat{\boldsymbol{\Pi}}\right)$ depends on the null hypothesis at step $s$, $H_{0,s}$. Particularly, at the first step $H_{0,1}:\,\operatorname{\text{rank}}(\boldsymbol{\Pi})=0$, i.e., under joint Gaussianity of the structural shocks, the population spectrum of order $k$ is zero for $k=3,4$. This result implies that the Jacobian of $\operatorname{\text{vec}}\left(\widehat{\boldsymbol{\Pi}}\right)$ at population values is null, making it not feasible to use the standard Delta method for finding the asymptotic distribution of $\operatorname{\text{vec}}\left(\widehat{\boldsymbol{\Pi}}\right)$. Proposition ((ref)) states the asymptotic distribution for $\operatorname{\text{vec}}\left(\widehat{\boldsymbol{\Pi}}\right)$ at different stages in the sequential testing procedure.
In Proposition (ref), $\boldsymbol{Q}\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\frac{1}{2}\left(I_{d^{2}}\otimes{\mathbf{X}^{\prime}}\right)\boldsymbol{\mathcal{H}}_{\hat{G}^{\Hat{u},2,(T)}_{k}}$ and $\boldsymbol{\mathcal{H}}_{\hat{G}^{\Hat{u},2,(T)}_{k}}$ is the vectorized Hessian \oldFootnote{The Hessian of a vector-valued function $\mathbb{R}^{d}\mapsto{\mathbb{R}^{q}}$ is a $3$-tensor. Intuitively, the Hessian matrix is $3$-dimensional array of dimensions $d\times{d}\times{q}$.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi of $\operatorname{\text{vec}}\left( \text{Re}\left\{\hat{G}_{k}^{\Hat{u},2,(T)}(\lambda)\right\} \right)$ evaluated at zero. $\mathbf{X}\mathrel{\vcenter{\baselineskip0.5ex \lineskiplimit0pt \hbox{\scriptsize.}\hbox{\scriptsize.}}} =\operatorname{\text{vec}}\left[{\mathcal{W}_{1}(\boldsymbol{\mathcal{V}}(\boldsymbol{\lambda}))}\right]$ with ${\mathcal{W}_{1}(\boldsymbol{\mathcal{V}}(\boldsymbol{\lambda}))}$ denotes a Wishart distribution with parameters $\boldsymbol{\mathcal{V}}(\boldsymbol{\lambda})$ and $1$ degree of freedom. Besides, $\boldsymbol{J}_{\hat{G}^{2}_{k}}(\boldsymbol{\lambda})$ is the Jacobian of $\operatorname{\text{vec}}\left( \text{Re}\left\{\hat{G}_{k}^{2,(T)}(\lambda)\right\} \right)$ evaluated at population values of cumulant spectrum of order $k$.
The results above are obtained for unobserved RF errors, $\boldsymbol{u}_{t}$. However, we employ the estimated RF residuals using a finite sample of size $T$ using the estimated RF parameters, $\hat{\boldsymbol{\vartheta}}_{f}$. In consequence, it is essential to show that asymptotic results remain valid when using estimated RF residuals, $\hat{\boldsymbol{u}}_{t}$. The spectrum of order $k$ for the estimated RF residuals can be obtained easily replacing $\boldsymbol{u}_{t}$ by $\hat{\boldsymbol{u}}_{t}$ in ((ref)), denoted by $\hat{g}_{\boldsymbol{c},k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$. The rest of the estimators can be obtained as explained above using $\hat{g}_{\boldsymbol{c},k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$. This analysis is omitted in maxand2020identification and guay2021identification works, maybe explained because the structural linear model is fundamental.
Proposition (ref) states the asymptotic equivalence between estimators of higher-order cumulant spectrum using unobserved RF errors and their sample counterparts. This implies that the discrepancy generated by using $\hat{g}_{\boldsymbol{c},k}^{\Hat{u},(T)}(\boldsymbol{\lambda})$ decreases to zero at the rate of $C H_{T}^{\frac{k-1}{2}}$, which is slower than $\sqrt{T}$. This convergence rate comes from the fact that we are using a non-parametric estimation for the higher order spectrum, which incorporates a kernel with band-with $H_{T}$.
In Corollary (ref), $\hat{\boldsymbol{\kappa}}^{\boldsymbol{u},(T)}_{k}$ and $\hat{\boldsymbol{\kappa}}^{\hat{\boldsymbol{u}},(T)}_{k}$ are the sample estimators for cumulant of order $k$ based on unobserved RF errors and RF residuals, respectively.
This corollary states that in the case of working with a more restricted structural model, since the researcher imposes the location of roots, Proposition ((ref)) still holds. The asymptotic equivalence remains if contemporaneous higher-order cumulants were employed instead of the higher-order spectrum. However, the convergence rate is faster than in the general setting. This latter result is because the estimation does not use kernel smoothing when using contemporaneous higher-order cumulants.
Show the validity of the restricted bootstrap sampling analytically is quite intricate. We believe that following procedures that have been proven to be consistent may assure the effectiveness of our approach, although this is not equivalent to formal proof. Additionally, we present some evidence from different Montecarlo exercises as additional support for the validity of our bootstrap sampling and test.
The data-generating process or the true structural model for all the exercises is a non-causal SVAR($1$). We employ different distributions for the structural disturbances. Montecarlo and bootstrap repetitions are set to $m=250$ and $B=500$, respectively. The nominal significance level is set to $\alpha=5\%$. The number of points for the DFT is set to a minimum even integer greater than or equal to the sample size $T$, and the size of the window is set to $\floor{T/4}$.
Table (ref) shows the rejection rates of our sequential procedure in the bivariate case ($d=2$). Besides, the values in a box represent the size of the test. At panel ((ref)), we employ the zero frequency ($\lambda=0$). The first row is associated with the Gaussian case. Its size is $4\%$. When we consider a mixed case (second row of the panel ((ref))), the rejection rate of $H_{0,1}$ represents the power, $32\%$; while its size is $7\%$. Finally, when structural shocks are fully non-Gaussian distributed (third row in panel ((ref))), both rejection rates of $H_{0,1}$ and $H_{0,2}$ represent the power at each step. These values are $42\%$ and $26\%$, respectively.
Employ estimates at a single frequency of higher order cumulant spectrum may be noisy, leading to imprecise results, especially when the sample size is small. To make our procedure robust to this unpleasant characteristic, we use a grid of frequencies instead of only the zero frequency. Let $\mathcal{G}$ denote the finite grid of frequencies that are selected, instead of constructing our matrix of interest as $\boldsymbol{\Pi}=\sum_{{\lambda}_{k}\in\mathcal{F}}{\text{Re}\left(G_{k}^{2}({\lambda_{k}})\right)}$ where $\mathcal{F}$ denotes a finite grid of frequencies, because it may entail the risk of gaining rank spuriously; we use the following statistic $\overline{KP}_{r}^{(T)}=\max_{\lambda\in\mathcal{F}}{\left\{ KP_{r}^{(T)}(\lambda) \right\}}$, where $KP_{r}^{(T)}(\lambda)$ denotes the KP-statistic for a particular frequency $\lambda$.
In panel ((ref)) of Table (ref), a grid of length $11$ is selected. For this grid, the test size when all disturbances are Gaussian distributed is $4\%$; when only one non-Gaussian shock is present, the size is around $6\%$. The power is quite similar to when a single frequency is employed. In panel ((ref)) of Table (ref), a grid of $26$ frequencies was selected. The size test when $\boldsymbol{\varepsilon}_{t}$ is a Gaussian process is $5\%$; and when we have only non-Gaussian shock, the size is around $6\%$.
Regarding power, the levels remain close to the previous exercises. Furthermore, in Table (ref), we can observe the effect of increasing the sample size. Regarding the test size for each case, the results show that it remains close to the values obtained with $T=250$. Moreover, concerning the power, we observe a significant increment in the power, especially in the case where only one single frequency is employed.
Suppose we take a strategy for restricting the bootstrap sample as the one followed by guay2021identification, which employs the first $r_{s}$ columns of $\boldsymbol{R}_{1}$ for projecting the residuals into the non-Gaussian dimension. In our case, since our matrix of interest is squared and symmetric, $\boldsymbol{R}_{1}=\boldsymbol{R}_{2}$, thus we use the first $r_{s}$ columns of $\boldsymbol{R}_{2}$. The rest of the components $d-r_{s}$ are drawn from a multivariate normal distribution with mean zero and variance $\boldsymbol{I}_{d-r_{s}}$. In Table ((ref)), we compute the rejection rates using Guay's bootstrap sampling approach. The size of the test when having fully Gaussian shocks is, disregarding the length of the frequency grid, around $5\%$. The size when $\boldsymbol{\varepsilon}_{t}$ has one Gaussian component is around $7\%$ for the different grid lengths. Nonetheless, the power using nordhausen2017asymptotic bootstrap sampling is consistently higher than using Guay's approach.
For larger dimensional models, we choose $d=3,4$. We perform these exercises with sample size $T=250$ and the grid $\lambda_{3}=0$. The results can be observed in Table ((ref)). The distribution $MN_{1}$ is a mixture of two normal distributions ($\mathcal{N}(10,0.75)$ and $\mathcal{N}(-2,4)$) such that the distribution has a positive skewness coefficient. It can be noted that, while the sample size is fixed, it is more difficult to reject the null hypothesis for larger $r_{s}$.
We select two well-known data sets in empirical macroeconomics. The first data set is used in the seminal work of blanchard1989dynamic (hereafter BQ), where they introduced the long-run restrictions as an identification scheme for causal structural VAR models. This database contains two endogenous variables: the US GNP's growth and the unemployment rate. The second data set is taken from the work of blanchard2002empirical (hereafter BP). Their data includes three endogenous variables: tax revenues, government spending, and GDP (all in real terms). In both cases, we use the exact SVAR specification for each work, avoiding our conclusions from being affected by a specification bias.
Table (ref) shows the results of applying our proposal to the BQ database. We can identify a single asymmetric structural shock using only third-order cumulant spectral density. When information in the fourth-order spectral cumulant is employed, the proposed method cannot identify any component with excess kurtosis. This result can seem contradictory, though the limited sample size ($T=148$) might affect the precision of fourth-order estimates. On the other hand, Table (ref) shows the results of applying our proposed method to the BP dataset. Using third-order information, we detect at most two asymmetric structural shocks at a significance level of $10\%$. This result is consistently found disregarding if only a single frequency or a grid of frequencies is employed. When fourth-order information is used, at $10\%$ of significance, we can detect three non-mesokurtic structural shocks at zero frequency. When a grid of frequencies is employed instead, the method does detect two non-mesokurtic structural shocks. Besides, if these results are compared to those obtained in guay2021identification, our proposal can detect at least one extra non-mesokurtic structural shock.
From the previous results, provided that structural shocks are independent across time and components, if the researcher imposes fundamentalness, applying the SIS to both BQ and BP datasets is feasible. A causal and invertible SVARMA model can be identified without imposing external identification restrictions or using proxy variables. Moreover, our results suggest that for the BP dataset, it is feasible to identify a possibly non-fundamental SVARMA model. In the case of the BQ dataset, the most limiting factor is the quite small sample size.
This paper aims to design a procedure to determine the number of non-Gaussian shocks in a structural linear VARMA model which is robust to the type of dynamic representation, i.e., whether the structural model is fundamental or non-fundamental. This objective is mainly motivated because knowing the number of non-Gaussian shocks in the structural model allows the researcher to implement an estimation procedure of structural parameters based on the statistical identification approach. In the fundamentalness of the structural model, the requirement is to have at most one Gaussian structural shock; if the researcher does not want to impose the root location, the requirement is that all the structural errors are non-Gaussian distributed.
We generalize the procedure in guay2021identification by exploiting that the rank of a matrix constructed from the third-order cumulant spectrum of the RF errors reveals the number of skewed or asymmetric structural errors in the SVARMA model. Meanwhile, the rank of an array constructed from the fourth-order cumulant spectrum reveals the number of non-mesokurtic structural shocks in the system. Simulation results show that our procedure correctly estimates the non-Gaussian dimension. Additionally, from a practice point of view, our proposal is intensive computationally, especially if we employ the cumulant spectrum of fourth order or when the dimension of the structural model increases or the sample size is large.