A fully nonlinear structural vector autoregressive model identified via independent innovation analysis

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

82,407 characters

\begin{titlepage}
\begin{center}
\linespread{1.2}
\Large{\textbf{A fully nonlinear structural vector autoregressive model identified via independent innovation analysis}}\\

\vspace{0.5cm}
\Large{Savi Virolainen}\\
\large{University of Helsinki}\\
\vspace{1.0cm}


\begin{abstract}
\noindent We develop a fully nonlinear structural vector autoregressive framework in which the contemporaneous structural mapping may be nonlinear and non-additive. Identification is achieved by exploiting variation in the conditional distributions of the mutually independent structural shocks induced by an observed exogenous variable. Specifically, a general contrastive learning framework that makes use of this variation together with the assumed exponential-family structure is employed to recover the shocks. Existing independent innovation analysis results identify such shocks only up to arbitrary componentwise invertible transformations, which is generally insufficient for structural econometric analysis. We strengthen this result by imposing a structured exponential-family specification for the conditional shock distributions. With the imposed sufficient statistics, the remaining ambiguity is reduced to a one-parameter transformed-scale map for each shock. We then show that, under a logistic specification used for the natural parameters, the identification is further strengthened up to permutation and componentwise sign changes. Once the shocks have been recovered, the fully nonlinear structural vector autoregression can be estimated using feed-forward neural networks, motivated by their universal approximation capabilities. The empirical application studies asymmetries in the responses of U.S. industrial production to the real oil price shock. We find modest asymmetries with respect to the sign of the shock and state of the economy. The accompanying R package iiasvar implements the introduced methods.

\noindent\textbf{Keywords:} Nonlinear SVAR, nonlinear structural vector autoregressive model, independent innovation analysis, feed-forward neural network, statistical identification\\%[2.0cm]
\end{abstract}

\vfill


  \begingroup
  \renewcommand{}\footnote{The author thanks the Research Council of Finland (Grant 347986) and the Yrjö Jahnsson Foundation (Grant 20247849) for the financial support.}
  \addtocounter{footnote}{-1}
  \endgroup


  \begingroup
  \renewcommand{}\footnote{Contact address: Savi Virolainen, Faculty of Social Sciences, University of Helsinki, P. O. Box 17, FI–00014 University of Helsinki, Finland; e-mail: [email removed]. ORCiD ID: 0000-0002-5075-6821.}
  \addtocounter{footnote}{-1}
  \endgroup


\end{center}
\end{titlepage}


\section{Introduction}
Structural vector autoregressive (SVAR) models play a central role in empirical macroeconomics by providing a framework for interpreting economic fluctuations through latent structural shocks and their dynamic propagation. In conventional SVAR analysis, identification is typically achieved by imposing restrictions on the contemporaneous relationship between reduced-form innovations and structural shocks. Much of the literature has focused on linear specifications, where the reduced-form innovations are assumed to be linear combinations of the structural shocks and identification may be achieved, for example, through exclusion restrictions, external instruments, heteroskedasticity, or non-Gaussianity \citep[see, e.g.,][]{Kilian+Lutkepohl:2017}.

While linear SVAR models are useful in many applications, they are not able to capture nonlinear features of economic dynamics. A growing literature has therefore considered nonlinear extensions of SVAR models. In nonlinear systems, the effects of shocks may depend on the state of the economy, the signs or sizes of shocks, or other nonlinear features of the transmission mechanism. More generally, the relationship between structural shocks and observables may itself be nonlinear, allowing for nonlinear interactions between shocks and lagged variables as well as non-additive transmission of the effects of the shocks.

At the same time, the recent theoretical work of \cite{Kolesar+Plagborg-Moller:2025} has demonstrated that statistical identification strategies developed under linearity may fail when the structural relationship between the reduced-form innovations and shocks is nonlinear, even if the shocks themselves satisfy otherwise standard assumptions. In particular, both identification by heteroskedasticity (\citealp{Rigobon:2003}, \citealp{Lanne+Lutkepohl:2010}, \citealp{Lutkepohl+Netsunajev:2017}, \citealp{Lewis:2021}, \citealp{Virolainen:2025}, and others) and identification by non-Gaussianity (\citealp{Lanne+Meitz+Saikkonen:2017}, \citealp{Lanne+Luoto:2021}, \citealp{Lanne+Liu+Luoto:2023}, and others) can fail in such a case. While \cite{Virolainen2:2024} has recently extended identification by non-Gaussianity to smooth-transition SVAR models, also his framework assumes that the mapping from structural shocks to reduced-form innovations is linear in each point of time, with nonlinearity arising only through time-variation in the parameter values. This motivates the development of identification frameworks specifically designed for fully nonlinear SVAR models.

This paper develops a fully nonlinear SVAR framework in which the shocks are identified statistically by imposing structure on the shock process itself and by utilizing independent innovation analysis based on general contrastive learning (IIA-GCL). Our framework builds on the independent innovation analysis approach of \citet{Morioka+Halva+Hyvarinen:2021}, who study identification of structural shocks in highly general nonlinear vector autoregressive models. Their IIA-GCL framework exploits conditional mutual independence and exponential-family structure together with variation in the conditional distributions of the structural shocks that is induced by an observed auxiliary variable. Under suitable conditions, this auxiliary variable modulation facilitates recovery of the structural shocks even when the contemporaneous structural mapping is fully nonlinear and shocks do not enter additively. However, the IIA-GCL framework of \citet{Morioka+Halva+Hyvarinen:2021} recovers the shocks only up to permutation and componentwise invertible transformations. From the perspective of nonlinear independent component analysis, this constitutes a strong identification result \citep[see, e.g.,][]{Hyvarinen+Pajunen:1999}. However, the remaining ambiguity is generally too large for structural econometric analysis, as economically meaningful shock interpretations are not preserved under arbitrary componentwise invertible transformations.

We strengthen the identification result of \citet{Morioka+Halva+Hyvarinen:2021} substantially by imposing a structured conditional shock distribution based on a particular exponential-family specification whose natural parameters vary systematically with the observed auxiliary variable. This additional structure reduces the componentwise nonlinear ambiguity to permutation, sign, and a one-parameter class of transformed-scale maps. We then show that the remaining transformed-scale ambiguity is generically eliminated under our logistic natural-parameter specification. The resulting framework thus identifies structural shocks up to permutation and componentwise sign changes in a fully nonlinear SVAR model, where the contemporaneous structural mapping itself need not be linear or additive.

Once the structural shocks have been recovered, estimation of the nonlinear SVAR model becomes a nonlinear regression problem. As proposed by \citet{Morioka+Halva+Hyvarinen:2021}, the model can be implemented in practice as a feed-forward neural network, motivated by their universal approximation capabilities \citep[see][]{Hornik+Stinchcombe+White:1989}. This flexible approach allows the dynamic relationships between the shocks and included variables to be learned from the data without having to specify any particular functional form for this relationship. Following the statistical identification literature \citep[see, e.g.,][]{Lanne+Meitz+Saikkonen:2017}, we propose labelling the identified shocks based on their estimated effects on the variables, here represented by the generalized impulse response functions \citep{Koop+Pesaran+Potter:1996}.

We illustrate the use of our methods in an empirical application studying asymmetries in the effects of the real oil price shock. Specifically, we consider a monthly bivariate system of U.S. industrial production growth and real oil price growth, covering the period from 1974:2 to 2026:1. Since the periods of elevated macroeconomic uncertainty can plausibly feature changes in the variances and the shapes of macroeconomic shock distributions, we assume that shock distributions are modulated by the (lagged) macro uncertainty index of \citet{Jurado+Ludvigson+Ng:2015}. We find modest asymmetries in the effects of the oil price shock with respect to its sign and the state of the economy. In particular, our findings suggest stronger effects for negative than positive oil price shocks, differing from the post-1973 (aggregate industrial production) result of \citet{Herrera+Lagalo+Wada:2011} for shocks of typical magnitude. Moreover, we find the effects of the oil price shock stronger when the recent industrial production growth is low than when it is high.

The remainder of the paper is organized as follows. Section~\ref{sec:model} introduces the nonlinear SVAR framework and the conditional shock distribution assumptions. Section~\ref{sec:identification} presents the identification result based on our modified IIA-GCL and summarizes the related Monte Carlo experiment. Section~\ref{sec:FFNN-SVAR} discusses the feed-forward neural-network implementation of the SVAR model and structural analysis. Section~\ref{sec:empirical} presents the empirical application, and Section~\ref{sec:conclusion} concludes. The appendices provide the proofs of the identification results, computational details, and further details on the empirical application. The accompanying R package iiasvar \citep{iiasvar} implements the introduced methods.


\section{The model}\label{sec:model}
\subsection{The nonlinear SVAR model}\label{sec:NVAR}
Consider the $d$-dimensional time series $\{y_t\}_{t\in\mathbb{N}}$ of interest, and the general nonlinear SVAR model of autoregressive order $p$:
\begin{equation}\label{eq:genvar}
y_t = f(\boldsymbol{y}_{t-1}, e_t),
\end{equation}
where $f:\mathbb{R}^{d(p+1)}\rightarrow \mathbb{R}^d$ represents the (unknown) functional form of the model, $\boldsymbol{y}_{t-1}=(y_{t-1},\ldots,y_{t-p})$, and $e_t=(e_{1t},\ldots,e_{dt})$ is a vector of unobserved structural shocks. In contrast to linear SVAR models, the function $f$ is allowed to be fully nonlinear in both past observations and the contemporaneous shocks. This general formulation accommodates a wide range of nonlinear dynamics, including state-dependent impulse responses and nonlinear propagation mechanisms.

Our main goal is to recover the structural shocks $e_t$ from the observed process $y_t$. As is well known, without further restrictions, such recovery is not possible due to the inherent non-uniqueness of nonlinear decompositions. The identification strategy developed in this paper builds on the independent innovation analysis (IIA) framework of \cite{Morioka+Halva+Hyvarinen:2021}, which exploits mutual independence together with variation in the shock distributions to resolve this ambiguity.


