LEM
WORKING PAPER SERIES
Identification of Vector Autoregressive Models
with Nonlinear Contemporaneous Structure

Francesco Cordoni 1
Nicolas Doremus 2
Alessio Moneta 2

1

Department of Economics, Royal Holloway University of London, UK
2
Institute of Economics, Scuola Superiore Sant’Anna, Pisa, Italy

2023/07

January 2024
ISSN(ONLINE) 2284-0400

Identification of Vector Autoregressive Models with
Nonlinear Contemporaneous Structure
Francesco C ORDONI1

Nicolas D OR ÉMUS2

Alessio M ONETA3

February 5, 2024

Abstract
We propose a statistical identification procedure for recursive structural vector autoregressive(VAR) models that present a nonlinear dependence (at least) at the contemporaneous level. By
applying and adapting results from the literature on causal discovery with continuous additivenoise models, we show that, under certain conditions, a large class of structural VAR models
is identifiable. We spell out these specific conditions and propose a scheme for the estimationof structural impulse response functions in a nonlinear setting. We assess the performance of
this scheme in a simulation experiment. Finally, we apply it in a study on the effects of themacroeconomic shocks that propagate through the economy, allowing for asymmetry between
responses from positive and negative impulses.

Keywords: Structural VAR models, Impulse Response Functions, Causal Discovery, Nonlinearity,
Additive Noise Models.
JEL classification: C32, C52, E52.

1

Department of Economics, Royal Holloway University of London. Email: Francesco.Cordoni@rhul.ac.uk
Institute of Economics, Scuola Superiore Sant’Anna, Pisa. Email: nicolas.doremus@santannapisa.it
3
Institute of Economics, Scuola Superiore Sant’Anna, Pisa. Email: a.moneta@santannapisa.it
The authors acknowledge support from the project “How good is your model? Empirical evaluation and validation ofquantitative models in economics” funded by PRIN grant no. 20177FX2A7. We thank Mario Martinoli and Sebastiano
Michele Zema for comments on an earlier version of the paper. We are also grateful to the participants of the ItalianWorkshop of Econometrics and Empirical Economics (Rimini 2022), the Econometric Society European Meeting (Milan 2022), the Workshop “Model Evaluation and Causal Search” (Pisa 2022), and the Italian Congress of Econometrics
and Empirical Economics (Cagliari 2023). We also thank Nicoló Maffei Faccioli for providing us with the Excess Bond
Premium data series. The R and Matlab code and the data set to replicate the results reported in this article are available
at the following repository: https://github.com/f-cordoni/NLAM_time_series.git
21 Introduction
Since the seminal work of Reiersøl (1950), the econometric literature has made clear that linear
additive noise models with Gaussian and spherical errors are not identifiable unless a priori restrictions or external instruments are employed. Structural vector autoregressive (VAR) models,
widely used in empirical macroeconomic research since Sims (1980), share the same problem ofidentification. This has been addressed by relying on constraints derived from economic theory(i.e., short-run, long-run, and sign restrictions, as pioneered, respectively, by Sims, 1980, Blanchard and Quah, 1989, and Faust, 1998), by using external instruments (from Romer and Romer,
1989 to Montiel Olea et al., 2021), or, alternatively, by exploiting specific statistical properties ofthe data that are at odds with the standard model.
This paper contributes to the statistical identification approach. As regards this approach, one
finds two strands of the structural VAR literature. One of them exploits non-Gaussianity of theerrors, the other exploits heteroskedasticity. In common, they have the use of higher moments
for identification (Montiel Olea et al., 2022; Herwartz et al., 2022). This is possible thanks to theintroduction of assumptions that depart from the standard model. In this contribution, we want toshow that another swerve from the linear model with spherical-normal errors can be exploited forthe identification of structural VAR models, namely nonlinearity.
Identification by non-Gaussianity is achieved by the application of ideas and techniques developed in the literature on Independent Component Analysis (ICA), a research field that emerged at
the intersection of statistics with signal processing (Comon, 1994; Hyvärinen et al., 2001). Studies in the machine learning literature have shown that ICA can be used for causal inference, andmore specifically, for the identification of linear structural equation models, under the assumptionof noise independence and non-Gaussianity, as well as the presence of a recursive (i.e. acyclic)
structure among the variables (see, e.g., Shimizu et al., 2006, 2011; Hyvärinen, 2013). Moneta
et al. (2013) have applied some of these techniques to structural VAR analysis. The structural
VAR-ICA literature has been rapidly growing in the recent years, with contributions by Lanne andLütkepohl (2010); Lanne et al. (2017); Gouriéroux et al. (2017); Herwartz (2018); Fiorentini and
Sentana (2023), among others. It is important to note that the recursiveness assumption has beenabandoned in these studies so that the impact matrix of the shocks is identified up to the postmultiplication of a generalized permutation matrix, which means, in practice, that one is able tolabel the shocks only after involving ex post some form of economic reasoning, jointly with an inspection, for instance, of the obtained impulse response functions. Recent studies have also relaxed
the independence assumption but maintaining non-Gaussianity and use of higher moments (Lanne
and Luoto, 2021; Guay, 2021; Mesters and Zwiernik, 2022). Identification by heteroskedasticity
was pioneered by Sentana and Fiorentini (2001) and further pursued by Rigobon (2003); Lanne
and Lütkepohl (2008); Sims (2020); Brunnermeier et al. (2021); Lewis (2021), among others.1
1

Identification by heteroskedasticity is a “higher-moment” procedure in the sense of using information from conditional (instead of unconditional) second moments (Montiel Olea et al., 2021), but its underlying model is not, differently
from the ICA model, part of the general class of models we consider in this paper.

2

We claim that nonlinearity, analogously to non-Gaussianity and heteroskedasticity, can be exploited for (recursive) structural VAR identification. This idea has been much less explored in theeconometric literature and we believe that it is important to trigger a novel discussion on this. We
put forth this idea by — similarly to what has been done in the VAR-ICA approach — importing techniques from the machine learning literature on causal discovery (see the review by Peters
et al., 2017) and by presenting a class of autoregressive processes in which it is possible to recover
structural shocks (without any problem of labelling) and their dynamic effects on variables of interests directly from the data. In tune with the causal discovery literature, we frame the problemof identifying structural shocks and their contemporaneous impacts in a structural VAR model asa particular case of the general problem of learning a causal structure (which can be usefully represented by a graph) from an observational joint distribution. Hoyer et al. (2008) and Peters et al.
(2014) have shown that if the observational distribution follows a structural equation model withan additive noise structure, then the causal graph becomes identifiable from the distribution under
specific assumptions. Peters et al. (2014) have also provided practical algorithms to retrieve thecausal structure from finite samples.
It is noteworthy that, in the causal discovery framework we introduce, identification by nonGaussianity and identification by nonlinearity are nested within the general framework of causal
discovery with additive noise models. In the first case the key assumptions are, as noted above,shocks’ non-Gaussianity plus independence, while in the second case nonlinearity in the shocktransmission mechanism plus (again) independence. A related contribution of this paper is toshed light on the connection between these two related identification methods. Furthermore, we
spell out the conditions that allow identification where the underlying model is a VAR processwith nonlinear dependence at the contemporaneous level. Such a process may or may not show
nonlinearities in the autoregressive structure, but the functional autoregressive form is not key for
identification. Similarly to the first wave of applications of ICA to VAR analysis (Moneta et al.,2013; Guerini and Moneta, 2017), however, the recursiveness assumption is here required.
We also provide a practical scheme (two algorithms) to recover the structure linking shocks to(reduced-form) VAR innovations and, on the basis of this, to estimate nonlinear structural impulse
response functions from time series data. Apart from the role it plays in identification, there aremany empirical cases in which the presence of nonlinearity in the processes underlying observed
data should not be underestimated. As Ramey (2016) points out, “positive shocks might havedifferent effects from negative shocks, effects might not be proportional to the size of the shock,or the effect of a shock might depend on the state of the economy when the shock hits.”
In the next section we present the model set up and the identification methodology, which
includes both theoretical results and practical algorithms. In section 3 we present the results of asimulation study. Section 4 presents an empirical application. Section 5 concludes.

3

2 Identification methodology
In this section, we present our method to achieve statistical identification of a specific class of
VAR models, namely VAR models with generic time dependence and nonlinear but recursive
contemporaneous causal structure.
In section 2.1, we introduce this class of models, which we call generalized VAR models
with additive noise innovations. This class of VAR models is characterised by a recursive causal
structure among innovation terms, which can be represented by a directed acyclic graph (DAG).
Recovering such causal structure is key for identification, since it allows transforming the reducedform model into a structural model with independent shocks. In section 2.2, we define and discuss
the concept of identification used in this framework.
In section 2.3, we show that the contemporaneous recursive structure among the innovation
terms can be recovered by exploiting nonlinearity. In section 2.4, we present an algorithm which
infers such structure from estimated reduced-form VAR innovations. In section 2.5, we introduce
an algorithm to estimate structural impulse response functions. Theoretical results related to theconsistency of the procedure are referred to and discussed in Appendix B.

2.1

Model setup

We consider here a general class of VAR models (see, e.g., the nonlinear structural VAR models
in Kilian and Lütkepohl, 2017, ch. 18), in which a vector of K time series variables yt depends on
its lags with an additive vector of disturbances ut :yt = Ft (yt−1 , . . . , yt−L ) + ut ,(1)

where Ft (·) is a generic (possibly time dependent) nonlinear or linear function, and ut is a zeromean i.i.d. vector of innovations terms and depends only on contemporaneous structural shocks
εt,1 , . . . , εt,K :ut = G(εt )(2)
We assume that equation (2), for an appropriate ordering of the variables yt = (yt,1 , . . . , yt,K ),can be written as:(3)ut,k = ϕk (ut,k ) + εt,k , for k = 1, . . . , K
where ut,k is a subset of the variables ut,1 , . . . , ut,k−1 , and ϕk (∅) = 0. In other words, ut canbe arranged in a recursive order. We also assume that εt,1 , . . . , εt,K are cross-sectionally mutual
independent and that each εt,i (i = 1, . . . , K) is i.i.d.. Henceforth, we will refer to the model
as formalized in equations (1-3) as the generalized VAR with additive noise innovation model(GVAR-ANIM).
Equations (1) - (2) can be seen as a generalization of the standard-linear structural VAR model:

4

yt =

L
X

Aℓ yt−ℓ + ut ,

(4)

ℓ=1

where
u t = A0 ε t .

(5)

If εt,1 , . . . , εt,K are mutually independent and non-Gaussian, equation (5) is an ICA model andthe matrix A0 can be identified up to a post-multiplication of a generalized permutation matrix(Eriksson and Koivunen, 2004; Gouriéroux et al., 2017). If, for any permutation matrix P , the−1Tmatrix P A−1
0 P is lower triangular (i.e. A0 is essentially triangular), under this assumption−12
the matrix A−1
0 is uniquely identified (Shimizu et al., 2006) . We call the assumption on A0
recursiveness assumption, since it imposes a recursive causal structure on the contemporaneous
variables of the structural VAR model: if this assumption is true, one can re-order the variables
entering in yt in a “Wold causal chain” (Wold, 1960), so that each variable yt,i causes yt,j and novariable yt,j causes yt,i (i < j, for i, j in 1, . . . , K).
If A−10 is essentially triangular, a convenient way to represent the ut from equation (5) is witha directed acyclic graph (DAG) (see Spirtes et al., 2000; Pearl, 2009). Let us suppose, only for theTsake of illustration, that P A−1
0 P is lower triangular for P = I and K = 3. Then we have thissystem of structural equations between the innovation terms:u= εt,1 t,1(6)ut,2 = α ut,1 + εt,2ut,3 = β ut,1 + γ ut,2 + εt,3 ,which is represented by the DAG in Figure 1.

ut,1

ut,2

ut,3

εt,1

εt,2

εt,3

Figure 1: Example of DAG between innovation terms.
If the true DAG were known, one would be able to put (just-identifying or over-identifying)zero restrictions on A−1
0 and eventually to recover the independent shocks εt (Swanson andGranger, 1997; Demiralp and Hoover, 2003; Moneta, 2008).
Knowing the causal structure is key for identification also in the general nonlinear case of2

