Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,680 characters · 12 sections · 69 citation commands
Semiparametric Local Projections
\doublespacing
Impulse response analysis is a cornerstone of empirical macroeconomics. Local projections have become a popular method for estimating impulse response functions (IRFs). In their simplest form, local projections consist of a sequence of OLS regressions, one for each horizon of interest. The impulse response of interest may be recovered from the estimated regressions without further transformations of the model coefficients or the need for Monte Carlo integration methods.
A large empirical literature has used generalizations of linear local projections to evaluate state-dependent impulse responses and other nonlinear responses. However, as shown in goncalves2021,goncalves2024b, goncalves2024a, such state-dependent local projections fail to recover the population responses when the state is endogenous and the shock is large as in much of applied work (e.g., ramey2018).
In this paper, we propose a semiparametric local projection estimator of nonlinear impulse response functions that is valid across a broad range of nonlinear settings relevant for applied macroeconomists. These include processes with nonlinearly transformed regressors herrera2015, tenreyro2016, benzeev2023, caravello2024, state-dependent coefficients ramey2018, and nonlinear interactions between shocks and state variables caramp2026monetary, cloyne2020decomposing, cloyne2023statedependent.
As is common in macroeconomics, our object of interest is the impulse response function of an outcome variable $y_{t+h}$ with respect to the primitive structural shock $\varepsilon_{1t}$ in the equation for variable $x_t$. Specifically, we aim to identify and estimate the response of $y_{t+h}$ to a shock of size $\delta$ in $\varepsilon_{1t}$. We assume that $x_t$ is predetermined with respect to $y_t$, an exclusion restriction that encompasses situations in which $x_t = \varepsilon_{1t}$ is an observed i.i.d.\ shock, as in the narrative approach to identification. In this case, $x_t$ is unconditionally independent of all other shocks driving the system between $t$ and $t+h$. When $x_t$ is not an observed shock $\varepsilon_{1t}$, the exclusion restriction together with the i.i.d.\ assumption on the structural shocks implies that $x_t$ is conditionally independent of all other shocks between $t$ and $t+h$, given control variables $\mathbf{z}_{t-1}$ that include the history of the system up to $t-1$.
We show how this conditional independence condition, combined with the assumption that $\varepsilon_{1t}$ enters $x_t$ additively, can be used to identify the IRF of $y_{t+h}$ with respect to $\varepsilon_{1t}$ using only the observables $(y_{t+h},x_t,\mathbf{z}_{t-1})$. In particular, the additive structure ensures that a $\delta$-perturbation in $\varepsilon_{1t}$ translates into a $\delta$-perturbation in $x_t$, holding fixed the control variables. The structural equation for the outcome variable is left unrestricted, permitting arbitrary nonlinearities.
This identification strategy places our problem within the semiparametric literature on inference for linear functionals of regression functions (newey1994, chernozhukov2018, chernozhukov2022). The key estimation challenge is that the conditional mean function $g_{0,h}(x,z) = E(y_{t+h}|x_t=x, \mathbf{z}_{t-1}=z)$ must be estimated nonparametrically, which can induce bias in the plug-in estimator. We address this concern by using the doubly robust moment condition of chernozhukov2022, augmented by a density ratio reflecting the relative change in the conditional distribution of $x_t$ given $\mathbf{z}_{t-1}$ when $x_t$ is shifted by $\delta$.
To handle the serial dependence of the data, we combine this moment condition with the NLO (“neighbors-left-out”) cross-fitting approach of semenova2023, which ensures approximate independence between training and evaluation sets. We derive the asymptotic distribution of the resulting estimator and show that it is $\sqrt{T}$-consistent and asymptotically normal, with the preliminary estimation of the nuisance functions having no effect on the first-order asymptotic distribution.
Our paper is related to a recent and growing literature on semiparametric and nonparametric inference for IRFs. The problem of estimating the average effects of policy interventions nonparametrically dates back at least to stock1989, but our focus is on the causal effect of structural macroeconomic shocks. Several recent papers have proposed nonparametric methods for estimating nonlinear IRFs. For example, gourieroux2023 propose a nonparametric local projection estimator for nonstructural IRFs identified from Gaussian shocks within a Markov process framework. ballarin2024 proposes a sieve-based nonparametric estimator for IRFs in models with nonlinearly transformed regressors, as in our Example (ref) below, but does not cover the doubly robust approach, the more general nonlinear settings we consider, or inference. While our paper deals with shocks of finite magnitude $\delta$, kolesar2025 focus on infinitesimally small shocks. They discuss the causal content of linear local projections when the data generating process is nonlinear and highlight challenges in applying doubly robust methods in the small samples typical of macroeconomics. Our paper builds on this work as well as two recent studies that apply doubly robust methods in macroeconometrics. ballinari2025 develop semiparametric inference for IRFs using double/debiased machine learning in a time series context, but focus on a binary treatment, so the adjustment term in their orthogonal moment condition is based on the propensity score rather than the density ratio we employ for continuous treatments. huang2026 develop a two-step high-dimensional nonparametric local projection estimator combining Neyman-orthogonal pseudo-outcomes with cross-fitting. Because they focus on an IRF that shifts the policy variable from a fixed baseline, their second step requires a nonparametric regression and yields convergence rates slower than $\sqrt{T}$; in contrast, our estimand averages over the distribution of the shock and can be estimated at rate $\sqrt{T}$. Finally, and independently, Nikolaishvili2026, building on an earlier version of our paper goncalves2024b, proposes a closely related doubly robust estimator for nonparametric local projections under different assumptions and without allowing for covariates in the conditioning set.
While our main focus is on identifying unconditional responses, our analysis also extends to average response functions conditional on a state variable $\Omega_t$. We provide identification conditions for this object and show how our semiparametric local projections estimator can be applied to each subsample $\{t: \Omega_t = \omega\}$ when $\Omega_t$ is discrete, as in state-dependent models.
The paper is organized as follows. Section (ref) introduces the structural model and three leading examples. Section (ref) defines the population IRFs of interest and contrasts our definition with alternative definitions used in the literature. Sections (ref) and (ref) discuss identification and estimation, and inference, respectively. Section (ref) briefly discusses how to extend the analysis to conditional IRFs. The simulation results are presented in Section (ref). Section (ref) contains two empirical illustrations focusing on possible nonlinearities in the pass-through of gasoline price shocks to inflation and in the response of motor vehicle sales to real gasoline price shocks. We conclude in Section (ref). The proofs are relegated to Appendices (ref) and (ref).
Let $z_t=(x_t, y_t)'$ denote a vector of observed time series, where $y_t$ is the outcome of interest and $x_t$ is predetermined with respect to $y_t$. For example, $y_t$ could be real GDP and $x_t$ government spending. For simplicity, we assume that $y_t$ is univariate, but extensions to multivariate outcomes could be easily accommodated. A general structural model for $z_t$ is a triangular system of the form
where $\phi$ and $\mu$ are (potentially unknown) nonlinear functions, and $\varepsilon_t = (\varepsilon_{1t}, \varepsilon_{2t})'$ is a vector of mutually independent structural shocks. We assume $\varepsilon_t$ to be i.i.d.\ with mean zero and diagonal covariance matrix $\Sigma = \mathrm{diag}(\sigma_1^2, \sigma_2^2)$. The vector $\mathbf{z}_{t-1}= (z_{t-1}, z_{t-2}, \ldots, z_{t-p})'$ contains lags of $z_t$; other variables can be included in $\mathbf{z}_{t-1}$ provided they are predetermined with respect to $\varepsilon_{1t}$. The exclusion of $y_t$ from equation ((ref)) is equivalent to the assumption of block recursiveness in linear structural VAR identification, where $x_t=\phi'\mathbf{z}_{t-1}+\varepsilon_{1t}$. An important special case is $x_t=\varepsilon_{1t}$, as in the narrative approach to identification.
Consistent with standard impulse response analysis in macroeconomics, our goal is to trace the effect over time of surprise changes in the variable $x_t$, as captured by a one-time change in the structural shock $\varepsilon_{1t}$, rather than changes in $x_t$ itself. This is the main reason why we assume in ((ref)) that $x_t$ and $\varepsilon_{1t}$ are separable, as this is crucial for identifying the impulse response function of $y_{t+h}$ with respect to $\varepsilon_{1t}$ from the observables $(y_{t+h},x_t,\mathbf{z}_{t-1})$. A non-separable specification $x_t=\phi(\mathbf{z}_{t-1},\varepsilon_{1t})$ could be considered if instead we targeted the impulse response function of $y_{t+h}$ with respect to $x_t$, as in kolesar2025 and huang2026.
Our framework accommodates a range of models of interest in applied work. One example is a model with nonlinearly transformed regressors. This model allows for a sign nonlinearity in the responses with the magnitude of the response depending on the sign of $x_t$ (e.g., $f(x_t)=\max \{x_t,0\}$ ) or a size nonlinearity with the magnitude of the response depending on the size of $x_t$ (e.g., $f(x_t)=x^3_t$). Although the regression model is linear in the parameters, the impulse response function is nonlinear, requiring the use of nonstandard estimation methods (e.g., kilian2011, goncalves2021). Models with nonlinearly transformed regressors have been used extensively in applied macroeconomics. Examples include studies of the asymmetry in the responses to positive and negative oil price shocks (e.g., herrera2015) as well as nonlinearities in the response of GDP to monetary policy shocks (e.g., tenreyro2016, ascari2022), financial shocks (e.g., forni2024) and fiscal shocks (e.g., benzeev2023).
A second example is the state-dependent model examined in goncalves2024a in which the response is allowed to differ between two observed states (e.g., expansion and recession) based on a dummy variable indicator $ S_{t-1}$. Models of this type have been used extensively to study the magnitude of the fiscal multiplier, the effectiveness of monetary policy, and the impact of uncertainty shocks in expansions and recessions (e.g., ramey2018, cacciatore2021, falck2021).
A final example is inspired by cloyne2020decomposing, cloyne2023statedependent and caramp2026monetary who consider a model in which the responses of $y_{t+h}$ to $\varepsilon_{1t}$ are allowed to be heterogeneous, with the heterogeneity being captured by an observable variable, say, $r_t$. For instance, imagine a situation in which monetary policy shocks, $\varepsilon_{1t}$, have a heterogeneous effect on GDP growth, $y_t$, that depends on the level of government debt, $r_{t}$. The level of debt, in turn, is a function of monetary policy in the previous period ($x_{t-1}$) through its effect on interest rates. The interaction between the debt level and the shock of interest induces a nonlinearity that needs to be taken into account when estimating the IRF. Note that this specification differs from the state-dependent model discussed earlier in that the model coefficients do not depend on the state, but the impulse response does.
Our main analysis focuses on unconditional versions of the IRF. Although the estimation and inference methods presented in Section (ref) are tailored to estimating unconditional IRFs, they can also be applied to conditional IRFs in state-dependent models such as in Example (ref), where the conditioning set is discrete. Section (ref) discusses this application as well as the challenges one would face in estimating conditional IRFs in other examples such as Example (ref).
Following the standard approach in macroeconomics, we care about the response of $y_{t+h}$ with respect to the structural shock $\varepsilon_{1t}$. As in the recent macroeconometrics literature, we adopt a potential outcomes framework (see e.g., goncalves2021 and goncalves2024a). One implication of the structural model (ref)--(ref) is that $y_{t+h}$ can be written as $y_{t+h} = m_h(\varepsilon_{1t},\, U_{t+h})$, where $m_h$ is obtained by iterating (ref) forward $h$ steps and substituting (ref), and $U_{t+h} \equiv (\varepsilon_{2t},\, \varepsilon_{1,t+1},\, \varepsilon_{2,t+1},\, \ldots,\, \varepsilon_{1,t+h},\, \varepsilon_{2,t+h},\, \mathbf{z}'_{t-1})'$ collects all remaining determinants of $y_{t+h}$. Since $\varepsilon_{1t}$ is i.i.d.\ and independent of $\{\varepsilon_{2t}\}$, $\varepsilon_{1t}$ is independent of $U_{t+h}$, which we write as $\varepsilon_{1t} \perp U_{t+h}$. The potential outcome associated with fixing $\varepsilon_{1t} = e$ is $y_{t+h}(e) = m_h(e, U_{t+h})$, where $e$ is any fixed value in the support of $\varepsilon_{1t}$. The observed outcome satisfies $y_{t+h} = y_{t+h}(\varepsilon_{1t})$, implying that it is the value that we observe when $e$ takes the value $\varepsilon _{1t}$ that generated the observed data. The fact that $\varepsilon _{1t}$ and $U_{t+h}$ are mutually independent implies that the potential outcomes are independent of $\varepsilon _{1t}$.
To define the response function of $y_{t+h}$ with respect to $\varepsilon_{1t}$, we compare the (observed) baseline value $y_{t+h}(\varepsilon_{1t})$ with the counterfactual (unobserved) value of $y$ at $t+h$ that would have been observed if $ \varepsilon_{1t}$ had been subject to a shock of size $\delta$, denoted $ y_{t+h}(\varepsilon_{1t}+\delta)$ (e.g., potter2000). In particular, following goncalves2024a, we adopt the following definition:
$\mathrm{ARF}_h(\delta)$ corresponds to the unconditional average response used in goncalves2021. A conditional average response function can also be defined as in Definition (ref) in Section (ref), following goncalves2024b.\footnote{Alternatively, one could also define versions of these IRFs where $\delta\to 0$, as discussed in goncalves2024b. In this paper, we focus on responses to a shock of finite magnitude $\delta$, as is common in applied work (e.g., ramey2018).}
Definition (ref) is not the only possible definition of an unconditional IRF. Other studies such as koop1996, rambachan2021common, and huang2026, for example, have instead compared the two potential outcomes $y_{t+h}(e^\prime)$ and $y_{t+h}\left(e\right)$, often setting $ e^{\prime}=\delta$ and $e=0$, which yields the alternative definition:
Whereas Definition (ref) has been used widely in the literature, Definition (ref) is more recent (goncalves2021, goncalves2024b, goncalves2024a). These two definitions are equivalent when the potential outcome is linear in $e$ for all horizons, as would be the case for a linear model or in special cases of Example (ref) and Example (ref) when the conditioning sets ($S_{t-1}$ and $r_t$, respectively) are exogenous. However, in general, the two definitions differ.
Figure (ref) illustrates these differences by example. Consider the nonlinear DGP:
where $\varepsilon_{1t}$ and $\varepsilon_{2t}$ are independent and have a standard normal distribution. For illustrative purposes, let the magnitude of the shock be $\delta=2$ and the functional forms $f(x_t)=\max(x_t,0)$ and $f(x_t)=x^3_t$, respectively. The solid red line in Figure (ref) denotes the $\mathrm{ARF}$ obtained as the average over the $y_{t}$ obtained for different realizations of $ \varepsilon_{1t}$, whereas the dashed line denotes the value of $ \mathrm{ARF}^{*}$ obtained by setting $\varepsilon_{1t}=e=0$. It is readily apparent that in this example, the two definitions of the IRF imply quite different measures of the conditional expectation of $y_{t}$ in the absence of a perturbation.
Which approach is the more natural one? The only difference between these two approaches is the treatment of the impact period. The baseline in computing any impulse response is the conditional expectation of $y_{t}$ in the absence of a perturbation $\delta$ (e.g., potter2000, p. 1430). In other words, the baseline is what we would have expected $y_{t}$ to be in the absence of a perturbation, possibly conditional on the history of the data. For example, if $\mathcal{F}^{t-1}$ denotes the information available up to time $t-1$, a natural baseline is the conditional expectation $E_{t-1}\left(y_t\right) \equiv E\left(y_t \mid \mathcal{F}^{t-1}\right)$. In this example, $ E_{t-1}(y_{t}) = 0.5y_{t-1}+0.3x_{t-1}-0.4E_{t-1}(f(x_t))-0.3f(x_{t-1}),$ where the predetermined values are known and we imposed $E_{t-1}(\varepsilon _{1t})=E_{t-1}(\varepsilon _{2t})=0.$ This expectation can only be evaluated by integrating $f(x_t)$ over all possible realizations of $x_t$, as in Definition 1. In contrast, Definition 2 evaluates this expression as $f(E(x_t))=f(0)$. By Jensen's inequality, this will not yield the desired baseline for computing the population impulse response to a shock of magnitude $\delta$ because $ E(f(x_t))$ is not $f(E(x_t))$. Thus, we work with Definition 1 throughout this paper.
We discuss the identification of $\theta_{0,h}\equiv \mathrm{ARF}_h(\delta)$ in Definition (ref). Motivated by Section 6 of kolesar2025, we first discuss identification based on a regression-based approach, where $\theta_{0,h}$ is identified using the conditional expectation function $g_{0,h}(x,z)\equiv E(y_{t+h}|x_t=x, \mathbf{z}_{t-1}=z)$, and then show how to obtain identification using a doubly robust approach. We use the subscript “0” to indicate true parameters and functions throughout.
Starting with the regression-based approach, note that the independence between $\varepsilon_{1t}$ and $U_{t+h}$ (which holds by (ref) and (ref) under the i.i.d.\ assumption on $\varepsilon_t$ and the mutually independent shocks assumption) allows us to identify $\theta_{0,h}$ as $\theta_{0,h} = E\!\left[g_{\varepsilon,0,h}(\varepsilon_{1t}+\delta) - g_{\varepsilon,0,h}(\varepsilon_{1t})\right]$, where $g_{\varepsilon,0,h}(e)\equiv E(y_{t+h}|\varepsilon_{1t}=e)$ (e.g., kolesar2025). This representation is useful for estimation when $\varepsilon_{1t}$ is an observed shock, as in the narrative approach to identification. However, it does not directly apply when $x_t$ is an observed variable and $\varepsilon_{1t}$ is its underlying (unobserved) structural shock. In what follows, we show how the additive structure of (ref) can be exploited to identify the average response function of $y_{t+h}$ to an impulse in $\varepsilon_{1t}$ (i.e., $\mathrm{ARF}_h(\delta)$ in Definition (ref)) in this more general context: since $x_t=\phi(\mathbf{z}_{t-1})+\varepsilon_{1t}$, a $\delta$-shift in $\varepsilon_{1t}$ holding $\mathbf{z}_{t-1}$ fixed is identical to a $\delta$-shift in $x_t$ holding $\mathbf{z}_{t-1}$ fixed, yielding an identification result for $\theta_{0,h}$ in terms of the observables $(y_{t+h},x_t,\mathbf{z}_{t-1})$ only.
The representation of $\mathrm{ARF}_h(\delta)$ given in (ref) is identified from observables because $g_{0,h}(x,z)$ is the conditional mean of $y_{t+h}$ given $(x_t,\mathbf{z}_{t-1})$, a functional of the joint distribution of $(y_{t+h}, x_t, \mathbf{z}_{t-1})$, and the outer expectation is taken over the marginal distribution of $(x_t, \mathbf{z}_{t-1})$, both of which are observable. No knowledge of $\phi$ or $\mu$ is required. Identification requires that $g_{0,h}$ be defined on the support of $(x_t+\delta,\mathbf{z}_{t-1})$ as well as on the support of $(x_t,\mathbf{z}_{t-1})$. When the support of $\varepsilon_{1t}$ (and hence of $x_t$ given $\mathbf{z}_{t-1}$) is bounded, the shifted support may fall outside the original one, in which case $g_{0,h}(x_t+\delta,\mathbf{z}_{t-1})$ is not identified for some $(x_t,\mathbf{z}_{t-1})$ pairs. For this reason, we do not impose bounded support conditions.
Proposition (ref) leads to a moment condition of the form
which identifies $\theta_{0,h}$ when $g_h=g_{0,h}$. A natural approach is to replace $g_{0,h}$ by a first-step estimator $\hat{g}_h$, yielding $\widehat{\mathrm{ARF}}_h(\delta)=T^{-1}\sum_{t=1}^{T}\left[\hat{g}_h(x_t+\delta, \mathbf{z}_{t-1})-\hat{g}_h(x_t,\mathbf{z}_{t-1})\right]$. When $\hat{g}_h$ is estimated by machine learning or nonparametric methods, this regression-based estimator can suffer from first-order bias, and inference requires adjusting for estimation uncertainty in $g_{0,h}$. This motivates using the double/debiased machine learning approach of chernozhukov2018, chernozhukov2022, chernozhukov2024covariate to augment the moment condition (ref) using Neyman orthogonality. When combined with a form of cross-fitting that handles time series dependence, inference based on this estimator can proceed as if the nuisance functions were fully observed.
To describe the doubly robust approach, let $f_{0,x|z}(x|z)$ denote the conditional density of $x_t$ given $\mathbf{z}_{t-1}$. It can be easily shown that for any function $g_h$,
where $\alpha_0(x,z)\equiv (f_{0,x|z}(x-\delta|z)-f_{0,x|z}(x|z))/ f_{0,x|z}(x|z)$ is the Riesz representer (we omit the dependence on $\delta$ throughout). This Riesz representer is a density ratio that captures the relative change in the conditional density of $x_t$ given $\mathbf{z}_{t-1}$ when $x_t$ is shifted by $\delta$ (see Section 6 of kolesar2025, who report the form of the Riesz representer when $\delta\to 0$).
As it turns out, Equation (ref) yields a doubly robust moment equation for $\theta_{0,h}$. Defining $$\psi(y_{t+h},x_t,\mathbf{z}_{t-1},g_h,\alpha,\theta_{h})=g_{h}(x_t+\delta,\mathbf{z}_{t-1})-g_{h}(x_t,\mathbf{z}_{t-1}) -\theta_{h}+\alpha(x_t,\mathbf{z}_{t-1})(y_{t+h} -g_{h}(x_t,\mathbf{z}_{t-1})),$$ we can show that $E[\psi(y_{t+h},x_t,\mathbf{z}_{t-1},g_h,\alpha,\theta_{0,h})]=0$ holds for any $g_h$ when $\alpha=\alpha_0$ and for any $\alpha$ when $g_h=g_{0,h}$. This follows by an application of Theorem 5 of chernozhukov2022. Hence, this moment condition yields an estimator that is insensitive to misspecification of $g_h$ provided $\alpha=\alpha_0$, and insensitive to misspecification of $\alpha$ provided $g_h=g_{0,h}$.
We estimate $\theta_{0, h}$ by combining the doubly robust moment condition in (ref) with cross-fitting. In the standard cross-fitting approach, the data are partitioned into non-overlapping blocks. For each block, the nuisance functions are estimated on the complement of that block and evaluated on the heldout observations. These out-of-sample nuisance estimates are then used to construct the moment condition for estimating $\theta_{0, h}$; see, for example, chernozhukov2022. A crucial assumption that justifies this approach is random sampling, which implies that the blocks are mutually independent (and the data within blocks are i.i.d.).
In a time series context, the presence of serial dependence in $z_{t}=(x_t,y_t)'$ violates this assumption even when $x_t$ is i.i.d., creating dependence among the blocks. Hence, we follow semenova2023 and rely on NLO (“neighbors-left-out”) cross-fitting.\footnote{This approach leaves out not only the target block but also its immediate neighbors when estimating the nuisance functions. The resulting training and evaluation sets are then approximately independent, with the approximation error controlled by the speed of mixing of the underlying time series. Recent applications of the NLO cross-fitting approach to inference on impulse response functions in time series include ballinari2025 and huang2026. Because their estimands are different than ours, their assumptions and asymptotic results also differ from ours.} We partition the index set $\{1,\ldots,T\}$ into $K\ge 4$ non-overlapping blocks $I_1,\ldots,I_K$ of contiguous time indices, each of size $T_\ell\equiv T/K$: $\{1,\ldots,T\} = I_1\cup\cdots\cup I_K$.\footnote{For simplicity we assume that $T$ is divisible by $K$.} We keep $K$ fixed as $T\to \infty$ when deriving the asymptotic theory below. For each $\ell\in\{1,\ldots,K\}$, let $\mathcal{N}(\ell)$ denote the set containing $\ell$ and its immediate neighbors in $\{1,\ldots,K\}$, i.e.\ $\mathcal{N}(\ell)=\{\ell-1,\ell,\ell+1\}\cap\{1,\ldots,K\}$, and define the quasi-complement of $I_\ell$ as $I^{\mathrm{qc}}_\ell = \bigcup_{j\notin\mathcal{N}(\ell)} I_j,$ so that $I^{\mathrm{qc}}_\ell$ is obtained from the full complement $I_{-\ell}=\bigcup_{j\ne\ell}I_j$ by additionally removing the two blocks adjacent to $I_\ell$. The key feature of NLO cross-fitting is that $I_\ell$ and $I^{\mathrm{qc}}_\ell$ are separated by at least $T_\ell$ time periods. Since $K$ is fixed, as $T\to\infty$ we have $T_\ell\to\infty$, so that under mixing-type conditions $I_\ell$ is approximately independent of $I^{\mathrm{qc}}_\ell$.
Given horizon $h$, the NLO cross-fitting estimator of $\theta_{0,h}$, which we refer to as DR-NLO, is computed as follows. For notational simplicity, we drop the horizon index $h$ in the nuisance functions and estimator that correspond to block $\ell$, with the understanding that all objects depend on $h$. For each $\ell=1,\ldots,K$:
We next state the regularity conditions used to derive the asymptotic distribution of $\hat{\theta}_h$. Following semenova2023 and huang2026, we impose geometric $\beta$-mixing on $\{z_t\}$. $\beta$-mixing is strictly stronger than the more standard strong mixing assumption, but allows us to use the Strassen coupling result that underlies the NLO theory of semenova2023. We define the $\beta$-mixing coefficients as $\beta(j)\equiv \sup_{t}\,\beta\!\left(\sigma(z_s, s\le t),\,\sigma(z_s,s\ge t+j)\right)$, where $\beta(\mathcal{A},\mathcal{B}) = E[\sup_{B\in\mathcal{B}}\left|P(B\mid \mathcal{A})-P(B)\right|]$ for any $\sigma$-algebras $\mathcal{A}$ and $\mathcal{B}$. The function $\psi(y_{t+h},x_t,\mathbf{z}_{t-1},g_0,\alpha_0,\theta_{0,h})$ in Assumption (ref) below is the doubly robust moment function defined in (ref) (see Proposition (ref)). We let $\mathcal{X}$ and $\mathcal{Z}$ denote the supports of $x_t$ and $\mathbf{z}_{t-1}$.
Assumption (ref)(i) requires the Riesz representer $\alpha_0(x,z)$ to be uniformly bounded over $\mathcal{X}\times\mathcal{Z}$, the support of $(x_t,\mathbf{z}_{t-1})$. This rules out distributions with thin-tails such as the Gaussian, but allows for Student-$t$ distributions. Assumption (ref)(ii) imposes a conditional $q$th moment bound on the regression residual $e_{t+h}$ for $q>2$, which is needed to control autocovariance terms arising from serial dependence. Assumption (ref)(iii) is a standard moment condition on the influence function that, together with Assumption (ref), ensures that a central limit theorem applies to the oracle estimator based on the true nuisance functions.
Assumption (ref)(i) imposes $L_q$ consistency on both nuisance estimators, which implies $L_2$ consistency since $q>2$. Assumption (ref)(ii) is the product rate condition standard in the double machine learning literature; see Assumption 2 of chernozhukov2024covariate. The $L_q$ rates in (i) are stronger than the $L_2$ rates typically assumed in the i.i.d.\ case and are needed here to control autocovariance terms arising from serial dependence in $e_{t+h}$ and $(x_t,\mathbf{z}_{t-1})$; see Remarks (ref) and (ref) below.
The proof of Theorem (ref) is in Appendix (ref). Letting $\tilde{\theta}_{h}$ denote the oracle estimator which assumes we know $g_0$ and $\alpha_0$, we decompose \[ \sqrt{T}(\hat{\theta}_h-\theta_{0,h}) = \sqrt{T}(\tilde{\theta}_h-\theta_{0,h})+ \sqrt{T}(\hat{\theta}_h-\tilde{\theta}_h), \]and show that $\sqrt{T}(\hat{\theta}_h-\tilde{\theta}_h)\equiv \sqrt{T}(R_1+R_2+R_3)=o_p(1)$, where $R_1$, $R_2$ and $R_3$ are remainder terms defined in Appendix (ref). Assumptions (ref), (ref) and (ref) suffice for proving that these remainders are $o_p(T^{-1/2})$.
The main implication of Theorem (ref) is that the preliminary estimation of $g_{0,h}$ and $\alpha_0$ does not affect the first-order asymptotic distribution of $\hat{\theta}_h$: the estimator is asymptotically equivalent to the infeasible oracle estimator that uses the true nuisance functions. This is a consequence of the Neyman orthogonality of the moment condition (ref) combined with the approximate independence between $I_\ell$ and $I^{\mathrm{qc}}_\ell$ under NLO cross-fitting. A feasible confidence interval for $\theta_{0,h}$ can be constructed using a standard HAC estimator of $V_h$ based on the estimated influence function $\psi(y_{t+h},x_t,\mathbf{z}_{t-1},\hat{g}_\ell,\hat{\alpha}_\ell,\hat{\theta}_h)$.
Conditional impulse response functions are often of interest in applications such as in Examples (ref) or (ref). A generalization of Definition (ref) to conditional IRFs is as follows.
Since $\varepsilon_{1t}$ is random, the conditional expectation in Definition (ref) averages over all possible realizations of $\varepsilon_{1t}$ (in addition to the other sources of randomness that enter into the potential outcomes through $U_{t+h}$), conditionally on $\Omega_{t}=\omega$. The choice of $\omega$ in $CAR_h(\delta,\omega)$ is context-dependent. For instance, in Example (ref) the conditioning set $\Omega_{t}$ is the state variable at time $t-1$, i.e.\ $\Omega_{t}=S_{t-1}$, so $\omega$ is either 0 or 1, while in Example (ref) the conditioning set is $\Omega_{t}=r_{t}$, so $\omega$ can take on any value in the support of $r_{t}$.
The following proposition provides conditions under which $CAR_h(\delta,\omega)$ is identified.
The key condition in Proposition (ref) is $x_t \perp U_{t+h} \mid (\mathbf{z}_{t-1}, \Omega_t)$, which requires that conditionally on $(\mathbf{z}_{t-1}, \Omega_t)$, all remaining variation in $x_t$ comes from $\varepsilon_{1t}$ alone. This condition holds in two leading cases. First, if $\Omega_t$ is a function of $\mathbf{z}_{t-1}$ alone, then conditioning on $(\mathbf{z}_{t-1}, \Omega_t)$ is the same as conditioning on $\mathbf{z}_{t-1}$ alone, and conditional independence follows directly from (ref) and (ref), as in Proposition (ref). This covers Example (ref), where $\Omega_t = S_{t-1}$ is a function of $\mathbf{z}_{t-1}$.\footnote{This corresponds to the case where the state of the economy depends on lags of the outcome of interest as it is the case when studying state-dependent government spending multipliers.} Second, if $\Omega_t$ contains variables not in $\mathbf{z}_{t-1}$, conditional independence is satisfied if $\Omega_t \perp \varepsilon_{1t}$. This covers Example (ref), where $\Omega_t = r_t = f(x_{t-1}) + \varepsilon_{3t}$: since $x_{t-1}$ is dated $t-1$ and $\varepsilon_{3t} \perp \varepsilon_{1t}$ by assumption, $r_t \perp \varepsilon_{1t}$ and conditional independence holds. However, this would fail if $r_t$ depended on $x_t$ (and hence on $\varepsilon_{1t}$), for instance if $r_t = f(x_t) + \varepsilon_{3t}$.
In the special case where $\Omega_t$ is binary, as in Example (ref) where $\Omega_t=S_{t-1}\in\{0,1\}$, $CAR_h(\delta,\omega)$ takes two values $CAR_h(\delta,0)$ and $CAR_h(\delta,1)$, corresponding to the impulse responses in each state. Each can be estimated by applying the NLO cross-fitting estimator of Section (ref) separately to each subsample $\{t:\Omega_t=\omega\}$, $\omega\in\{0,1\}$.
When $\Omega_t$ is continuous, as in Example (ref) where $\Omega_t=r_t$, $CAR_h(\delta,\omega)$ is a function of a continuous argument $\omega$. From (ref), it equals the conditional expectation of $g_{0,h}(x_t+\delta,\mathbf{z}_{t-1},\omega)-g_{0,h}(x_t, \mathbf{z}_{t-1},\omega)$ given $\Omega_t=\omega$, which must be estimated nonparametrically as a function of $\omega$. This introduces an additional nonparametric estimation step beyond what is required for $\mathrm{ARF}_h(\delta)$. We leave a formal treatment of this case for future work.
This section evaluates the finite sample performance of the semiparametric local projection estimator developed in Section (ref). We focus on the nonlinear regressors model of Example (ref) and the state-dependent design of Example (ref). In both cases, the structural shocks $\varepsilon_{1t}$ and $\varepsilon_{2t}$ are mutually independent, each drawn i.i.d.\ from a Student-$t$ distribution with $\nu=10$ degrees of freedom rescaled to unit variance, so that the size of the shock can be interpreted as a one-standard-deviation shock.\footnote{With $\delta=1$ and the use of a $t$-distribution for generating $\varepsilon_{1t}$, the population Riesz representer $\alpha_0(x)$ is uniformly bounded, satisfying Assumption (ref)(i). } Horizons range from $h = 0, 1, \ldots, 6$ for the first DGP where the IRFs are less persistent, whereas we set $h = 0, 1, \ldots, 8$ for the second DGP. The true impulse responses are computed by simulation from the structural model using $5{,}000$ counterfactual paths after a burn-in period of $T_0 = 500$. The DR-NLO estimator uses $K=10$ ($K=20$) contiguous blocks of equal length $\lfloor T/K \rfloor$, as described in Section (ref) for the state-dependent (nonlinear regressor) model. In both cases, $x_t$ is assumed predetermined and $x_{t-1}$ and $y_{t-1}$ are used as control variables. To estimate the conditional mean $g_{0,h}$, we use a nonparametric series estimator based on a Hermite polynomial with the total degree selected by the AIC up to a maximum value of $2$.\footnote{While Theorem 5.1 is stated for general nonparametric estimators whose complexity grows to satisfy the rate conditions in Assumption 3, our numerical implementation utilizes a low-complexity approximation to manage the bias-variance tradeoff.} The Riesz representer $\alpha_0$ is estimated by the LASSO minimum-distance procedure of chernozhukov2022, using the same Hermite dictionary at a fixed total degree of $2$. HAC standard errors for the DR-NLO estimator are computed from the cross-fitted influence function as described in Section (ref), using a Bartlett kernel and the andrews1991 automatic plug-in bandwidth. To prevent a small number of rare explosive Monte Carlo replications from dominating the results, we symmetrically trim the top and bottom 1% of the simulated IRF estimates.
This section summarizes the simulation results for a design similar to Example (ref) given by
For the sake of brevity and because simulations evaluating the performance of the plug-in and LP estimators in models with nonlinear regressors can be found in goncalves2021, goncalves2024b, we focus on the DR-NLO estimator. We use $20,000$ replications per sample size, letting $T\in\{250,500,1000,2000\}$. For each sample size, we report the median absolute bias, root mean squared error (RMSE) and nominal coverage of the 95% confidence intervals. The simulation results reported in Figure (ref) indicate that DR-NLO performs well when the object of interest is the ARF in models with nonlinear regressors. The bias is small and close to zero across horizons and sample sizes, indicating that the estimator is approximately unbiased. RMSE decreases with $T$, as expected for a consistent estimator. Although coverage falls below the nominal 95%, it improves with increasing sample size. Overall, the results suggest that the estimator performs reasonably well in modestly large samples.
We report simulation results for the state dependent model in Example (ref), where $x_t$ is predetermined so that $x_t=\rho_x x_{t-1}+\psi_x y_{t-1}+\varepsilon_{1 t}$ and $S_{t-1} = \mathbf{1}\{y_{t-1} > 0\}$ classifies the previous period as an expansion ($S_{t-1}=1$) or recession ($S_{t-1}=0$). The parameters are set to $\rho_x=0.3$, $\psi_x=-0.1$, $\beta_E=2.5$, $\beta_R=3.5$,$\gamma_E=0.9$, and $\gamma_R=-0.1$. For each replication, we apply each estimator separately to the expansion ($S_{t-1}=1$) and recession ($S_{t-1}=0$) sub-samples. We compare our estimator to the widely used state-dependent local projection (SD-LP), implemented by OLS with state-by-covariate interactions and one lag of $(x,y)$ as controls. The Monte Carlo design uses $30{,}000$ replications per sample size. For each estimator, sample size, state (expansion and recession) and horizon, we report the empirical median bias, RMSE and the coverage of nominal $95\%$ confidence intervals across replications.\footnote{In this example, we report median bias instead of median absolute bias in order to highlight the direction of the bias, namely that the LP estimator is negatively biased.}
Figure (ref) illustrates the systematic downward bias of the traditional state-dependent LP (SD-LP) estimator in both expansions and recessions, which does not vanish as the sample size grows. As explained in goncalves2024a, when the state is endogenous, the SD-LP specification does not recover the conditional average response to a fixed size shock because state-dependent dynamics in $(\beta_{t-1},\gamma_{t-1})$ introduce nonlinearities that are not accounted for by the regression. This is not a small sample issue, but a feature of the SD-LP itself. By contrast, the semiparametric DR-NLO estimator has smaller median bias at every horizon and for every sample size considered, and its bias decreases with $T$, as predicted by Theorem (ref). The substantial bias of SD-LP at shorter horizons explains why its RMSE is greater than that of DR-NLO, despite the latter having slightly higher variance. At longer horizons, the bias of the SD-LP declines, resulting in a lower RMSE. As $T$ increases, differences in RMSE between the two estimators vanish at longer horizons, with DR-NLO dominating SD-LP at shorter horizons.
Figure (ref) shows the empirical coverage probabilities of 95% confidence intervals based on SD-LP and DR-NLO, both using HAC variance estimators with a Bartlett kernel and Andrews (1991)'s automatic bandwidth. The SD-LP intervals exhibit substantial undercoverage for all sample sizes, which worsens as $T$ increases since the variance shrinks while the bias persists. DR-NLO slightly undercovers as well, but its coverage remains closer to the nominal 95% level and improves with $T$.
Semiparametric LP estimators such as the DR-NLO estimator are useful not only for capturing nonlinearities when the functional form of the nonlinearity is unknown, as illustrated by our simulation evidence, but also as a diagnostic tool for judging the adequacy of linear approximations. To illustrate how our semiparametric LP estimator may be used to assess the adequacy of the linear LP estimator, we apply both methods to study the pass-through from retail gasoline price shocks to inflation in the United States.\footnote{The reason the recent literature has focused on retail gasoline prices rather than the price of crude oil is that the relationship between crude oil and retail gasoline prices is unstable over time, reflecting large variation in the cost share of crude oil in the retail price of gasoline (see kilian2022oil).} This has been a question of continued policy interest, especially in recent years. There has been a proliferation of research addressing this question using linear VAR and distributed lag models (e.g., chudik2022estimation; kilian2022oil, kilian2025oil). The use of a semiparametric LP estimator is natural in this context since it has long been suspected that this pass-through may be nonlinear. A particular concern is that the pass-through may be stronger when the pre-existing level of inflation is higher. This concern has become particularly relevant in recent months, with the surge in headline inflation following the outbreak of the Iran War in late February 2026, but similar concerns already arose during the 2022 surge in inflation.
While this question has been addressed by a number of studies, such as clark2010time, grundler2024does or de2025energy, these studies are based on parametric nonlinear models such as threshold VAR models, regime-switching models, or time-varying coefficient VAR models. Two obvious concerns are that these specifications are mutually exclusive and that they do not exhaust the range of possible nonlinear specifications. Our semiparametric approach allows us to dispense with these parametric restrictions. We allow the effect of gasoline price shocks ($\varepsilon _{1t}$) on inflation to depend on lagged headline and core inflation, according to the model:
where $\mathbf{z}_{t-1}$ contains six lags of the percent change in retail gasoline prices ($x_t$), headline inflation ($y_{1t}$) and core CPI inflation for all urban consumers ($y_{2t}$). Core inflation is defined as inflation excluding food and energy. All data are monthly and seasonally adjusted. The data source is FRED. The model incorporates six lags consistent with other recent empirical studies. We are interested in the inflation responses at horizon $0,1,\ldots,12$. The estimation sample spans January 1974 through April 2026, which includes several high-inflation episodes.
We first compare the linear LP and DR-NLO estimates of the unconditional responses of headline and core inflation to a one percent gasoline price shock, using the same learners as in the Monte Carlo simulations and $K=10$. Evidence that the semiparametric LP estimates are very different from the linear LP estimates would cast doubt on prior linear estimates of the pass-through to inflation. Evidence that the two estimates are close, in contrast, would reassure policymakers that existing LP and VAR estimates based on linear approximations are informative. Figure (ref) illustrates that the choice of the estimator matters little. It shows point estimates of the unconditional impulse responses and 68% and 90% pointwise confidence bands. Both estimates indicate that headline inflation jumps by close to 0.05 percentage points (not annualized) in response to a one percentage point gasoline price shock, but the response declines quickly. By horizon 2, it is close to zero. There is no evidence that headline inflation rises persistently in response to gasoline price shocks. The response of core inflation is an order of magnitude smaller and only slightly positive. Both the linear LP and the DR-NLO estimate are marginally statistically significant at horizons 1 and 2. Overall, this evidence suggests that linear approximations are adequate for assessing the unconditional response of inflation to gasoline price shocks.\footnote{ Very similar results would have been obtained based on a bivariate model for gasoline prices and headline inflation.}
Next, we turn to the question of whether there is evidence of larger inflation responses conditional on the annualized inflation exceeding its historical average of 3.8% (i.e., 0.32% monthly) during the estimation period. Figure (ref) shows the unconditional response and the response conditional on inflation being above average based on the DR-NLO estimator. There is no material change in the headline and core inflation response estimates at horizons 0, 1 and 2. At longer horizons, the conditional responses tend to be larger, but also less precisely estimated. Even a user of the conditional response estimate would be unable to conclude that there are statistically significant increases in inflation in response to gasoline price shocks, however. Nor are the conditional point estimates consistent with gasoline price shocks causing persistent increases in headline or core inflation of a material magnitude.
We conclude with another example that illustrates that the impulse responses implied by the DR-NLO estimator may look materially different from those implied by the linear LP estimator and more economically plausible. The question of interest is how sales of U.S. motor vehicles respond to shocks to the real price of motor gasoline. This question is motivated by edelstein2009sensitive and ramey2006declining, ramey2011oil who document potential nonlinearities in the response of motor vehicle sales. We rely on equations ((ref)) and ((ref)) with $z_{t}$ containing the percent change in real U.S. retail gasoline prices and the percent change in U.S. sales of automobiles and light trucks. Retail gasoline prices are obtained from the CPI and deflated by the aggregate CPI. All data are monthly and seasonally adjusted. The data source is FRED. The model incorporates six lags. The estimation period is January 1974 through March 2026.
Figure (ref) focuses on the responses of the cumulative growth rate to a one standard deviation shock in the real price of gasoline.\footnote{The cumulative effect is estimated directly by applying the LP or DR-NLO to $\Delta^hy_{t+h}=y_{t+h}-y_{t-1}$.} The DR-NLO impact response is positive but statistically insignificant. There is no statistically significant response in sales for the first three months, according to the DR-NLO estimator, followed by a persistent decline in auto sales starting at horizon 3. The linear LP estimate, in contrast, suggests a statistically significant but economically counterintuitive increase in auto sales at horizons 0, 1 and 2. More generally, the persistently positive linear LP responses for the first four months are difficult to reconcile with economic reasoning. The numerical differences between the estimates compound at longer horizons. At horizon 12, the linear LP estimate shows a decline in auto sales that is only 28% of the DR-NLO response. These differences suggest that linear approximations may not adequately capture the dynamics in question.
This paper developed a semiparametric local projection estimator of unconditional nonlinear impulse response functions for a broad class of structural dynamic models that are widely used in applied macroeconomics, including models with nonlinearly transformed regressors, state-dependent coefficients, and nonlinear interactions between shocks and state variables. Under standard mixing and rate conditions on the nuisance estimators, the resulting estimator is $\sqrt{T}$-consistent and asymptotically normal, with the preliminary estimation of the two nuisance functions having no effect on the first-order asymptotic distribution. Inference is conducted via a HAC long-run variance estimator applied to the estimated influence function. We also showed how our framework accommodates conditional impulse responses in state-dependent models, where the conditioning variable takes discrete values.
Our Monte Carlo evidence indicates that the proposed estimator delivers substantially smaller bias than standard state-dependent local projection methods, while preserving competitive root mean squared errors and confidence-interval coverage close to the nominal level. We considered two empirical illustrations. The first one focused on the potentially nonlinear pass-through from gasoline price shocks to inflation. The second example examined whether linear LP estimators miss nonlinearities in the transmission of gasoline price shocks to motor vehicle sales.
There are several avenues for further research. First, our identification result relies on the additive separability of $\varepsilon_{1t}$ in equation (ref). One possible extension would be to allow the shock to interact nonlinearly with lagged state variables. Second, our treatment of conditional impulse responses assumes a discrete conditioning variable, leaving the continuous-state case as discussed in Example (ref) for future work.
\singlespacing \doublespacing \setcounter{equation}{0}