Following \cite{Morioka+Halva+Hyvarinen:2021}, it is convenient to consider the model in the following augmented form ("the mixing model"):
\begin{equation}\label{eq:ftilde}
\begin{bmatrix}
y_t \\
\boldsymbol{y}_{t-1}
\end{bmatrix}
= \tilde{f}\left(
\begin{bmatrix}
e_t \\
\boldsymbol{y}_{t-1}
\end{bmatrix}
\right)
=
\begin{bmatrix}
f(\boldsymbol{y}_{t-1},e_t) \\
\boldsymbol{y}_{t-1}
\end{bmatrix},
\end{equation}
where $\tilde{f}:\mathbb{R}^{d(p+1)}\rightarrow \mathbb{R}^{d(p+1)}$ is the augmented model, which incorporates~\eqref{eq:genvar} in the first $d$ entries and identity mapping for $\boldsymbol{y}_{t-1}$ in the remaining $dp$ entries. It is assumed that the augmented model is invertible (i.e., bijective) and twice continuously differentiable. In particular, the inverse of the augmented model $\tilde{f}$ ("the demixing model") can be expressed as:
\begin{equation}\label{eq:gtilde}
\begin{bmatrix}
e_t \\
\boldsymbol{y}_{t-1}
\end{bmatrix}
=
\tilde{g}\left(
\begin{bmatrix}
y_t \\
\boldsymbol{y}_{t-1}
\end{bmatrix}
\right)
=
\begin{bmatrix}
g(y_t,\boldsymbol{y}_{t-1}) \\
\boldsymbol{y}_{t-1}
\end{bmatrix},
\end{equation}
where $\tilde{g}:\mathbb{R}^{d(p+1)}\rightarrow \mathbb{R}^{d(p+1)}$ is the augmented demixing model of the (true) augmented model $\tilde{f}$, and $g(y_t,\boldsymbol{y}_{t-1})=(g_1(y_t,\boldsymbol{y}_{t-1}),\ldots,g_d(y_t,\boldsymbol{y}_{t-1}))\in\mathbb{R}^d$ represents a mapping from temporally consecutive observations to the shocks at each time period.


\subsection{Assumptions about the shocks}\label{sec:shockassumptions}

As described by \cite{Morioka+Halva+Hyvarinen:2021}, identification of the shocks requires assumptions on how their distributions vary over time. In the present framework, this variation is governed by an exogenous auxiliary variable $u_t$, which may be interpreted as describing the state of the shock process. The key idea is that, while mutual independence separates the shocks cross-sectionally, time variation in their conditional distributions provides the additional information needed for identification.

These requirements are formalized in the following assumption, which is close in spirit to Assumption~A1 of \cite{Morioka+Halva+Hyvarinen:2021}, but is stated here in a form adapted to the nonlinear SVAR setting. In particular, we make explicit that the auxiliary variable $u_t$ is observed and exogenous, is not included in the lag vector $\boldsymbol y_{t-1}$, and that the structural shocks are conditionally independent of $\boldsymbol y_{t-1}$ given $u_t$.
\begin{assumption}\label{as:shocks}
The following conditions hold:
\begin{enumerate}[label=(\roman*)]
    \item For each $i=1,\ldots,d$, the distribution of the shock $e_{it}$ is modulated by the observed exogenous auxiliary variable $u_t$, which is not included in the lag vector $\boldsymbol{y}_{t-1}$.

    \item Conditional on $u_t$, the shock vector $e_t$ is independent of the lag vector $\boldsymbol{y}_{t-1}$.

    \item For each $t$, the shocks $e_{1t},\ldots,e_{dt}$ are mutually conditionally independent given $u_t$.

    \item Conditional on $u_t$, the shock vector $e_t$ has density
    \begin{equation}\label{eq:shockdist}
        p(e_t\mid u_t) = \prod_{i=1}^d \frac{Q_i(e_{it})}{Z_i(u_t)} \exp\left\{\sum_{j=1}^k q_{ij}(e_{it})\lambda_{ij}(u_t)\right\},
    \end{equation}
    where $Q_i$ are the base measures, $Z_i$ are the normalizing constants, $k\geq 2$ is the number of sufficient statistics, $q_{ij}(e_{it})$ are the sufficient statistics, and $\lambda_{ij}(u_t)$ are the corresponding natural parameters.
\end{enumerate}
\end{assumption}

Assumption~\ref{as:shocks} gives the general exponential-family setting underlying IIA. As in \cite{Morioka+Halva+Hyvarinen:2021}, this assumption allows the conditional shock distribution to be learned in a flexible manner. However, the resulting identification is then only up to permutation and component-wise invertible transformations. In this paper, we sharpen the identification result further by imposing additional structure with restrictions on the sufficient statistics.

More specifically, we assume a particular exponential-family form with two sufficient statistics. The family is chosen to satisfy several criteria simultaneously. First, the shocks should have zero conditional mean, so that they retain the usual interpretation as structural shocks. Second, the family should fit into the exponential-family framework of Assumption~\ref{as:shocks}. Third, it should involve at least two sufficient statistics, which is essential for obtaining enough variation to facilitate identification. Fourth, the sufficient statistics should be smooth enough for the subsequent theoretical analysis. Finally, the family should allow the auxiliary variable $u_t$ to induce economically meaningful changes in the conditional shock distributions.
These requirements lead to the following assumption.
\begin{assumption}\label{as:suffstats}
Assumption~\ref{as:shocks} is satisfied with $k=2$, $Q_i(e_{it})=1$, $q_{i1}(e_{it})=e_{it}^2$, $q_{i2}(e_{it})=\sqrt{e_{it}^2+\varepsilon}-\sqrt{\varepsilon}$, $\varepsilon>0$, and $\lambda_{i1}(u_t)<0$ for all $i=1,\ldots,d$ and $u_t$.
\end{assumption}

Assumption~\ref{as:suffstats} strengthens the generic IIA framework of \cite{Morioka+Halva+Hyvarinen:2021} by fixing the sufficient statistics, thereby restricting the admissible family of conditional shock distributions. The first sufficient statistic is quadratic as in the Gaussian distribution, while the second introduces an additional nonlinear term that allows the conditional density to deviate from Gaussianity. The smoothing constant $\varepsilon>0$ ensures that $q_{i2}(e_{it})$ is differentiable at zero. While the resulting family is somewhat nonstandard, it is deliberately constructed to satisfy the set of above-described requirements that are difficult to meet simultaneously within standard parametric families. The practical appeal of the imposed distribution is that, through the natural parameters, it can accommodate changes in both scale and shape based on the level of the auxiliary variable $u_t$, and thereby allows the probability of large shocks to differ across economic states. Our specifications for the natural parameters are introduced next, and for intuition, we also illustrate how the resulting density can change with $u_t$.


\subsection{Specification of the natural parameters}\label{sec:natpar}

The density~\eqref{eq:shockdist} together with Assumption~\ref{as:suffstats} specifies the functional form of the conditional shock distribution, but leaves the natural parameters $\lambda_{ij}(u_t)$ unspecified. These parameters can be modeled as functions of the auxiliary variable $u_t$ as described below.

\begin{figure}[!t]
\centering
\includegraphics[width=\textwidth - 2cm]{Figures/Density_logistic.png}
\caption{Conditional densities implied by the logistic specification. The panels correspond, from left to right, to the parameter value configurations
$(a_{i1},b_{i1},\mathfrak{c}_{i1},\kappa_{i1},a_{i2},b_{i2},\mathfrak{c}_{i2},\kappa_{i2})=(-1.10,1.20,1.4,1.8,0.25,-0.12,2.4,1.4)$,
$(-1.40,1.85,2.5,2.0,0.20,-0.18,3.0,1.5)$,
$(-1.05,1.45,3.0,1.7,0.30,-0.22,1.7,1.7)$, and
$(-1.20,1.55,2.0,7.5,0.30,-0.20,3.6,5.0)$.
Within each panel, the density is plotted for $u_t\in\{0,2,4\}$ with a green solid line, a blue dashed line, and a red dotted line, respectively.}
\label{fig:density_logistic}
\end{figure}


Under Assumption~\ref{as:suffstats}, the conditional density of $e_{it}$ can be written as
\begin{equation}\label{eq:shock_dist_det}
p(e_{it}\mid u_t)\propto \exp\!\left\{ -\tau_i(u_t)^{-1} e_{it}^2 -\tau_i(u_t)^{-1/2}\rho_i(u_t)^{-1} \left(\sqrt{e_{it}^2+ \varepsilon}-\sqrt{\varepsilon}\right) \right\},
\end{equation}
with the natural parameter specifications
\begin{align}
\lambda_{i1}(u_t) &= -\tau_i(u_t)^{-1} \ \ \text{and}\label{eq:lambda1} \\
\lambda_{i2}(u_t) &= -\tau_i(u_t)^{-1/2}\rho_i(u_t)^{-1},\label{eq:lambda2}
\end{align}
where $\tau_i(u_t)>0$ governs the overall scale of the distribution and $\rho_i(u_t)>0$ governs its shape. The quadratic term determines the asymptotic tail behavior, which remains Gaussian, while the second term affects the shape of the distribution closer to the center. In particular, it allows the probability of large (but not asymptotically extreme) realizations to vary with $u_t$, thereby modifying the ``shoulders'' of the distribution, which becomes closer to Gaussian  the larger the value of $\rho_i(u_t)$ is.

We consider the following two-regime logistic  specifications for $\tau_i(u_t)$ and $\rho_i(u_t)$:
\begin{align}
\Lambda_{ij}(u_t) &= \left[1 + \exp\lbrace -\kappa_{ij}(u_t - \mathfrak{c}_{ij})\rbrace\right]^{-1}, \ \ j=1,2, \\
\tau_i(u_t) &= \exp\left\{(1 - \Lambda_{i1}(u_t))a_{i1} + \Lambda_{i1}(u_t)b_{i1}\right\},\label{eq:tau} \\
\rho_i(u_t) &= \exp\left\{(1 - \Lambda_{i2}(u_t))a_{i2} + \Lambda_{i2}(u_t)b_{i2}\right\},\label{eq:rho}
\end{align}
where $a_{ij},b_{ij},\mathfrak{c}_{ij}\in\mathbb{R}$ and $\kappa_{ij}>0$, $j=1,2$, are parameters and the exponential transformation forces positivity. This specification maps the auxiliary variable into a bounded range for the scale and shape parameters, capturing transitions between different states of the shock distribution.