Shimizu et al. (2006) prove identifiability of linear, non-Gaussian and recursive model (LiNGAM in their terminology) using ICA and, specifically, the Darmois-Skitovich theorem (Comon, 1994, Theorem 11) (see also Peters et al.,2017, Theorem 7.6).

5equations (1) - (2). Let us consider the following instance of equation (2), which we have assumed
that can be arranged in a recursive order with noises entering in an additive fashion:u= εt,1 t,1(7)ut,2 = ϕ2 (ut,1 ) + εt,2ut,3 = ϕ3 (ut,1 , ut,2 ) + εt,3The structure of equation (7) can also be represented by the DAG of Figure 1. Knowing the DAG,
or having a method that infers it, allows us to recover the shocks and the manner they impact
innovations from observed realizations of ut . This can be done in a recursive way. Suppose the
DAG in Figure 1 is true. Then, following equation (7), one can (i) assign εt,1 to ut,1 ; (ii) run
an adequate (e.g., nonparametric) regression of ut,2 on ut,1 and get εt,2 as residual; (iii) regress
ut,3 on ut,1 and ut,2 and get εt,3 as residual. In the next subsections, we will present a method torecover the DAG that represents the causal structure at stake from estimates of ut in the GVARANIM framework, where nonlinearity (with additive noise) is the general case and linearity (withnon-Gaussian noise) is only a special case.

2.2

A note on identification

The concept of identification of linear structural VAR models hinges on the theory of identification in parametric models (Rothenberg, 1971). In this framework, criteria for structural VAR
identification have been provided and discussed by Rubio-Ramirez et al. (2010) and Bacchiocchi
and Kitagawa (2022). The notion of identification involved in this paper relies on a general formulation, encompassing both parametric and nonparametric models, that has been proposed byKoopmans and Reiersøl (1950) and further developed by Roehrig (1988) (see Matzkin, 2007, fora definition of identification in a purely nonparametric setting). In this general formulation, thetarget of identification is not (necessarily) a set of parameter points, but, rather, a set of “characteristics” C of a structure S = (f, Φ). A structure consists of a set f of relationships of interest(“structural relationships”) between observed variables Y and latent variables V (f (Y, V ) = 0),
and a probability distribution on the latent variables Φ(V ). A structure implies a unique probability distribution Ψ on the observed variable Y . The set of possible structures are restricteda priori by a model M, which is defined as a set of structures which share certain pre-defined
characteristics.
Definition (observational equivalence): Two structures in M are said to be observational equivalent if they imply the same probability distribution Ψ for the observable variables Y (cf. Koopmans
and Reiersøl, 1950; Rothenberg, 1971).
Definition (identification): A structure S is identified in M if and only if there is no other S ∗ ∈M that is observational equivalent to S. A set of characteristics C(S) is identified in M if andonly if every structure S ∗ ∈ M that is observational equivalent to S has the same characteristics
C(S) (cf. Roehrig, 1988, p.435).

6

Thus, in our case, GVAR-ANIM, as defined in subsection 2.1 limits the set of possible structures to restricted characteristics that will be specified in subsection 2.3. The model is nonparametric in not imposing a parameterization of the functions Ft and G (equations 1 and 2). Therefore,
the characteristics that turn out to be identifiable in GVAR-ANIM are not parameter points, butrather the recursive causal structure among the innovations, like the one displayed in equation (7).

2.3

Exploiting nonlinearity

The causal discovery principle that we exploit here was first introduced in the literature by Hoyer
et al. (2008). This can be seen as a generalization of the principle underlying causal inference inthe linear non-Gaussian case. Let us first show this case, starting from a two-variables model.
Consider a bivariate VAR model, whose vector of reduced-form residuals is ut = (ut,1 , ut,2 ).
In the following, to improve readability, we drop the time-subscript t when it is evident from thecontext and is not relevant, so that, for instance, u ≡ ut and ui ≡ ut,i . The following result wasproven by Peters et al. (2017, Theorem 4.2).
Theorem 13 Consider the following linear model:
u1 = αu2 + ε1 ,ε1 ⊥⊥ u2 .

(8)

Then there exists β ∈ R and a random variable ε2 such thatu2 = βu1 + ε2 ,

ε2 ⊥
⊥ u1

(9)

if and only if ε1 and u2 are Gaussian (where ⊥⊥ denotes statistical independence).
As a corollary of this theorem, it follows that in a bivariate linear structural VAR, as the one inequation (4) but with K = 2 and non-Gaussian shocks, if ε1 is independent of u2 , then equation(8) is the true structural model and the causal structure can be represented by the DAG u2 −→ u1 .
Shimizu et al. (2011) extend this result to a multivariate framework. But let us focus here on
another possible extension, namely the generalization of this result to the nonlinear case, under
the condition that we preserve the noise-additivity assumption. It turns out that, in the strictly
nonlinear case, Gaussianity does not preclude identification any longer. Following the terminology
by Peters et al. (2017, Definition 4.4), we introduce the additive noise model (ANM) property.
Definition (bivariate ANM): The joint distribution P (u) is said to admit a (bivariate) ANM fromui to uj if, for any measurable function ϕj and a variable εj we have:
uj = ϕj (ui ) + εj ,

εj ⊥
⊥ ui

(10)

Identifiability of a bivariate ANM is based on the following theorem proven by Hoyer et al. (2008,
Theorem 1) (see also Peters et al., 2014, 2017):
3

The proof of this theorem hinges heavily on the Darmois-Skitovich Theorem (see Peters et al., 2017, Theorem 4.3).

7

Theorem 2 Let us assume that P (u) admits an ANM from ui to uj and that εj and ui have strictly
positive densities P (εj ), and P (ui ), with ϕj , P (εj ), and P (ui ) three times differentiable. To
simplify notation, let f := ϕj , and ξ := log P (ui ) and ν := log P (εj ), skipping the argumentsuj − f (ui ), ui and ui for ν, ξ and f and their derivatives, respectively. Consider the following
condition (see Peters et al., 2014, Condition 19).
Condition (C1): the triple (ϕj , P (ui ), P (εj )) does not satisfy the following differential equation for all ui , uj with ν ′′ (uj − f (ui ))f ′ (ui ) 6= 0:
ξ ′′′ = ξ ′′ (−ν ′′′ f ′ f ′′
ν ′ ν ′′′ f ′′ f ′ ν ′ (f ′′ )2′′ ′′ ′
′ ′′′
+
)
−
2ν
f
f
+
ν
f
+
−
,
ν ′′
f′
ν ′′
f′

(11)

If condition C1 is satisfied, then a backward ANM from uj to ui , i.e. ui = ϕi (uj ) + εi withεi ⊥⊥ uj , is not admitted.
Remark. As shown by Hoyer et al. (2008), a “generic” triple (ϕj , P (ui ), P (εj )) is expectedto satisfy condition C1. Zhang and Hyvärinen (2009) actually provide an exhaustive list of fivesettings that admit an ANM in both the forward and backward direction, the most remarkable casebeing the linear Gaussian case. Notice that if a joint distribution admits a (bivariate) ANM in one
direction (say from uj to ui ), but not in the backward direction, the DAG representing the ANM
is identifiable (uj −→ ui ).
Figure 2 shows results from an illustrative simulation (see also Peters et al., 2017), in whichdata x and y are generated from a structural model y = βx + Ny , with x and Ny following auniform distribution, with x ⊥⊥ Ny . In the top left panel, we show the results of a linear regressionof y on x, whose residuals are plotted over x in the top right panel. It is easy to notice the lack ofdependence between such residuals and x. The bottom left panel show the regression with exactly
the same data, but performed in the opposite direction, namely we regress x on y. Residuals ofsuch regression are plotted over y (the covariate in the new regression) in the bottom right panel.
It is easy to note a (higher-order) statistical dependence. This asymmetry can also be detected in
the nonlinear generic case. Figure 3 shows results from another illustrative simulation in whichdata x and y are generated from a structural model y = x3 + Ny , with x and Ny following anormal distribution, with x ⊥⊥ Ny . In the top left panel, we show the results of kernel regressionof y on x, whose residuals are plotted over x in the top right panel. The lack of dependence
between such residuals and x is confirmed by a nonparametric independence test, specificallythe Hilbert-Schmidt Independence Criterion proposed by Gretton et al. (2007). The p-value forthe null hypothesis of independence between residuals and covariate is reported inside the plot.
Bottom left panel show the kernel regression with exactly the same data, but performed in the
opposite direction, namely we regress x on y. Residuals of such kernel regression are plotted
over y (the covariate in the new regression) in the bottom right panel. The p-value for the null
hypothesis of independence between residuals and covariate is reported inside the plot and clearlysuggests to reject independence.
Peters et al. (2014) proved that the identifiability result stated in Theorem 2 can be extended
from the bivariate to the multivariate case. First of all, the definition of ANM given above can be
8

−0.2

0.0

0.2

0.0 0.2 0.4−0.4

residuals of regr. y on x

1.0
0.5
0.0
−1.0

−0.5

y

−0.4

0.4

−0.4

−0.2

0.4

0.2

0.4

1.0
0.5
−1.0

−0.5

y

0.0

0.5
0.0
−1.0

−0.5

y

0.2

x

1.0

x

0.0

−0.4

−0.2

0.0

0.2

0.4

−0.4

x

−0.2

0.0

residuals of regr. x on y

Figure 2: Illustration of the ANM principle with a sample generated from the linear modely = βx + Ny with x and Ny independent and drawn from a uniform distribution. Left-side
panels: scatter plot of y on x with OLS regression line (red) corresponding to the ‘forward’ (true)
specification (top panel) and to the ‘backward’ specification (bottom panel). Right-side panels:
corresponding regression residuals are plotted with respect to their regressor values, either x (toppanel) or y (bottom panel).
straightforwardly extended.
Definition (multivariate ANM). We call a system of K structural recursive equations a (multivariate) ANM if it can be written as:
uk = ϕk (P a(uk )) + εk

for k in 1, . . . , K,

(12)

where ε1 , . . . , εK are mutually independent. The set of variables P a(uk ), called parents of uk , isdefined as P a(uk ) ⊆ {u1 , . . . , uK }\{uk }. We denote by G the DAG representing the structural
relations between the u’s, i.e. ui −→ uj iff ui ∈ P a(uj ).
Notice that the causal structure G is a DAG because the system is recursive. Identifiability of
9

2
1
0
−1
−3

−2

residuals of regr. y on x

15
10
5
−10

0

y

dhsic.test − p.value= 0.77

−2

−1

0

1

2

−2

−1

1

2

5

10

dhsic.test − p.value= 0.006

−10

0

y
−10

0

y

5

10

15

x

15

x

0

−2

−1

0

1

2

−1.5 −1.0 −0.5

x

0.0

0.5

1.0

residuals of regr. x on y

Figure 3: Illustration of the ANM principle with a sample generated from the nonlinear model y =x3 +Ny with x and Ny independent and drawn from a normal distribution. Left-side panels: scatterplot of y on x with kernel regression line (red) corresponding to the ‘forward’ (true) specification(top panel) and to the ‘backward’ specification (bottom panel). Right-side panels: corresponding
regression residuals are plotted with respect to their regressor values, either x (top panel) or y(bottom panel).
the multivariate ANM is based on the following theorem:
Theorem 3. Consider a multivariate ANM
uk = ϕk (P a(uk )) + εk

for k in 1, . . . , K,(13)

associated to DAG G and entailing joint distribution P (u), with ϕk three times differentiable andP (εk ) strictly positive for all k. An alternative ANM with DAG G ′ , where G ′ 6= G, also entailing
P (u), is not admitted under the following conditions:
Condition (C2) The functions ϕ1 (·), . . . , ϕK (·) are not constant in any of their arguments.
Condition (C3) Let us denote with ND(uj ) the set of graphical “non-descendant” of uj , i.e.
10

