Identification of Impulse Response Functions for Nonlinear Dynamic Models

Identification of Impulse Response Functions for Nonlinear Dynamic Models

\setstretch{1}

Introduction

It is well known that in a linear dynamic framework, such as the Gaussian Structural Vector Autoregressive (SVAR) model, there exists an identification issue for the structural innovations and respectively, the Impulse Response Function (IRF). However, when a nonlinear dynamic is introduced for the autoregressive process, it is possible for some of these problems to be mitigated. The objective in this paper is to highlight the key identification issues in this context and to discuss some possible remedies. \\

We begin by establishing the nonlinear autoregressive representation for a Markov process. This leads to the definition of nonlinear innovations and their potential interpretation in terms of shocks. In general however, the multivariate nonlinear autoregressive representation, the nonlinear innovations, the associated shocks and IRF are not defined in a unique way without additional parametric or nonparametric restrictions. We characterize the identification issue in this setting by describing the identified set (IS), and demonstrating that underidentification arises from transformations which preserve either the normal or uniform distribution on $[0,1]^n$. \\

There are a number of ways to resolve such identification issues in a nonlinear multivariate dynamic model. This arises if the innovations are non Gaussian, which facilitates identification based on the different power moments of the errors, or if the underlying dynamic independent sources have different nonlinear dynamics. In both cases, identification is conditional on the uniqueness of deconvolution. For instance, we can impose independence restrictions on non-Gaussian sources, leading to identifications by means of linear ICA, which facilitates identification based on the different power moments of the errors. Under partial identification, an appropriate use of algorithmic estimation (learning) approaches and of simulations can allow the recovery of identifiable (functional) parameters, such as the term structure of some Pseudo Impulse Response Functions (PIRF). Identification can also be facilitated in nonlinear dynamic factor models, when the underlying factors are independent with different dynamics. Howeverm it is important to highlight that these notions and key results ultimately depend on the universe, that is, on the set of series under consideration.\\

This paper is organized as follows. Section 2 describes the nonlinear autoregressive representation of a Markov process and presents the notion of nonlinear impulse response function in the nonlinear dynamic framework. Section 3 explains the identification issue arising from the multiplicity of nonlinear innovations. Section 4 provides methods in which the identification issue can be reduced. Section 5 discusses identification in the context of nonlinear dynamic factor models where the underlying factors feature different dynamics. Section 6 discusses the role of the universe in the context of identification, and Section 7 concludes. Proofs are gathered in the appendices and addition results are provided in the online appendices.

Nonlinear Framework and Impulse Response Function

Nonlinear Autoregressive Representation of a Markov Process

Let us first introduce the notion of the nonlinear autoregressive representation for a Markov process. We consider an $n$-dimensional Markov process $(y_t$) of order 1 with values in $\mathcal{Y}=\mathbb{R}^n$. If its distribution is continuous, the Markov condition can be expressed in terms of its transition density, that is $f(y_t|\underline{y_{t-1}})=f(y_t|y_{t-1})$, $\forall y_t,\underline{y_{t-1}}$, where $\underline{y_{t-1}}=(y_{t-1},y_{t-2},...)$. Equivalently the Markov condition can be represented by the proposition below.

Proof: See Appendix A.1. \\

Nonlinear dynamic features can be introduced in this representation, since $y_t$ is a nonlinear function of its $y_{t-1}$ for a given $\varepsilon_t$, and/or a nonlinear function of $\varepsilon_t$ for a given $y_{t-1}$, and/or by nonlinear cross-effects of $y_{t-1}$ and $\varepsilon_{t}$. Moreover, such a representation exists even if $y_t$ has marginal and/or conditional fat tails. \\

The one-dimensional case has been discussed in Gouri\'eroux and Jasiak (2005, Section 2.1). Let us denote $F(\cdot|y_{t-1})$ as the conditional cumulative distribution function (c.d.f.) associated with $f(\cdot|y_{t-1})$. Under the assumption of Proposition 1, this function is invertible. We denote $Q(\cdot|y_{t-1})=F^{-1}(\cdot|y_{t-1})$ as its inverse, that is the conditional quantile function. Then we can choose:

where $\Phi$ denotes the cdf of the standard normal. This is the only solution, if we impose that the one-to-one function $g$ is also increasing in $\varepsilon$ \footnote{Another solution with $g$ decreasing is obtained by changing $\varepsilon_t$ into $-\varepsilon_t$.}. Proposition 1 shows that this representation corresponds to the inversion method used for drawing in a distribution [Gourieroux and Monfort (1997)]. This approach can be extended to the multidimensional framework. In particular, we have the following corollary:

Proposition 1 and Corollary 1 are in particular valid when the data themselves are i.i.d., that is when we look for independent sources $\varepsilon_t$ or $u_t$ with effects on $y_t$. This also says that the conditions of Gaussianity or of uniform distributions are simply normalization conditions. More generally, we could have introduced a nonlinear innovation with a distribution other than $N(0,1)$ or $U[0,1]$ for the strict white noise components.\\

This corresponds to nonlinear Independent Component Analysis (ICA). Corollary 2 shows that there is always a solution to the nonlinear ICA problem for continuous random vectors [see e.g. Hyvarinen and Pajunen (1999), Theorem 1].

Effect of Time Unit

The definition of the Gaussian nonlinear innovation in Proposition 1 (resp. of the uniform nonlinear innovation in Corollary 1) depends on the time unit. Let us discuss this by considering a time unit of 2. With this new time unit, the process $(y_{2\tau})$ is still a Markov process of order 1. Then we can write:

where the $\tilde{\varepsilon}_{2\tau}$'s are i.i.d. $N(0,Id)$ variables, defining the Gaussian nonlinear innovations at step 2. Simultaneously, we can also write:

Nonlinear Impulse Response Function

Our interest is to study the notion of shock and characterize its propagation mechanism in a nonlinear dynamic framework. It is important to discuss in this context the main assumptions in order to facilitate the introduction of shocks and the construction of their associated Impulse Response Function (IRF). For exposition, consider a bivariate framework with structural innovation $\varepsilon_t=(\varepsilon_{1,t},\varepsilon_{2,t}$). We consider “shocks" of magnitude $\delta_1$ on the first component $\varepsilon_{1,t}$ (or respectively of $\delta_2$ on $\varepsilon_{2,t}$). To get reasonable interpretations, we assume that:

Remark 1: The literature has suggested an alternative definition of nonlinear innovation (in the one dimensional framework) as:

[see Blanchard and Quah (1989), Koop, Pesaran and Potter (1996)]. It is easily checked that $\varepsilon_t^*$ is not independent of $y_{t-1}$ in general and thus this definition is not appropriate for shocking $\varepsilon_t^*$ in the construction of IRFs. \\

The nonlinear IRF corresponding to a given nonlinear autoregressive representation is defined as a comparison between the perturbed $(y^\delta_{t+h})$ and baseline $(y_{t+h})$ future trajectories of the process. At horizon 0, a transitory shock of magnitude $\delta$ hits the perturbed path, which yields:

Then, in all subsequent horizons, the perturbed and baseline paths can be defined recursively as follows:

for all $h=1,2,...$. The corresponding IRF is given by:

This is a (nonlinear) function of the horizon $h$ (therefore, we have a term structure of the IRF), the magnitude of the shock $\delta$, the information set summarized by $y_{t-1}$ and also the sequence of future innovations $\varepsilon_t,...,\varepsilon_{t+h}$. Moreover, the IRF is stochastic conditional on $y_{t-1}$, due to the presence of future innovations. Therefore, we may partly summarize its distribution conditional on $y_{t-1}$ through the lens of a conditional expectation:

$\forall \ h,\delta$, or a conditional covariance:

$\forall \ h,k,\delta$. \\

Remark 2: The IRF has a particular form in the special case of the linear Gaussian VAR framework\footnote{The Gaussian assumption is often implicit or explicit in the analysis of VAR(1) models [see e.g. Assumption 2 in Plagborg-Moller and Wolf (2021)].}:

Then we get:

which can be simplified to $(Id-\Phi)^{-1}(Id-\Phi^{k+1})D\delta$, if the autoregressive matrix $\Phi$ is such that $Id-\Phi$ is invertible. It is clear that the IRF in this context exhibits the following properties:

Indeed, the Gaussian VAR(1) is a very special case of dynamic models and such properties are in general, not indicative of what arises in a nonlinear dynamic framework.

Examples

Example 1: Conditionally Gaussian Model\\

The model can be written as:

where $\varepsilon_t\sim IIN(0,Id)$ and is assumed independent of $y_{t-1}$. This specification is used in Koop, Pesaran and Potter (1996) for instance. Except in the Gaussian linear VAR(1) framework: $m(y_{t-1})=A y_{t-1}, D(y_{t-1})=D$, independent of the past, the process is not conditionally Gaussian at horizons larger or equal to 2. At horizon 2, we get:

where the conditional distribution of $y_{t+1}$ given $y_{t-1}$ involves cross effects of $\varepsilon_t$, $\varepsilon_{t+1}$. Therefore, it is not Gaussian and it cannot be rewritten in the form:

where $\zeta_t$ is Gaussian. \\

Equivalently, this specification is not invariant by a change of time unit. In the one dimensional case, the specification (ref) includes the Double Autoregressive (DAR) model of order one [Weiss (1984), Borkovec and Kluppelberg (2001), Ling (2007)], introduced to account for conditional heteroscedasticity and given by:

where $\varepsilon_t$ is $IIN(0,1)$. The DAR model has a strictly stationary solution if the Lyapunov coefficient is negative:

The conditional heteroscedasticity can create fat tails for the stationary distribution of process $(y_t)$. This process has second-order moments if morever:

When $ \gamma^2 + \beta > 1$ and $\mathbb{E}\log |\gamma + \sqrt{\beta} \ \varepsilon|Z<0$, the process $(y_t)$ is stationary with infinite marginal variance (but finite conditional variance). \\

Its multivariate extension can be written as:

with $h_t=a + B(y^2_{1,t-1},...,y^2_{n,t-1})$ [see Zhu et al. (2017)]. \\

Example 2: Qualitative Threshold AR(1) Model \\

Such a model is defined by:

with $\varepsilon_t$ i.i.d. $N(0,\sigma^2)$. Then at horizon 2 we get:

Example 3: Time Discretized Diffusion Process\\

Due to the lack of coherency of the definitions with respect to the time unit, it is sometimes proposed to write the structural definition of innovations and/or impulse response functions in the infinitesimal time unit even if the observations are available in discrete time. This is a usual practice in mathematical finance or in engineering. The discrete time process is $y_t=y(t),t=1,2,...$, where the underlying process $y(\tau)$ is defined in continuous time $\tau \in (0,\infty)$ by a multivariate diffusion equation:

where $W$ is a multivariate Brownian motion with $\mathbb{V}[dW(\tau)]=Id. d\tau$. The diffusion equation (ref) is the analogue of the conditionally Gaussian model written in an infinitesimal time unit. For $n=1$, this class of time discretized diffusion contains the Gaussian AR(1) process, the autoregressive gamma (ARG) process, that is the time discretized Cox, Ingersoll, Ross process [Cox, Ingersoll and Ross (1985)], or the time discretized Jacobi process [Karlin and Taylor (1981), Gour\'ieroux, Jasiak and Sufana (2009) for its extension to the multivariate framework], with values in $(-\infty,+\infty),(0,+\infty),(0,1)$, respectively.

The Identification Issue

Except in irregular cases, the normalized nonlinear autoregressive representation is unique in the one-dimensional framework $n=1$ [see Appendix A.2. for an irregular case]. However, the uniqueness disappears in a multivariate setting. More precisely, let us consider the nonlinear autoregressive representation with uniform normalization and assume a true representation with true $g_0$, $u_t^0$. Then the identified set for $g$, $u_t$, is:

where $\overset{d}{=}$ means that the two variables have the same distribution. Since the functions $g_0$ and $g$ are invertible with respect to $u$, it follows that the identification issue arises from the possibility of having nonlinear transformations that preserve the uniform distribution on $[0,1]^n$.\footnote{The same reasoning can be applied with another normalization as the Gaussian one (see the proof in Section 3.2.)} This identification issue is linked to the identification issue in nonlinear Independent Component Analysis (ICA), where the literature considers the uniqueness of sources $u$ leading to a given output $y: \ y=g(u)$ [see Darmois (1953), Hyvarinen and Oja (2001), Roberts and Everson (2001) for the principles and practices of ICA.]. We are in the same framework with identical sources with varying information, that is, values of $y_{t-1}$.

Transformations Preserving the (Multivariate) Uniform Distribution

Proof: See Appendix A.3. \\

To summarize, for $n=1$, both the nonlinear autoregressive representation $g$ and the Gaussian nonlinear innovation $\varepsilon_t$ are identifiable. This is no longer the case for $n\geq2$, where there is a problem of nonparametric partial identification. Propositions 1 and 2 are valid under a property of null recurrence in the Markov process [see Tweedie (1975) for the definition of null recurrence], where, for any $y_{t-1}$, the support of the transition is $\mathcal{Y}$. Therefore, it is possible for a “return" to any region of $\mathcal{Y}$ from any value $y_{t-1}$. This recurrence property is compatible with both stationary and nonstationary features of process $(y_t)$ (see the discussion in Section 3). We also make no assumptions about the observability of $(y_t)$. Thus, the propositions are valid for instance in a nonlinear state space model with latent (Markov) factors. Moreover, the assumption $\mathcal{Y}=\mathbb{R}^n$ is not restrictive. Indeed, the Markov property in Proposition 1 can be equivalently written on $y=(y_t)$, or a transform $\tilde{y}=(\tilde{y}_t)$ with $\tilde{y}_t = g^*(y_t)$, where $g^*(y_t)$ is a diffeomorphism from $\mathbb{R}^n$ to another set $\mathcal{Y}^*=g^*(\mathbb{R}^n)$. After this transformation, we get:

with the same shock $\varepsilon_t$. \\

The Multiplicity of Nonlinear Innovations

Let us now illustrate the multiplicity of Gaussian (resp. uniform) nonlinear innovations. It is well known that the linear transformation $D\varepsilon_t$ of the standard normal is also standard normal if and only if $DD'=Id$, that is, if $D$ is an orthogonal matrix. However, there exists also nonlinear transformations of $\varepsilon_t$ leaving invariant the standard normal distribution.

Preliminary Lemmas

We see in Proposition 2 the key roles of the Jacobian matrices and their determinants. To apply the Jacobian formula, the Jacobian matrices must be invertible and matrices with strictly positive determinants form a multiplicative group. This group can be decomposed by subgroups with interesting interpretations. In particular, we have the following subgroups:

Then any matrix with positive determinant can be decomposed by considering successive quotient groups. We have the polar decomposition [see Hall (2015), Section 2.5]:

This decomposition can be also written by considering the matrix logarithm transforms. We denote $\text{Exp} \ A = \sum_{k=0}^\infty \frac{A^k}{k!}$. Then, we have the following lemma:

Note that the “dimension" of the group generated by the different components are $1$ for $\lambda$, $\frac{n(n-1)}{2}$ for $Q$ (or $B$), $\frac{n(n+1)}{2}-1$ for $\Omega$ (or $V$). These dimensions sum up to $n^2$, that is, the dimension for the group of matrices $A$ with positive determinants. These decompositions could be applied to the transformations of nonlinear Gaussian innovations $\varepsilon$ as well as of nonlinear uniform innovations $u$.

Transformations in the Identified Set

The definition of the identified set depends on the selected normalization of the nonlinear innovations either Gaussian, or uniform, say. In this section, we denote $u_t$, $v_t$, some potential uniform innovations and $\varepsilon_t$, $\eta_t$, potential Gaussian innovations, respectively. To make a link between these normalizations, we write:

and similarily:

To partly parameterize the identified set, it is easier to work with the Gaussian innovations and their “polar" coordinates. Indeed, the support restrictions are much weaker than with the uniform distribution. More precisely, the vector $\varepsilon_t$ can be reparameterized as $||\varepsilon_t||^2 = \rho_t^2$, and $\frac{\varepsilon_t}{||\varepsilon_t||} = \zeta_t$. The first component $\rho^2_t$ is positive and the second component is an element of the $n-1$ hypersphere $S^{n-1}$. It is well-known that the standard Gaussian distribution for $\varepsilon_t$ becomes in polar coordinates $(\rho^2,\zeta_t)$ a continuous distribution with respect to the joint Lebesgue measure on $\mathbb{R}^+ \times$ the uniform probability measure on the sphere, with a joint density function of the type $\ell(\rho)$. Therefore, in polar coordinates any local transformation that keeps $\rho$ invariant and is orthogonal in $\zeta$ given $\rho$ will lead to the same distribution. This can be rewritten in terms of $\varepsilon_t$, $\eta_t$.

These transformations are invertible with $\varepsilon_t = Q'(||\varepsilon_t||^2)\eta_t=\text{Exp} \left[-B(||\varepsilon||^2)\right] \eta_t$ and have their analogue in the uniform normalization using (ref).

By construction, the transformations in Proposition 3 and Corollary 3 are nonlinear and their determinant is equal to 1. \\

The transformations above are easily understood in dimension 2, where it is easier to work in polar coordinates: $(\rho, \theta)$, instead of initial coordinates $(\varepsilon_{1,t},\varepsilon_{2,t})$. Then, the density of the standard normal distribution becomes $f(\rho,\theta)=\frac{1}{2\pi}\rho\exp(-\rho^2/2)$, on the domain $(0,\infty)\times(0,2\pi)$ (up to modulo $2\pi$ for $\theta$). Let us now consider the transformation $(\rho,\theta) \rightarrow (\rho^*,\theta^*) = (\rho,\theta+a(\rho))$, where $a(\cdot)$ is a continuous differentiable function of $\rho$. Since the Jacobian is $\left|\det\frac{\partial(\rho^*,\theta^*)}{\partial(\rho,\theta)'}\right|=\left|\det

\right|=1$, we see that the density of $(\rho^*,\theta^*)$ is also of the form $\frac{1}{2\pi}\rho^2 \exp(-\rho^2/2)$. For each fixed distance to the origin, the transformation is a rotation. When the rotation is with an angle $a$ independent of $\rho$, we get the standard linear orthogonal transform. But in general, we can make the rotation depending on the distance to the origin. In general we have [see Appendix B.3 for details]:

To illustrate these transformations and their effect on the uniform normalization, we consider the bidimensional case and the following transformations $\theta+a(\rho)$ on angle $\theta$:

For each case, we represent in $

\in [0,1]^2$, the curves transformed of the segments $(u_1=\frac{k}{75}, 0\leq u_2 <1), k=0,1,...,75$ (resp. $0 \leq u_1 \leq 1, u_2 = \frac{k}{75}), k=0,1,...,75$. \\

They correspond to homotopic deformations of the square [Smale (1959)]. These transformations can be equivalently represented for the Gaussian normalization. They lead to the corresponding (truncated) diffeomorphisms in $\mathbb{R}^2$ which are plotted together in Figure 1 above. \\

The set of transformations above shows the high degree of underidentification of the set of nonlinear innovations. Do alternative transformations exist, other than local orthogonal transforms that preserve the standard Gaussian? A result in this direction can be given in the bivariate case $n=2$, under the additional condition of zero preserving, known as the Schwarz Lemma, where the identified set is a set of functions. Without loss of generality, we can restrict the functions on $[0,1]^2$ to the harmonic functions that admit a series expansion as:

and analyze the identification issues by means of the coefficients $(b_{1,h,k}, b_{2,h,k}, h,k=0,1,...)$. We get the following result [see Gourieroux and Jasiak (2022)]:

Proof: See Appendix A.4. \\

Therefore, the degree of underidentification is rather large and the identified set contains transformations other than the nonlinear diffeomorphisms given as examples above.

Identifying Restrictions and Identifiable Functions

When the dimension $n$ is larger than 2, the nonlinear AR(1) process, that is function $g$ in (ref), is not identifiable. This is a problem of nonparametric identification with the identified set described in Proposition 2. As a consequence, the nonlinear innovations and the impulse response functions are not identifiable either. There are different ways to solve this identification issue and make shocks of interest more identifiable. [1] We can introduce parametric or nonparametric identifying restrictions making function $g$ identifiable. For instance, zero or sign restrictions can be considered to identify or partially identify a specific type of shock, such as a monetary policy shock [Uhlig (2005)], a fiscal policy shock [Ramey (2016), Section 4], a technology shock, or an oil shock [Hamilton (2003), Kilian and Murphy (2012)]. [2] We can use narrative records to identify shocks, as suggested in Romer and Romer (1997) [see also Antolin-Diaz and Rubio-Ramirez (2018)]. [3] Even if some partial identification is still present, we can look for transforms of the function $g$ (resp. the IRF) that are identifiable, even if they themselves are not. For instance, in the standard linear VAR framework, the horizonal Forecast Error Variance Decomposition (FEVD) is identifiable [see Gorodinchenko and Lee (2019) for the derivation of the FEVD by local projection and Gourieroux and Lee (2025) for extensions of the FEVD in the nonlinear dynamic framework]. We first briefly review the example of the Gaussian VAR(1) model, and then models with i.i.d. non-Gaussian error terms when the observations are preliminary conditionally demeaned and reduced. Then we explain how learning algorithms, simulations and pseudo IRFs can be used in an underidentified framework.

Identifiable Function of the IRF

Let us consider a SVAR(1) model and focus on the maximum possible IRF under a normalization on the magnitude of shock $\delta$. We consider the IRF on a linear combination of the variable $a'y$, say, that becomes $IRF(h,\delta,a)=a'A^hD\delta$, and on the solution of the constrained optimization problem:

The Lagrangean $IRF(h,\delta,a)-\lambda \delta'\delta$ is optimized for:

where $\lambda$ denotes the Langrange multiplier and the normalization condition implies:

We deduce the value of the objective function at its optimum.

This approach can easily be extended to a specific type of shocks with sign restrictions on the IRF [see e.g. Plagborg-Moller and Wolf (2021), Example 3]. \\

An alternative is to introduce identifying restrictions linking the autoregressive parameter $A$ and the mixing matrix $D$. This is usually done on a “structural form" of the model:

or equivalently:

which is a VAR model with $A=\Phi^{-1}_0\Phi_1$, $D=\Phi_0^{-1}\text{diag}\sigma$. These can be exclusion (i.e. sparsity) restrictions with zero components in $\Phi_0$,$\Phi_1$, sign restrictions on IRF leading to partial identification only [Uhlig (2005), Antolin-Diaz and Rubio-Ramirez (2018)], or long run restrictions when technology shocks are defined as the only shocks that affect labour productivity in the long-run [Blanchard and Quah (1989)]. \\

Conditionally Demeaned and Reduced Non Gaussian Models with Independence Restrictions

They are written as:

where the components of $w_t$ are assumed to be independent. The specification (ref) does not assume that the errors $w_t$ are either uniformly, or normally distributed. However, we can easily recover the nonlinear autoregressive specification with uniform innovations. More precisely, model (ref) can be rewritten as:

where $u_t$ follows $U[0,1]^n$ and $F_i$ denotes the c.d.f. of $w_{it}$, for $i=1,...,n$. The representation (ref) shows the nonparametric restrictions introduced on the general autoregressive model (ref). The function $g(y_{t-1},u_t)$ is assumed linear affine in the specific transformations $F_i^{-1}$ of $u_{it}$, with coefficient functions of $y_{t-1}$. These restrictions allow for nonlinearities both in $y_{t-1}$ and $u_t$, but also on cross-effects. Nevertheless, these restrictions imply a significant reduction of the dimensions in terms of functions. In model (ref), function $g$ is from $\mathbb{R}^{2n}$ to $\mathbb{R}^n$ where $m(\cdot)$, $D(\cdot)$ depend on an argument of dimension $n$, and the $F_i^{-1}$'s depend on a one-dimensional argument. \\

Model (ref) includes the strong linear VAR(1) model where $m(y_{t-1})=Ay_{t-1}$ and $D(y_{t-1})=Id$. Let us assume that $m(y_{t-1})$ and $D(y_{t-1})$ are identified. Then we know asymptotically $ D^*w_t$. The problem of identification becomes a problem of multivariate deconvolution, also known as Independent Component Analysis (ICA). Can we identify $D^*$ and the distributions of $w_{1,t},...,w_{n,t}$ from the joint distribution of $D^*w_t$? The answer is given in the following proposition:

Proof: See Comon (1994), Erikson and Koivunen (2004), Gouri\'eroux, Monfort and Renne (2017), Velasco (2022). \\

Thus, the linearity in $w_t$ in (ref) and the independence assumptions make the model structural and facilitate the identification. In general the identified matrix $D^*$ is not triangular as assumed in the recursive definition of shocks by Sims (1980). Thus, the primitive shocks are not associated with an equation, but may enter in several equations, or in a combination of such equations. Proposition 6 shows also that, in such a model (ref), the degree of partial identification will depend on the number of Gaussian sources [See Guay (2021), Gouri\'eroux and Jasiak (2025)]. \\

Identifiable Parameters

As seen in Section 2, the distribution of the Markov process can be equivalently defined by means of its transition $\ell(y_{t}|y_{t-1})$, or by means of a nonlinear autoregressive representation. Usually the transition density is identifiable from a sequence of observations $y_t$, $t=1,...,T$, for instance if the process is stationary ergodic. Therefore, any function of this transition density is also identifiable. However, even under additional restrictions, such as a parametric form $\ell(y_t|y_{t-1};\theta)$, the nonlinear autoregressive model (that is $\theta$) is not always identifiable. Let us now discuss in detail the practice that consists of estimating nonparametrically the function $g$ (resp. the parameter $\theta$) by some learning algorithms, and the use of these estimates in the generative model (ref) for simulations, predictions, or computations of the IRF where the identification issue still exists.

Learning Algorithms

Learning algorithms to approximate the function $g$ (resp. the parameter $\theta$) are largely available. We focus on algorithms that require a starting value and then provide a well defined sequence of approximations $\hat{g}_{p,T}$ (resp. $\hat{\theta}_{p,T}$), where $p$ denotes the iteration step. They include gradient algorithms [Barzilai and Borwein (1988)], Average Stochastic Gradient Descent (ASGD) algorithms [Polyak and Juditsky (1992)], some Alternating Optimization algorithms and various algorithms used for ICA [Almeida (2003), Dinh, Krueger and Bengio (2015)]. Since we are in a framework, where $g$ (resp $\theta$) is not identifiable, we can just expect that some of these algorithms will provide for large $T$ and stopping rule $p_T$ an approximation close to the identified set. The approximation $\hat{g}_T$, say, depends in general on the selected starting value. It is out of the scope of this paper to derive theoretically the properties of these algorithms. We make the following assumptions:

Can we use such an estimate in practice and for which purpose?

Simulations

Let us first consider the simulation issues. Indeed, simulated trajectories are usually based on the nonlinear autoregressive specification in (ref). If $g_0$ and $g_1$ are two functions in the identified set, then we can simulate the series $y_{0t}^s$ and $y_{1t}^s$ from a same set of simulated errors $\varepsilon_t^s$, for $t=1,...,T$, with the same starting value $y_0$ by applying recursively:

The two simulated trajectories will significantly differ since they are based on different transformations $g_0$ and $g_1$. However, since $g_0,g_1 \ \in \ IS$, their historical distributions will be the same for large $T$, and their simulated analogues based on a consistent $\hat{g}_T$, i.e.:

will have the same asymptotic properties. To summarize, we can estimate consistently by simulations based on $\hat{g}_{p_T,T}$ any (functional) summary of the true distribution of the process, such as for instance the predictive distribution of $y_{t+h}$ given $y_t=y$. But the speed of convergence and asymptotic distribution of these estimators are not known in general.

Pseudo IRF

There exists a debate in Econometrics and Time Series on how to define shocks. Have they to be defined on an innovation $\varepsilon_t$, or directly on the series of interest $y_t$? [see e.g. Gallant et al. (1993) for this debate]. We have retained in Section 2 the definition based on innovations due to its more reasonable interpretation. Nevertheless, let us now change $y_t$ into $y_t + \Delta$, say, and let $\varepsilon_t$ be unchanged. The system (ref) - (ref) becomes:

where $y_t^\Delta=y_t + \Delta$ and $y^\Delta_{t-1} = y_{t-1}$. This leads to the definition of the so-called Pseudo IRF (PIRF):

Then we get the following result:

Therefore, the identification issue on the IRF is just due to the different effects at horizon 1 when applying a change on $y_t$ instead of a change on $\varepsilon_t$. The PIRF is in general stochastic and its expectation conditional on $y_t$ is equal to:

This is a function of the transition density $\ell(y_t|y_{t-1})$ that can easily be consistently approximated by simulations as in Section 4.2.2. \\

Remark 4: The difference between the IRF and the Pseudo IRF is easily seen in the example of a structural Gaussian VAR model:

where $\varepsilon_t \sim IIN(0,Id)$. We have:

The autoregressive coefficient $\Phi$ is identifiable, but not the rotation matrix $D$. This example is rather misleading since it gives the impression of similar term structures between the IRF and PIRF, independent of their environments. Proposition 6 shows that the link between them is much more complex in a nonlinear dynamic framework.

Nonlinear Dynamic Factor Models

The identification issues could also be solved when the multivariate observed time series is a transformation of independent dynamic factors (or sources in the operation research literature). This is the identification issue in dynamic nonlinear ICA, also called Blind Source Separation (BSS) to emphasize the difference with the i.i.d. static case [Jutten and Herault (1991), Hyvarinen and Morioka (2017)]. We consider below the bivariate case for expository purpose with linear and nonlinear transformations of the sources, respectively.

Linear Measurement Equation

Let us assume that the bivariate observable process $(y_t)$ is such that:

where the sources $(x_{1,t})$ and $(x_{2,t})$ are independent Markov processes:

$(\varepsilon_{1,t})$, $(\varepsilon_{2,t})$ are independent Gaussian noises, and matrix $A$ is invertible. The nonparametric identification of $A$, $g_1$, $g_2$, is equivalent to the nonparametric identification of matrix $A$. Note that we get the same type of dynamic model if $(x_{1,t})$ (or $(x_{2,t})$) is multiplied by a (signed) scale factor, and also if we permute the one-dimensional processes $(x_{1,t})$ and $(x_{2,t})$. This leads to the following definition: \\

The lack of identification due to scale effects can be solved by introducing an appropriate identification restriction as: “the diagonal elements of matrix $A$ are equal to 1". Then, we have the following result:

Proof: See Appendix A.5. \\

We have seen in Proposition 6 that there is an identification issue when $(x_{1,t})$ and $(x_{2,t})$ are independent Gaussian white noises. This lack of identification disappears, if the underlying independent factors (sources) have different second-order dynamics. This result can be extended to any dimension where the observations are combinations of independent AR(1) processes with different autoregressive coefficients in [Gouri\'eroux and Jasiak (2023b)]. Autocovariance based methods to separate sources in the linear framework have been proposed in Tong et al. (1991) and Belouchani et al. (1997). The condition in Proposition 8 is sufficient, but not necessary. Weaker identification conditions can sometimes be obtained by considering zero autocovariance restrictions or nonlinear transformations of $x_{1,t}$ (resp. $x_{2,t}$) [Gouri\'eroux and Jasiak (2022b)]. The idea is the following. Let us assume that the dynamics of the latent factors are parameterized, that is:

where $\theta_1$, $\theta_2$ are parameters. This system can be inverted as:

Let us denote $C=A^{-1}$ as the demixing matrix. Then, the model can be rewritten as:

with $\theta=(\text{vec} C,\theta^*)$. Since the $\varepsilon_t's$ are serially independent, it has been suggested in the literature to try to identify the parameters $C$ and $\theta$ by means of a set of nonlinear autocovariance restrictions of the form:

with a selected set of pairs of nonlinear transformations $a,\tilde{a}$ and lags. Such an identification approach is for instance introduced for the identification of mixed causal-noncausal models when using the Generalized Covariance (GCov) estimation approach.

Nonlinear Measurement Equation

The model in the previous subsection can be extended to account for nonlinear transformations of the sources, where (ref) is replaced by:

Then, the definition of essential identification has to be extended to account for the nonlinearity of the transformation $A$. Typically, $A$, $g_1$ and $g_2$ can be modified for the latent variables $(x_{1,t})$ and $(x_{2,t})$ to be marginally uniform on $[0,1]$, for $(y_t)$ to be marginally uniform on $[0,1]^2$, and up to permutation of indices. \\

By fixing the margins, the analysis of identification will be through pairwise copulas at different lags instead of covariances at different lags as in Proposition 8. We have:

Proof: See Appendix A.6. \\

Therefore, there is no identification issue if the sources are independent with different dynamics, even in the Gaussian framework. As an illustration, let us consider Gaussian AR(1) sources: $x_{j,t}=\rho_jx_{j,t-1}+\sigma_j\varepsilon_{j,t}$. The conditional transition at horizon $h$ is Gaussian $N(\rho^hx_{j,t},\eta^2_{j,h})$, where $\eta^2_{j,h}=\sigma_j^2\frac{1-\rho^{2h}_j}{1-\rho^2_j}$. Then we have:

and the condition of Proposition 8 is satisfied if and only if $\rho_1 \neq \rho_2$.

Impulse Response Functions

Under the condition of Proposition 8, $\varepsilon_t =

$ is a Gaussian nonlinear innovation defined in a unique way. Therefore, the associated IRF for $y_t$ are also identifiable. They are easily deduced from the IRF performed separately on $x_{1,t}$ and $x_{2,t}$. More precisely we define the shocked and unshocked trajectories of the sources with respective shocks $\delta_1$ and $\delta_2$. These are $x_{1,t+h}$, $x_{2,t+h}$, $x^{(\delta_1)}_{1,t+h}$ and $x^{(\delta_2)}_{2,t+h}$. Then the shocked and unshocked trajectories of $y_t$ are:

Due to the nonlinear transformation $A$, the $IRF(h,\delta,y_{t-1})$ is not a function of the IRF's of the two sources only. In other words, in a nonlinear dynamic framework, the notion of IRF is not coherent with respect to nonlinear aggregation. This approach is using implicitly the assumption of equal dimensions of the processes $(y_t)$ and $(x_t)$. If $(x_t)$ were with a dimension strictly larger than the dimension of $(y_t)$, the information in $\underline{y_t}$ and $\underline{x_t}$ would not be the same.

Augmented Nonlinear Autoregressive Models

It is known that the notion of causal inference, of Granger causality, and likely also the notion of Markov process, nonlinear innovation and IRF can depend on the universe, that is on the selected set of variables. For instance, narrative shocks are often found to be predictable, suggesting the possibility of endogeneity [Leeper (1997), Romer and Romer (1997); see also Ganies et al. (2019) for an IV adjustment for endogeneity]. In particular, starting from a given universe, they could be modified if the universe is either decreased, or increased. This leads for instance to the introduction of the Factor Augmented Vector Autoregressive (FAVAR) models in the linear dynamic VAR literature [Bernanke, Boivin and Eliasz (2002)]\footnote{These augmented vector autoregressive models concern the true dynamic of the variable. They have to be distinguished from the lag augmented local projections, that are instrumental models introduced for robust inference on standard errors of the IRF [Montiel-Olea and Plagborg-Moller (2021), (2022)]}.

The Markov Assumption

In fact there is a trade-off between the dimension of the system and the order of the Markov process\footnote{Assumed equal to 1 in our paper, but the discussion is the same for any fixed order $p$ [see e.g. Florens et al. (1993)].}. For instance, it is well known that a one-dimensional AR(2) process can be written as a bidimensional VAR(1) process. It is also known that the first component of a bivariate VAR(1) process satisfies an ARMA model (i.e. an AR($\infty$) model), when its dynamic is marginalized. This highlights the advantage of multivariate modelling to avoid the increased autoregressive order. \\

To understand why the Markov assumption depends on the universe, let us consider a (multivariate) process $(y_t)$, and the process augmented by another (multivariate) process ($\zeta_t$). We denote $\underline{y}_t$ (resp. $\underline{\zeta}_t$) the set of current and past values of the process $(y_t)$(resp. $(\zeta_t)$). There are three notions of Markov processes with conditions written in terms of transition densities. \\

(i) The augmented process $(y_t,\zeta_t)$ is Markov with respect to the augmented universe: \\

(ii) The initial process $(y_t)$ is Markov with respect to its own universe: \\

(iii) The initial process $(y_t)$ is Markov with respect to the augmented universe: \\

All these conditions can be rewritten in terms of conditional independence, denoted $\boldsymbol{\cdot} \perp \!\!\! \perp \boldsymbol{\cdot} \ | \ \boldsymbol{\cdot}$ [Florens et al. (1993)] as:

(i) $(y_t,\zeta_t)$ is Markov in the augmented universe, iff $(y_t,\zeta_t) \perp \!\!\! \perp (\underline{y}_{t-2},\underline{\zeta}_{t-2}) | \ y_{t-1},\zeta_{t-1}$. \\ (ii) $(y_t)$ is Markov in its universe, iff $y_t \perp \!\!\! \perp \underline{y}_{t-2} | \ y_{t-1}$. \\ (iii) $(y_t)$ is Markov in its augmented universe, iff $y_t \perp \!\!\! \perp (\underline{y}_{t-2},\underline{\zeta}_{t-1}) | \ y_{t-1}$. \\

It is easily seen that, under (ref), the condition (ref) means that process $(\zeta_t)$ does not (Granger) cause process $(y_t)$, that is a condition of dynamic exogeneity of process $(y_t)$. \\

In the previous sections, to define the nonlinear innovations, to define the IRF corresponding to the shocks on their innovation and to discuss the identification issues, we have assumed that the nonlinear VAR(1) model is well-specified, that is the Markov assumption is satisfied. This predictive assumption has now to be discussed. In particular, we have to consider the following question - in which case do we have simultaneously the joint and marginal Markov properties (ref),(ref) satisfied?\\

We have the following result:

This noncausality condition is sufficient, but not necessary [see Florens et al. (1993) for more detailed analysis]. This type of condition will also appear when the IRFs are obtained by “state dependent" local projection [Goncalves et al. (2022), (2023)].

Testing the Markov Assumption

The discussion in Section 5.3.1 shows that the Markov assumption is at the core of the definitions of IRFs and their identification. Therefore, in practice, this assumption would have to be tested (which is rarely the case, even in linear dynamic VARs). In the identified case, it is equivalent to test that the nonlinear innovations are independent white noises. This can be done by applying tests of independence and of strong white noise to the series of nonlinear residuals. These tests have to check for the absence of any type of cross-sectional and serial nonlinear correlation (not only linear ones). When there is an identification issue, tests of conditional independence can be used based on the interpretation of the Markov assumption in terms of conditional independence. This condition can be written in terms of appropriate moments. For instance, the Markov condition [ii] is equivalent to the covariance conditions:

for any integrable functions $a,b,c$. Such conditions can be checked by using pormanteau test statistics based on the sample covariance counterparts [Gourieroux and Jasiak (2023a),(2025), Velasco (2023)].

Concluding Remarks

This paper has extended the notions of nonlinear autoregressive representation, nonlinear innovations and impulse response functions to the multivariate nonlinear dynamic framework. We analyzed the identification issue for these innovations and IRFs. In this respect, the linear Gaussian VAR models and the unconstrained Markov model appear as special cases in which the identification issues are especially significant. The main identification issues can disappear when nonlinearities and/or independent latent sources with different dynamics exist. Such analysis can be applied to any nonlinear dynamic system. Can we use them as tools for economic policy? Following Bernanke (1986), p. 52-55, “Structures of shocks are primitive exogenous forces (i.e. the innovations), that are independent of each other and economically meaningful". Therefore, in practice, the estimation of the model, of the innovation, of the IRFs have to be followed by a structural step to interpret economically some of these innovations, not necessarily all of them. All the analyses in this paper have been developed for time series with multivariate observations indexed by time $t$. It would be interesting to extend this analysis to series that are doubly indexed: [1] By time and maturity when considering the term structure of interest rates or of volatilities, [2] indexed by the localization for spatial processes, or [3] indexed by a pair of individuals as in networks. These extensions are left for further research.

References

\setstretch{1}

Almeida, L. (2003). MISEP—Linear and Nonlinear ICA Based on Mutual Information. The Journal of Machine Learning Research, 4, 1297-1318.\\

Antolín-Díaz, J., and J., Rubio-Ramírez. (2018). Narrative Sign Restrictions for SVARs. American Economic Review, 108, 2802-2829.\\

Axler, S., Bourdon, P. and W., Ramey (2001): “Harmonic Function Theory", 2nd Ed., Graduate Texts in Mathematics, Springer.\\

Barzilai, J., and M., Borwein (1988). Two-Point Step Size Gradient Methods. IMA Journal of Mathematical Analysis, 8, 141-148. \\

Belouchani, A., Meraim, K., Cardoso, S., and E., Moulines (1997): A Blind Source Separation Technique Based on Second-Order Statistics", IEEE Trans on Signal Processing, 45, 434-444. \\

Bernanke, B. (1986). Alternative Explanations of the Money-Income Correlation. Carnegie Rochester Conference Series on Public Policy. 25, 49-99. \\

Bernanke, B., Boivin, J., and P., Eliasz. (2002). Measuring the Effects of Monetary Policy: A Factor Augmented Vector Autoregressive (FAVAR) Approach. Quarterly Journal of Economics, 120, 387-422. \\

Blanchard, O., and D., Quah. (1989). The Dynamic Effects of Aggregate Demand and Supply Disturbances. American Economic Review, 73, 655-673. \\

Borkovec, M., and C., Kluppelberg. (2001). The Tail of the Stationary Distribution of an Autoregressive Process with ARCH(1) Errors. Annals of Applied Probability, 11, 1220-1241. \\

Christiano, L. (2012). Christopher A. Sims and Vector Autoregressions. The Scandinavian Journal of Economics, 114, 1082-1104. \\

Comon, P. (1994). Independent Component Analysis, a New Concept?. Signal Processing, 36, 287-314. \\

Cox, J., Ingersoll, J., and S., Ross. (2005). A Theory of the Term Structure of Interest Rates. Theory of Valuation,129-164. \\

Darmois, G. (1953). Analyse Generale de Liaisons Stochastiques. Rev Inst Internat Stat, 21, 2-8. \\

Dinh, L., Krueger, D., and Y., Bengio (2015): “NICE: Nonlinear Independent Components Estimation", ICLR. \\

Eriksson, J., and V., Koivunen. (2004). Identifiability, Separability, and Uniqueness of Linear ICA models. IEEE Signal Processing Letters, 11, 601-604. \\

Florens, J.P., Mouchart, M. and J.M., Rolin. (1993). Noncausality and Marginalization of Markov Processes. Econometric Theory, 9, 239-260. \\

Gallant, A., Rossi, P., and G., Tauchen. (1993). Nonlinear Dynamic Structures. Econometrica, 61, 871-907.\\

Ganics, G., Inoue, A., and B., Rossi. (2021). Confidence Intervals for Bias and Size Distortion in IV and Local Projections-IV Models. Journal of Business & Economic Statistics, 39(1), 307-324. \\

Goncalves, S., Herrera, A., Kilian, L., and E., Pesavento. (2021). Impulse Response Analysis for Structural Dynamic Models with Nonlinear Regressors. Journal of Econometrics, 225, 107-130. \\

Goncalves, S., Herrera, A., Kilian, L., and E., Pesavento. (2022). When Do State Dependent Local Projections Work?, DP McGill University. \\

Gorodnichenko, Y., and B., Lee. (2020). Forecast Error Variance Decompositions with Local Projections. Journal of Business and Economic Statistics, 38(4), 921-933.\\

Gouri\'eroux, C., and J., Jasiak. (2005). Nonlinear Innovations and Impulse Responses with Application to VaR Sensitivity. Annales d'Economie et de Statistique, 78, 1-31. \\

Gouri\'eroux, C., and J., Jasiak. (2022). Nonlinear Forecasts and Impulse Responses for Causal-Noncausal (S)VAR Models. arXiv preprint arXiv:2205.09922. \\

Gouri\'eroux, C., and J., Jasiak. (2023a). Generalized Covariance Estimator. Journal of Business and Economic Statistics, 41, 1315-1327. \\

Gouri\'eroux, C., and J., Jasiak. (2023b). Dynamic Deconvolution of Independent Autoregressive Sources. Journal of Time Series Analysis, 44, 151-180. \\

Gouri\'eroux, C., and J. Jasiak. (2025). Generalized Covariance Based Inference for Models Partially Identified from Independence Restrictions, Journal of Time Series Analysis, 46, 300-324. \\

Gouriéroux, C., Jasiak, J., and R., Sufana. (2009). The Wishart Autoregressive Process of Multivariate Stochastic Volatility. Journal of Econometrics, 150(2), 167-181.\\

Gouri\'eroux, C., and Q., Lee. (2025). Forecast Relative Error Decompositions with Application to Cyber Risk. arXiv:2406.17708.\\

Gouri\'eroux, C., and A., Monfort. (1997). Simulation Based Econometric Methods, Oxford University Press. \\

Gouri\'eroux, C., Monfort, A., and J.P., Renne. (2017). Statistical Inference for Independent Component Analysis: Application to Structural VAR Models. Journal of Econometrics, 196, 111-126. \\

Guay, A. (2021). Identification of Structural Vector Autoregressions Through Higher Unconditional Moments. Journal of Econometrics, 225, 27-46.\\

Hall, B. (2015). Lie Groups, Lie Algebras and Representations: An Elementary Introduction", Graduate Texts in Mathematics, 222, (2nd Edition), Springer. \\

Hamilton, J. (2003). What is an Oil Shock?. Journal of Econometrics, 113, 363-398. \\

Hyvarinen, A., and H., Morioka. (2017). Nonlinear ICA of Temporally Dependent Stationary Sources. iin Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, 54, 460-469. \\

Hyvarinen, A., and P., Pajunen. (1999). Nonlinear Independent Component Analysis: Existence and Uniqueness Results. Neural Networks, 12, 429-439.\\

Hyvärinen, A., and E., Oja (2001). Independent Component Analysis: Algorithms and Applications. Neural Networks, 13(4-5), 411-430.\\

Hyvarinen, A., Sasaki, H., and R., Turner (2019). Nonlinear ICA Using Auxiliary Variables and Generalized Contrastive Learning. in Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, 89, 859 -868. \\

Isakin, M., and P., Ngo. (2020). Variance Decomposition Analysis for Nonlinear Economic Models. Oxford Bulletin of Economics and Statistics, 82, 1362-1374. \\

Jutten, C., and J., Herault. (1991). Blind Separation of Sources Part 1: An Adaptive Algorithm Based on Neuromimetic Architecture. Signal Processing. 94. 1-10. \\

Karlin, S., and H., Taylor. (1981). A Second Course in Stochastic Processes. Elsevier. \\

Karlsen, M. and D., Tjostheim. (2001). Nonparametric Estimation on Null Recurrent Time Series. Annals of Statistics, 29, 372-416. \\

Kilian, L., and D., Murphy. (2012). Why Agnostic Sign Restrictions are Not Enough: Understanding the Dynamics of Oil Market VAR Models. Journal of the European Economic Association, 10, 1166-1188.\\

Koop, G., Pesaran, H., and S., Potter. (1996). Impulse Response Analysis in Nonlinear Multivariate Models. Journal of Econometrics, 74, 119-147. \\

Leeper, E. (1997). Narrative and VAR Approaches to Monetary Policy: Common Identification Problems. Journal of Monetary Economics, 40, 641-657. \\

Ling, S. (2007). A Double AR(p) Model: Structure and Estimation. Statistica Sinica, 17, 161-175. \\

Montiel Olea, J., and M., Plagborg‐Møller. (2021). Local Projection Inference is Simpler and More Robust Than You Think. Econometrica, 89, 1789-1823. \\

Montiel Olea, J., and M., Plagborg‐Møller. (2022). Corrigendum: Local Projection Inference is Simpler and More Robust Than You Think. Online Manuscript. \\

Nakamura, E., and J., Steinsson. (2018). Identification in Macroeconomics. Journal of Economic Perspectives, 32, 59-86.\\

Plagborg‐Møller, M., and C., Wolf. (2021). Local Projections and VARs Estimate the Same Impulse Responses. Econometrica, 89, 955-980. \\

Polyak, B., and A., Juditsky (1992). Acceleration of Stochastic Approximation by Averaging. SIAM Journal on Control and Optimization, 30(4), 838-855.\\

Ramey, V. (2016). Macroeconomic Shocks and their Propagation. Handbook of Macroeconomics, 2, 71-162. \\

Roberts, S., and R., Everson (2001). Independent Component Analysis: Principles and Practice. Cambridge University Press. \\

Romer, C., and D., Romer. (1994). Monetary Policy Matters. Journal of Monetary Economics, 34, 75-88.\\

Romer, C., and D., Romer. (1997). Identification and the Narrative Approach: A Reply to Leeper. Journal of Monetary Economics, 40, 659-665.\\

Rosenblatt, M. (1952). Remarks on a Multivariate Transformation. The Annals of Mathematical Statistics, 23, 470-472.\\

Rubio-Ramirez, J. , Waggoner, D., and T., Zha. (2010). Structural Vector Autoregressions: Theory of Identification and Algorithms for Inference. The Review of Economic Studies, 77, 665-696.\\

Sims, C. (1980). Macroeconomics and Reality. Econometrica, 48, 1-48. \\

Smale, S. (1959). Diffeomorphisms of the 2-Sphere. Proceedings of the American Mathematical Society, 10, 621-626. \\

Tong, L., Liu, R., Soon, V., and Y., Huang. (1991). Indeterminacy and Identifiability of Blind Identification. IEEE Transactions on Circuits and Systems, 38(5), 499-509.\\

Tweedie, R. (1975). Sufficient Conditions for Ergodicity and Recurrence of Markov Chains on a General State Space. Stochastic Processes and their Applications, 3, 385-403. \\

Uhlig, H. (2005). What are the Effects of Monetary Policy on Output? Results from an Agnostic Identification Procedure. Journal of Monetary Economics, 52, 381-419.\\

Velasco, C. (2023). Identification and Estimation of Structural VARMA Models Using Higher Order Dynamics. Journal of Business & Economic Statistics, 41, 815-832.\\

Weiss, A. (1984). ARMA Models with ARCH Errors. Journal of Time Series Analysis, 5, 129-143. \\

Zhu, H., Zhang, X., Liang, X., and Y., Li. (2017). On a Vector Double Autoregressive Model. Statistics & Probability Letters, 129, 86-95.\\