Figure~\ref{fig:density_logistic} illustrates the logistic specification by depicting the implied conditional density $p(e_{it}\mid u_t)$ for four different parameter value configurations and for $u_t\in\{0,2,4\}$. The panels are chosen to represent distinct qualitative scenarios. In particular, from left to right, the figure includes a balanced benchmark configuration, a case with a stronger and later scale transition, a configuration in which the shape transition begins before the scale transition, and a case with much sharper logistic transitions. Across all four panels, the specification accommodates the interpretation that larger values of $u_t$ may be associated with a larger dispersion of the conditional shock distribution, although the timing and magnitude of the change differ across parameter configurations.



\section{Identification via independent innovation analysis}\label{sec:identification}

This section introduces our approach to identifying structural shocks based on the independent innovation analysis (IIA) of \cite{Morioka+Halva+Hyvarinen:2021}. The key idea of IIA is to exploit variation in the distributions of the shocks induced by an auxiliary variable $u_t$, and to recover the shocks by learning a representation that captures this variation. Depending on the nature of the auxiliary variable $u_t$, \cite{Morioka+Halva+Hyvarinen:2021} propose three different learning frameworks: general contrastive learning (IIA-GCL) for observed and possibly continuous $u_t$, time-contrastive learning (IIA-TCL) for observed discrete $u_t$, and a hidden Markov model-based approach (IIA-HMM) for unobserved $u_t$. We focus on the IIA-GCL framework, which allows the auxiliary variable to be specified directly and interpreted as an observed exogenous driver of distributional variation in the shocks.

Our contribution builds on the IIA-GCL framework of \cite{Morioka+Halva+Hyvarinen:2021}, but imposes the exponential-family structure introduced in Section~\ref{sec:shockassumptions} on the conditional distribution of the shocks within the learning problem. This restriction plays a central role in our identification strategy. While the original IIA-GCL framework identifies the shocks up to permutation and componentwise invertible transformations, the imposed exponential-family structure allows us to reduce this ambiguity substantially, and the identification results in Section~\ref{sec:ident_res} show that the remaining transformed-scale ambiguity is generically eliminated under the logistic natural-parameter specification introduced in Section~\ref{sec:natpar}.

\subsection{General contrastive learning framework (IIA-GCL)}\label{sec:IIA-GCL}

For an observed and possibly continuous exogenous variable $u_t$, \cite{Morioka+Halva+Hyvarinen:2021} propose applying general contrastive learning (GCL) to independent innovation analysis. In the IIA-GCL framework, a feature extractor and a logistic classifier are trained jointly to distinguish the real data $(y_t,\boldsymbol{y}_{t-1},u_t)$ from a version where randomization is performed on $u_t$. The objective is to learn a representation of $(y_t,\boldsymbol{y}_{t-1})$ that captures the dependence of the shock distribution on $u_t$.

More specifically, define:
\begin{equation}\label{eq:y_tilde}
\tilde{\boldsymbol{y}}_t = (y_t,\boldsymbol{y}_{t-1},u_t)
\quad \text{and} \quad
\tilde{\boldsymbol{y}}_t^\ast = (y_t,\boldsymbol{y}_{t-1},u_t^\ast),
\end{equation}
where $u_t^\ast$ is drawn from the marginal distribution of $u_t$ independently of $(y_t,\boldsymbol{y}_{t-1})$. In practice, $u_t^\ast$ is obtained by randomly permuting the sample of $u_t$, $t=1,\ldots,T$. This preserves the marginal distribution of $u_t$ while breaking its dependence on $(y_t,\boldsymbol{y}_{t-1})$. Consequently, $\tilde{\boldsymbol{y}}_t$ follows the joint density $p(y_t,\boldsymbol{y}_{t-1},u_t)$, whereas $\tilde{\boldsymbol{y}}_t^\ast$ follows $p(y_t,\boldsymbol{y}_{t-1})p(u_t^*)$.

A binary classifier is then trained to discriminate between $\tilde{\boldsymbol{y}}_t$ and $\tilde{\boldsymbol{y}}_t^\ast$. Let $D_t=1$ indicate that an observation is real and $D_t=0$ that it is randomized. The classifier is based on a scalar-valued nonlinear regression function $r(\tilde{\boldsymbol{y}}_t)$, whose logistic transformation gives the posterior probability that the observation is real:
\begin{equation}
p(D_t=1 \mid y_t,\boldsymbol{y}_{t-1},u_t)
 = (1 + \exp\{-r(y_t,\boldsymbol{y}_{t-1},u_t)\})^{-1}.
\end{equation}
At the population level, the classifier maximizes the objective
\begin{equation}\label{eq:popobj}
\mathcal{L}(r) = \mathbb{E}_{\tilde{\boldsymbol{y}}_t}\big[\log \Lambda^{\text{std}}(r(\tilde{\boldsymbol{y}}_t))\big] + \mathbb{E}_{\tilde{\boldsymbol{y}}_t^\ast}\big[\log (1-\Lambda^{\text{std}}(r(\tilde{\boldsymbol{y}}_t^\ast)))\big],
\end{equation}
where $\Lambda^{\text{std}}(x)=(1+e^{-x})^{-1}$ is the standard logistic function and $\mathbb{E}_x[\cdot]$ denotes expectation with respect to the distribution of $x$. In practice, this corresponds to minimizing the empirical binary cross-entropy loss:
\begin{equation}\label{eq:lossbin}
\text{LOSS}_{\text{bin}} = -\sum_{t=1}^T \left[\log \Lambda^{\text{std}}(r(y_t,\boldsymbol{y}_{t-1},u_t)) + \log (1 - \Lambda^{\text{std}}(r(y_t,\boldsymbol{y}_{t-1},u_t^\ast)))\right].
\end{equation}

A key property of this classification problem is that, when $\mathcal{L}(r)$ is maximized over all measurable functions, the maximizer is given by the log-density ratio \citep[see, e.g.,][Appendix~A, and the references therein]{Morioka+Halva+Hyvarinen:2021}:
\begin{equation}\label{eq:ldr}
r^\circ(\tilde{\boldsymbol{y}}_t) = \log\frac{p(y_t,\boldsymbol{y}_{t-1},u_t)}{p(y_t,\boldsymbol{y}_{t-1})\,p(u_t)}.
\end{equation}
Thus, the population-optimal classifier measures the extent to which the joint distribution of $(y_t,\boldsymbol{y}_{t-1},u_t)$ departs from the product distribution that would prevail under independence between $u_t$ and $(y_t,\boldsymbol{y}_{t-1})$. The role of the learned representation is therefore to capture the features of the data through which this dependence operates.

To exploit the assumed structure of the shocks, we restrict the regression function $r$ to a class of functions that is consistent with the exponential family specification introduced in Section~\ref{sec:shockassumptions}. Specifically, we consider functions of the form:
\begin{equation}\label{eq:r}
r(\tilde{\boldsymbol{y}}_t) = \sum_{i=1}^d\sum_{j=1}^2q_{ij}(h_i(y_t,\boldsymbol{y}_{t-1}))\lambda_{ij}(u_t) - \sum_{i=1}^d \log Z_i(u_t) + \phi(\boldsymbol{y}_{t-1},u_t) + \eta(\boldsymbol{h}(y_t,\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1}),
\end{equation}
where $q_{ij}(h_i(y_t,\boldsymbol{y}_{t-1}))$ are the sufficient statistics given by Assumption~\ref{as:suffstats} and $\lambda_{ij}(u_t)$ are the corresponding natural parameter functions; $Z_i(u_t)=\int \exp\{\sum_{j=1}^2 q_{ij}(x)\lambda_{ij}(u_t)\}\operatorname{d}x$ are the normalizing constants; $\boldsymbol{h}(y_t,\boldsymbol{y}_{t-1})=(h_1(y_t,\boldsymbol{y}_{t-1}),\ldots,h_d(y_t,\boldsymbol{y}_{t-1}))\in\mathbb{R}^d$ is a candidate demixing map; $\phi(\boldsymbol{y}_{t-1},u_t)$ is a function of $(\boldsymbol{y}_{t-1},u_t)$; and $\eta(\boldsymbol{h}(y_t,\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1})$ is a function of $(\boldsymbol{h}(y_t,\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1})$ that does not depend on $u_t$.

The first term in~\eqref{eq:r} captures the interaction between the learned features and the auxiliary variable through the natural parameters of the conditional shock distributions, and it is this term that carries the identifying information. The remaining terms absorb contributions to the log-density ratio that are not directly informative about the shocks. In particular, the normalizing constants ensure a proper conditional density, $\phi(\boldsymbol{y}_{t-1},u_t)$ captures dependence involving only $(\boldsymbol{y}_{t-1},u_t)$, and $\eta(\boldsymbol{h}(y_t,\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1})$ collects terms that may depend on the learned representation and lagged variables but not on $u_t$.

In the theoretical analysis, the functions $\boldsymbol{h}$, $\phi$, and $\eta$ are treated as flexible population functions of their arguments, subject to the regularity conditions imposed in the identification results below. In practical implementations, they can be parametrized using feed-forward neural networks, motivated by their universal approximation properties \citep{Hornik+Stinchcombe+White:1989}. The important point is that the restricted regression function~\eqref{eq:r} is constructed to match the separable part of the population log-density ratio~\eqref{eq:ldr} induced by the conditional shock distribution. This restriction is what makes the learned representation informative about the structural shocks. The next subsection formalizes this intuition and states the corresponding population identification results.


\subsection{The identification results}\label{sec:ident_res}