there is no directed path from uj to any variable in ND(uj ) (see Spirtes et al., 2000). For all
uj ∈ u, and for all ui ∈ P a(uj ), and all uS ∈ u such that P a(uj ) \ ui ⊆ uS ⊆ ND(uj ) \ {ui , uj },there is a value uS of uS with p(uS ) > 0 such that the triple
(14)
ϕj (P a(uj ) \ {ui }, ui ), p(ui |uS ), p(εj )satisfies condition C1. In equation (14) the upper bar indicates a specific value that a randomvariable or vector takes.
The conditions related to Theorem 3 and its proof can be found in (Peters et al., 2014, seesection 3.2 and appendix A therein). Thorem 3 allows us to recover the contemporaneous causal
structure among the innovation terms in the GVAR-ANIM. In the next subsection, we will present
an algorithm to learn such causal structure from estimated reduced-form VAR residuals.
Topological orderAs shown above, knowing the causal structure among innovation terms in the form of a DAG allows specifying the structural model, since the set of graphical parents determines the argumentsof ϕk (·) in equation (12). Structural shocks can be recovered by regression (parametric or nonparametric) methods applied to the same equation.
Let us now introduce the notion of topological order. Given a DAG over vertices v1 , . . . , vK ,it is possible to associate a total order on its vertices such that if vi is a parent of vj (i.e. there
is a directed edge from vi to vj ), then vi comes before vj in the ordering. Notice that there is
only one topological order associated to a fully connected DAG (a DAG in which each pair ofvertices is connected). But, in the case of a sparse DAG, there are possibly multiple topologicalorders associated to it. Also notice that removing an edge from a fully connected DAG does notinvalidate the original topological order, but there could be alternative ordering compatible with
it.4
Suppose that a set of structural recursive equations, associated to a graph G0 , has generated
some data, and suppose that there are multiple topological orders associated to G0 . One can deducea unique fully connected DAG, which we call G max , for each of these topological orders. It may
be the case G0 6= G max . But G max differs from G0 only in displaying extra edges, i.e. G0 is a
subgraph of G max . Therefore, if we write down a set of structural equations (as in equation 12)
associated to G max , this set of equations will differ from the “true” one, i.e. the one associatedto G0 , only in that it includes extra variables in right-hand side. These extra-variables are notin a parent set in (the structural model associated to) G0 . But since the original causal ordering(associated to G0 ) is maintained, no endogeneity issues arise when one uses G max as a template fora regression model. Indeed, if one applies regression methods to the set of equations associated to
G max , which may therefore include extra variables (with respect to the “true” structural model) on
4

Consider the fully connected DAG on v1 , v2 , v3 associated to the topological order hv1 , v2 , v3 i. From removing
v1 −→ v3 , we get v1 −→ v2 −→ v3 , and the original topological order hv1 , v2 , v3 i is still the only one valid. Fromremoving v2 −→ v3 , we get v2 ←− v1 −→ v3 , and there are two possible topological orders: hv1 , v2 , v3 i, hv1 , v3 , v2 i.

11the right hand side, one obtains consistent estimates, since these variables are at most irrelevant,but do not entail reverse causality problems.
In practice, if the aim of the study is to estimate the structural shocks (e.g., in impulse response
function analysis), one may aim at recovering a topological order instead of a DAG, so that to minimize inference errors in finite samples. Instead, if causal discovery is the goal, then a topologicalorder is not enough. This is why in the next section we will present a search algorithm that can be
stopped at the topological order step.
A simple case
While the paper, including the algorithms presented in the next subsections, covers a wide classof models that we have labelled GVAR-ANIM (see equations 1-3), in both the simulation andthe empirical analysis, we are going to focus on a simpler case, namely the case of linear lagged
relationships with nonlinear contemporaneous relationships between innovations (with additive
structural shocks). In terms of equations:

P
y
= L
t
ℓ=1 Aℓ yt−ℓ + ut ,(15)
u
= ϕ (P a(u )) + ε
for k in 1, . . . , K.
t,k

k

t,k

t,k

We decide to focus on this simple case because the crucial issue in VAR identification (see, e.g.
Stock and Watson, 2001) is the solution of the contemporaneous causality problem, which allows
to estimate the structural shocks. This holds both in the linear and nonlinear VAR case. Since
the method we propose exploits contemporaneous nonlinearity, equation (15) represents the most
simple time series model in which we can apply it.

2.4

Algorithm for the contemporaneous structure

Previous sections showed that, under the assumption of additive noise, the (graphical) causal
structure of a recursive structural equation model is generically identifiable. There are alternative strategies that exploit this identifying principle using sample data (see Peters et al., 2014 andPeters et al., 2017 for an overview). Here we are presenting the RESIT (REgression with Subsequent Independence Test) algorithm proposed in Peters et al. (2014), which incorporates theprinciples of identification in the nonlinear additive noise setting directly and straightforwardly. It
also has the advantage of segmenting the causal search in, first, a topological order search phase
and, subsequently, a DAG search phase. This modular architecture fits well with the possibility
that for shock recovery (and impulse response analysis) a topological order can be sufficient forthe scope (in finite samples), as noted above. The original RESIT algorithm of Peters et al., 2014
operates in two phases (see Appendix A). We add a preliminary phase, which we call Phase 0, in
which we estimate a reduced-form GVAR-ANIM (see equation 1) and use its estimated residuals
as input for the next phase.

12Phase 1 aims at ranking the input variables in a topological order. The topological order iselicited in a backward fashion by determining which variable is the latest (i.e. a sink node) among
the whole set of variables to be ordered, removing it from the set and repeating the search over the
resulting subset until only one element remains, which is given the first position in the ranking.
At each iteration, the search is performed by regressing each variable jointly on all the others andmeasuring the dependence (more on this below) between the regressors and the obtained residual.
The regressand yielding the weakest dependence is regarded as the sink node among the current
set of candidates; it is removed from the set of candidates and is assigned the first position among
the set of variables that have already been removed. The topological order simply corresponds
to the reversed order of elimination. The topological order determines, for each variable, a set of
potential parents. This set coincides with all the variables ranked higher.
Phase 2 (which we call pruning phase) further narrows down the sets of potential parents
by performing variable selection (see Appendix A). In other words, it prunes edges from thefully connected DAG which one can draw from the topological order of Phase 1. It essentiallyperforms a variable selection task. It regresses each variable on its putative parents (according
to the topological order derived above) but omitting iteratively one of them. If the estimatedresidual of the corresponding nested model is found independent of all preceding variables in the
topological order, the omitted regressor is removed from the set of parents of the regressand. The
search stops when no regressor can be removed without creating dependence.
Any regression method and dependence measure can be embedded in the RESIT algorithm, aslong as they are adequate with respect to the data-generating process. Not knowing the functional
dependence, i.e. ϕk (·) in equation (12), an appropriate strategy would be to use a nonparametric kernel estimator (Fan and Gijbels, 2018). In our implementation, in tune with Hoyer et al.
(2008), we opt for the Gaussian process approach (see Williams and Rasmussen, 2006). An alternative, more restricted, nonparametric model is the generalized additive noise model (Hastie andTibshirani, 2017). At the extreme, the functional form can be “restricted” down to linearity, as
long as one assumes non-Gaussianity of the noises. In this case, the ANM collapses to the Linear
Non-Gaussian Acyclic Model studied by Shimizu et al. (2006) (see also Shimizu et al., 2011).
As regards the dependence measure, any measure that accounts for higher order statistics(rather than linear dependence only) can in principle work. This feature is essential because regression techniques that minimize quadratic errors including linear components yield estimatednoises that are orthogonal to the covariates. Peters et al. (2014) propose to use the p-value of a
nonparametric independence test as dependence score, so that the dependence is minimized when
the p-value, under the null hypothesis of independence, is maximized. This follows the principleof Hodges-Lehmann estimation, which was proposed by Herwartz (2018) for independent component analysis. The specific independence test that we perform in our implementation, following(Peters et al., 2014), is the kernel independence test proposed by Gretton et al. (2007) with theHilbert-Schmidt Independence Criterion (HSIC). The asymptotic results presented by these authors guarantee a convergence of the p-value to zero only under dependence.
We mentioned above that the RESIT algorithm can be stopped at the topological order step13

(i.e. Phase 1) if the aim is to estimate structural shocks and impulse response functions. The
second phase should instead be applied if one is interested in uncovering the exact causal structure
among the innovation terms. Since both in our simulation study and in the empirical application
we focus on structural impulse response analysis, the algorithm shown in the main body of thetext skips the edge pruning phase (see Algorithm 1 here below). In Appendix A, one can find theoriginal RESIT algorithm embedded in a more general procedure that deals with time series data.
Notice that the motivation to skip Phase 2 of RESIT is only of practical nature, when one deals
with finite sample. Contrary to Phase 1, Phase 2 operates on the basis of statistical hypothesis
testing (independence test), in which one has to decide a rejection threshold. In finite sample,
this decision may lead to mistakenly classify weak dependence as independence and therefore to
mistakenly cutting edges between directly causally related variables. For the sake of retrieving
structural shocks from reduced-form VAR residuals, one has to run a battery of regressions of theform expressed in equation (13), in which shocks will be the residuals of these regressions. If weomit any variable from the parents’ set (i.e. covariates of these regressions) because we have madea mistake in Phase 2, we introduce an omitted variable bias. But if we keep an irrelevant variablein the parents’ set, we may add some variance in the estimates but we do not introduce any bias.
Given this trade-off, we opt for not using Phase 2 for the sake of estimating impulse response
functions.
Algorithm 1 VAR + RESIT (Phase 1 is from Peters et al., 2014)
1: Input: A K-dimensional time series vector (yt,1 , . . . , yt,K )′2: PHASE 0: Estimate the reduced-form model.
3: Estimate a reduced-form time series model of the class GVAR-ANIM, see eq. (1) and extract

residuals u
bt = (but,1 , . . . , ubt,K )′ . To simplify the notation, let us call u
bt,k ≡ vk , for k =1, . . . , K.
4: S := 1, ..., K, π := [ ]5: PHASE 1: Determine topological order and potential parent set
6: repeatfor k ∈ S doRegress vk on {vi }i∈S\{k} and obtain residuals ek
9:Measure dependence between ek and {vi }i∈S\{k}10:
end for11:
Let k ∗ be the k with the weakest dependence
12:
S := S \ k ∗13:pa(k ∗ ) := S
14:
π := [k ∗ , π]
15: until #S = 0
16: Output: π
17: Output: (pa(1), . . . , pa(K))
7:8:Note: with pa(k) we refer to the set of indices associated to the variables in P a(vk ), for any k.

14

A simple illustration
Let us suppose that a 3-variable GVAR-ANIM process of the type of equations (1-3) has generated
some data. Specifically, suppose that equation (3) corresponds to:u = εt,1

 t,1
(16)
ut,2 = ϕ2 (ut,1 , ut,3 ) + εt,2ut,3 = εt,3
The DAG associated to these equations is ut,1 −→ ut,2 ←− ut,3 . Now, consider to have estimatedut,1 , ut,2 , ut,3 as reduced-form residuals of a VAR of the form of equation (1), following Phase 0 of
Algorithm 1. A mistake-free output of Phase 1 will be either the topological order hut,1 , ut,3 , ut,2 i
or hut,3 , ut,1 , ut,2 i. A mistake-free output of Phase 2 (see Appendix A) will be ut,1 −→ ut,2 ←−ut,3 . In order to estimate the structural shocks ε1 , ε2 , ε3 one can run a set of regressions on thebasis of equations (16). But suppose one stops at Phase 1. Then one would estimate one of thefollowing regression models instead:u = εt,1
 t,1
ut,2 = ϕ2 (ut,1 , ut,3 ) + εt,2ut,3 = ϕ3 (ut,1 ) + εt,3

(17)



u = ϕ1 (ut,3 ) + εt,1
 t,1ut,2 = ϕ2 (ut,1 , ut,3 ) + εt,2ut,3 = εt,3(18)

Both systems of equations (17) and (18) differ from (16) in including irrelevant variables in
the regression systems, but not involving any reverse causality issue.

2.5

