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.
78,472 characters
Impulse Response Analysis of Structural Nonlinear Time Series Models
\onehalfspacing
\title{
Impulse Response Analysis of Structural Nonlinear Time Series Models
}
\author{Giovanni Ballarin
\thanks{
E-mail: \texttt{[email removed]}.
I thank Otilia Boldea, Timo Dimitriadis, Juan Carlos Escanciano, Lyudmila Grigoryeva, Klodiana Istrefi, Marina Khismatullina, So Jin Lee, Yuiching Li, Sarah Mouabbi, Andrey Ramirez, Christoph Rothe, Carsten Trenkler and Mengshan Xu, as well as the participants of the Econometrics Seminar at the University of Mannheim, the 2023 ENTER Jamboree, the 10th HKMEtrics Workshop, the GSS Weekly Seminar at Tilburg University, the Internal Econometrics Seminar at Vrije Universiteit Amsterdam, the Brown Bag Seminar at the University of St. Gallen and the 2024 ESAM for their comments, suggestions and feedback.
A significant part of this work was developed at the University of Mannheim thanks to the support of the Center for Doctoral Studies in Economics and the Chair of Empirical Economics.
}\\
University of St. Gallen
}
\date{
\today
}
\makeatletter
\let\thetitle\@title
\let\theauthor\@author
\let\thedate\@date
\makeatother
\maketitle
\vspace{-2em}
\begin{abstract}
\textbf{Abstract:}
This paper proposes a semiparametric sieve approach to estimate impulse response functions of
nonlinear time series within a general class of structural autoregressive models. We prove that a two-step procedure can flexibly accommodate nonlinear specifications while avoiding the need to choose fixed parametric forms. Sieve impulse responses are proven to be consistent by deriving uniform estimation guarantees, and an iterative algorithm
makes it straightforward to compute them in practice. With simulations, we show that the proposed
semiparametric approach proves effective against misspecification while suffering only from minor efficiency
losses. In a U.S. monetary policy application, the pointwise sieve GDP response associated
with an interest rate increase is larger than that of a linear model. Finally, in an analysis of interest rate uncertainty shocks, sieve responses indicate more substantial
contractionary effects on production and inflation.
\end{abstract}
\noindent\textit{Keywords:} macroeconometrics, semiparametric, sieve estimation, physical dependence
\noindent\textit{JEL:} C14, C22, C54, E52, F40
\pagebreak
\doublespacing
\section{Introduction}
Linearity is a foundational assumption in structural time series modeling.
For example, large classes of macroeconomic models in modern New Keynesian theory can be linearized, justifying the use of the linear time series toolbox from a theoretical point of view. The seminal work of \cite{simsMacroeconomicsReality1980} on vector autoregressive (VAR) models brought the study of dynamic economic relationships into focus within the macro-econometric literature, for which the estimation and analysis of impulse response functions (IRFs) is key \citep{hamilton1994TimeSeriesAnalysis,lutkepohlNewIntroductionMultiple2005,kilianStructuralVectorAutoregressive2017}.
The local projection (LP) approach of \cite{jordaEstimationInferenceImpulse2005} has also gained popularity as a flexible and easy-to-implement alternative.
Linear econometrics models are, however, limited in the kind of effects that they can describe.
In nonlinear DGPs, linear VARs as well as standard LP methods can only reconstruct the best linear impulse responses approximation \citep{plagborg-mollerLocalProjectionsVARs2021}.
And even though asymmetries in monetary policy and non-proportional shock effects are now commonly studied, most works still rely on parametric specifications. For example, \cite{tenreyro2016pushing} study both sign and size effects of monetary policy (MP) shocks using censoring and cubic transformations, respectively.
\cite{caggianoEconomicPolicyUncertainty2017,pellegrinoUncertaintyMonetaryPolicy2021} and \cite{caggianoUncertaintyShocksGreat2021} use multiplicative interacted VAR models to estimate the effects of uncertainty and MP shocks.
From a macro-finance perspective, \cite{forniNonlinearTransmissionFinancial2023,forniAsymmetricEffectsNews2023} study the economic effects of financial shocks following a quadratic VMA specification \citep{debortoliAsymmetricEffectsMonetary2020}.
\cite{gambettiBadNewsGood2022} study news shocks asymmetries by imposing that news changes enter an autoregressive model through a threshold map. Parametric nonlinear specifications are also common prescriptions in time-varying models \citep{auerbachMeasuringOutputResponses2012,caggianoEstimatingFiscalMultipliers2015} and state-dependent models \citep{ramey2018government,goncalvesStatedependentLocalProjections2024}.
In this paper, we aim to design a semiparametric, structural nonlinear time series modeling and estimation framework with explicit theoretical properties.
Our primary contribution is the extension and combination of the block-recursive structural framework of \cite{goncalvesImpulseResponseAnalysis2021} with the uniform sieve estimation theory of \cite{chenOptimalUniformConvergence2015} within a general physical dependence setup \citep{wuNonlinearSystemTheory2005}.
Under classical nonparametric assumptions, we show that a two-step semiparametric series estimation procedure can consistently recover the structural model in a uniform sense.
In order to be able to relax the assumptions on the linearity and fixed parametric model specifications, we restrict our study to the case of compactly supported, weakly dependent data.
We emphasize that, even in this constrained setting, to the best of our knowledge, our work is the first to offer a formal combination of these approaches.
We provide explicit guarantees for semiparametrically-estimated nonlinear IRFs: Nonlinear impulse response function estimates are asymptotically consistent and, thanks to an iterative algorithm, can also be straightforwardly computed.
To illustrate the validity of our proposed methodology, we first evaluate its performance with several simulations. With realistic sample sizes, the efficiency costs of our semiparametric procedure are small compared to correctly specified parametric responses. A second set of simulations provides a simple setup where the nonlinear parametric model is mildly misspecified. Still, the large-sample bias is considerable, while for semiparametric estimates it is negligible.
We then evaluate how the IRFs computed using the new method compare with those from two empirical exercises studied in the literature.
In a small, quarterly model of the U.S. macroeconomy based on \cite{tenreyro2016pushing}, we find that point estimates of linear and parametric nonlinear IRFs may underestimate in intensity the GDP responses by up to $13\%$ and $16\%$, respectively, after a large exogenous monetary policy shock. Moreover, sieve responses achieve maximum impact a year before their linear counterparts.
Then, we evaluate the effects of interest rate uncertainty on US output, prices, and unemployment following \cite{istrefiSubjectiveInterestRate2018}. In this exercise, the impact on industrial production of a one-deviation increase in uncertainty is $54\%$ stronger according to semiparametric IRFs than the comparable linear specification. These findings suggest that structural responses based on linear specifications can significantly underestimate the effects of shocks.
\paragraph{Literature Review.}
Let us mention some key references directly related to our discussion.
On the one hand, \cite{jordaEstimationInferenceImpulse2005} already proposed a ``{flexible local projection}'' approach based on the Volterra expansion.
The flexible LP proposal is effectively equivalent to adding polynomial terms to a linear regression, meaning it is a semiparametric method, and it should be analyzed as such. Yet, the Volterra expansion is not formally justified, nor is its truncation, which is key in studying its properties \citep{SirotkoSibirskaya2020volterraBootstrap,Movahedifar2023closedloopVolterra}.
On the other hand, smoothed LP methods, see e.g. \cite{plagborg-moller2016essays,barnichonImpulseResponseEstimation2019}, can address exclusively concerns of regularity in the shape of estimated impulse responses, but not any potential underlying nonlinearities in the DGP.
Recently, \cite{goncalvesStatedependentLocalProjections2024} outlined a general, nonparametric LP estimation procedure for nonlinear IRFs, which was later studied in \cite{goncalvesNonparametricLocalProjections2024} under high-level conditions on the functional form of the IRF itself.
\cite{gourieroux2023nonlinear} also devise a framework for nonparametric kernel estimation and inference of IRFs via local projections, although they mostly work in the one-dimensional, single lag case.
Finally, following the Generalized IRF (GIRF) approach \citep{Koop1996,potterNonlinearImpulseResponse2000,gourierouxNonlinearInnovationsImpulse2005,terasvirtaModellingNonlinearEconomic2010}, \cite{kanazawaRadialBasisFunctions2020} proposed to use radial basis function neural networks to estimate nonlinear reduced-form GIRFs for the U.S. economy. While GIRFs can be essentially characterized as impulse responses with more sophisticated conditioning sets, they are lacking in that they do not inherently address the problem of {structural} identification \citep{kilianStructuralVectorAutoregressive2017}.
\paragraph{}
The remainder of this paper is organized as follows. Section~\ref{section:framework} provides the general framework for the structural model.
Section~\ref{section:estimation} describes the two-step semiparametric estimation strategy, and Section~\ref{section:impulse_response_analysis} discusses nonlinear impulse response function computation, validity, and consistency.
In Section~\ref{section:simulations} we give a brief overview of simulation results, while Section~\ref{section:applications} contains the empirical analyses. Finally, Section~\ref{section:conclusion} concludes. All proofs and additional content can be found in the Appendix.
Concerning notation: scalar and vector random variables are denoted in capital or Greek letters, e.g. $Y_t$ or $\epsilon_t$, while realizations are shown in lowercase Latin letters, e.g. $y_t$. For a process $\{Y_t\}_{t\in \mathbb{Z}}$, we write $Y_{t:s} = (Y_t, Y_{t+1}, \ldots, Y_{s-1}, Y_s)$, as well as $Y_{*:t} = (\ldots, Y_{t-2}, Y_{t-1}, Y_t)$ for the left-infinite history and $Y_{t:*} = (Y_t, Y_{t+1}, Y_{t+2}, \ldots)$ for its right-infinite history. The same notation is also used for random variable realizations.
For a matrix $A \in \mathbb{R}^{d \times d}$ where $d \geq 1$, $\lVert A \rVert$ is the spectral norm, $\lVert A \lVert_\infty$ is the supremum norm and $\lVert A \rVert_r$ for $0 < r < \infty$ is the $r$-operator norm. For a random vector or matrix, we will use $\lVert \,\cdot\, \rVert_{L^r}$ to denote the associated $L^r$ norm.
\section{Model Framework}\label{section:framework}
In this section, we introduce the general nonlinear time series model, which is a generalization of the one developed in \cite{goncalvesImpulseResponseAnalysis2021}.
In terms of structural shocks identification, the idea is straightforward: A scalar series, $X_t$, is chosen to be the \textit{structural variable} identifying shocks, and it explicitly determines the dynamic effects on the remaining data, vector $Y_t$. This will enable the derivation of economically meaningful (structural) impulse responses due to an exogenous shock impacting $X_t$.
\subsection{General Model}\label{section:model_setup}
This paper focuses on the family of nonlinear autoregressive models of the form
\begin{equation}\label{eq_main:general_nonlin_model}
\begin{split}
X_t & = \mu_1 + A_{12}(L) Y_{t-1} + A_{11}(L) X_{t-1} + u_{1t}, \\
Y_t & = \mu_2 + G_2(Y_{t-1}, \ldots, Y_{t-p}, X_t, X_{t-1}, \ldots, X_{t-p}) + u_{2t} .
\end{split}
\end{equation}
where $X_t \in \mathcal{X} \subseteq \mathbb{R}$ and $Y_t \in \mathcal{Y} \subseteq \mathbb{R}^{d_Y}$ are scalar and $d_Y$-dimensional time series, respectively, $u_t = (u_{1t}, u_{2t}')' \in \mathcal{U} \subseteq \mathbb{R}^d$ are innovations, $d = 1 + d_Y$, $G_2 : \mathbb{R}^{1 + pd} \to \mathbb{R}$ is a generic nonlinear map, and $A_{12}(L)$ and $A_{11}(L)$ are lag polynomials \citep{lutkepohlNewIntroductionMultiple2005}. We let $Z_t := (X_t, Y_t')' \in \mathbb{R}^{d}$ be the full data vector.
Let us provide some examples for the model classes nested by \eqref{eq_main:general_nonlin_model}.
\begin{example}[Linear VAR]
In the simplest case, $G_2(Y_{t-1}, \ldots, Y_{t-p}, \allowbreak X_t, X_{t-1}, \ldots, \allowbreak X_{t-p}) = A_{22}(L) Y_{t-1} + A_{21}(L) X_{t}$, and we recover the class of linear vector autoregressive models.
\end{example}
\begin{example}[Additively separable model]
When $G_2(Y_{t-1}, \ldots, Y_{t-p}, \allowbreak X_t, X_{t-1}, \ldots, \allowbreak X_{t-p}) = \sum_{i=1}^p G_{i,22}(Y_{t-i}) + \sum_{j=0}^p G_{j,21}(X_{t-j})$, model \eqref{eq_main:general_nonlin_model} is additively separable \citep{fan2003nonlinear}.
\end{example}
\begin{example}[Nonlinear impact model]\label{example:nonlin_impact}
A parsimonious semiparametric class, which may be informally termed the ``{nonlinear impact model class}'', involves specification
\begin{equation*}
\begin{split}
Y_t & = \mu_2 + A_{22}(L) Y_{t-1} + \sum_{j=0}^p G_{j,21}(X_{t-j}) + u_{2t} ,
\end{split}
\end{equation*}
see e.g. \cite{goncalvesImpulseResponseAnalysis2021}.
An equivalent representation for $Y_t$ is
\begin{equation*}
Y_t = \mu_2 + A_{22}(L) Y_{t-1} + A_{21}(L) X_{t-1} + \sum_{j=0}^p \mathpalette\arc@arc{G}_{j,21}(X_{t-j}) + u_{2t} ,
\end{equation*}
where now to identify nonlinear functions $\mathpalette\arc@arc{G}_{j,21} : \mathbb{R} \to \mathbb{R}^{d_Y}$, $0 \leq j \leq p$, we require that constant and linear factors be not included at indices $j \geq 1$. To make this more compact, write $Z_t = \mu + A(L) Z_{t-1} + \mathpalette\arc@arc{G}(L) X_t + u_t$, where
\begin{equation*}
A(L) :=
\begin{bmatrix}
A_{11}(L) & A_{12}(L) \\
A_{21}(L) & A_{22}(L)
\end{bmatrix}
\quad\textnormal{and}\quad
\mathpalette\arc@arc{G}(L) :=
\begin{bmatrix}
0 \\
\mathpalette\arc@arc{G}_{0,21} + \mathpalette\arc@arc{G}_{1,21} L + \ldots + \mathpalette\arc@arc{G}_{p,21} L^p
\end{bmatrix} ,
\end{equation*}
with the minor abuse of notation that $\mathpalette\arc@arc{G}_2(L) := \mathpalette\arc@arc{G}_{0,21} + \allowbreak \ldots + \mathpalette\arc@arc{G}_{p,21} L^p$ is now intended as a \textit{functional} lag polynomial, meaning $\mathpalette\arc@arc{G}_2(L) X_t \equiv \sum_{j=0}^p \mathpalette\arc@arc{G}_{j,21}(X_{t-j})$.\footnote{The choice to use a functional matrix notation is due to the ease of writing multivariate additive nonlinear models such as \eqref{eq_main:structural_model} in a manner consistent with standard formalisms of linear VAR models, following again e.g. \cite{lutkepohlNewIntroductionMultiple2005}.}
\end{example}
\subsection{Structural Framework}
Model \eqref{eq_main:general_nonlin_model} involves only reduced-form innovations $u_{1t}$ and $u_{2t}$, and additional assumptions are necessary to provide a structural interpretation. Many such assumptions have been devised in the macroeconomic literature, but few can be directly applied to nonlinear models \citep{kilianStructuralVectorAutoregressive2017}. We follow the block-recursive identification strategy outlined in \cite{goncalvesImpulseResponseAnalysis2021} and originally due to \cite{kilian2011responses}.
From \eqref{eq_main:general_nonlin_model} we derive
\begin{equation*}
\begin{split}
X_t & = \mu_1 + A_{12}(L) Y_{t-1} + A_{11}(L) X_{t-1} + u_{1t}, \\
Y_t & = \mu_2 + A_{22}(L) Y_{t-1} + A_{21}(L) X_{t-1} + \mathpalette\arc@arc{G}_2(Y_{t-1:t-p}, X_{t:t-p}) + u_{2t} ,
\end{split}
\end{equation*}
where, without loss of generality, we have assumed (as in Example~\ref{example:nonlin_impact}) that we can separate the linear and nonlinear ($\mathpalette\arc@arc{G}_2$) components from $G_2$. In general, it can be the case that $\mu_2 = 0$, $A_{22}(L) = 0$ or $A_{21}(L) = 0$ if e.g. $G_2$ is strictly nonlinear.
In vector form:
\begin{equation}\label{eq_main:pseudo_reduced_form_single_eq}
Z_t = \mu + A(L) Z_{t-1} + \mathpalette\arc@arc{G}(Z_{t:t-p}) + u_t ,
\quad\textnormal{where}\quad
\mathpalette\arc@arc{G}(Z_{t:t-p}) :=
\begin{bmatrix}
0 \\
\mathpalette\arc@arc{G}_2(Y_{t-1:t-p}, X_{t:t-p})
\end{bmatrix}
.
\end{equation}
We can now formalize the structural specification of our model.
\begin{assumption}\label{assumption:structural_model}
There exist (i) a vector $B_0^{21} \in \mathbb{R}^{d_Y}$ and a matrix $B_0^{22} \in \mathbb{R}^{d_Y \times d_Y}$ such that
\begin{equation*}
\begin{bmatrix}
1 & 0 \\
B_0^{21} & B_0^{22}
\end{bmatrix}
=:
B_0^{-1}
\end{equation*}
is invertible and has unit diagonal, and (ii) mutually independent innovations sequences $\{\epsilon_{1t}\}_{t \in \mathbb{Z}}$, $\epsilon_{1t} \in \mathcal{E}_1 \subseteq \mathbb{R}$, and $\{\epsilon_{2t}\}_{t \in \mathbb{Z}}$, $\epsilon_{2t} \in \mathcal{E}_2 \subseteq \mathbb{R}^{d_Y}$, such that
\begin{equation*}
\begin{bmatrix}
\epsilon_{1t} \\
\epsilon_{2t}
\end{bmatrix}
\:\overset{\text{i.i.d.}}{\sim}\:
\left(
\begin{bmatrix}
0 \\
0
\end{bmatrix} ,
\begin{bmatrix}
\Sigma_1 & 0 \\
0 & \Sigma_2
\end{bmatrix}
\right) ,
\end{equation*}
where $\Sigma_1 > 0$ and $\Sigma_2$ is a diagonal positive definite matrix so that
\begin{equation}\label{eq_main:pseudo_reduced_form_structural_model}
\begin{split}
X_t & = \mu_1 + A_{12}(L) Y_{t-1} + A_{11}(L) X_{t-1} + \epsilon_{1t}, \\
Y_t & = \mu_2 + A_{22}(L) Y_{t-1} + A_{21}(L) X_{t-1} + \mathpalette\arc@arc{G}_2(Y_{t-1:t-p}, X_{t:t-p}) + B_0^{21} \epsilon_{1t} + B_0^{22} \epsilon_{2t} ,
\end{split}
\end{equation}
where $u_{1t} \equiv \epsilon_{1t}$, $u_{2t} := B_0^{21} \epsilon_{1t} + B_0^{22} \epsilon_{2t}$ and thus $u_t = B_0^{-1}\epsilon_t$ for $\epsilon_t = (\epsilon_{1t}, \epsilon_{2t}')' \in \mathcal{E} \subseteq \mathbb{R}^d$.
\end{assumption}
\begin{remark}
Assumption~\ref{assumption:structural_model} follows \cite{goncalvesImpulseResponseAnalysis2021} closely.
By design, one does not need to identify the model fully, meaning that fewer assumptions on $Z_t$ and $\epsilon_t$ are needed to estimate the individual structural effects of $\epsilon_{1t}$ on $Y_t$. This comes at the price of not being able to simultaneously study structural effects for shocks impacting $\epsilon_{2t}$.
\end{remark}
Note that inverting $B_0^{-1}$ gives
\begin{equation*}
B_0 = \begin{bmatrix}
1 & 0 \\
-B_{0,12} & B_{0,22}
\end{bmatrix} ,
\end{equation*}
and, multiplying both sides of \eqref{eq_main:pseudo_reduced_form_single_eq} by $B_0$, we find
\begin{equation}\label{eq_main:structural_model}
B_0 Z_t = b + B(L) Z_{t-1} + \mathpalette\arc@arc{F}(Z_{t:t-p}) + \epsilon_t ,
\end{equation}
where $b = (b_1, b_2')' \in \mathbb{R}^d$ and $\mathpalette\arc@arc{F}(Z_{t:t-p}) = (0, \mathpalette\arc@arc{F}_2(Z_{t:t-p}))'$ for $\mathpalette\arc@arc{F}_2 : \mathbb{R}^{1+pd_Y} \to \mathbb{R}^{d_Y}$, $F_2 = B_{0,22} \mathpalette\arc@arc{G}_2$.
In practice, to estimate the model's coefficients, we will leverage \eqref{eq_main:pseudo_reduced_form_structural_model}. This latter form was termed the \textit{pseudo-reduced form} by \cite{goncalvesImpulseResponseAnalysis2021}.
Observe that $\mathpalette\arc@arc{G}_2(Y_{t-1:t-p}, X_{t:t-p})$ is correlated with $u_{2t}$ through $B_0^{21} \epsilon_{1t}$. As $X_t$ depends linearly on $\epsilon_{1t}$, if $B_0^{21} \not= 0$ and $\mathpalette\arc@arc{G}_2(Y_{t-1:t-p}, X_{t:t-p})$ is not independent of $X_t$, there is an endogeneity problem.
\cite{goncalvesImpulseResponseAnalysis2021} address the issue by proposing a two-step estimation procedure wherein one proxies for $\epsilon_{1t}$ with residual $\widehat{\epsilon}_{1t}$. As we prove in Section~\ref{section:estimation} below, this approach also generally allows for consistent semiparametric estimation.
\begin{remark}{(Moving Average Identification).}
\cite{forniNonlinearTransmissionFinancial2023,forniAsymmetricEffectsNews2023} work with an alternative nonlinear structural identification framework to the block-recursive form. Their approach follows \cite{debortoliAsymmetricEffectsMonetary2020} and is based on a vector MA representation. Under appropriate assumptions, the structural model studied by \cite{forniNonlinearTransmissionFinancial2023} is
\begin{equation}\label{eq_main:debortoli_forni_structural_form}
Z_t = \mu + A(L) Z_t + Q_0 F(\epsilon_{1t}) + B_0 \epsilon_t ,
\end{equation}
where $\epsilon_t$ are independent innovations with zero mean and identity covariance, and $\epsilon_{1t}$ identifies the shocks of interest. $Q(L)$ and $B(L)$ are both linear lag polynomials, and $F(x) = x^2$ in their baseline specification.
For \eqref{eq_main:debortoli_forni_structural_form} to overlap with \eqref{eq_main:pseudo_reduced_form_structural_model} one must impose that (i) $X_t$ is exogenous and independently distributed and (ii) only $\epsilon_{1t}$ has nonlinear effects.
We emphasize that, if the innovation sequence $\epsilon_{1t}$ is assumed to be observable, applying our results to the framework of \cite{debortoliAsymmetricEffectsMonetary2020} is straightforward.
\end{remark}
\subsection{Structural Nonlinear Impulse Responses}\label{section:nonlin_irfs}
Starting from pseudo-reduced equations \eqref{eq_main:pseudo_reduced_form_structural_model}, we begin by assuming that the linear autoregressive component is stable.
\begin{assumption}\label{assumption:roots_linear_part_model}
The roots of $\det(I_{d} - A(L) L) = 0$ are outside the complex unit circle.
\end{assumption}
This standard stability assumption enables us to write impulse responses in a manner that can yield useful simplifications for additively separable models.\footnote{Stability of the linear VAR component is neither necessary nor sufficient for ensuring stability and stationarity of the entire nonlinear process, cf. Assumption~\ref{assumption:physical_dep} in Section~\ref{section:estimation} below.}
Then, letting $\Psi(L) = (I_{d} - A(L) L)^{-1}$, one can write
\begin{equation}
Z_t = \eta + \Theta(L) \epsilon_t + \Gamma(Z_{t:*}) ,
\end{equation}
where
$\eta := \Psi(1) (\mu_1, \mu_2')'$,
$\Theta(L) := \Psi(L) B_0^{-1}$ and
$\Gamma(Z_{t:*}) := \Psi(L)(0, \mathpalette\arc@arc{G}_2(Y_{t-1:t-p}, X_{t:t-p})')'$.
We emphasize that the nonlinear term $\Gamma(Z_{t:*})$ generally depends on the entire past of the process $Z_t$, as $\Psi(L)$ is an infinite-order MA polynomial.
To define impulse responses, we partition the polynomial $\Theta(L)$ according to
$\Theta(L) := [\Theta_{\cdot 1}(L) \:\vert\: \Theta_{\cdot 2}(L)]$ ,
where $\Theta_{\cdot 1}(L)$ represents the first column of matrices in $\Theta(L)$, and $\Theta_{\cdot 2}(L)$ the remaining $d_Y$ columns.
Given impulse $\delta \in \mathbb{R}$ at time $t$, define the shocked innovation process as $\epsilon_{1 s}(\delta) = \epsilon_s$ for $s \not= t$ and $\epsilon_{1 t}(\delta) = \epsilon_{1 t} + \delta$, as well as the shocked structural variable as $Z_s(\delta) = Z_s$ for $s < t$ and $Z_s(\delta) = X_s(\epsilon_{s:t+1}, \epsilon_t + \delta, \epsilon_{t-1:*})$ for $s \geq t$. Further, for a given horizon $h \geq 0$, let
\begin{align*}
Z_{t+h} & := \eta + \Theta_{\cdot 1}(L) \epsilon_{1 t+h} + \Theta_{\cdot 2}(L) \epsilon_{2 t+h} + \Gamma(Z_{t:*}) , \\[2pt]
Z_{t+h}(\delta) & := \eta + \Theta_{\cdot 1}(L) \epsilon_{1 t+h}(\delta) + \Theta_{\cdot 2}(L) \epsilon_{2 t+h} + \Gamma(Z_{t:*}(\delta)) ,
\end{align*}
be the time-$t$ baseline and shocked series, respectively.
Then,
\begin{equation}\label{eq_main:def_irf}
\textnormal{IRF}_{h}(\delta) = \mathbb{E}\left[ Z_{t+h}(\delta) - Z_{t+h} \right]
\end{equation}
is the unconditional impulse response at horizon $h$ due to shock $\delta$. The difference between series is
$ Z_{t+h}(\delta) - Z_{t+h}
= \Theta_{h,\cdot 1} \delta + \Gamma(Z_{t:*}(\delta)) - \Gamma(Z_{t:*}) $,
hence
\begin{equation}\label{eq_main:irf_h}
\textnormal{IRF}_{h}(\delta)
= \Theta_{h,\cdot 1} \delta + \mathbb{E}\left[ \Gamma(Z_{t:*}(\delta)) - \Gamma(Z_{t:*}) \right] .
\end{equation}
\begin{remark}
In additively separable models, $\Gamma(Z_{t:*})$ is also additively separable over lags of $Z_t$. Accordingly, the baseline and shock series have an additive form,
as terms with time indices $s < t$ remain unaffected by the shock.
Therefore, \eqref{eq_main:irf_h} reduces to
\begin{equation}\label{eq_main:irf_h_separable}
\textnormal{IRF}_{h}(\delta)
= \Theta_{h,\cdot 1} \delta + \mathbb{E}\left[ \Gamma_0(Z_{t+h}(\delta)) - \Gamma_0(Z_{t+h}) \right] + \ldots + \mathbb{E}\left[ \Gamma_h(Z_{t}(\delta)) - \Gamma_h(Z_{t}) \right] .
\end{equation}
Coefficients $\Gamma_j$ are still functional, and cannot be collected across $X_{t+j}(\delta)$ and $X_{t+j}$.
\end{remark}
Closed-form computation of nonlinear IRFs is highly non-trivial. Even in the separable case \eqref{eq_main:irf_h_separable}, while one can linearly separate expectations in the impulse response formula, terms $\mathbb{E}\left[ \Gamma_j(Z_{t+j}(\delta)) - \Gamma_j(Z_{t+j}) \right]$ for $0 \leq j \leq h$ cannot be meaningfully simplified further. Moreover, these expectations involve nonlinear functions of lags of $Z_t$ and are impractical to derive explicitly.
To avoid working with $\Theta(L)$ and $\Gamma(L)$, we now present an iterative algorithm which allows one to easily and efficiently compute nonlinear IRFs.\footnote{The algorithm we propose is a natural counterpart to the one in Proposition 3.1 of \cite{goncalvesImpulseResponseAnalysis2021}, wherein they suggest to estimate the MA form coefficients recursively. Our approach instead relies on directly iterating forward the model's equations, which is more computationally straightforward.}
\begin{proposition}\label{prop:irf_iterate_algorithm}
For any $h = 0, 1, \ldots, H$, with $H \geq 1$ fixed, if impulse response $\textnormal{IRF}_{h}(\delta)$ is finite and well-defined, it can be computed with the following steps:
\begin{description}
\item[($\text{i}$)] For $j = 0$, let $X_t(\delta) = X_t + \delta$ and
$Y_{t}(\delta) = \mu_2 + G_2(Y_{t-1}, \ldots, Y_{t-p}, X_t(\delta), X_{t-1}, \ldots, X_{t-p}) + B_0^{21} (\epsilon_{1t} + \delta) + \xi_{2t}$, where $\xi_{2t} = B_0^{22} \epsilon_{2t}$.
\item[($\text{ii}$)] For $j = 1, \ldots, h$, let
\begin{equation*}
\begin{split}
X_{t+j}(\delta) & = \mu_1 + A_{12}(L) Y_{t+j-1}(\delta) + A_{11}(L) X_{t+j-1}(\delta) + \epsilon_{1t+j} , \\
Y_{t+j}(\delta) & = \mu_2 + G_2(Y_{t-1}(\delta), \ldots, Y_{t-p}(\delta), X_t(\delta), X_{t-1}(\delta), \ldots, X_{t-p}(\delta)) + B_0^{21} \epsilon_{1t+j} + \xi_{2t+j} .
\end{split}
\end{equation*}
where $X_{t}(\delta)$ and $Y_{t}(\delta)$ are the shocked sequences determined by forward iteration after time $t$, equaling baseline sequences $X_{t}$ and $Y_{t}$ at lags before $t$, respectively.
\end{description}
Setting $Z_{t+j}(\delta) = ( X_t(\delta), Y_t(\delta) )'$, it holds $\textnormal{IRF}_h(\delta) = \mathbb{E}[ Z_{t+j}(\delta) - Z_{t+j} ]$.
\end{proposition}
Proposition~\ref{prop:irf_iterate_algorithm} follows directly from the definition of the unconditional impulse response \eqref{eq_main:def_irf} combined with a direct forward iteration of \eqref{eq_main:pseudo_reduced_form_structural_model}, sidestepping the explicit MA($\infty$) formulation in \eqref{eq_main:irf_h}. This approach dispenses from the need to simulate innovations $\{\epsilon_{t+j}\}_{j=1}^{h-1}$, as the joint distribution of $\{X_{t+h-1}, X_{t+j-1}, \ldots, X_{t}\}$ contains all relevant path information.
When the model is estimated from data, for residuals $\widehat{\epsilon}_{1t}$ and $\widehat{\xi}_{2t}$ it trivially holds
\begin{equation*}
\begin{split}
X_{t} & = \widehat{\mu}_1 + \widehat{A}_{12}(L) Y_{t-1} + \widehat{A}_{11}(L) X_{t-1} + \widehat{\epsilon}_{1t}, \\
Y_{t} & = \widehat{\mu}_2 + \widehat{G}_2(Y_{t-1}, \ldots, Y_{t-p}, X_t, X_{t-1}, \ldots, X_{t-p}) + \widehat{B}_0^{21} \widehat{\epsilon}_{1t} + \widehat{\xi}_{2t} .
\end{split}
\end{equation*}
In practice, this means that one can numerically construct the shocked sequence as
\begin{equation*}
\begin{split}
\widehat{X}_{t+j}(\delta)
& = \widehat{\mu}_1 + \widehat{A}_{12}(L) \widehat{Y}_{t+j-1}(\delta) + \widehat{A}_{11}(L) \widehat{X}_{t+j-1}(\delta) + \widehat{\epsilon}_{1t+j}, \\
\widehat{Y}_{t+j}(\delta)
& = \widehat{\mu}_2 + \widehat{G}_2(\widehat{Y}_{t-1}(\delta), \ldots, \widehat{Y}_{t-p}(\delta), \widehat{X}_t(\delta), \widehat{X}_{t-1}(\delta), \ldots, \widehat{X}_{t-p}(\delta)) + \widehat{B}_0^{21} \widehat{\epsilon}_{1t+j} + \widehat{\xi}_{2t+j} ,
\end{split}
\end{equation*}
for $j = 1, \ldots, h$ where $\widehat{X}_t(\delta) = X_t + \delta$, $\widehat{X}_{t-s} = X_{t-s}$ for all $s \geq 1$, and similarly for $\widehat{Y}_t(\delta)$.
\section{Estimation}\label{section:estimation}
To discuss estimation, we will rewrite the equations in \eqref{eq_main:pseudo_reduced_form_structural_model} with some minor reordering as
\begin{equation}\label{eq_main:regression_model}
\begin{split}
X_t & = \Pi_1' W_{1t} + \epsilon_{1t} , \\
Y_t & = \Pi_2' W_{2t} + \xi_{2t} ,
\end{split}
\end{equation}
where
$\xi_{2t} = B_0^{22} \epsilon_{2t}$,
$\Pi_1 := ( \eta_1, A_{1,11}, \cdots, A_{p,11}, A_{1,12}', \cdots, A_{p,12}' )' \in \mathbb{R}^{1 + p d}$,
\begin{equation*}
\Pi_{2}
:=
\begin{bmatrix}[c|c|c]
& G_{1,2}(\cdot) & \\
\mu_2 & \cdots & B_0^{21}\\
& G_{d_Y,2}(\cdot) &
\end{bmatrix}' : \mathbb{R}^{2 + p d} \to \mathbb{R}^{d_Y} ,
\end{equation*}
$W_{1t} := ( 1, X_{t-1}, \ldots, X_{t-p}, Y_{t-1}', \ldots, Y_{t-p}' )' \in \mathbb{R}^{1 + p d}$, and
$W_{2t} := (1, X_t, \allowbreak X_{t-1}, \ldots, X_{t-p}, \allowbreak Y_{t-1}', \ldots, \allowbreak Y_{t-p}', \allowbreak \epsilon_{1t} )' \in \mathbb{R}^{2 + p d}$.
With a slight abuse of notation, similar that of Example~\ref{example:nonlin_impact}, we have written the functional terms in $\Pi_2$ as a ``vector product'', $G_2 \cdot (X_{t:t-p}', Y_{t-1:t-p}')' \equiv G_2(X_{t:t-p}, Y_{t-1:t-p})$, where $G_2$ is a vector of functions,
one for each component of $Y_t$.
Whenever $\Pi_1 \not= 0$, $W_{2t}$ is an infeasible vector of regressors due to term ${\epsilon}_{1t}$. To estimate $\Pi_2$, one can use $\widehat{W}_{2t} = (1, X_t, \allowbreak X_{t-1}, \ldots, X_{t-p}, \allowbreak Y_{t-1}', \ldots, Y_{t-p}', \allowbreak \widehat{\epsilon}_{1t})'$ instead, which contains generated regressors in the form of residual $\widehat{\epsilon}_{1t}$.
A valid two-step estimation procedure \citep{goncalvesImpulseResponseAnalysis2021} is:
(I) Regress $X_t$ on $W_{1t}$ to get estimate $\widehat{\Pi}_1$ and residuals $\widehat{\epsilon}_{1t} = X_t - \widehat{\Pi}_1' W_{1t}$;
(II) Semiparametrically regress $Y_t$ on $\widehat{W}_{2t}$ to get estimate $\widehat{\Pi}_2$.
There are many ways to implement step (II), given that the literature on non- and semiparametric regression is mature.
We rely on the sieve framework of \cite{chenOptimalUniformConvergence2015} as the workhorse to derive the main theoretical results. The sieve framework is quite rich, encompassing e.g. neural networks \citep{chenImprovedRatesAsymptotic1999,shenAsymptoticPropertiesNeural2023}.
\subsection{Semiparametric Series Estimation}
The semiparametric regression step we require is more readily analyzed by working on each component of $Y_t$. For $i \in \{1, \ldots, d_Y\}$, consider
\begin{equation}\label{eq_main:regression_eq_i}
Y_{t,i} = \mu_{2,i} + G_{2,i}(Y_{t-1}, \ldots, Y_{t-p}, X_t, X_{t-1}, \ldots, X_{t-p}) + B^{21}_{0,i} \epsilon_{1t} + \xi_{2t,i} .
\end{equation}
Let then $\pi_{2,i} := [ \mu_{2,i}, \: G_{2,i}, \: B^{21}_{0,i} ]'$. The regression equation for $\pi_{2,i}$ is thus $Y_{i} = \pi_{2,i}' W_{2} + \xi_{2i}$, where $Y_{i} = (Y_{1,i}, \ldots, Y_{n,i})'$ and $\xi_{2i} = (\xi_{2t,1}, \ldots, \xi_{2t,n})'$. The estimation target is the conditional expectation $\pi_{2,i}(w) = \mathbb{E}[ Y_{t,i} \:\vert\: W_{2t} = w ]$ under the assumption $\mathbb{E}[ \xi_{2t,i} \:\vert\: W_{2t} ] = 0$.
Assume that $G_{2,i} \in \Lambda$, where $\Lambda$ is a sufficiently regular function class to be specified in the following.
Given a collection $b_{1\kappa}, \ldots, b_{\kappa\kappa}$ of $\kappa \geq 1$ basis functions belonging to sieve $\mathcal{B}_\kappa$, define $b^\kappa(\cdot) := \left( b_{1\kappa}(\cdot), \ldots, b_{\kappa\kappa}(\cdot) \right)'$ and
$
B_\kappa := \big(
b^\kappa(Y_{0:1-p}, X_{1:1-p}),
\ldots,
b^\kappa(Y_{n-1:n-p}, X_{n:n-p})
\big)'
$.
For univariate functions, one can directly apply spline, wavelet and Fourier sieves; in the multivariate case, tensor-product sieves are straightforward generalizations \citep{chenOptimalUniformConvergence2015}.
To construct the final semiparametric sieve for $\pi_{2,i}$, indicated by $\mathcal{B}_\pi$, let $b_{\pi,1K}, \ldots, b_{\pi,KK}$ be the sieve basis in $\mathbb{R} \times \mathcal{B}_\kappa \times \mathbb{R}$ for $\kappa \geq 1$ and $K = 2 + \kappa$ given by
$b_{\pi,1 K}(W_{2t}) = 1$,
$b_{\pi,\ell K}(W_{2t}) = b_{\ell \kappa}(Y_{t-1:t-p}, X_{t:t-p})$, for $2 \leq \ell \leq \kappa+1$, and
$b_{\pi,K K}(W_{2t}) = \epsilon_{1t}$.
Note that $K$, the overall size of the sieve, grows linearly in $\kappa$, which itself controls the effective dimension of the nonparametric component of the sieve, $b_{\pi,2 K}, \ldots, b_{\pi,(\kappa+1) K}$.
Introducing $b^K_\pi(w) := ( b_{\pi,1K}(w), \ldots, b_{\pi,KK}(w) )'$ and $B_\pi := ( b^K_\pi(W_{21}), \ldots, b^K_\pi(W_{2n}) )'$, the generally \textit{infeasible} least squares series estimator $\widehat{\pi}_{2,i}^*(w)$ is given by $\widehat{\pi}^*_{2,i}(w) = b^K_\pi(w)' ({B}_\pi' {B}_\pi)^{-1} {B}_K' Y_i$.
Similarly, the feasible series regression matrix $\widehat{B}_\pi := ( b^K_\pi(\widehat{W}_{21}), \ldots, b^K_\pi(\widehat{W}_{2n}) )'$ yields the \textit{feasible} least squares series estimator, $\widehat{\pi}_{2,i}(w) = b^K_\pi(w)' (\widehat{B}_\pi' \widehat{B}_\pi)^{-1} \widehat{B}_K' Y_i$.
To further streamline notation, wherever it does not lead to confusion, we will let $\pi_2$ be a generic coefficient vector belonging to $\{\pi_{2,i}\}_{i=1}^p$, as well as define $\widehat{\pi}_{2}$, $Y$ and $u_2$ associated to the same regression equation.
\subsection{Distributional and Sieve Assumptions}
To derive asymptotic consistency results, we begin by stating conditions on the basic probability structure of the model.
\begin{assumption}\label{assumption:stationarity}
$\{Z_{t}\}_{t \in \mathbb{Z}}$ is a strictly stationary and ergodic time series.
\end{assumption}
\begin{assumption}\label{assumption:compactness}
$X_{t} \in \mathcal{X} \subset \mathbb{R}$, $Y_{t} \in \mathcal{Y} \subset \mathbb{R}^{d_Y}$ and $\epsilon_t \in \mathcal{E} \subset \mathbb{R}^{d}$ for all $t \in \mathbb{Z}$, where $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{E}$ are compact, convex sets with nonempty interior.
\end{assumption}
Assumption~\ref{assumption:stationarity} follows both \cite{goncalvesImpulseResponseAnalysis2021} and \cite{chenOptimalUniformConvergence2015}.
Note that, as $W_{2t}$ depends only on $X_{t:t-p}$, $Y_{t-1:t-p}$ and $\epsilon_{1t}$, the entries of $\xi_{2t}$ in \eqref{eq_main:regression_model} are independent of $W_{2t}$, so that $\mathbb{E}[ u_{2it} \:\vert\: W_{2t} ] = 0$.
Assumption~\ref{assumption:compactness} implies that $X_t$, $Y_t$, as well as $\epsilon_t$ are bounded random variables. In (semi-)nonparametric estimation, imposing that $X_t$ is bounded almost surely is a standard assumption. Since lags of $Y_t$ and innovations $\epsilon_t$ contribute linearly to all components of $Z_t$, it follows that they too must be bounded.
In practice, Assumption~\ref{assumption:compactness} is not particularly restrictive, as many credibly stationary economic series often have reasonable implicit (e.g., inflation) or explicit bounds (e.g., employment rate).
The analysis of impulse responses on compact domains is, however, non-trivial. In Section~\ref{section:relaxed_shocks} below, we provide an IRF shock relaxation framework that can accommodate this setting.
\begin{remark}\label{remark:compactness_assumption}
Bounded support assumptions are uncommon in time series econometrics, as boundedness is not necessary in the analysis of linear models \citep{hamilton1994state,lutkepohlNewIntroductionMultiple2005,kilianStructuralVectorAutoregressive2017,stock2016dynamic}.
Unbounded regressors are significantly more complex to handle when working in the nonparametric setting.
\cite{chenOptimalUniformConvergence2015} do work in weighted sup-norms, but their uniform results are stated only under a compact domain assumption.
Avoiding Assumption~\ref{assumption:compactness} can be achieved by changing the model's equations -- e.g., the lags of $Y_t$ only affect $X_t$ via bounded functions -- but this further restricts the model.
Establishing a general (uniform) theory of nonparametric regressions with unbounded data domains, on the other hand, is a complex question. For kernel, partitioning and nearest-neighbor methods and i.i.d. data, a handful of papers develop results in $L^1$ and $L^2$ norms, see \cite{kohlerRatesConvergencePartitioning2006,kohlerOptimalGlobalRates2009} and \cite{kohlerOptimalGlobalRates2013}. For wavelet estimators in the i.i.d. regression setting, \cite{zhouUniformConvergenceRates2022} provided the first sup-norm result in Besov spaces with suboptimal rates.
\cite{hansenUNIFORMCONVERGENCERATES2008} is, to the best of our knowledge, the only work providing convergence rates for local constant and local linear regression estimators in a dependent data setting without bounded support restrictions. Yet, even in the nonparametric LP setup of \cite{goncalvesNonparametricLocalProjections2024}, the authors argue that it is not clear if these results allow for IRF estimation guarantees over $\mathbb{R}$.
Construction of a comprehensive nonparametric framework to handle non-independent, unbounded data should thus be considered an important objective of future research.
\end{remark}
Without loss of generality, let $\mathcal{Y} = [0,1]^{d_Y}$ and $\mathcal{X} = [0,1]$.
\begin{assumption}\label{assumption:regressor_density}
The unconditional densities of $Y_t$ and $X_{t}$ are uniformly bounded away from zero and infinity over $\mathcal{Y}$ and $\mathcal{X}$, respectively.
\end{assumption}
\begin{assumption}\label{assumption:function_class}
For all $1 \leq i \leq d_Y$ the restriction of $G_{2,i}$ to $\mathcal{Y}^{\, p} \times \mathcal{X}^{1+p} \equiv [0,1]^{1+pd}$ belongs to the Hölder class $\Lambda^s([0,1]^{1+pd})$ of smoothness $s \geq 1$.
\end{assumption}
Assumptions \ref{assumption:regressor_density} and \ref{assumption:function_class} are classical in the nonparametric regression literature.
Let then $\mathcal{W}_2 \subset \mathbb{R}^d$ be the domain of $W_{2t}$. By assumption, $\mathcal{W}_2$ is compact and convex and is given by the direct product
$
\mathcal{W}_2 = \{1\} \times \mathcal{Y}^{\, p} \times \mathcal{X}^{1+p} \times \mathcal{E}_1
$,
where $\mathcal{E}_1$ is the domain of structural innovations $\epsilon_{1t}$ i.e. $\mathcal{E} \equiv \mathcal{E}_1 \times \mathcal{E}_2$.
\begin{assumption}\label{assumption:sieve_regularity}
Define $\zeta_{K,n} := \sup_{w \in \mathcal{W}_2} \lVert b^K_\pi(w) \rVert$ and $\lambda_{K,n} := [ \lambda_{\min}(\mathbb{E}[\, b^K_\pi(W_{2t}) b^K_\pi(W_{2t})' \,]) ]^{-1/2}$.
It holds: (i) there exist $\omega_1, \omega_2 \geq 0$ s.t. $ \sup_{w \in \mathcal{W}_2} \lVert \nabla b^K_\pi(w) \rVert \lesssim n^{\omega_1} K^{\omega_2}$;
(ii) there exist $\overline{\omega}_1 \geq 0$, $ \overline{\omega}_2 > 0$ s.t. $ \zeta_{K,n} \lesssim n^{\overline{\omega}_1} K^{\overline{\omega}_2} $;
(iii) $\lambda_{\min}(\mathbb{E}[\, b^K(W_{2t}) b^K(W_{2t})' \,]) > 0$ for all $K$ and $n$.
\end{assumption}
Assumption \ref{assumption:sieve_regularity} provides mild regularity conditions on the families of sieves that can be used for the series estimator. More generally, letting $\mathcal{W}_2$ be compact and rectangular makes Assumptions~\ref{assumption:sieve_regularity}(i)-(ii) hold for commonly used basis functions \citep{chenOptimalUniformConvergence2015}.
The approximation properties of these sieves are well understood \citep{chenChapter76Large2007}.\footnote{See also \cite{chenPenalizedSieve2013,belloniNewAsymptoticTheory2015} for additional discussion and examples of sieve families.}
In particular, Assumption \ref{assumption:sieve_regularity}(i) holds with $\omega_1 = 0$ since the domain is fixed over the sample size.
What is also needed is that the nonparametric components of the sieve given by $b_{\pi,1K}, \ldots, b_{\pi,KK}$ are able to approximate $G_{2,i}$ with an error that decays sufficiently fast with $K$.
Lastly, Assumption~\ref{assumption:sieve_regularity}(iii) is a mild assumption on the conditioning of the semiparametric sieve.
\begin{assumption}\label{assumption:sieve_type}
Sieve $\mathcal{B}_\kappa$ belongs to $\text{BSpl}(\kappa, \mathcal{W}_2, r)$ or $\text{Wav}(\kappa, \mathcal{W}_2, r)$, the tensor B-spline and tensor wavelet sieve, respectively, of degree $r$ over $\mathcal{W}_2$,
with $r \geq \max\{ s, 1 \}$.
\end{assumption}
We define
$\widetilde{b}^K_\pi(w) := \mathbb{E}[\, {b}^K_\pi(W_{2t}) {b}^K_\pi(W_{2t})' \,]^{-1/2}\, {b}^K_\pi(w)$ and
$\widetilde{B}_\pi := \big( \widetilde{b}^K_\pi(W_{21}), \allowbreak \ldots, \allowbreak \widetilde{b}^K_\pi(W_{2n}) \big)'$
to be the orthonormalized vector of basis functions and the orthonormalized regression matrix, respectively.
To derive uniform converges rates under dependence, we require that the Gram matrix of orthonormalized sieve converges to the identity matrix.
\begin{assumption}\label{assumption:series_gram_matrix_convergence}
It holds that $\lVert (\widetilde{B}_\pi' \widetilde{B}_\pi / n) - I_K \rVert = o_P(1)$.
\end{assumption}
\cite{chenOptimalUniformConvergence2015} introduced Assumption~\ref{assumption:series_gram_matrix_convergence} as a key ingredient for their proofs, while also showing that it holds whenever $\{W_{2t}\}_{t\in\mathbb{Z}}$ is either an exponential or algebraic $\beta$-mixing process.
Unfortunately, mixing conditions are difficult to verify or test with respect to model specification, as they rely on bounding the worst-case ``independence gap'' between probability events (see Appendix~\ref{appendix:dependence}).
We extend their approach to the case of geometrically decaying physical dependence, a metric proposed by \cite{wuNonlinearSystemTheory2005}. This is a setting where many estimation and inference results have been derived, see for example \cite{wuKernelEstimationTime2010,wuAsymptoticTheoryStationary2011a,chenSelfnormalizedCramertypeModerate2016} and references within.
\begin{assumptionp}{\ref*{assumption:series_gram_matrix_convergence}$\,'$}\label{assumption:physical_dep}
Let $\{Z_t\}_{t \in \mathbb{Z}}$ be such that we can write $Z_{t+h} = \Phi^{(h)}(Z_t,\allowbreak \epsilon_{t+1:t+h})$ for some nonlinear maps $\Phi^{(h)}$ and innovations $\{\epsilon_t\}_{t \in \mathbb{Z}}$ over all $h \geq 1$. Then, for $r \geq 2$, there exists constants $a_1 > 0$, $a_2 > 0$ and $\tau \in (0,1]$ such that it holds
\begin{equation*}
\sup_t \big\lVert\, Z_{t+h} - \Phi^{(h)}(Z_t, \epsilon_{t+1:t+h}) \,\big\rVert_{L^r}
\leq
a_1 \exp(- a_2 \, h^\tau) .
\end{equation*}
\end{assumptionp}
Assumption~\ref{assumption:series_gram_matrix_convergence} is subsumed by Assumption~\ref{assumption:physical_dep}. Using a physical dependence measure, we argue that it is also possible to swap mixing conditions with more explicit, primitive conditions derived exclusively in terms of model specification \eqref{eq_main:general_nonlin_model}. In particular, for specific semiparametric model specifications, it is possible to verify Assumption~\ref{assumption:physical_dep} directly by leveraging stability/contractivity theory of dynamic systems.
We refer the reader to Appendix~\ref{appendix:dependence} for an additional, in-depth discussion of dependence and physical conditions.
\subsection{Uniform Convergence and Consistency}
We can now state our main result, which shows that the two-step estimation procedure for \eqref{eq_main:regression_model} provides consistent estimates.
\begin{theorem}\label{theorem:twostep_estimator_consistency}
Let $\{Z_t\}_{t \in \mathbb{Z}}$ be determined by structural model \eqref{eq_main:structural_model}. Under Assumptions \ref{assumption:structural_model}, \ref{assumption:stationarity}, \ref{assumption:compactness}, \ref{assumption:regressor_density}, \ref{assumption:function_class}, \ref{assumption:sieve_regularity}, \ref{assumption:sieve_type} and \ref{assumption:physical_dep}, let $\widehat{\Pi}_1$ and $\widehat{\Pi}_2$ be the least squares and two-step semiparametric series estimators for $\Pi_1$ and $\Pi_2$, respectively. Then,
$ \lVert \widehat{\Pi}_1 - \Pi_1 \rVert_\infty = O_P(n^{-1/2\,}) $
and
\begin{equation*}
\lVert \widehat{\Pi}_2 - \Pi_2 \rVert_\infty
\leq
O_P\left( \zeta_{K,n} \lambda_{K,n} \, \frac{K}{\sqrt{n}} \right) + \lVert \widehat{\Pi}^*_2 - \Pi_2 \rVert_\infty ,
\end{equation*}
where $\widehat{\Pi}^*_2$ is the infeasible series estimator involving $\epsilon_{1t}$.
\end{theorem}
Sup-norm bounds for $\lVert \widehat{\Pi}^*_2 - \Pi_2 \rVert_\infty$ may be obtained from Lemma~2.3 and Lemma~2.4 in \cite{chenOptimalUniformConvergence2015}.
Assuming $s \geq 1$ and $d = 1$, such as in the setting of the additively separable model in Example~\ref{example:nonlin_impact} and in our empirical applications, it is possible to show that, if the optimal nonparametric rate for $K$ is used and the additive sieve inherits the conditioning of the underlying sieve bases, then $\widehat{\Pi}_2$ is sup-norm consistent.
\begin{corollary}\label{corollary:twostep_estimator_op1}
Under the same assumptions as Theorem~\ref{theorem:twostep_estimator_consistency}, further assume that $s \geq 1$, $G_{2,i}$, $1 \leq i \leq d_Y$, in \eqref{eq_main:regression_eq_i} is additively separable in all its components and $\lambda_{K,n} \lesssim 1$. Then for the choice $K \asymp (n / \log(n))^{1/(2 s + 1)}$ it holds that
\begin{equation*}
\lVert \widehat{\Pi}_2 - \Pi_2 \rVert_\infty
=
O_P \left(
n^{- \frac{s-1}{2s+1}}
\log(n)^{-\frac{3}{2(2s+1)}}
\right)
\end{equation*}
and, in particular, $\lVert \widehat{\Pi}_2 - \Pi_2 \rVert_\infty = o_P(1)$.
\end{corollary}
The requirement $\lambda_{K,n} \lesssim 1$ for additively separable sieves is mild given the known properties of B-spline and wavelet sieves, although nontrivial.
Since one cannot exploit sparsity in the case of non-locally supported bases, as is the case with linearly separable sieves, we assume $\lambda_{K,n}$ is upper bounded by a constant to streamline the analysis of the empirical sieve projection operator (see also the discussion in \citealp{huangLocalAsymptoticsPolynomial2003a}, Section 7).
\begin{remark}
Several methods can be used to select $K$ in practice: Cross-validation, generalized cross-validation, Mallow's criterion, and others \citep{li2009nonparametric}. In the case of piecewise splines, once size is selected, knots can be chosen to be the $K$ uniform quantiles of the data. In our simulations and applications, for simplicity, we select sieve sizes manually, while knots are located following empirical quantiles.
\end{remark}
\section{Impulse Response Analysis}\label{section:impulse_response_analysis}
After discussing the estimation of the structural model's coefficients, we can now address the derivation of nonlinear impulse responses.
To ensure compatibility with bounded support assumptions, we introduce an extension of the classical IRF definition, termed \textit{relaxed impulse response function}, which differs only in the form of the shock applied to the model.
We then show that nonlinear relaxed IRFs can be consistently estimated, and uniformly so for shocks picked within a compact range.
\subsection{Relaxed Shocks}\label{section:relaxed_shocks}
Under Assumptions~\ref{assumption:compactness} and \ref{assumption:regressor_density}, the standard construction of impulse responses following Section~\ref{section:nonlin_irfs} is, unfortunately, improper. This is immediately seen by noticing that, at impact, $X_t(\delta) = X_{t} + \delta$, meaning that $\mathbb{P}( X_t(\delta) \not\in \mathcal{X} ) > 0$ since there is a translation of size $\delta$ in the support of $X_t$.
To address this problem, we introduce an extension to the standard additive shock that is used to define impulse responses.
We begin by defining mean-shift shocks, that is, shocks such that the distribution of time $t$ innovations is shifted to have mean $\delta$, while retaining compact support almost surely.
\begin{definition}
Let $\mathcal{E}_1 \subseteq \mathbb{R}$ and $\mathbb{P}(\epsilon_{1t} \in \mathcal{E}_1) = 1$. A {mean-shift structural shock} $\epsilon_{1t}(\delta)$ is an appropriately chosen transformation of $\epsilon_{1t}$ such that $\mathbb{P}(\epsilon_{1t}(\delta) \in \mathcal{E}_1) = 1$ and $\mathbb{E}[\epsilon_{1t}(\delta)] = \delta$.
\end{definition}
With a mean-shift shock, at impact it holds $X_{t}(\delta) = X_{t} + (\epsilon_{1t}(\delta) - \epsilon_{1t})$. In the standard setting, where $\mathbb{E}[\epsilon_t] = 0$ and $\mathcal{E}_1 \equiv \mathbb{R}$, $\epsilon_{1t}(\delta) = \epsilon_{1t} + \delta$ is clearly valid. More generally, however, imposing $\mathbb{E}[\epsilon_{1t}(\delta)] = \delta$ requires the distribution of $\epsilon_{1t}$ to be known.
If instead one is willing to assume only that $\mathbb{E}[\epsilon_{1t}(\delta)] \approx \delta$, it is possible to sidestep this need by introducing a \textit{shock relaxation function}.
\begin{definition}
Assume $\mathcal{E}_1 = [a, b]$. A shock relaxation function is a map $\rho : \mathcal{E}_1 \to [0, 1]$ such that $\rho(e) = 0$ for all $e \in \mathbb{R} \,\setminus\, \mathcal{E}_1$, $\rho(e) \geq 0$ for all $e \in \mathcal{E}_1$ and there exists $e_0 \in \mathcal{E}_1$ for which $\rho(e_0) = 1$. Moreover, for a given shock $\delta \in \mathbb{R}$,
\begin{itemize}
\item[(i)] If $\delta > 0$, $\rho$ is said to be right-compatible with $\delta$ if $e + \rho(e)\delta \leq b$ for all $e \in \mathcal{E}_1$.
\item[(ii)] If $\delta < 0$, $\rho$ is said to be left-compatible with $\delta$ if $e + \rho(e)\delta \geq a$ for all $e \in \mathcal{E}_1$.
\item [(iii)] $\rho$ is compatible with shock magnitude $|\delta| > 0$ if it is both right- and left-compatible.
\end{itemize}
\end{definition}
By setting $\epsilon_{1t}(\delta) = \epsilon_{1t} + \delta \rho(\epsilon_{1t})$ for a $\rho$ compatible with $\delta$, it follows that $ X_{t}(\delta) = X_{t} + \delta \rho(\epsilon_{1t})$ and $\lvert \mathbb{E}[\epsilon_{1t}(\delta)] \rvert = \lvert \delta \mathbb{E}[\rho(\epsilon_{1t})] \rvert \leq \lvert \delta \rvert$ since $\mathbb{E}[\rho(\epsilon_{1t})] \in [0, 1]$ by definition of $\rho$. If $\rho$ is a bump function, a relaxed shock is a structural shock that has been mitigated proportionally to the density of innovations at the edges of $\mathcal{E}_1$ and the squareness of $\rho$.
\begin{remark}
When studying impulse responses, a researcher is primarily interested in shock $\delta$ itself, not in $\rho$, and the latter plays the role of a tuning parameter.
From a practical perspective, given a choice of $\delta$ (or a range $\mathcal{D}$) of interest, the researcher should explicitly select $\rho$ to \textit{minimize distortions} implied by using $\epsilon_{1t} + \delta \rho(\epsilon_{1t})$ instead of a pure shift $\epsilon_{1t} + \delta$.
If $\delta$ is sufficiently small and $\epsilon_{1t}$ is sufficiently concentrated, negligible distortions can be achieved.
Importantly, one can empirically check the impact of $\rho$ on the nonparametric IRFs by comparing relaxed and non-relaxed response estimates.
For example, in Appendix~\ref{appendix:robustness}, we provide robustness checks showing that, for both applied examples in Section~\ref{section:applications}, our chosen relaxation functions introduce negligible distortions.
\end{remark}
It is important to emphasize that shock relaxation is a generalization of standard shock designs. Indeed, when $\mathcal{X} = \mathbb{R}$ and $\mathcal{E}_1 = \mathbb{R}$, $\rho = 1$ is a relaxation function compatible with all $\delta \in \mathbb{R}$.
Nonetheless, we may also wonder how much information on the nonlinear term $G_2$ we can recover at the ``boundary'' of a finite sample.
If $X_t$ is unbounded but well-concentrated, even under strong smoothness conditions and strictly positive density, little can be learned about the \textit{local} structure of regression functions in regions of low density.\footnote{In our regression setting, for example, Theorem 1 in \cite{kohlerOptimalGlobalRates2009} on $L_2$ kernel regression error, assuming $\mathbb{E}[\lvert X_t \rvert^{\beta}] \leq M < \infty$ for some constant $\beta > 2s$, would require the bandwidth to grow over $\mathcal{X}$ faster than $\lvert X_t \rvert$. This question is also linked to issues in kernel density estimation over sets with boundary, see e.g. \cite{karunamuni2005on,malec2014nonparametric,berry2017density} and references therein.}
\begin{remark}\label{remark:shock_relax_choice}
In this paper, and more specifically in Sections~\ref{section:simulations} and \ref{section:applications}, we choose $\rho$ to be a symmetric exponential bump function, $\rho \in \{ x \mapsto \mathbb{I}\{x \leq c\} \exp(1 + (|x/c|^\alpha - 1)^{-1} ) \:|\: \alpha > 0 \}$ for some constant $c > 0$.
This $\mathcal{C}^\infty$ bump class is widely studied in both functional \citep{mitrovic1997fundamentals} and Fourier analysis \citep{stein2011fourier}.\footnote{For generic shock distributions, one can for also consider the class $\{ x \mapsto \mathbb{I}\{a \leq x \leq b\} \exp(1 + (|2(x-b)/(b-a) + 1|^\alpha - 1)^{-1} ) \:|\: \alpha > 0 \}$ of exponential bump functions with domain $[a,b] \subset \mathbb{R}$.}
We aim to set $\alpha$ to be as large as possible to minimize distortions from a linear shift, while retaining compatibility with $\delta \in \mathcal{D}$, where $\mathcal{D}$ is a set of shocks of empirical interest.
\end{remark}
\subsection{Relaxed Impulse Response Consistency}
We will now study relaxed impulse responses in the setting of additively separable models. Additive separability is a common assumption in applied work, as we shall impose it in the empirical analyses of Section~\ref{section:simulations} and \ref{section:applications}. Further, collecting nonlinear terms over lags significantly streamlines notation and analysis, and aligns with the setup of Corollary~\ref{corollary:twostep_estimator_op1}. It would be straightforward, if tedious, to extend our derivations below to the more general setting of Theorem~\ref{theorem:twostep_estimator_consistency}.
Given $\delta \in \mathbb{R}$ and compatible shock relaxation function $\rho$, let $\widetilde{\delta}_t := \delta \rho({\epsilon}_{1t})$.
Starting from a path $X_{t+j:t}$ and \eqref{eq_main:irf_h_separable}, the relaxed shock path is
\begin{equation*}
X_{t+j}(\widetilde{\delta}_t)
= X_{t+j} + \Theta_{j,11} \widetilde{\delta}_t + \sum_{k=1}^j \left[ \Gamma_{k,11} X_{t+j-k}(\widetilde{\delta}_t) - \Gamma_{k,11} X_{t+j-k} \right]
=: \gamma_{j}(X_{t+j:t}; {\widetilde{\delta}}_t) .
\end{equation*}
The relaxed-shock impulse response is thus given by
\begin{align*}
{\widetilde{\textnormal{IRF}}}_{h}(\delta)
&
:= \mathbb{E}[Z_{t+j}(\widetilde{\delta}_t) - Z_{t+j}]
= \Theta_{h,\cdot 1} \delta \, \mathbb{E}\left[ \rho(\epsilon_{1t}) \right] + \sum_{k=1}^j \mathbb{E}\left[ \Gamma_{k} X_{t+j-k}(\widetilde{\delta}_t) - \Gamma_{k} X_{t+j-k} \right] .
\end{align*}
For $1 \leq \ell \leq d$, we define ${V}_{j,\ell}(\delta)$ to be the sample analog of the horizon $j$ nonlinear effect on the $\ell$th variable,
\begin{equation*}
{V}_{j,\ell}(\delta)
:=
\frac{1}{n-j} \sum_{t=1}^{n-j} \left[ {\Gamma}_{j,\ell} {\gamma}_{j}(X_{t+j:t}; {\widetilde{\delta}}_t) - {\Gamma}_{j,\ell} X_{t+j} \right]
=
\frac{1}{n-j} \sum_{t=1}^{n-j} v_{j,\ell}(X_{t+j:t}; {\widetilde{\delta}}_t) ,
\end{equation*}
where ${\Gamma}_{j,\ell}$ is the $\ell$th component of functional vector ${\Gamma}_{j}$.
As $\epsilon_{1t}$ is not universally observable, we introduce its residual counterpart, $\widehat{\widetilde{\delta}}_t = \delta \rho(\widehat{\epsilon}_{1t})$.
The associated plug-in sample estimates are
$\widehat{V}_{j,\ell}(\delta)
=
({n-j})^{-1} \sum_{t=1}^{n-j} \widehat{v}_{j,\ell}\big( X_{t+j:t}; \widehat{\widetilde{\delta}}_t \big)$,
$\widehat{v}_{j,\ell}(X_{t+j:t}; \widehat{\widetilde{\delta}}_t)
=
\widehat{\Gamma}_{j,\ell} \widehat{\gamma}_{j}(X_{t+j:t}; \widehat{\widetilde{\delta}}_t) - \widehat{\Gamma}_{j,\ell} X_{t+j}$,
and
\begin{equation*}
\widehat{\widetilde{\textnormal{IRF}}}_{h,\ell}(\delta)
=
\widehat{\Theta}_{h,\cdot 1} \delta \left( \frac{1}{n} \sum_{t=1}^{n} \rho(\widehat{\epsilon}_{1t}) \right)
+
\sum_{j=0}^h \widehat{V}_{j,\ell}(\delta) .
\end{equation*}
Our next theorem proves the consistency of the relaxed impulse responses estimator based on semiparametric series estimates. We leverage the sup-norm bounds of Theorem~\ref{theorem:twostep_estimator_consistency} to derive a result that is uniform in $\delta$ over a compact interval $[-\mathcal{D}, \mathcal{D}]$, $\mathcal{D} > 0$. This allows us to make valid comparisons between IRFs due to shocks of different sizes.
\begin{theorem}\label{theorem:consistent_irf}
Let $\widehat{\widetilde{\textnormal{IRF}}}_{h,\ell}(\delta)$ be the semiparametric estimate for the horizon $h$ relaxed shock IRF of variable $\ell$ based on relaxation function $\rho$ with compatibility range $[-\mathcal{D}, \mathcal{D}]$. Under the assumptions in Theorem \ref{theorem:twostep_estimator_consistency} and Assumption~\ref{assumption:roots_linear_part_model},
\begin{equation*}
\sup_{\delta \in [-\mathcal{D}, \mathcal{D}]} \left\lvert \widehat{\widetilde{\textnormal{IRF}}}_{h,\ell}(\delta)
-
\widetilde{\textnormal{IRF}}_{h,\ell}(\delta)
\right\rvert
=
o_P(1)
\end{equation*}
for any fixed integers $0 \leq h < \infty$ and $1 \leq \ell \leq d$.
\end{theorem}
\begin{remark}
By construction of $\widehat{\widetilde{\textnormal{IRF}}}_{h,\ell}(\delta)$, Proposition~\ref{prop:irf_iterate_algorithm} remains valid when computing ${\widetilde{\textnormal{IRF}}}_{h}(\delta)$ instead of $\textnormal{IRF}_{h}(\delta)$. The only adjustment to be made is that in step (i) one must set $X_{t}(\delta) = X_{t} + \delta \rho(\epsilon_{1t})$ and iterate forward accordingly. Assumptions \ref{assumption:structural_model}, \ref{assumption:stationarity}, and \ref{assumption:physical_dep} ensure that the IRFs of interest are well-defined.
\end{remark}
\begin{remark}
Our definition of a compatible relaxation function is \textit{static}, as it considers only the impact effect of a shock. Nonetheless, $X_t(\delta) \in \mathcal{X}$ for all $t$ must hold to properly define ${\widetilde{\textnormal{IRF}}}_{h}(\delta)$. In theory, given $\delta$, one can always either expand $\mathcal{X}$ or strengthen $\rho$ so that compatibility is enforced at all horizons $1 \leq h \leq H$. In simulations, the choice of domains and relaxation functions can be done transparently. When working with empirical data, unless $X_t$ is exogenous or strictly autoregressive, more care has to be taken to check that there is no dynamic domain violation. In Section~\ref{subsection:app_istrefi}, where $X_t$ is an endogenous series, we discuss such a robustness check.
\end{remark}
\section{Simulations}\label{section:simulations}
To analyze the performance of the two-step semiparametric estimation strategy discussed above, we begin by considering the two simulation setups employed by \cite{goncalvesImpulseResponseAnalysis2021}. We compare the bias and MSE of the estimated relaxed shocked impulse response functions for different methods. The population responses we consider in this section are also constructed using the same shock relaxation scheme; therefore, both relaxed IRF estimators are correctly specified. Appendix~\ref{appendix:robustness} includes a robustness analysis wherein non-relaxed nonlinear population IRFs are targeted.
We also provide simulations under a misspecified design, which highlight how, in larger samples, the nonparametric sieve estimator consistently recovers impulse responses, whereas a least-squares estimator constructed with a pre-specified nonlinear transform may not.\footnote{Population impulse responses are estimated with $10^5$ replications, while MSE and bias of both semiparametric and parametric IRFs are computed with $10^4$ Monte Carlo replications. In all setups, a cubic B-spline sieve is used.}
\paragraph*{Benchmarks.}
Like in \cite{goncalvesImpulseResponseAnalysis2021}, we consider two simulation setups: A bivariate design with identified shocks (DGPs 1-3) and a three-variable design with partial block-recursive identification (DGPs 4-6).
In both, we set a sample size of $n = 240$, which is realistic for most macroeconomic data settings: this is approximately equivalent to 20 years of monthly data or 60 years of quarterly data \citep{goncalvesImpulseResponseAnalysis2021}. We discuss here only the bivariate simulation design with shock $\delta = +1$, and refer the reader to Appendix~\ref{appendix:sim_details} for the block-recursive setup.
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{figures/plot_mse_bias_DGP2__n=240_deg=3_B=100000_M=10000.pdf}
\caption{Simulation results for DGP 2 with $\delta = +1$.}
\label{fig:mse_bias_DGP_2}
\end{figure}
We set either $X_t = \epsilon_{1t}$ (DGP 1), $X_t = 0.5 X_{t-1} + \epsilon_{1t}$ (DGP 2) or $X_t = 0.5 X_{t-1} + 0.2 Y_{t-1} + \epsilon_{1t}$ (DGP 3), and
\begin{equation*}
Y_t = 0.5 Y_{t-1} + 0.5 X_{t} + 0.3 X_{t-1} - 0.4 \max(0, X_t) + 0.3 \max(0, X_{t-1}) + \epsilon_{2t} .
\end{equation*}
Innovations $\epsilon_{1t}$ and $\epsilon_{2t}$ are drawn as independent, truncated standard Gaussian variables over $[-3, 3]$. The shock relaxation function is $\rho(z) = \mathbb{I}\{|z| \leq 3 \}\exp\left( 1 + ( \lvert {z}/{3} \rvert^4 - 1 )^{-1} \right)$, cf. Remark~\ref{remark:shock_relax_choice}.
In Figure~\ref{fig:mse_bias_DGP_2} we show MSE and bias curves for the IRF on $Y_t$ in DGP 2, where $X_t$ is an exogenous AR(1) process. One can see that the sieve IRF leads only to a minor increase in mean squared error at short horizons compared to directly estimating the parameter of the true specification. This marginal increase in MSE is consistent across DGPs 1 through 3.
These simulations show that there is negligible loss of efficiency in terms of either MSE or bias when implementing the fully flexible semiparametric estimates at realistic sample sizes. We confirm these results when studying DGPs 4-6, where estimation of the structural matrix $B_0$ is included in the regression problem. Detailed results can be found in Appendix~\ref{appendix:sim_details}.
\paragraph*{Misspecified Model.}
To assess the robustness of the proposed semiparametric approach versus the parametric nonlinear model, we consider a modified process (DGP 7):
\begin{equation}\label{sim_eg:DGP_7}
\begin{split}
X_t & = 0.8 X_{t-1} + \epsilon_{1t} , \\
Y_t & = 0.5 Y_{t-1} + 0.9 \varphi(X_t) + 0.5 \varphi(X_{t-1}) + \epsilon_{2t} ,
\end{split}
\end{equation}
where $\varphi(x) := (x - 1)(0.5 + \tanh(x - 1)/2)$. In this design, we assume that the researcher's prior is $\varphi(x) = \max(0, x)$, as in the benchmark simulations. To emphasize the difference in estimated IRFs, in this setup we focus on $|\delta| = 2$ and $n = 2400$; innovations $\epsilon_{1t}$ and $\epsilon_{2t}$ are drawn from a standard Gaussian distribution truncated over $[-5, 5]$, and $\rho(z) = \mathbb{I}\{|z| \leq 5 \}\exp( 1 + ( \lvert {z}/{5} \rvert^{3.9} - 1 )^{-1} )$.
As Figure~\ref{fig:mse_bias_DGP_2_plus} shows, positive-shock parametric nonlinear IRF estimates are severely biased, while semiparametric sieve IRFs have comparatively negligible error: This yields an up to 4 times reduction of overall MSE at short horizons. Appendix~\ref{appendix:sim_details} provides additional simulation results showing that the same improvements hold when $\delta = - 2$. There, we also discuss the setting where $\varphi(x)$ is replaced with map $\widetilde{\varphi}(x) = \varphi(x+1)$, which agrees closely with $\max(0, x)$. In this last setting, we find that parametric nonlinear regression dominates in MSE and bias terms. As one might expect, therefore, parametric modeling is reliable only in cases where a sufficiently good model prior is available.
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{figures/plot_mse_bias_DGP2-smooth__n=2400_deg=3_B=100000_M=10000_delta=2.pdf}
\caption{Simulation results for DGP 7 with shock $\delta = +2$.}
\label{fig:mse_bias_DGP_2_plus}
\end{figure}
\section{Empirical Applications}\label{section:applications}
In this section, we showcase the practical utility of the proposed semiparametric sieve estimator by considering two applied exercises.
In line with previous work applying nonlinear IRF methods to macroeconomic data, such as e.g. \cite{kilian2011responses,goncalvesImpulseResponseAnalysis2021,goncalvesStatedependentLocalProjections2024} and \cite{goncalvesNonparametricLocalProjections2024}, our discussion is focused on point impulse response estimates.
In both cases, we will consider sieve IRFs constructed with shock relaxation: Appendix~\ref{appendix:robustness} shows that our analysis remains valid also when evaluating non-relaxed semiparametric responses.
\subsection{Monetary Policy Shocks}\label{subsection:app_goncalves}
We first consider a four-variable model identical to the one analyzed by \cite{goncalvesImpulseResponseAnalysis2021}, and based on \cite{tenreyro2016pushing}. Let $Z_t = (X_t, \textnormal{FFR}_t, \textnormal{GDP}_t, \textnormal{PCE}_t)'$, where $X_t$ is the series of narrative U.S. monetary policy shocks, $\textnormal{FFR}_t$ is the federal funds rate, $\textnormal{GDP}_t$ is log-real GDP and $\textnormal{PCE}_t$ is PCE inflation.\footnote{In \cite{goncalvesImpulseResponseAnalysis2021} p.~122, it is mentioned that CPI inflation is included in the model, but both in the replication package made available by one the authors (\url{https://sites.google.com/site/lkilian2019/research/code}) from which we source the data, and in \cite{tenreyro2016pushing}, PCE inflation is used instead. Moreover, the authors say that both the FFR and PCE enter the model in first differences, yet, in their code, these variables are kept in levels. We thus consider a model in levels to allow for a proper comparison between estimation methods, although the series are highly persistent.}
As a pre-processing step, GDP is transformed to log GDP and then linearly detrended. The data is available quarterly and spans from 1969:Q1 to 2007:Q4.
As in \cite{tenreyro2016pushing}, we use a model with one lag, $p=1$. The narrative shock $X_t$ is considered to be an i.i.d. sequence, i.e. $X_t = \epsilon_{1t}$, therefore we assume no dependence on lagged variables when implementing the pseudo-reduced form \eqref{eq_main:pseudo_reduced_form_structural_model}. Like in \cite{goncalvesImpulseResponseAnalysis2021}, we consider positive and negative shocks of size $|\delta| = 1$ and choose
$\rho(z) = \mathbb{I}\{ |z| \leq 4 \} \exp( 1 + ( \lvert {z}/{4} \rvert^{6} - 1 )^{-1} )$
to be the shock relaxation function. Figure \ref{fig:plot_app_gonc2021_meanshift} in the Online Appendix provides a check for the compatibility of $\rho$ given the sample distribution of $X_t$. Knots for sieve estimation are located at $\{-1, 0, 1\}$. The model is block-recursive, and U.S. monetary policy shocks are identified without the need to impose additional assumptions on the remaining shocks.
\cite{goncalvesImpulseResponseAnalysis2021}, like \cite{tenreyro2016pushing}, use two nonlinear transformations, $F(x) = \max(0, x)$ and $F(x) = x^3$, to try to gauge how negative versus positive and large versus small shocks, respectively, affect the U.S. macroeconomy. They find that the two maps yield very similar responses, so we focus on comparing the IRFs estimated via sieve regression with the ones obtained by setting $F(x) = \max(0, x)$, as well as linear IRFs.
\begin{figure}[t!]
\centering
\begin{subfigure}[b]{\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/plot_app_cub_gonc2021_irfs__delta=1.pdf}
\caption{$\delta = +1$}
\end{subfigure}
\\[15pt]
\begin{subfigure}[b]{\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/plot_app_cub_gonc2021_irfs__delta=-1.pdf}
\caption{$\delta = -1$}
\end{subfigure}
\\[5pt]
\caption{Effect of an unexpected U.S. monetary policy shock on federal funds rate, GDP and inflation. Linear (gray, dashed), parametric nonlinear with $F(x) = \max(0, x)$ (red, point-dashed) and sieve (blue, solid) structural impulse responses. For $\delta = +1$, the lowest point of the GDP response is marked with a dot. Note that $\delta = 1 \approx 1.7 \times \sigma_{\epsilon,1}$.}
\label{fig:app_gonc2021_irfs}
\end{figure}
Figure \ref{fig:app_gonc2021_irfs} plots estimated impulse responses to both positive and negative monetary policy shocks.
The impact on the federal funds rate is consistent across all three procedures. The semiparametric nonlinear response for GDP, unlike in the case of linear and parametric nonlinear IRFs, is nearly zero at impact and has a monotonic decrease until around 10 quarters ahead. The change in shape is meaningful, as the procedure of \cite{goncalvesImpulseResponseAnalysis2021} still yields a small short-term upward jump in GDP when a monetary tightening shock hits. Moreover, after the positive shock, the sieve GDP responses reaches its lowest value 4 and 2 quarters before the linear and parametric nonlinear responses, while its size is 13\% and 16\% larger, respectively.\footnote{The strength of this effect changes across different shock sizes, as Figure \ref{fig:plot_app_gonc2021_scale} in Appendix \ref{appendix:additional_plots} proves. As shock sizes get smaller, nonlinear IRFs, both parametric and sieve, show decreasing negative effects.} Finally, the sieve PCE response is positive for a shorter interval, but looks to be more persistent once it turns negative, also 10 months after impact.
When the shock is expansionary, one sees that the semiparametric FFR response is marginally mitigated compared to the alternative estimates. An important puzzle is due to the negative impact on GDP: Both types of nonlinear responses show a drop in output in the first 5 quarters. Such a quick change seems unrealistic, as one does not expect inflation to suddenly reverse sign, but, as \cite{goncalvesImpulseResponseAnalysis2021} also remark, the overall impact on inflation of both shocks is small when compared to the change in federal funds rate.
\subsection{Uncertainty Shocks}\label{subsection:app_istrefi}
Traditional central bank policymaking is heavily guided by the principle that a central bank can and should influence expectations. Therefore, controlling the (perceived) level of ambiguity in current and future commitments is key. \cite{istrefiSubjectiveInterestRate2018} provide an analysis of the impact of unforeseen changes in the level of subjective interest rate uncertainty on the macroeconomy.
For the sake of simplicity, our evaluation will focus only on their 3-month-ahead uncertainty measure for short-term interest rate maturities (3M3M) and the U.S. economy.
Like in \cite{istrefiSubjectiveInterestRate2018}, let $Z_t = (X_t, \textnormal{IP}_t, \textnormal{CPI}_t, \textnormal{PPI}_t, \textnormal{RT}_t, \textnormal{UR}_t)'$ be a vector where $X_t$ is the chosen uncertainty measure, $\textnormal{IP}_t$ is the (log) industrial production index, $\textnormal{CPI}_t$ is the CPI inflation rate, $\textnormal{PPI}_t$ is the producer price inflation rate, $\textnormal{RT}_t$ is (log) retail sales and $\textnormal{UR}_t$ is the unemployment rate. The nonlinear model specification is given by
$Z_t = \mu + A_1 Z_{t-1} + A_2 Z_{t-1} + F_1(X_{t-1}) + F_2(X_{t-2}) + D W_t + u_t ,$
where $W_t$ includes a linear time trend and oil price $\textnormal{OIL}_t$.\footnote{Inclusion of linear exogenous variables in the semiparametric theoretical framework in Section~\ref{section:estimation} is straightforward as long as one can assume that they are stationary and weakly dependent. The choice of using $p=2$ is identical to that of the original authors, based on BIC.}
The data has a monthly frequency and spans the period between May 1993 and July 2015.\footnote{We utilize the original data employed by the authors, who kindly shared it upon request. However, we rescale retail sales ($\textnormal{RT}_t$) so that the level in January 2000 equals 100.}
Note here that nonlinear functions $F_1$ and $F_2$ are assumed not to affect $X_t$, which is the structural variable. The linear VAR specification of \cite{istrefiSubjectiveInterestRate2018} is recovered by simply assuming $F_1 = F_2 = 0$ prior to estimation. Since they use recursive identification and order the uncertainty measure first, this model too is block-recursive.
We consider a positive shock with intensity $\delta = \sigma_{\epsilon,1}$, where $\sigma_{\epsilon,1}$ is the standard deviation of structural innovations. In this empirical exercise, the relaxation function is
$\rho(z) = \mathbb{I}\left\{ |z| \leq 1/4 \right\} \exp( 1 + ( \lvert 4 x \rvert^{8} - 1 )^{-1} ) $
and we set $\{0.1, 0.3\}$ to be the cubic spline knots.
As 3M3M is a non-negative measure of uncertainty, some care must be taken to make sure that the shocked paths for $X_t$ do not reach negative values. Figure \ref{fig:plot_app_istrefi2018_meanshift} in Appendix \ref{appendix:additional_plots} shows that the relaxation function is compatible, and also that the shocked nonlinear paths of $X_t$ with impulse $\delta$ and $\delta'$ all do not cross below zero.
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{figures/plot_app_cub_istref2015_alt_irfs__delta=1.pdf}
\\[5pt]
\caption{Effect of an unexpected, one-standard-deviation uncertainty shock to U.S. macroeconomic variables. Linear (gray, dashed) and sieve (blue, solid) structural impulse responses. The extreme points of the responses are marked with a dot.}
\label{fig:app_istrefi2018_irfs}
\end{figure}
Figure \ref{fig:app_istrefi2018_irfs} presents both the linear and nonlinear structural impulse responses obtained. Importantly, even though \cite{istrefiSubjectiveInterestRate2018} estimate a Bayesian VAR model and here we consider a frequentist vector autoregressive benchmark, the shape of the IRFs is retained, cf. the median response in the top row of their Figure 4. When uncertainty increases, industrial production drops, and the size and extent of this decrease are intensified in the nonlinear responses. The sieve IP response reaches a value that is $54\%$ lower than that of the respective linear IRF.\footnote{Figure \ref{fig:plot_app_istrefi2018_scale_IP} in Appendix \ref{appendix:additional_plots} confirms that this difference is consistent over a range of shock sizes, too.} A similar behavior holds true for retail sales ($38\%$ lower) and unemployment ($23\%$ higher), proving that this shock is more profoundly contractionary than suggested by the linear VAR model.
Further, CPI and PP inflation both display short-term fluctuations, which strengthen the short- and medium-term impact of the shock. CPI and PP nonlinear inflation responses are $76\%$ and $41\%$ stronger than their linear counterpart, respectively. These differences show that linear IRFs might be both under-estimating the short-term intensity and misrepresenting the long-term persistence of inflation reactions.
Given the strength of nonlinear IRFs, this discrepancy may also suggest that the 3M3M uncertainty measure partially captures the financial channel, too.
Hence, we believe our analysis provides some evidence that the linear VAR used by \cite{istrefiSubjectiveInterestRate2018} may miss some key impulse response features.\footnote{See also Figure \ref{fig:plot_app_istrefi2018_regfuns}, which plots the estimated functions of the endogenous variables.}
\section{Conclusion}\label{section:conclusion}
This paper studies the application of semiparametric series estimation to the problem of structural impulse response analysis for time series. After first discussing partial block-recursive identification, we have shown that, for models with moderate physical dependence, series estimation can be employed and structural IRFs are consistently estimated. Simulations showcase that this approach is both valid in moderate samples and has the added benefit of being robust to misspecification of the nonlinear model components. Finally, two empirical applications showcase the potential insights gained by departing from either linear or parametric nonlinear specifications when estimating structural responses.
A key aspect that we have not touched upon is inference in the form of confidence intervals. The development of inferential theory appears feasible in light of the uniform inference results obtained by e.g. \cite{belloniNewAsymptoticTheory2015} in the i.i.d. setting and \cite{liUniformNonparametricInference2020} for time series data, and it is an important direction for future research.
Studying other sieve spaces, such as neural networks \citep{chenImprovedRatesAsymptotic1999,farrellDeepNeuralNetworks2021a} or shape-preserving sieves \citep{chenChapter76Large2007}, would also be highly desirable. Finally, in the spirit of \cite{kangINFERENCENONPARAMETRICSERIES2021}, deriving new inference results that are uniform in the selection of series terms is important, as, in practice, the sieve should be tuned in a data-driven way.
\pagebreak
\bibliography{nonlin_irf_bib}
\newpage