We now establish that, under suitable conditions, the IIA-GCL framework described above identifies the structural shocks. The argument proceeds in three steps. Theorem~\ref{thm:IIA-GCL} first shows that the restricted IIA-GCL representation rules out the arbitrary componentwise nonlinear indeterminacies present in the original IIA-GCL result of \citet{Morioka+Halva+Hyvarinen:2021}. With the fixed sufficient statistics in Assumption~\ref{as:suffstats}, the remaining ambiguity is reduced to permutation, componentwise sign, and a one-parameter family of componentwise transformed-scale maps. Corollary~\ref{cor:nonclosure_ident} then gives a general condition on the admissible natural-parameter class under which this remaining transformed-scale ambiguity is resolved and the shocks are identified up to permutation and componentwise sign. Finally, Corollary~\ref{cor:logistic_ident} verifies this condition generically for the logistic natural-parameter specification introduced in Section~\ref{sec:natpar}. Together, these results imply that, under the logistic specification, the shocks are generically identified up to permutation and componentwise sign.
\begin{theorem}\label{thm:IIA-GCL}
Let
\begin{equation}\label{eq:qfixed}
\boldsymbol{q}(x) = \big(q_1(x),q_2(x)\big) = \big(x^2,\, \sqrt{x^2+\varepsilon}-\sqrt{\varepsilon}\big),
\end{equation}
where $\varepsilon>0$ is fixed and known. Assume the following:
\begin{enumerate}[label=(\alph*)]
    \item Observations $(y_t,\boldsymbol{y}_{t-1},u_t)$ are generated by the nonlinear SVAR~\eqref{eq:genvar}. The augmented model $\tilde{f}$ defined in~\eqref{eq:ftilde} is bijective, twice continuously differentiable, and its inverse $\tilde{g}$ defined in~\eqref{eq:gtilde} is twice continuously differentiable.\label{cond:smooth}

    \item The shocks $e_t$ satisfy Assumptions~\ref{as:shocks} and~\ref{as:suffstats}, with the fixed sufficient statistics $q_1$ and $q_2$ in~\eqref{eq:qfixed}.\label{cond:sufstat}

    \item Let $\nu$ denote the marginal law of $u_t$ and let $\mathcal U_0$ denote its support. There exists a nonempty open connected set $\mathcal Y_0\subset\mathbb R^{dp}$ and conditional densities specified for every $u\in\mathcal U_0$ such that the conditional distribution of $(e_t,\boldsymbol y_{t-1})$ given $u_t=u$ assigns probability one to $\mathbb R^d\times\mathcal Y_0$ and has Lebesgue density
    \begin{equation}\label{eq:pointwise_conditional_factorization}
    p(e_t,\boldsymbol y_{t-1}\mid u)
    =p(e_t\mid u)p(\boldsymbol y_{t-1}\mid u)
    \end{equation}
    on $\mathbb R^d\times\mathcal Y_0$. Here $p(e_t\mid u)$ is the exponential-family density in Assumptions~\ref{as:shocks} and~\ref{as:suffstats}, and both factors in~\eqref{eq:pointwise_conditional_factorization} are finite and strictly positive on their respective domains. Define
    \begin{equation}\label{eq:common_supports}
    \mathcal Z_0=\mathbb R^d\times\mathcal Y_0
    \quad\text{and}\quad
    \mathcal X_0=\tilde f(\mathcal Z_0).
    \end{equation}
    For $\boldsymbol x=(y_t,\boldsymbol y_{t-1})\in\mathcal X_0$, let $p_{\mathcal X}(\boldsymbol x\mid u)$ denote the conditional Lebesgue density induced by~\eqref{eq:pointwise_conditional_factorization} and the augmented model $\tilde f$, and let $p_{\mathcal X}(\boldsymbol x)$ denote its marginal density. Both densities are finite and strictly positive for every $(\boldsymbol x,u)\in\mathcal X_0\times\mathcal U_0$. The population log-density ratio~\eqref{eq:ldr} is then represented pointwise on this domain by
    \begin{equation}\label{eq:pointwise_ldr_theorem}
    r^\circ(\boldsymbol x,u)
    =\log p_{\mathcal X}(\boldsymbol x\mid u)-\log p_{\mathcal X}(\boldsymbol x).
    \end{equation}
    This function is jointly continuous on $\mathcal X_0\times\mathcal U_0$. Finally, $\mathcal U_0$ contains $2d+1$ distinct values $u^{(0)},u^{(1)},\ldots,u^{(2d)}$ for which the true variability matrix
    \begin{equation}\label{eq:L}
    L^0 = \left[\boldsymbol{\lambda}^0(u^{(1)})-\boldsymbol{\lambda}^0(u^{(0)}):\cdots : \boldsymbol{\lambda}^0(u^{(2d)})-\boldsymbol{\lambda}^0(u^{(0)})\right]
    \end{equation}
    is invertible, where
    \begin{equation*}
    \boldsymbol{\lambda}^0(u)
    =
    \big(\lambda_{11}^0(u),\lambda_{12}^0(u), \ldots,
    \lambda_{d1}^0(u),\lambda_{d2}^0(u)\big)
    \end{equation*}
    is the true vector of natural-parameter functions.\label{cond:variability}

    \item The population objective $\mathcal{L}(r)$ in~\eqref{eq:popobj} has at least one maximizer over the class of functions
    of the form~\eqref{eq:r}, with the fixed sufficient statistics $q_1$ and $q_2$ in~\eqref{eq:qfixed} and with learned natural-parameter functions. Every such maximizer satisfies the regularity requirements stated below. For an arbitrary maximizer, denote its feature extractor by
    \begin{equation*}
    \boldsymbol{h}^\ast(y_t,\boldsymbol{y}_{t-1}) = \big(h_1^\ast(y_t,\boldsymbol{y}_{t-1}),\ldots, h_d^\ast(y_t,\boldsymbol{y}_{t-1})\big),
    \end{equation*}
    and its learned natural-parameter vector by
    \begin{equation*}
    \boldsymbol{\lambda}^\ast(u)
    =
    \big(\lambda_{11}^\ast(u),\lambda_{12}^\ast(u), \ldots,
    \lambda_{d1}^\ast(u),\lambda_{d2}^\ast(u)\big).
    \end{equation*}
    Let $r^\ast$ denote the scalar function generated pointwise by these learned constituents through~\eqref{eq:r}. The function $r^\ast$ is jointly continuous on $\mathcal X_0\times\mathcal U_0$. The augmented mapping
    \begin{equation*}
    (y_t,\boldsymbol{y}_{t-1})
    \mapsto
    \big(\boldsymbol{h}^\ast(y_t,\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1}\big)
    \end{equation*}
    is a $C^2$ diffeomorphism from $\mathcal{X}_0$ onto its image.\label{cond:regulatory}
\end{enumerate}
Then there exist a permutation $\pi$ of $\{1,\ldots,d\}$, constants $\iota_1,\ldots,\iota_d\in\{-1,+1\}$, and constants
$c_1,\ldots,c_d>0$ such that, for all $(y_t,\boldsymbol{y}_{t-1})\in\mathcal{X}_0$,
\begin{align}
h_i^\ast(y_t,\boldsymbol{y}_{t-1}) &= \iota_i\,T_{c_i}\!\left(g_{\pi(i)}(y_t,\boldsymbol{y}_{t-1})\right),\ \ i=1,\ldots,d, \ \ \text{where}\\
T_c(x) &= \operatorname{sgn}(x)\sqrt{\left(\sqrt{\varepsilon} + c\left(\sqrt{x^2+\varepsilon}-\sqrt{\varepsilon}\right)\right)^2 - \varepsilon}, \ \ \operatorname{sgn}(0)=0.\label{eq:Tc}
\end{align}
Consequently, without further restrictions on the learned natural-parameter family, the structural shocks are identified up to permutation, componentwise sign, and the componentwise transformations $T_{c_i}$.
\end{theorem}
Theorem~\ref{thm:IIA-GCL} is proven in Appendix~\ref{sec:proof_thm1}. It is a population identification result for the restricted IIA-GCL classification problem. Condition~\ref{cond:smooth} ensures that the mapping between observables and shocks is invertible and sufficiently smooth, whereas Condition~\ref{cond:sufstat} imposes the distributional structure that restricts the admissible componentwise transformations. Condition~\ref{cond:variability} supplies the density and common-support regularity needed to compare the population log-density-ratio representation~\eqref{eq:ldr} pointwise over the auxiliary-variable support, and requires sufficiently rich variation in the natural parameters. Finally, Condition~\ref{cond:regulatory} ensures that each learned population maximizer is sufficiently regular for the almost-sure equality of the optimal regression functions to extend to the full common support.

Relative to the original IIA-GCL result of \citet[][Theorem~1]{Morioka+Halva+Hyvarinen:2021}, where the shocks are identified only up to arbitrary componentwise invertible transformations, Theorem~\ref{thm:IIA-GCL} substantially strengthens the identification by fixing the sufficient statistics~\eqref{eq:qfixed}. However, it does not give identification up to permutation and sign alone but leaves the componentwise transformed-scale maps $T_{c_i}$ in~\eqref{eq:Tc}. This remaining ambiguity arises because the chosen sufficient statistics have an algebraic closure property. Applying $T_c$ to a shock changes the sufficient-statistic vector in a way that can, in principle, be offset by a corresponding transformation of the natural-parameter path. Whether this compensation is admissible depends on the class of natural-parameter functions imposed on the true and learned conditional shock densities. The next corollary gives a general non-closure condition under which such compensation is impossible, so that the transformed-scale ambiguity collapses to componentwise sign.
\begin{corollary}\label{cor:nonclosure_ident}
Suppose the conditions of Theorem~\ref{thm:IIA-GCL} hold. Let $\mathfrak L$ be a single admissible class of two-dimensional componentwise natural-parameter paths, imposed on both the true and learned conditional shock densities and common to all component labels. Thus, for each $i=1,\ldots,d$, the paths
\begin{equation}
\boldsymbol\lambda_i^0(u)=\big(\lambda_{i1}^0(u),\lambda_{i2}^0(u)\big)
\quad\text{and}\quad
\boldsymbol\lambda_i^\ast(u)=\big(\lambda_{i1}^\ast(u),\lambda_{i2}^\ast(u)\big)
\end{equation}
belong to $\mathfrak L$.