Algorithm for the structural impulse response functions

The advantage of identifying a structural model is that one can estimate the dynamic effects ofstructural shocks. But, under the generic GVAR-ANIM framework, this cannot be done by using the conventional tools of impulse response analysis. Indeed in such a framework the variables’
responses to a specific shock at a certain time depend on the system’s past history, sign and magnitude of all contemporaneous shocks and subsequent ones. Following Koop et al. (1996) and Kilianand Lütkepohl (2017, ch.18), we define the structural (nonlinear) impulse response function (IRF)as:

IRF(h, δ, Ωt−1 ) = E(yt+h |εt,i = δ, Ωt−1 ) − E(yt+h |Ωt−1 ),(19)

where δ is the (positive or negative) magnitude of the shock εt,i , whose effects one wants to study,Ωt−1 is the history of the model data up to time t − 1, and h is the horizon point up to which
the impulse response function is studied. To estimate equation (19), we can use the Monte Carlointegration approach suggested by Kilian and Lütkepohl (2017).

15

The idea of this approach, whose detailed scheme is presented in Algorithm 2, is to take pairsof simulation runs: Ỹ δ,n and Ỹ n (for n = 1, . . . , N ), in which Ỹ δ,n denotes the “treated” timeseries and Ỹ n the “control” one, where N is the number of Monte Carlo runs. Both elements
share a common history Ωt−1 , but diverge at date t, in which a structural shock of interest is set
to δ in one run (treatment), while no such restriction applies to the other (control). Apart fromthis restriction, all shocks are free to fluctuate in both groups of series from time t onwards. The
difference is then taken for each pair. Finally, estimated impulse response functions are calculated
by taking average of these simulated differences.
Algorithm 2 requires as input estimates (from the observed data) of the functional relationshipϕbk (·) between innovations in ut,k = ϕk (P a(ut,k )) +εt,k (for k in 1, . . . , K); estimates of thelagged functional relationship fbk (·) in yk,t = fk (yt−1 , . . . , yt−p )+ uk,t (for k in 1, . . . , K);5
estimates of the structural shocks εt,1 , . . . , εt,K . These can be estimated through the following
procedure. For each k = 1, . . . , K, if P a(ut,k ) = ∅ then set εbt,k := ubt,k ; else (non-parametrically)
regress u
bt,k on {but,i }i∈pa(k) and set εbt,k := ubt,k − ϕbk ({but,i }i∈pa(k) ).6

Each simulated7 path (Ỹ δ,n and Ỹ n ) is computed iteratively at each time period, starting fromthe time point of the “treatment”. At each step, new values are generated summing (i) lagged
effects fk (yt−1 , . . . , yt−L ); with (ii) innovation terms ut,k = ϕk (P a(ut,k )) + εt,k , for every k.
New values of εt,k are independently drawn from their marginal empirical distribution at eachtime horizon (except for the one set at δ at time t in the treatment run). Since the innovation
model is recursive, each ũt,k is computed iteratively at each time step by adding its structural
shock εt,k to the contemporaneous effect of its parents (P a(ut,k )). Parents’ sets can be assignedin two alternative ways: (i) by simply taking the parents’ sets as given by the full version of
RESIT as it appears in Appendix A; or (ii) in a more conservative fashion (see discussion above),by considering the sole topological order (see Algorithm 1 in section 2.4) and by regarding any
variable ut,k as a parent of those it precedes in that ordering. Either way, the functional forms ϕbk (·)of the contemporaneous effects are estimated via (nonparametric) regression of each variable overthe parent set considered.

3 Simulation analysis
The simulation study documented in this section aims at evaluating the performance of our approach in estimating structural impulse response functions. We simulate a simple VAR with different contemporaneous causal structures and different parametrizations. We then estimate thestructural IRFs using Algorithm 2, but under different identification approaches, and, we systematically compare them. Our simulation results, as reported below, demonstrate that, in case of
nonlinearity in the data generating process, our procedure is able to satisfactorily estimate the5

For simplicity, we dropped the time index from the lagged effects fk,t (·). A time-varying function can be incorporated into the algorithm we are discussing in this subsection by updating the fk,t (·) at each step.
6
With pa(k) we refer to the set of indices associated to the variables in P a(ut,k ), for any k.
7
Simulated variables are marked with a tilde sign (∼).

16

Algorithm 2 Structural IRFs at time t∗
1: Setting Values:• t∗ : date of the shock to be studied• H: max horizon• k ∗ : index of the shocked variable• δ: shock magnitude• N : number of Monte Carlo runs
2: Input data

• history (Ωt∗ −1 ) at date t∗b• estimated structural shock matrix E• estimated fbk (·), ϕbk (·) (for each k)• topological order: π = (π1 , . . . , πK ); parents set: (pa(1), . . . , pa(K))

3: PHASE 1: Simulate N couples of time paths and get the difference
4: for n ∈ {1, . . . , N } do
5:
6:
7:

Create a (t∗ + H) × K empty matrix Ỹ n . Set 1 : (t∗ − 1) rows := Ωt∗ −1Ỹ δ,n := Ỹ nb (columnwise)
Create a (H + 1) × K matrix Ẽ, randomly sampling from E
Note: denote any (t, k) entry of any matrix E with Et,k

8:
9:
10:
11:
12:
13:
14:

Create a (H + 1) × K matrix Ẽ δ , as in line 7, but set Ẽtδ∗ ,k∗ := δCreate a (H + 1) × K matrices Ũ and Ũ δ (initially empty)∗
Create a (H + 1) × K matrices I˜n;k ,δ (initially empty)
for h ∈ {0, . . . , H} dofor k ∈ {π1 . . . πK } do
δ := ϕδ }δ
Ũh,kck ({Ũh,ii∈pa(k) ) + Ẽh,k (with ϕk (∅) = 0)
) + Ũ δỸ δ,n:= fbk (Ỹ δ,n, . . . , Ỹ δ,n∗∗∗t +h,kt +h−1t +h−Lh,k

Ũh,k := ϕck ({Ũh,i }i∈pa(k) ) + Ẽh,kn
16:Ỹt∗ +h,k := fbk (Ỹtn∗ +h−1 , . . . , Ỹtn∗ +h−L ) + Ũh,kn,k∗ ,δn
17:
I˜h,k:= Ỹtδ,n∗ +h,k − Ỹt∗ +h,k18:
end for19:
end for
20: end for
21: PHASE 2: Average values∗∗
1 PN
˜n,k∗ ,δ,t∗22: I˜k ,δ,t := N
n=1 I
∗
∗
23: Output: I˜k ,δ,t15:

17

structural IRFs and that a linear approximation may lead astray even under correctness of therecursive order.

3.1

Data generating processes

We simulate a 3-variables VAR model with one lag, which can be written as
yt = A1 yt−1 + ut ,(20)

where A1 , to control for the persistence of the process, is a lower triangular matrix with all elements below and on the main diagonal equal to 0.5. The ut terms are generated following three
types of causal structures, which we call the “causal chain” (21), “common cause” (22), and “vstructure” (23):Causal chain
u = εt,1
 t,1
ut,2 = ϕ2 (ut,1 ) + εt,2ut,3 = ϕ3 (ut,2 ) + εt,3(21)Common causeu = εt,1
 t,1
ut,2 = ϕ2 (ut,1 ) + εt,2ut,3 = ϕ3 (ut,1 ) + εt,3(22)v-structureu = εt,1
 t,1
ut,2 = εt,2ut,3 = ϕ3 (ut,1 , ut,2 ) + εt,3(23)Furthermore, each generic model will be given three alternative parametrizations to (ϕ• (·), ε• ).
We first define a “linear Gaussian” setting in which each ϕ• (x) = x and εt,1 , εt,2 , εt,3 are mutually independent and (standard) normally distributed. Second, the “linear non-Gaussian” setting
is characterized by the same identity function but the εt,1 , εt,2 , εt,3 are drawn from a Laplace distribution with scales parameters 1, 2, 4, respectively. Finally, we define a “nonlinear Gaussian”
setting in which shocks are drawn in the same fashion as in the “linear Gaussian” setting and thefunctional forms taken by ϕ2 (·) and ϕ3 (·) in each causal structure are reported in Table 1. We
generate samples {yt }Tt=1 of size T = 250, 500, 1000.
Table 1: Functional forms of nonlinear models. The parameters α and β are independently drawn
from a uniform distribution with support [1, 4].
ϕ2 ()
ϕ3 ()

3.2

Causal chain
sign(ut,1 )|ut,1 |αsin(sign(ut,2 )|ut,2 |β )Common causesign(ut,1 )|ut,1 |αsin(sign(ut,2 )|ut,2 |β )v-structuresign(ut,1 )|ut,1 |αsign(ut,1 )|ut,1 |α + sin(sign(ut,2 )|ut,2 |β )Identification schemes

For each selected T , type of structure and parametrization we generate 200 artificial (3-variate)
time series, from which we estimate IRFs. These are constructed following Algorithm 2, according
to five different schemes, which we call: (i) “theoretical” IRFs, (ii) Cholesky-based (linear) IRFs,(iii) true-graph-based nonlinear IRFs, (iv) true-topological-order-based nonlinear IRFs, and (v)Algorithm-1-based nonlinear IRFs.
18

Scheme (i) builds IRFs by applying Algorithm 2 to each artificial time series, except thatinstead of the estimated structural shock matrix (see “Input data” in the pseudo code) we use thetrue DGP’s one, instead of the estimated functional forms we use the autoregressive matrix A1 ,and the true generating functions ϕk (·), and we use the correct set of parents (and topologicalorder) implied by the true causal structure. We calculate IRFs using this scheme because we wantto have a ground truth against which we can assess our procedure. In practice, scheme (i) is theapplication of Algorithm 2 to the true model, bypassing the model estimation step. For this reason,
we call the produced IRFs “theoretical”.
In contrast, the remaining four procedures are based on estimates of the model. Scheme (ii)is meant to provide a benchmark to the approaches that, differently from this scheme, allow for
nonlinearity. Here structural shocks are computed via Cholesky factorization of the reduced-form
residuals’ covariance matrix, where the recursive ordering of the innovation structure is the true
one, instead of being determined in a data driven fashion as Algorithm 1 does. Moreover, for thesake of comparison with the other schemes, we do not apply here the standard linear approach thatderives IRFs from MA coefficients, but we feed Algorithm 2 with the estimated linear coefficients.
We label IRFs from this scheme “CHOL” (from Cholesky).
Scheme (iii) is similar to scheme (i) except that the functional (autoregressive and contemporaneous) dependencies are estimated from the data, on the basis of the true causal graph (which
automatically implies a correct topological order). Functional forms are estimated via OLS as
regards the autoregressive components and through Gaussian process regression as regards thecontemporaneous structure. In practice, in this scheme we are bypassing the data-driven causal(order) search by using the correct specification of the model. The aim here is to isolate the performance of Algorithm 2, avoiding potentially compounding mistakes of Algorithm 1. We label
IRFs from this scheme “TDAG” (i.e. true DAG).
Scheme (iv) is a slight variant of scheme (iii): we feed Algorithm 2 with a correct topologicalorder and a fully connected DAG deduced from it, which may have extra edges with respect to thetrue DAG. Thus, we also isolate here the performance of Algorithm 2 (as in (iii)), but we allow
the introduction of potentially irrelevant variables in its inputs, to see if they make a differencein terms of performance. Notice that Algorithm 1, which is involved in the next scheme, mayalso introduce these extra variables, so that scheme (iv) is a closer benchmark to scheme (v).
Furthermore, scheme (iv) is a nonlinear counterpart of the Cholesky scheme (ii) because they areboth based on the true topological order, without any further restrictions (i.e., they may include
irrelevant covariates). We label IRFs from this scheme “TTOP” (i.e. true topological order).
Scheme (v) implements the data-driven approach proposed in Section 2 so that its assessment
is the focus of this analysis, the previous schemes ultimately being benchmarks for this one. Here,from each generated time-series sample we estimate the topological order, from which we deducea fully connected DAG and a corresponding set of parents, by Algorithm 1 and subsequently weapply Algorithm 2. As input for the functional forms we use OLS for the autoregressive part (asin schemes (ii-iii-iv)) and Gaussian process regression for the contemporaneous structure (as inschemes (iii-iv)). Recall from Section 2.5 that the structural shocks (further input of Algorithm 2)
19are recovered from the residuals of the regression model at the contemporaneous level. We label
IRFs from this scheme “ALG1” (i.e. based on Algorithm 1).

3.3

Simulation results

Let us call IRF•k,h the IRF referring to the responses of variable yt,k to shock εt,1 , at horizon h,using scheme •, where k ∈ {2, 3}, h ∈ {0, 1, 2, 3, 4}, • ∈ {CHOL, TDAG, TTOP, ALG1}. For agiven data generating process, we calculate the following average mean squared error:
4

200

X XX
•
2
d •k,h − IRF
d theor.
\ = 1 1(IRF
AMSEk,h ) ,
10 200(24)

k=2,3 h=0 l=1

where l indicates a simulation run. We focus on the effect of the first shock to other variables since
the first variable is always exogenous in all the causal structures we consider.
Table 2 shows estimated AMSE, averaged across 200 simulation experiments, for different
values of T and parametrizations, along with standard errors, for the causal-chain structure. Analogous results for the common-cause (Table 3) and v-structure (Table 4) are reported in Appendix
C since they are qualitatively similar.
Looking at the top panels (of all these tables), we note that the ALG1 scheme (v) outperforms
the Cholesky one (scheme ii) whatever the structure and sample size when the DGP is nonlinear(second row). Recall that scheme (v) has the advantage, over scheme (ii), of allowing for nonlinearity in the estimation of impulse response functions (since it uses Gaussian process regressions
via Algorithm 2), when the data are generated by a nonlinear DGP, as in this case. But it undertakes the risky task of inferring the causal order from data.8 On the other hand, scheme (ii) runsin this case under a misspecified model (the linear one). But it has the advantage of overriding thesearch for a topological order from the data: it is given a correct one.
We observe a reversal in the relative accuracy of scheme (ii) and (v) when the model becomeslinear (see first and third row in Table 2, subtables CHOL and ALG1). The poor performance ofscheme (v) within the linear Gaussian setting does not come as a surprise since this class of ANMs
is not identified. In contrast, the linear non-Gaussian model is identified, but scheme (ii) has theadvantage of being endowed with the true topological order. Scheme (v) has the disadvantage
of learning the causal order from data and having an estimation method which is not the mostefficient for the linear case. Scheme (iv) (TTOP) fills the gap between scheme (v) and scheme(ii) in being a nonlinear approach — like the former — and hacking the order search — like thelatter. As expected, scheme (iv) systematically beats scheme (v) in terms of AMSE. To assess in
8Peters et al. (2014) report simulation results about the accuracy of the RESIT algorithm, showing that the performance (in a nonlinear setting) in terms of distance between the estimated structure and the correct DAG is good both inoutperforming other established causal discovery algorithms, but also in absolute terms. We have replicated this analysis by applying RESIT to VAR residuals and we found similar results. We also studied the accuracy of Algorithm 1 inreturning a correct topological order. Both these results are available upon request. The accuracy rate is satisfactory inmost cases, but we also found that small samples (i.e. T < 500) can be problematic for some causal structure (specifically for the common cause). Note that, however, this problem seems to have only a limited impact on the estimates ofthe structural IRFs (which is the main object of this study), as shown in the subsequent analysis reported in Figure 4.

20

Table 2: Average MSE, across 200 simulations, between IRFs estimated using different schemes(i.e. CHOL, TDAG, TTOP, ALG1) and the theoretical IRFs. MSE are calculated over the responses of yt,2 and yt,3 to the first shock (of one standard-deviation magnitude) and over the first
five periods. The contemporaneous structure is the “causal chain”. We consider different values of T (sample size) and three different model classes (Linear with Gaussian noises, non-linearwith Gaussian noises and linear with non-Gaussian noises). Standard errors are reported in roundbrackets.
CHOL
T

250

500

ALG1
1000

250

500

1000

LinearGaussian

0.034

0.020

0.015

0.257

0.244

0.249

(1.42E-03)(7.51E-04)(5.20E-04)(8.26E-03)(7.81E-03)(7.97E-03)

Nonlinear
Gaussian

1.428

1.023

1.303

0.188

0.124

0.085

(1.05E-01)(7.10E-02)(8.90E-02)(1.20E-02)(1.00E-02)(5.00E-03)

0.078

0.051

0.038

0.293

0.179

0.094

(3.53E-03)(1.92E-03)(1.82E-03)(1.36E-02)(8.11E-03)(4.88E-03)

Linear
Non-Gaussian

TDAGTLinearGaussian

TTOP

250

500

1000

250

500

1000

0.029

0.023

0.018

0.036

0.022

0.019

(1.12E-03)(8.51E-04)(6.44E-04)(1.30E-03)(8.35E-04)(7.18E-04)

Nonlinear
Gaussian

0.142

0.086

0.091

0.135

0.086

0.081

(8.00E-03)(9.00E-03)(8.00E-03)(7.00E-03)(7.00E-03)(8.00E-03)

Linear
Non-Gaussian

0.152

0.085

0.060

0.230

0.123

0.079

(6.91E-03)(3.30E-03)(2.70E-03)(1.35E-02)(6.67E-03)(4.09E-03)

isolation the effect of the estimation method, one can compare scheme (ii) with scheme (iv) (cf.
subtables CHOL and TTOP in Table 2). To assess in isolation the effect of causal learning, one
can compare scheme (v) with (iv) (cf. subtables ALG1 and TTOP in Table 2). One can noticethat the main source of error of scheme (v) depends on the distribution of the shocks. In the linear
non-Gaussian case, the performance of scheme (v) is mostly hindered by the estimation method(AMSEs of TTOP are much closer to those of ALG1 than of CHOL). Notice that, in this setting,the estimation method improves as the sample size increases, so that the performance of scheme(iii) and (iv) gets closer to scheme (ii). In the linear Gaussian case, the performance of scheme (v)is mostly hindered by the structure learning method (AMSEs of TTOP are much closer to those of
CHOL than of ALG1).
Comparisons between scheme (iv) and scheme (iii) (subtables TTOP and TDAG) do not
clearly discriminate between the two. It seems that within our sample range the less parsimonious
modeling strategy (scheme iv) often allows to better mimic the DGP than what can be achieved by21

faithfully sticking to the actual DGP specification (scheme iii). This finding is important because itsupports our decision to use Algorithm 1, i.e. a truncated version of the original RESIT by Peters
et al. (2014). As regards of how AMSE changes with respect to the sample size T , we notice thatwith scheme (ii) AMSE does not show any downward trend over T when the DGP is nonlinear,while with scheme (v) AMSE shows a clear tendency to decrease, at least over the values of T we
consider.
Figure 4 shows structural IRFs in the causal chain setting computed via Monte Carlo integration (Algorithm 2) under the theoretical scheme (i), the CHOL scheme (ii), and the ALG1 scheme(v), plotted from horizon 0 to 10. Dashed lines give the 68% confidence intervals of those averages.
Specifically, the top panels show the dynamic responses of variable 2 (yt,2 ) to setting structural
shock εt,1 to one at time t∗ = 2, while the bottom panels display their counterparts for variable 3(yt,3 ). For the sake of conciseness and readability, we do not plot IRFs obtained with TDAG and
TTOP since they are very similar to each other and tend to overlap with those obtained by ALG1.
For the same reasons, we display simulation results only for the smallest samples (T = 250, left
panels) and largest (T = 1000, right panels). Analogous figures for the common-cause and vstructure settings can be found in Appendix C (Figures 8 and 9, respectively). From these plots,
one can notice that the ALG1-IRFs are in most cases, even at horizon equal to zero, indistinguishable from the theoretical ones. Moreover, the failure of the CHOL scheme in recoveringcontemporaneous shocks’ impacts is quite evident. Since CHOL is built on the correct (imposed apriori) topological order, this suggests that knowing this order is not enough to accurately estimate
structural IRFs, but it is also crucial to allow for nonlinearities in the regression method, if present.
The nonparametric regression we use (GP regression) seems to successfully capture the nonlinear
contemporaneous effects.

4 Empirical application
Our empirical application studies the effects of macroeconomic shocks on nominal interest rate,inflation, output and financial conditions, using U.S. data. The study allows for nonlinearity at thecontemporaneous level. We compare our findings with the results from a linear recursive SVAR
model, highlighting the difference between the two results. This fact allows us to depart from atypical feature of linear SVAR models, namely the fact that the effect of a positive shock is byconstruction symmetric to the effect of a negative shock. Asymmetric effects of shocks have beenstudied, among others, by Lo and Piger (2005), Höppner et al. (2008), Kilian (2014), Hussain andMalik (2016). We should notice, however, that we are imposing a recursive structure with additivenoises on the relationship among innovations, in a setting where the reduced form VAR is linear.
This implies that the shock which is last in the topological order can only have a linear impact on
the variables. Specifically, this means that if the monetary policy instrument is last, as it turns out
to be the case in our application, its shocks have linear effects by construction.
We estimate a linear reduced-form VAR model with five variables: inflation (pit ), output gap(outt ), the federal funds rate (rt ), industrial production (ipt ), and excess bond premium (ebpt ).
22

IRF from shock 1 to variable 2 (T =1000)
THEOR
CHOL
ALG1

0

0

1

1

2

2

3

3

THEOR
CHOL
ALG1

4

IRF from shock 1 to variable 2 (T =250)

0

2

4

6

8

0

2

4

6

8

IRF from shock 1 to variable 3 (T =1000)
3.0

IRF from shock 1 to variable 3 (T =250)

2.5

THEOR
CHOL
ALG1

0.0

0.0

0.5

0.5

1.0

1.0

1.5

1.5

2.0

2.0

2.5

THEOR
CHOL
ALG1

0

2

4

6

8

0

2

4

6

8

Figure 4: DGP with “causal chain”. Theoretical (THEOR), linear (CHOL), and nonlinear (ALG1)
IRFs for T = 250 (left panels) and T = 1000 (right panels). Each solid line exhibits the average,
across 200 simulations, IRF from a unitary shock 1 (εt,1 ) to variable 2 (yt+h,2 ) (top panels) orvariable 3 (yt+h,3 ) (bottom panels), for h = 0, . . . , 9. Confidence interval at 68% are reported indashed lines.
We take the data from the Federal Reserve Economic Database (FRED), except ebpt , which we
took from Forni et al. (2024). The data are collected at the quarterly frequency in U.S. and cover
the period from 1973-Q1 to 2019-Q4 (188 observations). Inflation is computed as the change inthe logarithm of the GDP deflator (FRED mnemonic GDPDEF). Output gap is computed as thedeviation of the logarithm of real GDP (FRED mnemonic GDPC1) from the logarithm of potential
GDP (FRED mnemonic GDPPOT), which we take from the U.S. Congressional Budget Office’s
estimate. The series rt is set to the federal funds effective rate (FRED mnemonic FEDFUNDS).
The series ipt is taken from the index of industrial production (FRED mnemonic INDPRO). The
five time series are plotted in Figure 5. A three-lags VAR specification is selected according to the
Akaike information criterion. The “portmanteau” test for residual autocorrelation does not reject
the null hypothesis of zero autocorrelation at 0.01 level of significance.
We identify a SVAR model of the type of equation (15). As mentioned in Section 2, identification can be obtained by knowledge of the causal structure among innovation terms. As noted, for
23

8 10
6
0

2

4 0

2

4

pi

0

1

2

3

380

420

460 0

5

10

15

−8

−4

out
r
ip
ebp

1980

1990

2000

2010

2020

Figure 5: Employed time series. From the top, inflation (π), output gap (o) and nominal interest
rate (r). Data sample from 1953-Q3 to 2019-Q4.
the sake of estimating structural impulse response functions, knowledge of the topological orderis sufficient in our framework, given the recursiveness assumption. Therefore, as in the simulation
analysis, we apply Algorithm 1 to the data (i.e., first step of RESIT to the reduced-form VAR
residuals). Having obtained the topological order and a set of (potential) graphical parents from
Algorithm 1, we feed Algorithm 2 with these two pieces of information. Thus, we can estimate
the structural IRFs, setting the number of simulations for the Monte Carlo integration N = 100.
Both in Algorithm 1 and Algorithm 2, estimates of ϕk (·) (i.e., functional dependence between
innovation terms) are obtained through Gaussian process regressions.
To address uncertainty in the estimates of the structural IRFs, we perform a bootstrap analysis
with 500 iterations. For each bootstrap iteration, we estimate a new VAR model and, maintaining
the topological order that we got from the application of Algorithm 1 to the original data, we

24compute new contemporaneous effects and structural shocks.9 On the basis of this, we calculate
structural IRFs following Algorithm 2.
The topological order that we obtain by applying Algorithm 1 to our time series is hpit , ipt ,outt , ebpt , rt i, which supports the idea that the federal funds rate is set in reaction to movementsof the other variables, since the Fed is supposed to continuously monitor the macroeconomic andfinancial conditions of the system. Inflation, macroeconomic activity, financial conditions do notrespond to changes in the federal funds rate within the period, according to our findings. This
is in tune with typical recursiveness assumptions used in SVAR-based monetary analysis (seeChristiano et al., 1999).
For each nonlinear IRF we analyze, we consider the response of a variable y(t+h),k to a shock
εt,j of magnitude equal to one standard deviation of the variable yt,j .10 We do it for h = 0, . . . , 20.
We compare these nonlinear IRFs with linear IRFs estimated using the following scheme: (i) we
impose the Cholesky order derived from the topological order estimated by Algorithm 1 (i.e. hpit ,ipt , outt , ebpt , rt i); (ii) hence we get the structural shocks and the coefficients of the impact
matrix; (iii) using the latter terms as input, we apply Algorithm 2 and get IRFs.
Figure 6 (left panels) shows the median effects of a positive interest rate shock (i.e. a contractionary monetary policy shock) on the other variables, derived from bootstrapped nonlinear
and linear IRFs. Confidence intervals refer to the interquartile range (i.e. the 25th and 75th percentiles). Notice that the effects of an expansionary monetary policy shock (right panels) are by
construction completely symmetric to the contractionary ones, figuring r as last in the topologicalorder. Figure 7 shows the median effects of a positive (left panel) and negative (right panel) inflation shock. In Appendix D (see Figure 10 and Figure 11), we also show the effects of the output
and the excess bond premium shocks.
As regards the r (positive) shock (Figure 6, left panel), the reaction of inflation is positive
for two quarters and then becomes negative after one year. The reactions of industrial production
and output gap are decisively negative. The excess bond premium responds positively, but only
after 3 quarters. Notice that the nonlinear and linear scheme cannot depart from each other byconstruction as regards this shock.
Looking at the responses of a positive inflation shock (Figure 7, left panels), we see thatindustrial production has an instantaneous positive reaction. This holds both for the linear andnonlinear schemes. If we look at the responses of a negative inflation shock (Figure 7, rightpanels), the instantaneous effects are symmetric to the positive shock in the linear scheme byconstruction. But we find an asymmetry in the nonlinear scheme: the responses of ip to thenegative inflation shock are also positive at the impact. In the subsequent periods, for 2-3 quarters,
the effects of a positive pi shock on ip vanish and then backfire, under both the linear and nonlinear
9
We maintain the same topological order at each bootstrap iteration because we do not want to account for uncertainty about inference of the causal structure here. To do that, one can adapt the stability selection procedure by
Meinshausen and Bühlmann (2010).
10
As starting date, we simply choose t∗ = 4, considering that the number of lags is equal to 3. Note, however, that
due to the linearity of the reduced-form model, different choices would not affect our results.

25

schemes. Under the latter scheme, the effects of a negative pi shock on ip remain on the positivesign. Other asymmetries are detectable also as regards the reactions of out, ebp, and r to thesame shock: the instantaneous reactions of output gap and excess bond premium (in the nonlinear
scheme) are not symmetric between the positive and negative pi shock. In particular, although
output gap reacts instantaneously to the pi shock with the opposite sign, the absolute value ofsuch a reaction is greater and more evident (considering the position of confidence interval above
the zero line) in the case of a negative pi shock. The short-term response of the federal funds
rate (r) to positive and negative inflation shocks under the nonlinear scheme is noteworthy for its
asymmetry. Notably, when the system experiences a positive inflation shock, the federal funds
rate reacts much more strongly, nearly doubling in absolute terms compared to its response toa negative inflation shock, as regards the point estimates. In terms of confidence intervals, theimpact is significantly positive after a positive inflation shock, but not significant after a negative
inflation shock. This suggest that the FED is more reactive, in terms of changing the short terminterest rate, to a positive inflation shock than a negative one.
Looking at Figures 10 and 11 (Appendix D), we can detect some asymmetries in the contemporaneous responses of r to both the out and ebp shocks. In particular, looking at the nonlinear
scheme, the contemporaneous median response of r is positive to out shocks, no matter the sign,and, conversely, the contemporaneous median response of r is negative to ebp shocks, no matterthe sign.

5 Conclusions
There is a departure from the linear model with Gaussian disturbances that can be exploited for
identification, namely nonlinearity, as long as the disturbances are additive. We have shown that
this feature can be exploited to identify a structural VAR model that contains nonlinearities in the
cross-sectional structure. The identification criterion we have exploited in this paper is based on
the following simple idea: in case of nonlinearity it is admitted an additive noise model (i.e. amodel in which disturbances are independent of covariates) from the cause to the effect, but, inthe generic case (see conditions C1-C3 in section 2.3), not from the effect to the cause. Thus, by
iterative regressions and subsequent independence tests, it is possible to determine the contemporaneous causal structure underlying a VAR model, under the assumption of recursiveness andmutual independence of the shocks. On this basis, we can recover structural impulse response
functions, which we define as difference of conditional expectations and estimate through Monte
Carlo integration. It also turns out that, for the sake of estimating impulse response functions,the inferred causal structure can contain some redundant links, as long as the topological orderis correct. Simulation results have shown that the proposed search procedure and our scheme toestimate impulse response functions perform correctly under data generating processes satisfying
the theoretical conditions. Moreover, if the data generating process is nonlinear, they outperform
methods based on linearity assumptions, even under the correct zero restrictions. Our empiricalanalysis has shown that taking into account nonlinearity helps recovering the contemporaneous26

IRF of pi to a −1sd r shock

−1.0

−1.5

−0.5

0.5 1.0

0.0 0.5 1.0 1.5

IRF of pi to a +1sd r shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

16

18

20

16

18

20

16

18

20

16

18

20

IRF of ip to a −1sd r shock

−1

−4

0

−3

1

−2

2

−1

3

0

4

1

IRF of ip to a +1sd r shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of out to a −1sd r shock

−0.5

−1.5

0.0

−1.0

0.5

−0.5

0.0

1.0

0.5

IRF of out to a +1sd r shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of ebp to a −1sd r shock

−0.2

−0.2

0.0

0.0

0.1

0.1

0.2

0.2

IRF of ebp to a +1sd r shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of r to a −1sd r shock

−4

0

−3

1

−2

2

−1

3

0

4

IRF of r to a +1sd r shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

Figure 6: Structural IRFs to a positive (left panels) and negative (right panels) interest rate shockof one-standard deviation magnitude. Each panel shows nonlinear structural IRFs obtained from
Algorithm 2 with GP regressions (red lines) and linear structural IRFs obtained from Algorithm
2 but fed with coefficients derived from Cholesky factorization (blue lines). Solid lines refer tothe median effects while dashed lines refer to the interquartile range computed over 500 bootstrap
simulations.

27

IRF of pi to a −1sd pi shock

−2.5

0.0

−1.5

1.0

−0.5

2.0

IRF of pi to a +1sd pi shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

16

18

20

16

18

20

16

18

20

16

18

20

IRF of ip to a −1sd pi shock

−3

0

−2

1

−1

2

0

3

1

4

IRF of ip to a +1sd pi shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of out to a −1sd pi shock

−1.0

0.0

−0.6

0.4

−0.2

0.8

0.2

IRF of out to a +1sd pi shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of ebp to a −1sd pi shock

−0.2

−0.1

−0.1

0.0

0.0

0.1

0.1

0.2

IRF of ebp to a +1sd pi shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of r to a −1sd pi shock

−0.8

0.0

0.4

−0.4

0.8

0.0

1.2

IRF of r to a +1sd pi shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

Figure 7: Structural IRFs to a positive (left panels) and negative (right panels) inflation shockof one-standard deviation magnitude. Each panel shows nonlinear structural IRFs obtained from
Algorithm 2 with GP regressions (red lines) and linear structural IRFs obtained from Algorithm
2 but fed with coefficients derived from Cholesky factorization (blue lines). Solid lines refer tothe median effects while dashed lines refer to the interquartile range computed over 500 bootstrap
simulations.

28

structure of the SVAR system. It may also make a substantial difference in terms of predicting
the effect of macroeconomic shocks to the economy, particularly because it allows for asymmetriceffects between positive and negative structural shocks.
The idea of exploiting nonlinearity for identification has been put into practice in this paperby recovering a topological order through Algorithm 1, but one can, of course, conceive different
ways to implement it. Other procedures worth exploring are those in the family of the scorebased approach to causal discovery (see, e.g., Peters et al., 2014: section 4.2, and Bühlmann et al.,2014). A limitation of our procedure can be seen in the recursiveness assumption. This has beenrelaxed, in the bivariate case, by Mooij et al. (2011), but further research is needed for the case withmore than two variables and, more in general, for its adaptation to the structural VAR framework.
Another interesting path of research is the nonlinear ICA approach, which has been pioneered byHyvärinen and Pajunen (1999) and recently developed by Monti et al. (2020) and Gunsilius andSchennach (2023). This class of blind source separation methods departs from the additive noisemodel framework studied in this paper and its applicability to macroeconomic data is still to be
explored. Our hope is that the potentiality that nonlinearity has shown for identification will foster
the adoption of this idea within the field of econometrics.

References
Bacchiocchi, E. and T. Kitagawa (2022). Locally- but not globally-identified SVARs. Working
Paper, Bologna.
Blanchard, O. J. and D. Quah (1989). The dynamic effects of aggregate demand and supply
disturbances. The American Economic Review 79(4), 655–673.
Brunnermeier, M., D. Palia, K. A. Sastry, and C. A. Sims (2021). Feedbacks: financial markets
and economic activity. American Economic Review 111(6), 1845–79.
Bühlmann, P., J. Peters, and J. Ernest (2014). CAM: Causal additive models, high-dimensional
order search and penalized regression. The Annals of Statistics, 2526–2556.
Choi, T. and M. J. Schervish (2007). On posterior consistency in nonparametric regression problems. Journal of Multivariate Analysis 98(10), 1969–1987.
Christiano, L. J., M. Eichenbaum, and C. L. Evans (1999). Monetary policy shocks: What have
we learned and to what end? Volume 1 of Handbook of Macroeconomics, pp. 65–148. Elsevier.
Comon, P. (1994). Independent component analysis, a new concept? Signal Processing 36(3),287–314.
Deistler, M. and H.-G. Seifert (1978). Identifiability and consistent estimability in econometricmodels. Econometrica, 969–980.
Demiralp, S. and K. D. Hoover (2003). Searching for the causal structure of a vector autoregression. Oxford Bulletin of Economics and statistics 65, 745–767.
29

Eriksson, J. and V. Koivunen (2004). Identifiability, separability, and uniqueness of linear ICA
models. IEEE Signal Processing Letters 11(7), 601–604.
Fan, J. and I. Gijbels (2018). Local polynomial modelling and its applications. Routledge.
Faust, J. (1998). The robustness of identified VAR conclusions about money. In CarnegieRochester conference series on public policy, Volume 49, pp. 207–244. Elsevier.
Fiorentini, G. and E. Sentana (2023). Discrete mixtures of normals pseudo maximum likelihood
estimators of structural vector autoregressions. Journal of Econometrics 235(2), 643–665.
Forni, M., L. Gambetti, N. Maffei-Faccioli, and L. Sala (2024). Nonlinear transmission of financial
shocks: Some new evidence. Journal of Money, Credit and Banking 56(1), 5–33.
Gabrielsen, A. (1978). Consistency and identifiability. Journal of Econometrics 8(2), 261–263.
Gouriéroux, C., A. Monfort, and J.-P. Renne (2017). Statistical inference for independent component analysis: Application to structural VAR models. Journal of Econometrics 196(1), 111–126.
Gretton, A., K. Fukumizu, C. Teo, L. Song, B. Schölkopf, and A. Smola (2007). A kernel statistical
test of independence. Advances in neural information processing systems 20.
Guay, A. (2021). Identification of structural vector autoregressions through higher unconditional
moments. Journal of Econometrics 225(1), 27–46.
Guerini, M. and A. Moneta (2017). A method for agent-based models validation. Journal of
Economic Dynamics and Control 82, 125–141.
Gunsilius, F. and S. Schennach (2023). Independent nonlinear component analysis. Journal of theAmerican Statistical Association 118(542), 1305–1318.
Györfi, L., M. Kohler, A. Krzyzak, H. Walk, et al. (2002). A distribution-free theory of nonparametric regression, Volume 1. Springer.
Härdle, W., H. Lütkepohl, and R. Chen (1997). A review of nonparametric time series analysis.
International statistical review 65(1), 49–72.
Hastie, T. J. and R. J. Tibshirani (2017). Generalized additive models. Routledge.
Herwartz, H. (2018). Hodges–Lehmann detection of structural shocks–an analysis of macroeconomic dynamics in the Euro area. Oxford Bulletin of Economics and Statistics.
Herwartz, H., A. Lange, and S. Maxand (2022). Data-driven identification in SVARs—when
and how can statistical characteristics be used to unravel causal relationships? Economic Inquiry 60(2), 668–693.
Höppner, F., C. Melzer, and T. Neumann (2008). Changing effects of monetary policy in theus–evidence from a time-varying coefficient var. Applied Economics 40(18), 2353–2360.

30

Hoyer, P., D. Janzing, J. M. Mooij, J. Peters, and B. Schölkopf (2008). Nonlinear causal discovery
with additive noise models. Advances in neural information processing systems 21, 689–696.
Hussain, S. M. and S. Malik (2016). Asymmetric effects of exogenous tax changes. Journal of
Economic Dynamics and Control 69, 268–300.
Hyvärinen, A. (2013). Independent component analysis: recent advances. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371(1984),
20110534.
Hyvärinen, A., J. Karhunen, and E. Oja (2001). Independent component analysis. John Wiley &Sons.
Hyvärinen, A. and P. Pajunen (1999). Nonlinear independent component analysis: Existence anduniqueness results. Neural networks 12(3), 429–439.
Kilian, L. (2014). Oil price shocks: Causes and consequences. Annu. Rev. Resour. Econ. 6(1),133–154.
Kilian, L. and H. Lütkepohl (2017). Structural vector autoregressive analysis. Cambridge University Press.
Koop, G., M. H. Pesaran, and S. M. Potter (1996). Impulse response analysis in nonlinear multivariate models. Journal of econometrics 74(1), 119–147.
Koopmans, T. C. and O. Reiersøl (1950). The identification of structural characteristics. The
Annals of Mathematical Statistics 21(2), 165–181.
Lanne, M. and J. Luoto (2021). GMM estimation of non-gaussian structural vector autoregression.
Journal of Business & Economic Statistics 39(1), 69–81.
Lanne, M. and H. Lütkepohl (2008). Identifying monetary policy shocks via changes in volatility.
Journal of Money, Credit and Banking 40(6), 1131–1149.
Lanne, M. and H. Lütkepohl (2010). Structural vector autoregressions with nonnormal residuals.
Journal of Business & Economic Statistics 28(1), 159–168.
Lanne, M., M. Meitz, and P. Saikkonen (2017). Identification and estimation of non-Gaussian
structural vector autoregressions. Journal of Econometrics 196(2), 288–304.
Lewis, D. J. (2021). Identifying shocks via time-varying volatility. The Review of Economic
Studies.
Lo, M. C. and J. Piger (2005). Is the response of output to monetary policy asymmetric? Evidence
from a regime-switching coefficients model. Journal of Money, credit and Banking, 865–886.
Matzkin, R. L. (2007). Nonparametric identification. Handbook of econometrics 6, 5307–5368.

31

Meinshausen, N. and P. Bühlmann (2010). Stability selection. Journal of the Royal StatisticalSociety: Series B (Statistical Methodology) 72(4), 417–473.
Mesters, G. and P. Zwiernik (2022). Non-independent components analysis. arXiv preprint
arXiv:2206.13668.
Moneta, A. (2008). Graphical causal models and VARs: An empirical assessment of the realbusiness cycles hypothesis. Empirical Economics 35(2), 275–300.
Moneta, A., D. Entner, P. O. Hoyer, and A. Coad (2013). Causal inference by independent component analysis: Theory and applications. Oxford Bulletin of Economics and Statistics 75(5),705–730.
Monti, R. P., K. Zhang, and A. Hyvärinen (2020). Causal discovery with general non-linear
relationships using non-linear ica. In Uncertainty in artificial intelligence, pp. 186–195. PMLR.
Montiel Olea, J. L., M. Plagborg-Møller, and E. Qian (2022). SVAR identification from higher
moments: Has the simultaneous causality problem been solved? AEA Papers and Proceedings 112, 481–85.
Montiel Olea, J. L., J. H. Stock, M. W. Watson, et al. (2021). Inference in structural vector
autoregressions identified with an external instrument. Journal of Econometrics 225(1), 74–87.
Mooij, J. M., D. Janzing, T. Heskes, and B. Schölkopf (2011). On causal discovery with cyclic
additive noise models. Advances in neural information processing systems 24.
Pearl, J. (2009). Causality. Models, Reasoning, and Inference. Cambridge University Press.
Peters, J., D. Janzing, and B. Schölkopf (2017). Elements of causal inference: foundations andlearning algorithms. The MIT Press.
Peters, J., J. Mooij, D. Janzing, and B. Schölkopf (2014). Causal discovery with continuous
additive noise models. Journal of Machine Learning Research 15(1), 2009–2053.
Ramey, V. A. (2016). Macroeconomic shocks and their propagation. Handbook of Macroeconomics 2, 71–162.
Reiersøl, O. (1950). Identifiability of a linear relation between variables which are subject to error.
Econometrica, 375–389.
Rigobon, R. (2003). Identification through heteroskedasticity. Review of Economics and Statistics 85(4), 777–792.
Roehrig, C. S. (1988). Conditions for identification in nonparametric and parametric models.
Econometrica, 433–447.
Romer, C. D. and D. H. Romer (1989). Does monetary policy matter? a new test in the spirit offriedman and schwartz. NBER macroeconomics annual 4, 121–170.
32

Rothenberg, T. J. (1971). Identification in parametric models. Econometrica, 577–591.
Rubio-Ramirez, J. F., D. F. Waggoner, and T. Zha (2010). Structural vector autoregressions:
Theory of identification and algorithms for inference. The Review of Economic Studies 77(2),665–696.
Sentana, E. and G. Fiorentini (2001). Identification, estimation and testing of conditionally heteroskedastic factor models. Journal of Econometrics 102(2), 143–164.
Shimizu, S., P. O. Hoyer, A. Hyvärinen, and A. Kerminen (2006). A linear non-Gaussian acyclic
model for causal discovery. Journal of Machine Learning Research 7, 2003–2030.
Shimizu, S., T. Inazumi, Y. Sogawa, A. Hyvärinen, Y. Kawahara, T. Washio, P. O. Hoyer, andK. Bollen (2011). Directlingam: A direct method for learning a linear non-gaussian structural
equation model. The Journal of Machine Learning Research 12, 1225–1248.
Sims, C. A. (1980). Macroeconomics and reality. Econometrica, 1–48.
Sims, C. A. (2020). SVAR identification through heteroskedasticity with misspecified regimes.
Princeton University.
Spirtes, P., C. N. Glymour, R. Scheines, and D. Heckerman (2000). Causation, prediction, andsearch. MIT press.
Stock, J. H. and M. W. Watson (2001). Vector autoregressions. Journal of Economic perspectives 15(4), 101–115.
Swanson, N. R. and C. W. Granger (1997). Impulse response functions based on a causal approachto residual orthogonalization in vector autoregressions. Journal of the American Statistical
Association 92(437), 357–367.
Williams, C. K. and C. E. Rasmussen (2006). Gaussian processes for machine learning, Volume 2.
MIT press Cambridge, MA.
Wold, H. O. (1960). A generalization of causal chain models. Econometrica, 443–463.
Zhang, K. and A. Hyvärinen (2009). On the identifiability of the post-nonlinear causal model.
Proceedings of the 25th Annual Conference on Uncertainty in Artificial Intelligence (UAI), 647–
655.

33

Appendix
A

Complete RESIT algorithm

Algorithm 1full VAR + RESIT (Phases 1-2 are from Peters et al., 2014)
1: Input: A K-dimensional time series vector (yt,1 , . . . , yt,K )′2: PHASE 0: Estimate the reduced-form model.
3: Estimate a reduced-form time series model of the class GVAR-ANIM, see eq. (1) and extract

residuals u
bt = (but,1 , . . . , ubt,K )′ . To simplify the notation, let us call u
bt,k ≡ vk , for k =1, . . . , K.
4: S := 1, ..., K, π := [ ]5: PHASE 1: Determine topological order.
6: repeatfor k ∈ S doRegress vk on {vi }i∈S\{k} and obtain residuals ek
9:Measure dependence between ek and {vi }i∈S\{k}10:
end for11:
Let k ∗ be the k with the weakest dependence
12:
S := S \ k ∗13:pa(k ∗ ) := S
14:
π := [k ∗ , π]
15: until #S = 0
16: Output: π
7:

8:

17: PHASE 2 (pruning): Remove superfluous edges.
18: for k ∈ {2, . . . , K} dofor p ∈ pa(π(k)) do
Regress vπ(k) on {vi }i∈pa(π(k))\{p} .
21:
if residuals are independent of {vi }i∈{π(1),...,π(k−1)} then22:pa(π(k)) := pa(π(k)) \ {p}23:
end if24:
end for
25: end for
26: Output: (pa(1), . . . , pa(K))
19:

20:

Note: with pa(k) we refer to the set of indices associated to the variables in P a(vk ), for any k.

34

BDiscussion on consistency

The method proposed provides the topological order and causal structure among innovation terms
and allows us to estimate structural impulse response functions. As one gets more data, one would
hope that the recovered order/structure is correct and that the estimated IRFs would converge tothe true IRFs, under the assumptions that the underlying data generating process is GVAR-ANIMand that the relationships between innovation terms (equation 3) can be described by a multivariate
ANM under assumptions (C1)-(C3). While it is clear that an unidentified model cannot be consistently estimated, it is not necessarily true that identifiability implies consistency (see Deistler andSeifert, 1978; Gabrielsen, 1978). Studying consistency and other statistical properties (e.g. asymptotic normality, efficiency) of the entire procedure, considering the different options in terms ofmodel specification and estimation, is a task we leave for future research. In the following, however, we refer to some key results for such task.
1. Assume a joint distribution P (u) generated by a GVAR-ANIM, associated to a DAG G
representing the structural relations between innovation terms ut = (ut,1 , . . . , ut,K )′ . Assume
that: (i) conditions C1-C3 are satisfied; (ii) the ut are consistently estimated from the reducedform VAR model (Alg. 1, Phase 0); (iii) ϕk (·) is estimated with a consistent regression method,for each k; (iv) the dependence measures (Alg. 1, Phase 1) is correct and independence tests(Alg. 1, Phase 2, see Appendix A) are consistent. Then the output of the (full) Algorithm 1 isguaranteed to be correct. This statement generalizes Theorem 34 in Peters et al. (2014) to the caseof VAR-estimated input and presupposes the identifiability result of Theorem 3 (section 2.3).
2. The regression method to estimate ϕk (·) is user-specified in our procedure, with the caveatthat linear regression is not allowed, except in the case in which the structural shocks are nonnormal. In our simulation analysis and empirical application, we use Gaussian process regressions, whose consistency has been studied by Williams and Rasmussen (2006, Chapter 7) in thegeneral case and by Choi and Schervish (2007) in the univariate case. Following the latter authors, let us consider Yi = η(Xi ) + ǫi (i = 1, . . . , n), where we are interested in estimating
η(x) = E[Y |X = x]. Assume: (i) ǫi ∼ N (0, σ 2 ) or DE(0, σ) (DE stands for double exponential, i.e. Laplace, distribution); (ii) η(·) is a random process with a Gaussian process priordistribution: η(·) ∼ GP (µ(·), R(·, ·)); (iii) smoothness conditions on µ(·), R(·, ·)), and η(·) (inparticular, continuously differentiable mean function µ(x) and covariance function R(·, ·) having
continuous fourth partial derivatives). Let the posterior distribution pn,N (η|X, Y ) be, for eachneighborhood N of the true regression function η0 and each sample size n: pn,N = Pr({η ∈N }|X1 , . . . , Xn , Y1 , . . . , Yn ). Under the assumptions above, Choi and Schervish (2007) showthat the posterior distribution of η is almost surely consistent, i.e., for every N , limn→∞ pn,N = 1
a.s..
As an alternative nonparametric regression method to estimate ϕk (·), let us consider the NadarayaWatson estimator. Considering again the regression function η(x), this estimator computes η̂(x) =Pn
PnKh (Xi −x), being Kh a kernel function with bandwidth h. It has been
i=1 wi Yi where wi =j=1 Kh (Xj −x)established (see, e.g., Györfi et al., 2002) that if nh → ∞ and n → ∞ then the estimator is con-35