For $c>0$, define
\begin{equation}\label{eq:Bc}
B(c)=
\begin{bmatrix}
c^2 & 2\sqrt{\varepsilon}\,c(1-c)\\
0 & c
\end{bmatrix}.
\end{equation}
Assume that $\mathfrak L$ is not closed under a nontrivial transformed-scale representation of any true component. Specifically, for each $k=1,\ldots,d$, there do not exist a constant $c>0$ with $c\neq1$, an alternative admissible path $\bar{\boldsymbol\lambda}(u)=\big(\bar\lambda_1(u),\bar\lambda_2(u)\big)\in\mathfrak L$, and a constant vector $\boldsymbol a=(a_1,a_2)\in\mathbb R^2$ such that
\begin{equation}\label{eq:natpar_nonclosure}
\boldsymbol\lambda_k^0(u)=B(c)'\bar{\boldsymbol\lambda}(u)+\boldsymbol a
\end{equation}
for every $u\in\mathcal U_0$.

Then any population maximizer satisfies
\begin{equation}
h_i^\ast(y_t,\boldsymbol y_{t-1}) = \iota_i g_{\pi(i)}(y_t,\boldsymbol y_{t-1}), \quad i=1,\ldots,d,
\end{equation}
on the common support $\mathcal X_0$, where $\pi$ is a permutation of $\{1,\ldots,d\}$ and $\iota_i\in\{-1,+1\}$. Hence, under the natural-parameter non-closure condition~\eqref{eq:natpar_nonclosure}, the shocks are identified up to permutation and componentwise sign changes.
\end{corollary}
Corollary~\ref{cor:nonclosure_ident} is proven in Appendix~\ref{sec:proof_cor1}. The non-closure condition rules out the possibility that a nontrivial componentwise transformation $T_{c_i}$ can be absorbed into another admissible natural-parameter path in $\mathfrak L$. Under this condition, the constants $c_i$ in Theorem~\ref{thm:IIA-GCL} must equal one, so that the remaining transformed-scale ambiguity collapses to componentwise sign. Thus, the admissible natural-parameter class plays two roles. It determines how the conditional shock distribution is allowed to vary with the auxiliary variable, and it determines whether the transformed-scale ambiguity left by Theorem~\ref{thm:IIA-GCL} is admissible.

The non-closure condition is stated abstractly because Theorem~\ref{thm:IIA-GCL} itself does not require a parametric model for the natural-parameter paths. In empirical work, however, it is often useful to estimate the natural parameters within a chosen admissible class. The following corollary specializes the non-closure argument to the logistic specification introduced in Section~\ref{sec:natpar} and shows that, under this specification, sign identification holds generically.
\begin{corollary}\label{cor:logistic_ident}
Suppose the conditions of Theorem~\ref{thm:IIA-GCL} hold. Suppose, in addition, that the true and learned componentwise natural-parameter paths are restricted to the logistic class $\mathfrak L_{\mathrm{log}}$ defined by Equations~\eqref{eq:lambda1}--\eqref{eq:rho} in Section~\ref{sec:natpar}. Assume that the support $\mathcal U_0$ of $u_t$ contains a nonempty open interval.

Then, for Lebesgue-almost every value of the true logistic parameter vector in its full parameter space, any population maximizer satisfies
\begin{equation}
h_i^\ast(y_t,\boldsymbol y_{t-1}) = \iota_i g_{\pi(i)}(y_t,\boldsymbol y_{t-1}), \quad i=1,\ldots,d,
\end{equation}
on the common support $\mathcal X_0$, where $\pi$ is a permutation of $\{1,\ldots,d\}$ and $\iota_i\in\{-1,+1\}$. Hence, under the logistic natural-parameter specification, the shocks are generically identified up to permutation and componentwise sign changes.
\end{corollary}
Corollary~\ref{cor:logistic_ident} is proven in Appendix~\ref{sec:proof_cor2}. The proof shows that the set of true logistic parameter values for which the non-closure condition of Corollary~\ref{cor:nonclosure_ident} can fail is Lebesgue-null. Consequently, under the maintained assumptions, the modified IIA-GCL stage generically recovers the structural shocks up to permutation and componentwise sign.\footnote{This is a population identification statement. In finite samples, the transformed-scale ambiguity removed by Corollary~\ref{cor:logistic_ident} may still be weakly separated from ordinary rescaling, particularly when the smoothing constant $\varepsilon$ is very small. It is therefore often useful to normalize the estimated shocks before estimating the nonlinear SVAR mapping~\eqref{eq:genvar} (see Section~\ref{sec:FFNN-SVAR}).}



\subsection{A Monte Carlo experiment}\label{sec:mc_evidence}

We examine the finite-sample performance of our modified IIA-GCL method with a small-scale Monte Carlo (MC) experiment. The experiment focuses on a two-variable nonlinear SVAR in which the true structural map is a feed-forward neural network with two hidden layers, four neurons in each hidden layer, and sigmoid activation (see Section~\ref{sec:FFNN-SVAR-subsec}). The benchmark design is correctly specified in the sense that the structural shocks follow the exponential-family conditional distribution in Assumption~\ref{as:suffstats}, with logistic natural-parameter functions of the form~\eqref{eq:tau}--\eqref{eq:rho}. To examine robustness to distributional misspecification, we also consider two designs in which the same nonlinear structural map is used but the conditional shock distributions are replaced by time-varying Student-$t$ and skewed-$t$ distributions \citep[see][for the latter]{Hansen:1994}.

The estimator is evaluated by comparing the recovered shocks with the true shocks after aligning the remaining sign and permutation indeterminacies. Specifically, for each replication we compute the componentwise absolute correlations between the true and recovered shocks under the permutation and sign choices that maximize the mean absolute correlations. We then report the mean correlations (and their standard deviations) across Monte Carlo replications. Details on the IIA-GCL training procedure are provided in Appendix~\ref{sec:iia_gcl_training_details}, whereas further details on the Monte Carlo experiment are reported in Appendix~\ref{sec:mc_details}.


\begin{table}[!t]
\centering
\begin{tabular}{lcccc}
\hline\\[-1.3ex]
DGP & $T=250$ & $T=500$ & $T=1000$ & $T=2000$\\
\hline\\[-1.3ex]
Exponential-family shocks & 0.87 (0.04) & 0.88 (0.04) & 0.89 (0.01) & 0.90 (0.02)\\
Student-$t$ shocks & 0.87 (0.05) & 0.88 (0.03) & 0.88 (0.02) & 0.89 (0.01)\\
Skewed-$t$ shocks & 0.88 (0.06) & 0.88 (0.05) & 0.89 (0.03) & 0.90 (0.02)\\
\hline
\end{tabular}
\caption{Monte Carlo shock-recovery correlations under minimum-BCE selection based on $100$ Monte Carlo replications. The entries report the Monte Carlo mean of the average aligned componentwise absolute correlation between the true and recovered shocks, with the Monte Carlo standard deviation in parentheses.}
\label{tab:mc_main_results}
\end{table}


The results in Table~\ref{tab:mc_main_results} show a reasonable recovery already in small samples, but the recovery improves only slowly with the sample size. The MC standard deviations mainly decrease with sample size, but not uniformly in our results based on $100$ MC replications. The results are very similar between the benchmark specification (Exponential-family) and the two misspecified cases (Student-$t$ and Skewed-$t$), suggesting robustness to misspecification of the shock distribution. The fitted demixer $\boldsymbol{h}$ uses a parsimonious one-hidden-layer architecture with $12$ neurons, which we find to provide more stable optimization across random initializations and better shock recovery at the sample sizes considered than comparable deeper specifications. Recovery may therefore eventually level off as network approximation error becomes relatively more important. With larger samples, increasing the capacity of the demixer (and the nuisance networks $\phi$ and $\eta$) could potentially improve recovery.

Because the neural-network optimization problem is nonconvex, different optimization starts can converge to different local solutions. Different starts may recover nearly the same shocks while obtaining slightly different terminal binary cross-entropies (BCE) through different nuisance-network fits (see~\eqref{eq:lossbin}). Conversely, two starts with similar BCEs may recover materially different shocks. For this reason, we group the local solutions into basins according to the similarity of the recovered shocks, after which the shock basin needs to be selected. In the benchmark results reported in Table~\ref{tab:mc_main_results}, we select the shock basin containing the solution with the smallest BCE. Since the selected shock basin can affect the recovery, and in empirical applications it may be useful to include diagnostic criteria in the basin selection rule, we report the MC results in Appendix~\ref{sec:basin_selection} for various basin selection methods. In this MC experiment, the considered alternative basin selection methods produce results similar to those reported in Table~\ref{tab:mc_main_results}, i.e., a reasonable recovery in small samples that improves only slowly with the sample size.


\section{A neural network SVAR and structural analysis}\label{sec:FFNN-SVAR}


\subsection{A feed-forward neural network SVAR}\label{sec:FFNN-SVAR-subsec}

After identifying the structural shocks, the next step is to approximate and estimate the unknown function $f$ in the SVAR model~\eqref{eq:genvar}, which maps lagged observations and contemporaneous shocks to the observed variables. Following \citet{Morioka+Halva+Hyvarinen:2021}, we approximate this function using a feed-forward neural network (FFNN), which provides a flexible parametric representation capable of capturing complex nonlinear relationships. In particular, sufficiently large feed-forward networks can approximate a wide class of measurable functions arbitrarily well \citep{Hornik+Stinchcombe+White:1989}, making them an appealing choice when the functional form of $f$ is left unrestricted.

In practice, it is convenient to normalize the recovered shocks prior to estimating the structural mapping. Under the logistic natural-parameter specification and the generic conditions of Corollary~\ref{cor:logistic_ident}, the remaining population ambiguity is limited to permutation and componentwise sign. Nevertheless, finite-sample estimates may still differ in scale across components, and the absolute scaling of the recovered shocks is inconsequential for the subsequent neural-network approximation because it can be absorbed by the nonlinear function $f_\theta$. For these reasons, we rescale each recovered shock component by its sample standard deviation and denote the resulting normalized shocks by $\tilde e_t$. In addition, a sign convention is fixed for each component, which is without loss of generality.

We consider the nonlinear SVAR model
\begin{equation}\label{eq:nn_svar}
y_t = f_{\theta}(\boldsymbol{y}_{t-1}, \tilde{e}_t),
\end{equation}
where $f_{\theta}:\mathbb{R}^{d(p+1)}\to\mathbb{R}^d$ is an FFNN with parameter vector $\theta$.
An FFNN maps an input vector to an output vector through a sequence of layers. Each layer consists of a collection of nodes (or neurons), where every neuron in a given layer is connected to all neurons in the preceding layer via a system of weights and biases. Specifically, each layer computes an affine transformation of its inputs, followed by an elementwise nonlinear activation function.

To illustrate the structure of an FFNN, let $z_t = (\boldsymbol{y}_{t-1}, \tilde e_t) \in \mathbb{R}^{d(p+1)}$ denote the input vector, and consider a network with two hidden layers containing $n_1$ and $n_2$ neurons, respectively. Then the network output can be written as
\begin{equation}\label{eq:nn_example}
f_{\theta}(z_t) = W_3 \,\sigma_2\!\left(W_2 \,\sigma_1\!\left(W_1 z_t + b_1\right) + b_2\right) + b_3,
\end{equation}
where $W_1 \in \mathbb{R}^{n_1 \times d(p+1)}$ and $b_1 \in \mathbb{R}^{n_1}$ are the weight matrix and bias vector of the first hidden layer; $W_2 \in \mathbb{R}^{n_2 \times n_1}$ and $b_2 \in \mathbb{R}^{n_2}$ correspond to the second hidden layer; and $W_3 \in \mathbb{R}^{d \times n_2}$ and $b_3 \in \mathbb{R}^{d}$ correspond to the output layer. The functions $\sigma_1:\mathbb{R}^{n_1}\to\mathbb{R}^{n_1}$ and $\sigma_2:\mathbb{R}^{n_2}\to\mathbb{R}^{n_2}$ are activation functions applied elementwise.
The activation functions introduce nonlinearity into the model, allowing the network to represent complex nonlinear relationships between inputs and outputs. A common choice is the leaky rectified linear unit (leaky ReLU), defined elementwise as $\sigma(x) = \max\{\alpha x, x\}$ for some small $\alpha > 0$, which mitigates certain optimization issues associated with standard ReLU (defined as $\max\{0, x\}$) while retaining computational simplicity \citep{Maas+Hannun+Ng:2013}.

More generally, a feed-forward neural network with $\mathcal{H}$ hidden layers can be defined recursively as
\begin{equation}
h_0(z_t) = z_t, \ \ h_\ell(z_t) = \sigma_\ell\!\left(W_\ell h_{\ell-1}(z_t) + b_\ell\right), \ \ \ell=1,\ldots,\mathcal{H},
\end{equation}
with output
\begin{equation}
f_{\theta}(z_t) = W_{\mathcal{H}+1} h_H(z_t) + b_{\mathcal{H}+1},
\end{equation}
where each hidden layer $\ell$ contains $n_\ell$ neurons, so that $W_\ell \in \mathbb{R}^{n_\ell \times n_{\ell-1}}$, $b_\ell \in \mathbb{R}^{n_\ell}$, $n_0 = d(p+1)$, and $n_{\mathcal{H}+1} = d$. The parameter vector $\theta$ collects all weights and biases across layers.
The flexibility of the network is governed by its architecture, in particular the number of layers $H$ and the number of neurons $n_\ell$ in each layer. Increasing these quantities improves the network's flexibility, but also increases the number of parameters and the risk of overfitting.

The use of a feed-forward architecture is natural in the present framework, where the structural model is specified in terms of a finite set of lagged variables and contemporaneous shocks. While recurrent neural networks, such as long short-term memory networks, are often used in time series applications \citep[see, e.g.,][and references therein]{Altmeyer+Agusti+Costa:2022}, they are designed to learn latent state representations over long sequences of observations. In contrast, the present model explicitly specifies the relevant state vector as $\boldsymbol{y}_{t-1}$, so an FFNN provides a direct and transparent representation of the structural SVAR.

The network parameters $\theta$ are estimated by minimizing a prediction loss, typically using numerical gradient-based methods. In particular, gradients of the loss function with respect to the network parameters can be computed efficiently using the backpropagation algorithm, which applies the chain rule through the layered structure of the network. The resulting gradients are then used within iterative optimization schemes such as stochastic gradient descent or its variants. See, for example, \citealp{Hastie+Tibshirani+Friedman:2009}; \citealp{Goodfellow+Bengio+Courville:2016}, for a more detailed discussion on neural networks and their training in general. Appendix~\ref{sec:svar_training_details} provides details on the FFNN-SVAR training employed in our empirical application and implemented in the accompanying R package iiasvar \citep{iiasvar}.


\subsection{Generalized impulse response function}\label{sec:GIRF}

To analyze the dynamic effects of structural shocks in our nonlinear SVAR model, it needs to be taken into account that the response of the variables may depend on both the history of the process as well as on the sign and size of the shock. To accommodate these features, we consider the generalized impulse response function (GIRF) \citep{Koop+Pesaran+Potter:1996}, which is a nonlinear counterpart of the conventional impulse response function. The GIRF is defined as:
\begin{equation}\label{eq:girf}
\text{GIRF}(h,\delta_i,\mathcal{F}_{t-1}) = \text{E}[y_{t+h}|\delta_i,\mathcal{F}_{t-1}] - \text{E}[y_{t+h}|\mathcal{F}_{t-1}],
\end{equation}
where $h$ is the horizon. The first term on the right side of (\ref{eq:girf}) is the expected realization of the process at time $t+h$ conditionally on the $i$th structural shock of sign and size $\delta_i \in\mathbb{R}$ at time $t$, and the previous observations. The latter term is the expected realization of the process conditionally on the previous observations only. The GIRF thus expresses the expected difference in the future outcomes when the $i$th structural shock of sign and size $\delta_i$ arrives at time $t$ as opposed to all shocks being random.

The dynamics of our model depend only on the current shocks and the previous $p$ observations, so the conditioning set $\mathcal{F}_{t-1}$ can be replaced by the lag vector $\boldsymbol{y}_{t-1}=(y_{t-1},...,y_{t-p})$. Because the contemporaneous dynamics may depend on the entire structural shock vector, the response generally depends not only on the specified shock $e_{it}=\delta_i$ but also on the remaining contemporaneous shocks. Since (\ref{eq:girf}) conditions only on the $i$th shock, the conditional expectation averages over the distribution of the remaining contemporaneous shocks, yielding the average effect of the $i$th structural shock of sign and size $\delta_i$ rather than the response associated with any particular realization of the other contemporaneous shocks.\footnote{Unlike in conventional nonlinear SVAR models, in our SVAR model~\ref{eq:genvar}, the effects of a shock may generally depend on the full contemporaneous shock vector and not just on the sign and size of the shock of interest. Nevertheless, our GIRFs average over the distribution of the rest of the shocks to obtain impulse response with interpretations similar to those of conventional nonlinear SVAR models.}

In practice, we evaluate GIRFs for histories drawn from the observed sample \citep[cf.][]{Lanne+Virolainen:2025}. For a given history $\boldsymbol{y}_{t-1}$, the response to a structural shock is computed by simulating the model forward under the estimated dynamics. The impact shock is set to a given value $\delta_i$, and in particular we set $\delta_i$ equal to the recovered structural shock associated with the history, so that the responses reflect empirically relevant combinations of states and shock signs and sizes. Future shocks are generated from the fitted conditional distribution, which depends on the auxiliary variable $u_t$. For the initial period, we use the observed value of $u_t$ corresponding to the given history. Along the simulated paths, the observed future values $u_{t+h}$ are used whenever they are available in the sample. If the simulation horizon extends beyond the available observations, the missing values of $u_{t+h}$ are generated from a fitted univariate autoregressive model for $u_t$, conditional on the available past observations of $u_t$.

This procedure yields, for each history, a path-specific response that depends on both the initial state and the magnitude of the shock. Consequently, state-dependent effects can be analyzed by comparing GIRFs across different subsets of histories. Algorithm~\ref{algo:girf} \citep[adapted from the algorithm in Appendix~B of][]{Lanne+Virolainen:2025} formalizes the computation of the GIRF for a given history. The GIRFs are computed under the fitted model, combining the estimated dynamics $f_{\theta}$ and the fitted conditional distribution of shocks, and are therefore subject to the associated estimation uncertainty.\footnote{The conditional shock distribution is fitted in the IIA-GCL training for shocks that have not been normalized for sign and scale. To simulate observations from the (componentwise) normalized shock distributions, we first simulate a shock from the nonnormalized shock distribution, and then normalize it \textit{ex post}.}

\begin{algorithm}[!ht]
\caption{Generalized impulse response function}
\begin{algorithmic}[1]\label{algo:girf}
  \STATE Set the horizon $H$ and the number of replications $R$. For the $j$th Monte Carlo replication, let $y_{t+h}^{(j)}(\delta_i,\boldsymbol{y}_{t-1})$ denote a realization of the process at time $t+h$ conditional on the history $\boldsymbol{y}_{t-1}$ and the $i$th structural shock of sign and size $\delta_i \in\mathbb{R}$ arriving at time $t$, and let $y_{t+h}^{(j)}(\boldsymbol{y}_{t-1})$ denote an alternative realization conditional on the history $\boldsymbol{y}_{t-1}$ only.
  \FOR{each $j$ in $1,2,\ldots,R$}
    \FOR{each $h$ in $0,1,\ldots,H$}
      \STATE Take the observed value of $u_{t+h}$ from the data if available; otherwise, simulate $u_{t+h}$ from a fitted univariate autoregressive model for $u_t$ conditional on the previous observations.
\label{step:girf_shock_starts}
      \STATE Draw $e_{t+h}$ from its fitted conditional distribution given $u_{t+h}$.
      \IF{$h=0$}
        \STATE Impose the sign and size $\delta_i$ to the $i$th element of $e_{t+h}$ to obtain the modified structural shock vector $\tilde{e}_{t+h}$.
      \ENDIF \label{step:girf_shock_ends}
      \STATE Calculate $y_{t+h}^{(j)}(\delta_i,\boldsymbol{y}_{t-1})$ from the fitted model~\eqref{eq:nn_svar} using the previous observations $(y_{t}^{(j)}(\delta_i,\boldsymbol{y}_{t+h-1}),\ldots,y_{t}^{(j)}(\delta_i,\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1})$ and the structural shock vector $e_{t+h}$ (or $\tilde{e}_{t+h}$ if $h=0$) obtained from Steps~\ref{step:girf_shock_starts}--\ref{step:girf_shock_ends}.
      \STATE Calculate $y_{t+h}^{(j)}(\boldsymbol{y}_{t-1})$ from the fitted model~\eqref{eq:nn_svar} by using the previous observations $(y_{t+h-1}^{(j)}(\boldsymbol{y}_{t-1}),\ldots,y_{t}^{(j)}(\boldsymbol{y}_{t-1}),\boldsymbol{y}_{t-1})$ and the structural shock vector $e_{t+h}$ obtained from Steps~\ref{step:girf_shock_starts}--\ref{step:girf_shock_ends}.
    \ENDFOR
    \STATE Calculate $y_{t+h}^{(j)}(\delta_i,\boldsymbol{y}_{t-1}) - y_{t+h}^{(j)}(\boldsymbol{y}_{t-1})$.
  \ENDFOR