sistent (the result has also been extended to the multivariate case).
3. Conditions for a consistent estimation of a fully non-parametric (but time invariant) reducedform GVAR-ANIM (see equation 1) are discussed in Härdle et al. (1997). Nonparametric kernel estimation turns out to be consistent under the conditions discussed above for the NadarayaWatson estimator with the addition that the underlying process must be α- or φ-mixing. In our
simulation and empirical application, we actually estimate a linear reduced-form VAR. Here theleast squares estimator is consistent and asymptotic normal (even in the presence of unit roots witha number of lags strictly greater than one) under general conditions (see Kilian and Lütkepohl,
2017). Notice, however, that we cannot assume normality of the reduced-form residuals here,
since this would imply a linear form in the conditional expectation function. Conditions for consistency in the case of a linear non-Gaussian SVAR are studied by Gouriéroux et al. (2017).
4. The dependence measure and the independence test used in the Algorithm 1 (full version,
see Appendix A) are also user-specified in our procedure, under the condition that they accountfor order statistics. In our applications, we use the nonparametric measure of dependence andthe kernel independence test proposed by Gretton et al. (2007), who prove consistency of thetest by deriving distribution of the test both under the null hypothesis of independence and under
the alternative. The distribution of the test statistics, as the sample size approaches infinity, isparameterised in terms of kernels of the data.