\STATE Calculate the sample mean of $y_{t+h}^{(j)}(\delta_i,\boldsymbol{y}_{t-1}) - y_{t+h}^{(j)}(\boldsymbol{y}_{t-1})$ across the Monte Carlo replications $j=1,\ldots,R$ for all $h=0,1,\ldots,H$ to obtain the GIRF related to the history $\boldsymbol{y}_{t-1}$ for the $i$th structural shock of sign and size $\delta_i$.
\end{algorithmic}
\end{algorithm}


\subsection{Shock labelling}\label{sec:shocklabelling}

After normalizing the sign and scale of the structural shocks, they remain identified up to permutation and are subject to the usual labelling problem. In other words, in line with the statistical identification literature, labelling the identified shocks as economic shocks requires external information. To address this, we follow a procedure similar in spirit to \cite{Lanne+Meitz+Saikkonen:2017} and propose labelling the shocks based on their estimated impulse response functions.

Specifically, we compute GIRFs for histories drawn from the observed sample and use the recovered shocks associated with these histories. To make the GIRFs comparable across different signs and sizes of the shocks, we normalize each response by dividing by the corresponding recovered shock. This transformation rescales the responses to correspond to a positive unit shock. The normalized responses are then aggregated across histories by taking the median response in each time period, yielding a representative response profile for each shock. Finally, the shocks of interest are labelled based on the patterns of these median responses across variables and horizons. The interpretation of each shock depends on how it affects the variables in the system over time and therefore on the specific empirical application.


\section{Empirical application}\label{sec:empirical}

We illustrate the use of our methods to study the effects of real oil price shocks on U.S. industrial production. In particular, our empirical application is motivated by \cite{Herrera+Lagalo+Wada:2011}, who study whether the response of U.S. industrial production to real oil price shocks is asymmetric. Their results show that the evidence depends on the estimation period and the level of aggregation. In post-1973 data, evidence against symmetry is weak for aggregate industrial production, but it is stronger for some sector-specific series.

We consider a bivariate monthly U.S. dataset covering the period from February 1974 through January 2026 ($624$ observations). Industrial production is measured by the industrial production index (IPI), which is logarithmized and detrended by taking first differences. Following \cite{Herrera+Lagalo+Wada:2011}, the nominal oil price is measured by the U.S. crude-oil composite refiner acquisition cost, which is deflated by the consumer price index to obtain the real oil price (OIL), which is then likewise logarithmized and detrended by taking first differences.

\begin{figure}[!ht]
    \centerline{\includegraphics[width=\textwidth - 2cm]{Figures/Seriesplot.png}}
    \caption{Monthly U.S. time series from 1974:2 through 2026:1. The top panel shows the first difference of the logarithm of the seasonally adjusted industrial production index (multiplied by $100$). The middle panel shows the first difference of the logarithm of the real oil price, constructed by deflating the U.S. crude-oil composite refiner acquisition cost by the consumer price index (multiplied by $100$). The bottom panel shows the one-month-ahead macro uncertainty index of \citet{Jurado+Ludvigson+Ng:2015}, lagged by one month and reported in its original units. The shaded areas indicate the periods of U.S. recessions defined by the NBER.}
\label{fig:seriesplot}
\end{figure}