36

C

Further simulation results: common-cause and v- structure

Table 3: Average MSE, across 200 simulations, between IRFs estimated using different schemes(i.e. CHOL, TDAG, TTOP, ALG1) and the theoretical IRFs. MSE are calculated over the responses of yt,2 and yt,3 to the first shock (of one standard-deviation magnitude) and over the first
five periods. The contemporaneous structure is the “common cause”. We consider different values of T (sample size) and three different model classes (Linear with Gaussian noises, non-linearwith Gaussian noises and linear with non-Gaussian noises). Standard errors are reported in roundbrackets.
CHOL
T

ALG1

250

500

1000

250

500

1000

LinearGaussian

0.031

0.019

0.014

0.200

0.209

0.189

(1.29E-03)(7.05E-04)(4.74E-04)(7.27E-03)(7.40E-03)(7.22E-03)

Nonlinear
Gaussian

1.403

1.01

1.28

0.185

0.12

0.077

(1.04E-01)(7.00E-02)(8.80E-02)(1.20E-02)(9.00E-03)(4.00E-03)

Linear
Non-Gaussian

0.068

0.045

0.035

0.268

0.162

0.089

(3.09E-03)(1.74E-03)(1.58E-03)(1.22E-02)(7.12E-03)(4.51E-03)

TDAG
T