As the auxiliary variable $u_t$, we use the first lag of the one-month-ahead macro uncertainty index (MUI) of \citet{Jurado+Ludvigson+Ng:2015}. The MUI summarizes uncertainty about the unpredictable component of a broad set of macroeconomic indicators, and periods of elevated macroeconomic uncertainty can plausibly feature changes in the variances and the shapes of macroeconomic shock distributions. It therefore provides an economically meaningful source of variation in the conditional shock distributions. Using the first lag, in turn, makes the auxiliary variable predetermined relative to the time-$t$ shocks, facilitating plausibility of Assumption~\ref{as:shocks}. The series of the considered variables are depicted in Figure~\ref{fig:seriesplot}.\footnote{The IPI series is obtained from \url{https://fred.stlouisfed.org/series/INDPRO}, the consumer price index from \url{https://fred.stlouisfed.org/series/CPIAUCSL}, the nominal oil price series from \url{https://www.eia.gov/dnav/pet/pet_pri_rac2_dcu_nus_m.htm}, and the macro uncertainty series from \url{https://www.sydneyludvigson.com/macro-and-financial-uncertainty-indexes}. The availability of the oil price data determines the starting period and the availability of the macro uncertainty data the ending period of our sample.}


\subsection{Training and model selection}\label{sec:emp_training}

Training our IIA-SVAR model is conducted in two stages. In Stage~1, the structural shocks are estimated with the IIA-GCL method described in Section~\ref{sec:IIA-GCL}. Specifically, the demixing network $\boldsymbol{h}$ and the logistic classifier~\eqref{eq:r} are trained jointly by minimizing the binary cross-entropy loss~\eqref{eq:lossbin}. In Stage~2, the recovered shocks are first normalized with respect to permutation, sign, and scale. Then, the nonlinear SVAR model~\eqref{eq:nn_svar} is estimated by training the neural network $f_{\theta}$ to predict the observations $y_t$ from the lag vector $\boldsymbol{y}_{t-1}$ and normalized recovered shocks $\tilde{e}_t$.

Model specification consists of selecting the autoregressive order $p$ as well as the architectures of the involved feed-forward neural networks. In particular, for the demixing network $\boldsymbol{h}$, the nuisance neural networks $\phi$ and $\eta$ (see~\eqref{eq:r}), as well as the SVAR network $f_{\theta}$, one must choose the number of hidden layers $\mathcal{H}$, the number of neurons in each hidden layer $n_\ell$, and the activation function $\sigma_{\ell}$, $\ell=1,...,\mathcal{H}$. The Akaike information criterion selects $p=3$ for a conventional linear VAR, indicating that three lags could be sufficient to capture the majority of the autocorrelation structure of the data, so we also set $p=3$ for our nonlinear SVAR model.

Table~\ref{tab:nn_architectures} summarizes the selected network architectures. In particular, the demixing and SVAR networks are assigned greater capacity than the nuisance networks because they perform the primary shock-recovery and SVAR-modeling tasks. Relatively shallow architectures are employed to balance flexibility and parsimony. The nuisance and SVAR networks use leaky-ReLU activation, whereas the demixer uses sigmoid activation because, in our finite-sample experiments (not shown), it produced more stable shock estimates across optimizer initializations without materially weakening shock recovery.\footnote{In particular, the demixer and nuisance networks follow the architecture from our Monte Carlo experiment (see Section~\ref{sec:mc_evidence}).}


\begin{table}[!ht]
\centering
\begin{tabular}{lccccc}
\hline \\[-1.3ex]
 & Demixer $\boldsymbol{h}$ & Nuisance $\phi$ & Nuisance $\eta$ & SVAR $f_\theta$ \\
\hline \\[-1.3ex]
Hidden layers $\mathcal{H}$        & 1 & 1 & 1 & 2 \\
Neurons in layer $n_\ell$          & 12 & 8 & 8 & 12 \\
Activation function $\sigma_\ell$ & Sigmoid & Leaky ReLU & Leaky ReLU & Leaky ReLU \\
\hline
\end{tabular}
\caption{Selected neural network architectures. The sigmoid activation is defined as $\sigma_\ell(x)=(1+e^{-x})^{-1}$ and Leaky-ReLU as $\sigma_\ell(x) = \max\{\alpha x, x\}$ for some small $\alpha > 0$, where we use $\alpha=0.01$.}\label{tab:nn_architectures}
\end{table}


\begin{figure}[!ht]
    \centerline{\includegraphics[width=\textwidth - 1cm]{Figures/Shockplot.png}}
    \caption{Each column corresponds to one recovered structural shock component obtained from the Stage-1 IIA-GCL estimator. The top panels display the recovered shocks $\hat{e}_{it}$ (black solid line) together with the standard deviations of their fitted distributions (purple dotted line). The bottom panels display the (standardized) auxiliary variable $u_t$ (solid black line) together with the estimated modulation functions governing the conditional distribution of the shock component, namely the scale parameter function $\hat{\tau}_i(u_t)$ (blue dashed line) and the shape parameter function $\hat{\rho}_i(u_t)$ (orange dashed line). The modulation functions are plotted on a common rescaled axis for visualization, with their original scale indicated by the colored right-hand axis. The shaded areas indicate the U.S. recessions defined by the NBER.}
\label{fig:shockplot}
\end{figure}

Training the IIA-GCL model for our dataset, with further details on the training procedure provided in Appendix~\ref{sec:iia_gcl_training_details}, recovers the normalized shocks presented in the top panels of Figure~\ref{fig:shockplot}.  The bottom panels of Figure~\ref{fig:shockplot} depict the fitted scale (blue dashed line) and shape (orange dashed line) parameter values for each shock along with the auxiliary variable $u_t$, showing how their values vary with $u_t$. For both shocks, higher macroeconomic uncertainty is associated with larger values of the scale parameter. Interestingly, for Shock~1, elevated uncertainty is associated with larger values of the shape parameter, whereas for Shock~2 it is the other way around. An increase in the shape parameter implies that the shock distribution becomes more Gaussian (see Section~\ref{sec:natpar}), while both the scale and shape parameters jointly affect its dispersion. To provide a more direct measure of time-varying dispersion, the estimated standard deviations of the fitted shock distributions are shown in the top panels of Figure~\ref{fig:shockplot} (purple dotted line). These estimates indicate that the dispersion of both shock distributions tends to increase with macroeconomic uncertainty. The diagnostics presented in Appendix~\ref{sec:empapp_shockdiag} (Figure~\ref{fig:shock_acf}) show that there is not much auto- or crosscorrelation in the recovered shocks.\footnote{
We also check that Condition~\ref{cond:variability} of Theorem~\ref{thm:IIA-GCL} (the variability condition) is satisfied in the estimate by randomly sampling $100000$ length $2d+1$ subsets of observed $u_t$ and constructing empirical analogues of the matrix $L$ in~\eqref{eq:L}. In the subset that has the largest smallest singular value of $L$, the smallest singular value is $0.04$ and the ratio between the largest and smallest singular value is $29$, suggesting that the variability condition is well satisfied in the estimate.}


Using the recovered (normalized) shocks, we then train the nonlinear SVAR model~\eqref{eq:nn_svar} implemented as the FFNN described in Table~\ref{tab:nn_architectures}. Details on the training procedure are provided in Appendix~\ref{sec:svar_training_details}. The related prediction error diagnostics presented in Appendix~\ref{sec:svar_diagnostics} suggest a reasonable fit. In particular, there is no obvious systematic variation nor substantial auto- or crosscorrelation in the prediction errors, and apart from two particularly large outliers for oil prices (related to the Gulf War oil shock in August 1990 and the COVID-19 shock in March 2020), the prediction errors seem reasonable.


\subsection{Labelling the shocks}\label{sec:emp_shocklabelling}


Following the procedure described in Section~\ref{sec:shocklabelling}, we label the structural shocks based on their median GIRFs across all length $p$ histories in the data. These median GIRFs computed for our IIA-SVAR model fitted in Section~\ref{sec:emp_training} are presented in Figure~\ref{fig:labellingplot}. Shock~1 moves output and oil prices in the same direction, while Shock~2 moves them in opposite directions. Moreover, Shock~1 has a larger impact effect on production than Shock~2, while Shock~2 has a larger impact effect on oil prices. Hence, the median effects of Shock~2 seem to resemble those of a (real) oil price shock, so we deem it as the oil price shock.

\begin{figure}[!t]
    \centerline{\includegraphics[width=\textwidth - 1cm]{Figures/Girf_labelling.png}}
    \caption{The median GIRFs across all the length $p$ histories in the data for the recovered structural shocks, where the responses are normalized to correspond to a positive unit shock. Each panel shows the response of one variable to each shock, with the left panel for IPI and the right panel for OIL, covering the horizons of $h=0,1,\ldots,24$ months. The responses are in the approximate percentage growth rate scale on which the variables are presented in Figure~\ref{fig:seriesplot}.}
\label{fig:labellingplot}
\end{figure}

\subsection{Impulse response analysis}

Median GIRFs provide compact summary statistics on the effects of the shocks useful for labelling them, but in general the effects of the shocks may depend on the related history and sign and size of the shock. Here we are interested in studying asymmetries in the effects of the oil price shock.
The GIRFs are computed using Algorithm~\ref{algo:girf} with a procedure comparable to \cite{Lanne+Virolainen:2025} and \cite{Virolainen2:2024}. Specifically, for each length $p$ history in the data, we take the corresponding oil price shock recovered from the data and compute the corresponding GIRF. Each GIRF is assigned to a sign group according to the sign of the recovered shock. The GIRFs are then normalized to correspond to a positive unit impact response of OIL (i.e., roughly a one-percentage-point increase in the real oil price). Consequently, the comparisons below summarize model-implied heterogeneity across the observed history-shock pairs rather than ceteris-paribus effects of changing only the shock sign or history. The reported $90\%$ GIRF intervals are empirical quantile ranges across these pairs and do not quantify parameter uncertainty.

Firstly, motivated by \cite{Herrera+Lagalo+Wada:2011}, we study how the effects of the oil price shock vary between positive and negative shocks. Figure~\ref{fig:posnegplot} illustrates the distribution of the GIRFs for positive and negative shocks, where the median GIRF is depicted with a solid line (blue for negative and red for positive shocks), with the shaded areas quantifying the dispersion of the GIRFs across the related histories (and shocks). The median responses of both IPI and OIL are slightly stronger for negative than positive shocks. However, the GIRF intervals show that the dispersion of the responses to negative shocks is clearly skewed towards stronger responses, while for positive shocks it is skewed towards weaker responses. The effects of negative shocks thereby seem overall stronger than the effects of positive shocks. Our findings thus differ from the post-1973 (aggregate IPI) result of \citet{Herrera+Lagalo+Wada:2011} for shocks of typical magnitude. In particular, they find no evidence against symmetry in that case, although they do find evidence of asymmetry for large shocks.

\begin{figure}[!t]
    \centerline{\includegraphics[width=\textwidth - 1cm]{Figures/Posneg_girf.png}}
    \caption{Generalized impulse response functions to the identified oil price shock $h=0,1,\ldots,24$ months ahead, computed with Algorithm~\ref{algo:girf}. The left panel shows the responses of IPI and the right panel the responses of OIL, both accumulated to ($100\times$)log-levels. The blue solid line shows the responses for negative shocks and the red solid line shows the responses for positive shocks. There are $212$ negative and $409$ positive shocks. The shaded areas (also indicated by dashed lines) are the symmetric $90\%$ GIRF intervals that contain $90\%$ of the GIRFs across the histories and shocks of the given sign. All the GIRFs have been scaled so that the instantaneous increase in OIL is unity.}
\label{fig:posnegplot}
\end{figure}

\begin{figure}[!t]
    \centerline{\includegraphics[width=\textwidth - 1cm]{Figures/Momentum_girf.png}}
    \caption{Generalized impulse response functions to the identified oil price shock $h=0,1,\ldots,24$ months ahead, computed with Algorithm~\ref{algo:girf}. The left panel shows the responses of IPI and the right panel the responses of OIL, both accumulated to ($100\times$)log-levels. The blue solid line shows the responses for histories with low IPI momentum and the red solid line shows the responses for histories with high IPI momentum. There are $117$ histories with low IPI momentum and $122$ histories with high IPI momentum. The shaded areas (also indicated by dashed lines) are the symmetric $90\%$ GIRF intervals that contain $90\%$ of the GIRFs across the histories (and shocks) of the given momentum. All the GIRFs have been scaled so that the instantaneous increase in OIL is unity.}
\label{fig:momentumplot}
\end{figure}

Secondly, our IIA-SVAR model has the interesting property that it allows effects of the shocks to vary depending on the joint features of the length $p$ history of the included variables and the contemporaneous shocks, where this dependency is learned from the data during the Stage-2 neural network training. We find that the effects of the oil price shock seem to particularly depend on the momentum of the industrial production growth in the previous $p$ periods. Specifically, we define the momentum as
$\text{IPI}^{\text{mom}}_t = \sum_{i=1}^p \text{IPI}_{t-i}$,
where $\text{IPI}_{t-i}$ is the approximate percentage growth rate of the industrial production index at lag $i$. Then, we separately consider histories with low momentum, $\text{IPI}^{\text{mom}}_t<-0.5$, and histories with high momentum, $\text{IPI}^{\text{mom}}_t>1.5$.

Figure~\ref{fig:momentumplot} shows that even though the responses of oil prices are similar for low and high momentum histories, the responses of production are somewhat different. In particular, they are generally stronger when the recent IPI momentum is low than when it is high. Overall, we therefore find modest asymmetries in the effects of the oil price shock with respect to its sign and the state of the economy.


\section{Conclusion}\label{sec:conclusion}
We develop a fully nonlinear structural vector autoregressive framework identified through variation in the conditional distributions of the mutually independent structural shocks, induced by an observed exogenous variable. In particular, a general contrastive learning framework that makes use of this variation together with the assumed exponential-family structure is employed to recover the shocks. The framework allows the contemporaneous structural mapping to be fully nonlinear and non-additive, thereby substantially extending the class of admissible structural dynamics relative to conventional nonlinear SVAR models.

Our identification strategy builds on the independent innovation analysis framework of \citet{Morioka+Halva+Hyvarinen:2021} based on general contrastive learning. Their identification results recover the structural shocks only up to permutation and componentwise invertible transformations, which is generally insufficient for structural econometric analysis. By imposing a structured exponential-family specification for the conditional shock distributions, we reduce this ambiguity substantially. The general identification result leaves only a transformed-scale ambiguity, and we show that this remaining ambiguity is generically eliminated under our logistic natural-parameter specification, yielding identification up to permutation and componentwise sign changes.

Once the structural shocks have been recovered, the fully nonlinear SVAR model can be estimated flexibly using feed-forward neural networks, motivated by their universal approximation capabilities \citep[see][]{Hornik+Stinchcombe+White:1989}. The resulting framework permits structural analysis through generalized impulse response functions \citep{Koop+Pesaran+Potter:1996} that accommodate nonlinear and state-dependent dynamics. The accompanying R package iiasvar \citep{iiasvar} implements the introduced methods.

We illustrate the use of methods in an empirical application considering a monthly bivariate system of U.S. industrial production growth and real oil price growth, covering the period from 1974:2 to 2026:1. To facilitate identification, the shock distributions are assumed to be modulated by the (lagged) macro uncertainty index of \citet{Jurado+Ludvigson+Ng:2015}. We find modest asymmetries in the effects of the oil price shock with respect to its sign and the state of the economy. In particular, our findings suggest stronger effects for negative than positive oil price shocks, differing from the post-1973 (aggregate industrial production) result of \citet{Herrera+Lagalo+Wada:2011} for shocks of typical magnitude. Moreover, we find the effects of the oil price shock stronger when the recent industrial production growth is low than when it is high.




\section*{Declaration of the use of AI and AI-assisted technologies}
During the preparation of this work, the author used ChatGPT and Codex to improve the language, write the practical R implementation of the method, and brainstorm ideas, including the preparation of mathematical proofs. After using these tools, the author reviewed and edited the content as needed and takes full responsibility for the content of the publication.

\bibliography{/Users/savi/Documents/texrefs/masterrefs.bib}

\pagebreak