250

500

TTOP
1000

250

500

1000

LinearGaussian

0.028

0.023

0.018

0.031

0.019

0.017

(1.01E-03)(8.27E-04)(6.16E-04)(1.01E-03)(6.85E-04)(6.20E-04)

Nonlinear
Gaussian

0.155

0.093

0.094

0.14

0.089

0.081

(8.00E-03)(9.00E-03)(8.00E-03)(7.00E-03)(6.00E-03)(8.00E-03)

0.179

0.102

0.069

0.182

0.102

0.070

(9.27E-03)(4.80E-03)(3.30E-03)(1.04E-02)(5.37E-03)(3.52E-03)

Linear
Non-Gaussian37

Table 4: Average MSE, across 200 simulations, between IRFs estimated using different schemes(i.e. CHOL, TDAG, TTOP, ALG1) and the theoretical IRFs. MSE are calculated over the responses of yt,2 and yt,3 to the first shock (of one standard-deviation magnitude) and over the first
five periods. The contemporaneous structure is the “v-structure”. We consider different valuesof T (sample size) and three different model classes (Linear with Gaussian noises, non-linearwith Gaussian noises and linear with non-Gaussian noises). Standard errors are reported in roundbrackets.
CHOL
T

250

500

ALG1
1000

250

500

1000

LinearGaussian

0.019

0.010

0.008

0.084

0.089

0.094

(8.23E-04)(4.23E-04)(2.97E-04)(4.56E-03)(5.01E-03)(5.11E-03)

Nonlinear
Gaussian

0.936

0.678

0.861

0.118

0.064

0.048

(9.90E-02)(6.80E-02)(8.50E-02)(1.10E-02)(5.00E-03)(4.00E-03)

0.057

0.034

0.023

0.255

0.135

0.057

(2.79E-03)(1.42E-03)(1.26E-03)(1.40E-02)(7.07E-03)(3.44E-03)

Linear
Non-Gaussian

TDAGTLinearGaussian

TTOP

250

500

1000

250

500

1000

0.016

0.012

0.009

0.023

0.014

0.011

(7.14E-04)(5.00E-04)(3.65E-04)(9.47E-04)(5.73E-04)(4.37E-04)

Nonlinear
Gaussian

0.088

0.049

0.053

0.089

0.053

0.046

(7.00E-03)(6.00E-03)(7.00E-03)(7.00E-03)(6.00E-03)(5.00E-03)

Linear
Non-Gaussian

0.121

0.064

0.049

0.207

0.102

0.066

(8.11E-03)(3.94E-03)(3.19E-03)(1.26E-02)(5.60E-03)(3.95E-03)

38

IRF from shock 1 to variable 2 (T =1000)
THEOR
CHOL
ALG1

0

0

1

1

2

2

3

3

THEOR
CHOL
ALG1

4

IRF from shock 1 to variable 2 (T =250)

0

2

4

6

8

0

2

4

6

8

IRF from shock 1 to variable 3 (T =1000)
3.0

IRF from shock 1 to variable 3 (T =250)

2.5

THEOR
CHOL
ALG1

0.0

0.0

0.5

0.5

1.0

1.0

1.5

1.5

2.0

2.0

2.5

THEOR
CHOL
ALG1

0

2

4

6

8

0

2

4

6

8

Figure 8: DGP with “common cause”. Theoretical (THEOR), linear (CHOL), and nonlinear(ALG1) IRFs for T = 250 (left panels) and T = 1000 (right panels). Each solid line exhibits
the average, across 200 simulations, IRF from a unitary shock 1 (εt,1 ) to variable 2 (yt+h,2 ) (toppanels) or variable 3 (yt+h,3 ) (bottom panels), for h = 0, . . . , 9. Confidence interval at 68% arereported in dashed lines.

39

IRF from shock 1 to variable 2 (T =1000)

0.6

IRF from shock 1 to variable 2 (T =250)

0.5

THEOR
CHOL
ALG1

0.0

0.0

0.1

0.2

0.2

0.3

0.4

0.4

THEOR
CHOL
ALG1

0

2

4

6

8

0

IRF from shock 1 to variable 3 (T =250)

2

4

6

8

IRF from shock 1 to variable 3 (T =1000)
4

THEOR
CHOL
ALG1

0

0

1

1

2

2

3

3

THEOR
CHOL
ALG1

0

2

4

6

8

0

2

4

6

8

Figure 9: DGP with “v-structure”. Theoretical (THEOR), linear (CHOL), and nonlinear (ALG1)
IRFs for T = 250 (left panels) and T = 1000 (right panels). Each solid line exhibits the average,
across 200 simulations, IRF from a unitary shock 1 (εt,1 ) to variable 2 (yt+h,2 ) (top panels) orvariable 3 (yt+h,3 ) (bottom panels), for h = 0, . . . , 9. Confidence interval at 68% are reported indashed lines.

40

D

Empirical results: further IRFs
IRF of pi to a −1sd out shock

−0.5

−0.5

0.0

0.0

0.5

0.5

1.0

IRF of pi to a +1sd out shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

16

18

20

16

18

20

16

18

20

16

18

20

IRF of ip to a −1sd out shock

−2

−4

0

−2

1

2

0

3

1

4

2

IRF of ip to a +1sd out shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of out to a −1sd out shock

−0.5

−2.5

0.5

−1.5

1.5

−0.5

2.5

0.5

IRF of out to a +1sd out shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of ebp to a −1sd out shock

−0.4

−0.4

−0.2

−0.2

0.0

0.0

0.2

0.2

0.4

0.4

IRF of ebp to a +1sd out shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of r to a −1sd out shock

−0.5

−1.5

0.0

0.5

−0.5

1.0

0.5

1.5

IRF of r to a +1sd out shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

Figure 10: Structural IRFs to a positive (left panels) and negative (right panels) output gap shockof one-standard deviation magnitude. Each panel shows nonlinear structural IRFs obtained from
Algorithm 2 with GP regressions (red lines) and linear structural IRFs obtained from Algorithm
2 but fed with coefficients derived from Cholesky factorization (blue lines). Solid lines refer tothe median effects while dashed lines refer to the interquartile range computed over 500 bootstrap
simulations.

41

IRF of pi to a −1sd ebp shock

−0.2

−0.4

0.0

−0.2

0.2

0.0

0.4

0.2

IRF of pi to a +1sd ebp shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

16

18

20

16

18

20

16

18

20

16

18

20

IRF of ip to a −1sd ebp shock

−2.5

−0.5

0.5

−1.5

1.5

−0.5

0.5

2.5

IRF of ip to a +1sd ebp shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of out to a −1sd ebp shock

−1.0

−0.4

−0.6

0.0

−0.2

0.4

0.2

0.8

IRF of out to a +1sd ebp shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of ebp to a −1sd ebp shock

−0.6

0.0

−0.4

0.2

−0.2

0.4

0.0

0.6

IRF of ebp to a +1sd ebp shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

IRF of r to a −1sd ebp shock

−0.2

−1.0

0.2

−0.6

0.6

−0.2

1.0

0.2

IRF of r to a +1sd ebp shock

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

Figure 11: Structural IRFs to a positive (left panels) and negative (right panels) excess bond premium shock of one-standard deviation magnitude. Each panel shows nonlinear structural IRFs
obtained from Algorithm 2 with GP regressions (red lines) and linear structural IRFs obtained
from Algorithm 2 but fed with coefficients derived from Cholesky factorization (blue lines). Solidlines refer to the median effects while dashed lines refer to the interquartile range computed over
500 bootstrap simulations.

42