Dynamic Causal Effects in a Nonlinear World: the Good, the Bad, and the Ugly

Dynamic Causal Effects in a Nonlinear World: the Good, the Bad, and the Ugly

Keywords: dynamic treatment effect, impulse response, local projection, semiparametric identification, structural vector autoregression.

Introduction

Impulse response functions are key objects in macroeconomic analysis. Since they measure dynamic causal effects of surprise changes in policy or fundamentals on subsequent macroeconomic outcomes, they provide calibration targets for structural modeling and help validate model predictions. They also inform optimal economic policy questions, both directly and indirectly Christiano1999,McKay2023.

Applied researchers typically report impulse response estimators motivated by linear time series models, such as \acp{VAR} or local projections. Although there exists a wealth of nonlinear alternatives (Fan2003; Herbst2016; Kilian2017), linear methods are attractive due to their simplicity and the difficulty of clearly detecting nonlinear relationships in typical macroeconomic data. At the same time, both macroeconomic theorists and policymakers think nonlinearities are important: structural models with essential nonlinearities have become dominant in recent decades, and many economic policy debates concern state-dependence and asymmetries. How can we justify using linear methods if we think the world is a nonlinear place?

This paper studies the causal interpretation of impulse response estimators based on linear models when the data is generated by an essentially unrestricted nonparametric structural model. We first deliver good news for linear local projection or \ac{VAR} estimators that project directly on an observed shock or proxy: their estimand (i.e., probability limit) equals a weighted average of the true nonlinear causal effects, regardless of the extent of nonlinearities in the \ac{DGP}. By contrast, the news is bad or even ugly for estimators that identify latent shocks via heteroskedasticity or non-Gaussianity: they generally do not estimate a meaningful causal summary under departures from linearity. Thus, the hard work needed to directly measure shocks (or proxies) using historical or institutional data buys insurance against nonlinearities that other identification approaches lack.

Our good news are based on an extension of the results in Yitzhaki1996 and Rambachan2021: impulse response estimands from linear local projections and \acp{VAR} that project on observed shocks or proxies correspond to positively-weighted averages of marginal effects---causal effects of changing the shock variable from $x$ to $x+\delta$ for infinitesimally small $\delta$, averaged out over all (past, present, and future) shocks other than the contemporaneous shock of interest, and weighted over different baseline values $x$ for the shock variable. Thus, these estimands provide a scalar causal summary of the full richness of the nonlinear causal effects; the positive weights ensure that the researcher gets the sign right if the true marginal effects are uniformly positive or negative. Our assumptions drop restrictions imposed in the existing literature that ruled out models with kinks or discontinuous regime-switches or shocks with unbounded support.

In a nonlinear \ac{DGP}, both the sign and the magnitude of the causal effects can depend on the baseline shock value (e.g., whether it is positive or negative), so how these values are weighted can matter a lot. Fortunately, as we illustrate using several empirical examples, the weight function used by local projections and \acp{VAR} is straightforward to estimate and report. In many applications, the researcher does not directly observe the shock of interest but only a proxy, also known as an external instrument Stock2018. In this case, we show that an easily-interpretable monotonicity condition is required to guarantee a positive weight function. However, we also show that when control variables are needed to isolate a true shock (i.e., recursive or Cholesky identification), then positive weights can only be guaranteed if the linearly residualized shock is nonlinearly unpredictable by the controls, which may be a strong assumption in the absence of detailed institutional knowledge and high-quality data.

One implication of these results is that linearity-based estimators are useful even when economic theory predicts a nonlinear relationship between the shock and the outcome of interest. For example, if the outcome variable has limited support, such as when it is binary or censored (say, due to a zero lower bound), nonlinearities are inherently present. If one is interested in characterizing the nonlinearities, then it makes sense to model them, and it is of course always a good idea to plot the raw data regardless. However, if one is interested in an overall summary of marginal effects, then linear local projections and \acp{VAR} are theoretically coherent estimators, as discussed earlier. In fact, we show that directly modeling nonlinearities can be counterproductive unless the researcher is confident in their modeling: under functional form misspecification, local projections with higher-order terms still estimate a weighted average of marginal effects, but some of the weights may be negative, which risks getting the sign of the causal effects wrong. This echoes the message from an earlier JBES lecture by angrist_dummy_2001 that in a cross-section context with limited dependent variables, linear methods provide more robust estimates of treatment effects than nonlinear ones.

When there is a dearth of direct shock measures or proxies, applied researchers frequently resort to identification via heteroskedasticity Sentana2001,Rigobon2003,Lewbel2012. Unfortunately, we show that these estimation approaches are sensitive to the assumption that the structural model is linear: the estimand can easily be nonzero when there is no causal effect, or negative when the true shock has a uniformly positive effect on the outcome of interest. Fixing these issues while still delivering informative inference appears difficult, since a natural nonparametric generalization of the identification strategy yields very wide identified sets. The intuition for these negative results is that the identification exploits a source of exogenous variation that shifts the scale of the latent shock of interest but not its mean. Without strong functional form assumptions, this type of exogenous variation is uninformative about the effect of a location shift in the shock on the conditional mean of the outcome, i.e., the impulse response. However, a silver lining is that the linear model delivers testable restrictions.

The sensitivity to nonlinearity is even greater for identification via non-Gaussianity Comon1994,Gourieroux2017,Lanne2017. Also known as \ac{ICA}, this identification approach has recently increased significantly in popularity in the \ac{VAR} literature. We show that the nonparametric analogue of the identification assumptions yields an identified set so large that effectively any function of the data can be construed as a “shock”. Intuitively, the mere assumptions that the latent shocks are independent and non-Gaussian are vacuous in a nonparametric context: any collection of random variables can always be represented as some nonlinear function of independent uniformly distributed random variables. Moreover, we give examples of simple \acp{DGP} featuring slight nonlinearity for which any linearity-based \ac{ICA} procedure is highly biased asymptotically, yet in these \acp{DGP} one cannot reject the validity of the linear model.

The building block underlying most of the above findings is a set of results on the identification of weighted averages of marginal treatment effects using weighted regressions, which connects our analysis to a large literature in microeconometrics Yitzhaki1996,NeSt93,Angrist1999,GoldsmithPinkham2024. We extend existing results in this literature by unifying the treatment of continuous, discrete, and mixed regressors, and by substantially weakening the regularity conditions: we allow for regressors with unbounded support, impose minimal regularity on the regression function, and our moment conditions essentially only require the existence of the probability limit of the regression estimator.

An important limitation of our results is that they only concern identification. While we are motivated by the observation that full-fledged nonparametric estimation is challenging in realistic macroeconomic data sets, we do not explicitly analyze the precision or small-sample bias of the estimators we study. We refer to Herbst2024 for a discussion of finite-sample biases of local projections and \acp{VAR} in linear models. Another limitation is that we do not consider identification via non-recursive short-run restrictions, long-run restrictions, or sign restrictions.

\paragraph{Literature.} Pioneering work on semiparametric causal time series analysis includes Gallant1993, Potter2000, White2006, White2009, Angrist2011, and Angrist2018, see also Goncalves2021,goncalves24nonparametric,Goncalves2024, Gourieroux2023, and Kitagawa2023 for recent contributions. Our result on the causal interpretation of local projections with observed shocks is very closely related to Rambachan2021 and subsequent work by Caravello2024 and Casini2024.

As for identification via heteroskedasticity or non-Gaussianity, we are not aware of other work in a nonparametric vein. \Citet{MontielOlea2022} criticize linearity-based versions of these identification strategies for being seemingly sensitive to functional form assumptions, and potentially being subject to weak identification. The present analysis quantifies this sensitivity more precisely by deriving both the identified sets for the nonparametric analogues of these identification assumptions, and the estimands of linearity-based procedures.

\paragraph{Outline.} (ref) defines a nonparametric framework for identification of dynamic causal effects. (ref) argues that local projection and \ac{VAR} estimands based on observed shocks or proxies have a robust causal interpretation regardless of the extent of nonlinearities. (ref) show, on the other hand, that estimands based on identification through heteroskedasticity or non-Gaussianity are sensitive to the assumption that the structural function is linear. (ref) provides the theoretical basis for the results in the earlier parts of the paper by extending results from the microeconometric literature on the interpretation of regression estimators as weighted marginal treatment effects; this section may be of independent interest for readers outside macroeconomics. (ref) concludes. Technical details and proofs are relegated to the online supplement.

Nonparametric framework for dynamic causality

In this section we set up a nonparametric framework for dynamic causal identification.

Model

We are interested in the dynamic response of a scalar outcome variable $Y_{t}$ to an impulse in the scalar shock variable $X_t$. As a leading example, one may think of $X_t$ as a variable controlled by a policy-maker, such as a surprise change in the policy interest rate set by the central bank. For ease of exposition, we restrict attention to continuously distributed shocks $X_t$ for now, but (ref) shows that our results generalize to handle continuous, discrete, or mixed distributions in a unified manner.

The outcome variable is determined by an underlying dynamic structural model. Our causal framework doesn't restrict this model; we only assume that when evaluated $h$ periods after the realization of the shock $X_{t}$, the outcome admits the nonparametric structural representation

For each horizon $h$, $\psi_{h}(\cdot, \cdot)$ is an unknown measurable function that we call the structural function, while $\boldsymbol{\mathbf{U}}_{h, t+h}$ is a vector of all variables (dated before, on, and after time $t$) that causally affect $Y_{t+h}$, other than $X_t$. Without restrictions on $\boldsymbol{\mathbf{U}}_{h, t+h}$ or the structural function, the representation (ref) is without loss of generality. In typical recursive time-series models, however, $\boldsymbol{\mathbf{U}}_{h, t+h}$ will contain the vector $\boldsymbol{\mathbf{Y}}_{t-1}$ of observed data at time $t-1$ as well as shocks dated $t, t+1,\dotsc, t+h$, but exclude $X_t$ and shocks dated after $t+h$ White2006,White2009,Caravello2024,Goncalves2024. A leading special case is the linear structural \ac{VAR} model, which additionally implies that $\psi_{h}$ is linear in both $X_{t}$ and $\boldsymbol{\mathbf{U}}_{h, t+h}$ Kilian2017.

We assume throughout that it is meaningful to think of varying $X_t$ while keeping $\boldsymbol{\mathbf{U}}_{h, t+h}$ constant, so that the random function $Y_{t+h}(x)\equiv \psi_h(x, \boldsymbol{\mathbf{U}}_{h, t+h})$ defines a potential outcome function at horizon $h$. \Citet{Angrist2011}, Angrist2018, and Rambachan2021 work with this potential outcome notation $Y_{t+h}(x)$, and keep all other past, present, and future shocks (captured by $\boldsymbol{\mathbf{U}}_{h, t+h}$ in our model) implicit. Our structural function framework (ref) is mathematically equivalent, but facilitates comparisons with the linear structural \ac{VAR} literature.

For now, we focus on the case where $X_{t}$ is a “shock” (such as a surprise change in a policy instrument), and assume that it is independent of $\boldsymbol{\mathbf{U}}_{h, t+h}$:

This assumption is common in the literature, and essentially just normalizes the structural function $\psi_{h}$, so that its first argument captures the total causal effect of the shock $X_t$ on $Y_{t+h}$, including its direct effect and any indirect effects, both contemporaneous and dynamic. This is illustrated in the following simple example.

When $X_{t}$ is not a shock but corresponds to a policy instrument (such as the policy interest rate) or any other kind of serially correlated variable, it will typically be correlated with $\boldsymbol{\mathbf{U}}_{h, t+h}$, both because of the serial correlation and because it may be correlated with other determinants of $Y_{t+h}$---we discuss this case in (ref).

Causal effects

A familiar issue in nonlinear models is that there are multiple possible definitions of an impulse response, i.e., a dynamic causal effect. In a linear model, the effect of exogenously changing $X_{t}$ from $x_{0}$ to $x_{1}$ is a linear function of the difference $x_{1}-x_{0}$: it equals $\psi_h(x_1,\boldsymbol{\mathbf{U}}_{h, t+h})-\psi_h(x_0,\boldsymbol{\mathbf{U}}_{h, t+h})=\theta_h(x_{1}-x_{0})$ for some constant scalar $\theta_{h}$. By contrast, in a nonlinear model, the effect is in general a nonlinear function of both $x_{1}$ and $x_{0}$ (not just their difference), and it also depends on the past history and the current and future nuisance shocks via $\boldsymbol{\mathbf{U}}_{h, t+h}$.

In theoretical macroeconomic modeling, researchers often report the impulse responses with respect to a so-called “MIT shock”, which starts the economy at steady state, then hits the economy with a one-off impulse to $X_t$, and subsequently sets all other current and future shocks to zero: $\psi_h(X_{t}, \boldsymbol{\mathbf{0}})-\psi_h(0,\boldsymbol{\mathbf{0}})$, where we normalize the steady-state values of $X_{t}$ and $\boldsymbol{\mathbf{U}}_{h, t+h}$ to zero. While computationally convenient, this impulse response concept has no empirical counterpart; moreover, for the measure to be policy relevant, it is necessary that the model satisfy certainty equivalence.

We shall instead focus on impulse responses (causal effects) defined as expected counterfactual changes in the outcome of interest, averaging out over all other shocks. Specifically, define the average structural function

which corresponds to the expected potential outcome function. Here the expectation is taken over the marginal distribution of $\boldsymbol{\mathbf{U}}_{h, t+h}$. The expectation is implicitly assumed to exist for all $x$. The average structural function measures the counterfactual average value of the future outcome $Y_{t+h}$ that we would observe if the policy-maker engineered a particular fixed value $x$ for the policy variable at time $t$, averaging out over the randomness caused by all other factors that influence the outcome independently of the policy decision at time $t$. Even though in some nonlinear models the structural function $\psi_h$ is discontinuous in the policy variable $x$, the average structural function $\Psi_h$ will typically be a smoother function of $x$, as it averages out over the realizations of other shocks. For example, this is the case in the regime-switching model in (ref) if $\xi_t$ is continuously distributed. A certain amount of smoothness in $\Psi_h$ will be important for the identification of causal effects, as we discuss below.

Typical data samples in macroeconomics are too small to permit accurate nonparametric estimation of the entire average structural function $x \mapsto \Psi_h(x)$. A pragmatic alternative is to target weighted averages of the structural function---average causal effects---or its derivatives---average marginal effects Rambachan2021,Goncalves2021,Goncalves2024. This paper focuses on estimation of average marginal effects

where $\omega(\cdot)$ is a weight function averaging across the baseline values of the shock variable $X_{t}$. We reserve the term average marginal effect to weight functions that are convex, i.e., $\omega(x)$ is nonnegative for all $x$ and integrates to one, $\int \omega(x)\, dx=1$. This ensures that $\theta_{h}(\omega)$ is a meaningful causal summary of the average structural function $\Psi_{h}(x)$ in that it prevents what strlb17 call a sign-reversal: if $\Psi_{h}'(x)$ has the same sign for all $x$ ($+,0$ or $-$), then $\theta_{h}(x)$ will also have this sign.\footnote{blandhol2022tsls call estimands with convex $\omega$ “weakly causal”. Convex weighting schemes also satisfy what rslr07 call “boundedness”: $\theta_{h}(\omega)$ lies in the support of $\Psi_{h}'(x)$.} This property is particularly useful when qualitatively validating predictions of structural macroeconomic models.

Depending on the form of the weight function, $\theta_{h}(\omega)$ has two interpretations in terms of an average causal effect of a shock with magnitude $\delta>0$,

First, the average marginal effect corresponds to the average causal effect for infinitesimally small shocks: $\theta_{h}(\omega)=\lim_{\delta\to 0}\theta_{h}(\delta, \omega)$, provided we can pass the limit as $\delta\to 0$ under the integral sign in (ref). Second, if the weighting in (ref) admits the integral representation $\omega(x)=\frac{1}{\delta}\int_{x-\delta}^{x}\omega_{0}(x)\, dx$, substituting $\Psi_h(x+\delta)-\Psi_h(x)=\int_{x}^{x+\delta}\Psi'(\chi)\, d\chi$ into (ref) and changing the order of integration yields $\theta_{h}(\omega)=\theta_{h}(\delta, \omega_{0})$. For this reason, focusing on average marginal effects is without loss of generality.

In a linear model, the weighting does not matter, since $\Psi_{h}'(x)$ does not depend on $x$. But in nonlinear models, it could matter greatly whether we attach most weight to positive or negative shocks, or to shocks with small or large magnitude. Therefore, accounting for the form of the weighting $\omega$ is important when using estimates of $\theta_{h}(\omega)$ to calibrate or validate structural macroeconomic models. In the next section, we discuss identification approaches that deliver weighted averages of marginal effects under a particular weighting scheme that depends on the shock distribution. In (ref), we discuss estimation approaches that target any pre-specified weighting scheme.

The good: observed shocks and proxies

If the researcher directly observes the shock of interest, or at least a valid proxy for it, then there is good news: conventional local projections or structural \ac{VAR} impulse responses estimate average marginal effects with an interpretable weighting scheme, regardless of how nonlinear the underlying \ac{DGP} is. Moreover, the weights can be estimated from the data, and we give several empirical examples of how to interpret them. In contrast to linear estimators, we demonstrate using a simple example that nonlinear extensions of local projections or \acp{VAR} do not generally provide meaningful causal summaries under misspecification. Finally, we extend the analysis to shocks that are recursively identified, i.e., by controlling for covariates.

Identification with observed shocks

We start off by assuming that the researcher directly observes (or consistently estimates) the shock $X_t$ of interest. This would be the case, for example, if the shock is identified through a “narrative approach”. See Ramey2016 for several empirical examples.

Under the nonlinear structural model (ref) and the shock independence assumption (ref), the conditional expectation of the outcome given the shock,

nonparametrically identifies the average structural function:

Hence, in principle, we could estimate any weighted causal effect of interest by running a nonparametric regression of $Y_{t+h}$ on $X_t$ to obtain $g_h(\cdot)$ in the first step, and then averaging this function according to the desired weighting scheme in the second step, as suggested by goncalves24nonparametric,Goncalves2024. In (ref), we discuss a complementary strategy that identifies the same estimand via weighted averages of the observed outcomes, and how both strategies can be combined. However, as discussed in more detail in (ref), these strategies may yield noisy and sensitive estimates in the relatively small samples available in macroeconomics.

\paragraph{Interpretation of linear projection estimates.} We take a cue from Rambachan2021 and instead aim for a less ambitious goal. Rather than targeting a pre-specified weighted average of causal or marginal effects, we focus on simple local projection and \ac{VAR} estimators, which are relatively precise even with small sample sizes. We demonstrate that these simple estimators have an attractive robustness property: even though they are motivated by a linear model, when the \ac{DGP} is nonlinear, their estimand can still be interpreted as an average marginal effect with a particular weight function.

The local projection estimator of Jorda2005 estimates the impulse response of $Y_t$ with respect to $X_t$ at horizon $h$ as the coefficient $\hat{\beta}_h$ in the \ac{OLS} regression

where $\boldsymbol{\mathbf{W}}_t$ is a vector of control variables (typically including a constant and lagged outcomes and shocks). For now, we will assume that the shock $X_t$ is in fact a “shock”, so that it is linearly unpredictable using the controls: $\operatorname*{Cov}(X_t, \boldsymbol{\mathbf{W}}_t)=0$. Then the set of controls $\boldsymbol{\mathbf{W}}_t$ affects only the precision of $\hat{\beta}_{h}$, but not its probability limit. In particular, under standard stationarity and ergodicity assumptions, the local projection estimator $\hat{\beta}_h$ will converge in probability to the population projection coefficient

PMW2021 show that a \ac{VAR} which includes $X_t$ ordered first has the exact same population estimand (ref), provided that the number of lags in the \ac{VAR} is sufficiently large.\footnote{If $X_t$ is linearly unpredictable from lagged data, it is sufficient that the lag length weakly exceed $h$.} It is a textbook result that the linear function $\beta_{h}x$ provides the best linear approximation to the potentially nonlinear average structural function $g_{h}(x)=\Psi_{h}(x)$ Angrist2009, so that it approximates the average structural function in a prediction sense. However, this result is not directly informative about whether $\beta_{h}$ has a causal interpretation---whether it can be interpreted as an average marginal effect if $\Psi_{h}(x)$ is nonlinear.

The following proposition shows that the local projection and \ac{VAR} estimand (ref) achieves our goal: it has a causal interpretation as an average marginal effect (ref). The result is not new---it appeared previously in Yitzhaki1996 and Rambachan2021; as we discuss below, the novelty lies in substantively weakening the regularity conditions.

Combined with the identification result (ref) for the average marginal effect, (ref) shows that linear local projections and \acp{VAR} remain useful in a nonlinear world: they estimate an average causal effect $\theta_h(\omega_X) = \int \omega_X(x)\Psi_h'(x)\, dx$ for infinitesimal shocks, with a convex weighting scheme $\omega_{X}$. Furthermore, the scheme gives most weight to shocks close to the mean $E[X_t]$, with little weight given to extreme values. In the special case where $X_{t}$ is normally distributed, (ref) reduces to Stein's lemma stein81: the weight function $\omega_X$ reduces to the normal density function, so that $\beta_{h}$ equals the expected marginal effect, $E[\Psi_{h}'(X_{t})]$, as noted by Yitzhaki1996. The fact that the weighting scheme depends only on the marginal distribution of $X_{t}$ and not the particular outcome variable $Y_{t+h}$ or horizon $h$ allows for comparisons of average marginal effects for different outcomes or across different horizons $h$.\footnote{Since the weighting does not depend on the outcome variable or horizon, multipliers (i.e., ratios of cumulative impulse response estimands, see Jorda2023) can be expressed as a ratio of two average marginal effects with the same weight function, where we now additionally average across horizons.} If the true \ac{DGP} is in fact linear, then the weighting of course does not matter, and we recover the conventional linear impulse response.

While we focus here on interpreting the proposition in the context of the causal model in (ref), the result does not require the structural assumptions (ref)--(ref). This is relevant in settings in which the conditional mean $g_h(x)=E[Y_{t+h} \mid X_t=x]$ is a useful descriptive object even if it does not have a direct causal interpretation.

The assumption that $X_{t}$ is continuously distributed can be dropped without changing the result, as we show in (ref). In cases where there are gaps in the support of $X_{t}$, such as when the shock is discrete or mixed, one just needs to extend the definition of the conditional mean function $g_{h}$ to the whole interval $I$ by linear interpolation. To our knowledge, this unification of the treatment of continuous, discrete, and mixed distributions is novel.

Even in the case of a continuously distributed shock, the assumptions in (ref) are substantively weaker than those in the literature, and accommodate all textbook linear models as well as a wide range of nonlinear models. The assumption that $g_{h}$ is locally absolutely continuous is necessary to ensure that weighted marginal effects are well-defined. As discussed earlier, this assumption will typically hold even in models with discrete regimes or kinks, since the expectation (ref) averages over the distribution of the nuisance shocks. The moment conditions and integrability condition $\int_{I} \omega_X(x)\abs{g_h'(x)}\, dx<\infty$ just ensure that the estimand $\beta_{h}$ and the weighted marginal effect exist.\footnote{\Cref*{lemma:finite_integral_omega} in \Cref*{app:proofs} shows that for the integrability condition to hold, it is sufficient to assume the tails of $g_{h}(x)$ are monotone.} In contrast, the original work by Yitzhaki1996 does not provide a formal proof or regularity conditions on $g_{h}$ (neither does the discussion by Angrist2009). Analogous results in Rambachan2021, Graham2022, Caravello2024, and Casini2024 require the potential outcome function (not its expectation) to be smooth, which rules out models with kinks or discrete regimes, and require the interval $I$ to be bounded, which rules out the textbook case of normally distributed shocks.\footnote{Our result is not strictly more general than that of Casini2024, since they allow for non-stationary potential outcomes.} The restrictiveness of these conditions led Goncalves2024 to question the applied relevance of the causal interpretation of the estimand (ref), but our weaker conditions demonstrate that this concern is unfounded when interest centers on average marginal effects.

\paragraph{Estimating the weight function.} As argued by Angrist1999 for the case of discrete $X_t$, the weight function $\omega_X$ defined in (ref) can be estimated in the data. This allows the researcher to gauge which weighted causal effect is being estimated: does it attach most weight to negative or positive shocks, small or large shocks? Since the weight function depends only on the shock variable itself and not the outcome variable or the impulse response horizon, it is only necessary to estimate a single function. We therefore recommend that researchers always estimate and plot this function.

Estimation is simple: $\omega_X(x)$ equals the slope coefficient in a (population) regression of the indicator $\1{X_t \geq x}$ on $X_t$. This regression can be implemented in the data via \ac{OLS}, separately for each value $x$ of $X_t$ observed in the data.\footnote{Pointwise confidence intervals can be obtained with conventional heteroskedasticity-robust standard errors. One could also use autocorrelation robust standard errors to allow for time series dependence of $X_t$, but causal interpretation is more challenging if the shocks are not independent.} In applications, it may also be of interest to report an integral $\int_{\underline{x}}^{\overline{x}} \omega_X(x)\, dx$ of the weight function over an interval $x \in [\underline{x}, \overline{x}]$. \Cref*{app:weight_integr} shows that we can estimate this integral by the slope coefficient in an OLS regression of $M_t \equiv \max\lbrace \min\lbrace X_t, \overline{x} \rbrace, \underline{x} \rbrace$ on $X_t$. In particular, to estimate the total weight $\int_0^\infty \omega_X(x)\, dx$ given to positive shocks, we simply regress $M_t \equiv \max\lbrace X_t,0\rbrace$ on $X_t$.

To illustrate, we now empirically estimate the weight function $\omega_X$ for various macroeconomic shocks considered in the handbook chapter by Ramey2016. We use Ramey2016's replication code and data off the shelf. In particular, prior to computing weights, all shocks are residualized on the same control variables that she uses in her VARs and local projections. The estimates of the weight functions are obtained from OLS regression output, as described above. To demonstrate the ease of implementation, all steps of the computations are carried out in Stata, like Ramey2016's replication code.\footnote{Our code and data are available at \url{https://github.com/mikkelpm/nonlinear_dynamic_causal}}

(ref) shows the estimated weight functions for four identified government spending shocks from the applied literature. Note that the shocks are not entirely comparable due to differences in their precise definitions and sample periods. The Blanchard2002 and Fisher2010 shocks, which are intended to capture general government spending shocks, yield approximately symmetric weight functions. By contrast, the BenZeev2017 and Ramey2011 shocks, which capture news about future defense spending, generate weight functions that are skewed towards positive shocks. In fact, both these shocks exhibit a large positive outlier in the 3rd quarter of 1950 (the onset of the Korean War), reflected in the fat right tail of the weight functions. In other words, impulse responses from local projections or \acp{VAR} estimated off the latter two shocks will largely reflect the causal effects of sharp military buildups, rather than retrenchments. This is important to remember when using empirical impulse responses to discipline structural models that feature asymmetries (such as downward nominal wage rigidity or borrowing constraints), since then model-implied impulse responses with respect to positive government spending shocks will differ from those for negative shocks. \Cref*{app:weight_empir} gives further examples of weight functions for several identified tax, technology, and monetary policy shocks. As these examples illustrate, plotting estimates of the weights $\omega_{X}$ is useful in interpreting the results of any subsequent impulse response analysis and for comparing with prior studies.

\paragraph{Parametric nonlinear specifications.} In many cases, economic theory predicts that the average structural function $\Psi_{h}$ is likely nonlinear. For example, if the outcome variable has limited support, such as due to censoring or when it is discrete, the structural function must necessarily be nonlinear. In such cases, it seems natural to model the nonlinearity directly, rather than to stick to a linear specification as in (ref). For example, Jorda2005 and Jorda2025 suggest including powers of the shock in local projections. Similarly, there is a rich literature on nonlinear extensions of \ac{VAR} models, see for example Kilian2017. Such direct modeling of the nonlinearities is sensible if the goal of the analysis is to directly characterize the extent and types of nonlinearity present in the data, e.g., threshold effects or sign and size dependence Caravello2024.

However, for estimating average causal effects, simple linear local projections or \acp{VAR} appear more robust than nonlinear parametric specifications. As shown in (ref), the linear specification in (ref) is robust to misspecification in that it estimates a well-defined average marginal effect regardless of the form of nonlinearity in the structural function $\Psi_{h}$. By contrast, we now show that this is not the case for a local projection specification that includes a quadratic term, echoing similar results in angrist_dummy_2001 regarding robustness of parametric nonlinear limited dependent variable models. These results suggest that nonlinear specifications do not generally have such a robustness property.

Consider a quadratic local projection of $Y_{t+h}$ on $X_t$, $X_t^2$, and an intercept.\footnote{While we focus on the quadratic case for simplicity, (ref) below can be shown to generalize to a polynomial specification of any fixed order.} Assume for analytical simplicity that $X_t$ has a standard normal distribution, so in particular $X_t$ and $X_t^2$ are uncorrelated (our qualitative conclusions can be shown to go through without the normality assumption). Then the population version of the projection is

with implied derivative of the regression function at $X_t=x$ given by

and population regression coefficients

The first expression in (ref) shows that the estimated derivative $\bar{\beta}_h(x)$ equals a weighted average of the true derivative function $g_h'(\cdot)$, but with weights that are negative whenever $1+X_{t}x<0$.\footnote{It also follows from the proposition that any estimated weighted average derivative $\int \omega(x) \bar{\beta}(x)\, dx$ that is a nontrivial function of the coefficient $\beta_{2,h}$ (i.e., whenever $\int x\omega(x)\, dx \neq 0$) equals a weighted average of $g_h'(\cdot)$ with weights that are negative for some $x$.} If the true regression function $g_h$ is in fact quadratic, then $\bar{\beta}_{h}(x)$ is consistent for the marginal effect function $\Psi_h'(x)$. But if the regression function is misspecified, the negative weighting leads to a sign reversal: the second expression in (ref) implies that even if $g_h$ is monotonically increasing, the estimated derivative $\bar{\beta}_h(x)$ will be negative for sufficiently large $x$ whenever $E[g_h''(X)]<0$. Such sign reversal is not shared by the linear estimator (ref), for which the weighting scheme $\omega_{X}$ is convex. This lack of robustness of a quadratic (or more generally polynomial) specification of the regression function to functional form misspecification is related to the observation in white_using_1980 that polynomial approximations to the conditional mean function $g_{h}$ cannot be interpreted as providing a Taylor series approximation to $g_{h}$.

\paragraph{State-dependent specifications.} A particularly popular nonlinear local projection specification in applied work is a state-dependent specification that interacts the shock with a binary regime indicator $S_t \in \lbrace 0,1 \rbrace$ (see Cloyne2023, Goncalves2024, and references therein):

For example, $S_t$ may indicate whether the economy is in a recession. Assuming that the local projection is fully interacted as above (i.e., all control variables $\boldsymbol{\mathbf{W}}_t$ are interacted with $S_t$), then the procedure is tantamount to running separate regressions on the subsamples with $S_t=0$ and $S_t=1$, respectively.\footnote{If we instead omit the interaction terms from the regression and only control linearly for $S_t$, then we are in the case of (ref) below.} It follows that all the analysis surrounding (ref) above applies upon conditioning on $S_t=s \in \lbrace 0,1\rbrace$. In particular, the probability limit of the state-dependent impulse response estimate $\hat{\beta}_{s, h}$ equals a positively weighted average of conditional marginal effects $\partial E[Y_{t+h} \mid X_t=x,S_t=s]/\partial x$, which have a clear causal interpretation provided the shock independence assumption (ref) holds conditional on $S_t$ (i.e., within each regime). Thus, despite their apparent linearity conditional on regime, state-dependent local projections identify causal estimands even when the true \ac{DGP} has a nonlinear form, such as a model with smooth or discrete regime-switching. However, consistent with the discussion in (ref), it is important to interpret the impulse responses as averaging over all future shocks, including potential future regime switches. In other words, the local projection estimand does not hold the regime fixed within the impulse response horizon.

Identification with proxies

In many applications, observations of the shock are contaminated by measurement error, such as when accurate measurements are available only in a subset of the time periods. In such cases, researchers typically treat the measurements $Z_{t}$ as a proxy for the shock of interest $X_{t}$, or, equivalently, an instrument for the shock Stock2018. We now show that when the structural function is nonlinear, linear \acp{VAR} and local projections onto the proxy identify average marginal effects up to scale, provided that the conditional mean of the proxy given the shock is monotone in the shock.

We assume that the proxy $Z_{t}$ is valid, in the sense that it satisfies the exclusion restriction

formalizing the notion that if the shock $X_{t}$ were observed, the proxy $Z_{t}$ would not provide any further predictive power. It is implied by the standard assumption in the measurement error literature that the measurement error in $Z_{t}$ is non-differential, i.e., the whole conditional distribution of $Y_{t+h}$ given $(X_{t}, Z_{t})$ depends only on $X_{t}$ (or equivalently that $Z_{t}$ is independent of $\boldsymbol{\mathbf{U}}_{h, t+h}$) CaRuSt06.

We consider the “reduced-form” local projection of the outcome $Y_{t+h}$ on the proxy $Z_t$. Under (ref), the population version of this regression has slope coefficient

where

As shown by PMW2021, $\tilde{\beta}_{h}$ also corresponds to the probability limit of an impulse response from a structural \ac{VAR} where the proxy is ordered first, and the specification controls for sufficiently many lags.

Combining (ref) with the identification result (ref) implies that linear proxy regressions identify weighted averages of marginal effects, $\theta_h(\tilde{\omega}_Z) = \int \tilde{\omega}_Z(x)\Psi_h'(x)\, dx$, just as in the case of directly observed shocks. Unlike in the observed shocks case, the weights $\tilde{\omega}_{Z}$ will not be positive unless the proxy satisfies the condition in point ((ref)) of (ref)---this condition is slightly weaker than monotonicity of $\zeta(x)$.\footnote{For instance, the condition may still hold even if monotonicity is violated over a sufficiently small interval in the middle of the support of $X_{t}$.} However, monotonicity of $\zeta$ ensures not just that the weights are positive, but also that they have an intuitive hump-shape, giving most weight to shocks in the middle of the distribution.

Monotonicity of $\zeta$ is implied by, but much weaker than the continuous-treatment version of the ImAn94 monotonicity condition, needed for causal interpretation of two-stage least squares estimands under endogeneity. We defer the details to \Cref*{app:ident-with-instr}, where we generalize the identification results in AnGrIm00 by allowing for non-smooth potential outcome functions and non-binary $Z_{t}$; we also extend (ref) to the case with covariates. It follows from this identification result that monotonicity of $\zeta$ holds under much weaker conditions than those required for causal interpretation of $\tilde{\beta}_{h}$ under endogeneity. \Citet[Theorem 7]{Rambachan2021} derive an alternative characterization of the proxy estimand (ref) involving derivatives of the reduced-form potential outcome as a function of the proxy $Z_t$ (rather than of the shock $X_t$), and therefore the monotonicity assumption has no counterpart in their analysis.

A practical implication of (ref) is that applied researchers should seek to construct proxies that are credibly positively related to the unobserved latent shock of interest. However, it is not essential that the relationship is linear or indeed of any particular known functional form. Of course, constructing valid proxies may be challenging in practice: measurement error in their construction may fail to be non-differential, leading to violation of the exclusion restriction (ref), and lack of variability of $Z_{t}$ in the data may lead to large standard errors for $\tilde{\beta}_{h}$. The purpose of (ref) is to clarify the value of this hard work, if done well.

The weight function (ref) does not integrate to 1 due to attenuation bias, so that we only identify average marginal effects up to scale. However, since the weight function doesn't depend on the outcome, this is not an issue in practice: we can scale $\tilde{\beta}_{h}$ by the response of some normalization variable to the proxy (this is the so-called unit effect normalization) to identify a relative marginal effect. The local projection instrumental variable estimator of Stock2018, which is a two-stage least squares version of local projection, automatically performs this normalization.

Since the shock $X_t$ is not directly observed, we cannot generally estimate the weight function $\tilde{\omega}_Z$ in the data. Instead, it may be useful to plot the observed-shock weight function (ref) pretending that $Z_{t}$ is the actual shock of interest. If it happens that $Z_t \approx X_t$, then these weights will be close to the proxy weights $\tilde{\omega}_Z$, so the plot provides a “best-case” scenario.

Identification with control variables

In applications where it is challenging to isolate purely exogenous shifts in policy or fundamentals, researchers may be willing to assume that the observed variable $X_t$ (which could be a policy instrument) is exogenous conditional on some control variables $\boldsymbol{\mathbf{W}}_t$ (such as variables that comprise the policy-makers information set):

This is a selection on observables assumption as in Angrist2011 and Angrist2018. For example, the assumption holds if $X_t = \Upsilon(\varepsilon_t, \boldsymbol{\mathbf{W}}_t)$, where $\varepsilon_t$ is a shock that is independent of $(\boldsymbol{\mathbf{W}}_t', \boldsymbol{\mathbf{U}}_{h, t+h}')'$, a nonparametric version of the recursive (or Cholesky) assumption in linear structural \ac{VAR} identification Christiano1999. Then the conditional expectation function

equals the conditional average structural function in the causal model (ref):

where the expectation is taken with respect to the conditional distribution of the nuisance shocks $\boldsymbol{\mathbf{U}}_{h, t+h}$ given $\boldsymbol{\mathbf{W}}_{t}$.

Even under the selection on observables assumption (ref), the local projection with controls (ref) need not estimate an average marginal effect if the relationship between $X_t$ and the controls is nonlinear. This result extends to recursively identified structural \acp{VAR}, due to the nonparametric equivalence between these procedures PMW2021. (ref) below shows that the population local projection coefficient $\beta_h$ can still be written as a weighted average of the marginal effects $\partial g(x, \boldsymbol{\mathbf{w}})/\partial x=\partial \Psi_h(x, \boldsymbol{\mathbf{w}})/\partial x$, but the weights can be negative if the true “propensity score” $\pi^*(\boldsymbol{\mathbf{w}}) \equiv E[X_t \mid \boldsymbol{\mathbf{W}}_t = \boldsymbol{\mathbf{w}}]$ is nonlinear (or, in other words, if the purported “shock” that is the residual from a linear projection $X_t$ on the controls $\boldsymbol{\mathbf{W}}_t$ is predictable by nonlinear transformations of the controls). We leave the details, which extend the analysis of GoldsmithPinkham2024 to cases with non-discrete $X_t$, to (ref). As usual, negative weights are worrying, as they may lead to a sign-reversal. If $\boldsymbol{\mathbf{W}}_{t}$ just consists of a set of mutually exclusive dummies, then linearity of the propensity score comes for free, and the weights are guaranteed to be positive. However, outside this case, we recommend that researchers do careful sensitivity checks with respect to both the set of controls and the functional form for the controls (e.g., whether the addition of nonlinear transformations of the controls, such as interactions and polynomials, improves prediction of $X_{t}$). More research is warranted on best practices for such sensitivity checks, given the data limitations in applied macroeconomics.

The bad: identification via heteroskedasticity

Identification via heteroskedasticity has become a popular procedure for causal identification in applications where direct shock measures are unavailable, following Sentana2001, Rigobon2003, and Rigobon2004.\footnote{Similar identification approaches were developed in the signal-processing literature in the 1990s, see the review by Hyvarinen2001.} In a pair of highly-cited papers, Lewbel2012,Lewbel2018 exploits this idea to achieve identification in cross-sectional regressions with endogenous variables and no external instruments Klein2010. In stark contrast to (ref), the results in this section deliver bad news regarding the sensitivity of identification approaches via heteroskedasticity to the assumption that the underlying structural function is linear: the Rigobon-Sack-Lewbel estimator does not generally estimate average marginal effects; more generally, we show that the nonparametric analogue of the identification approach yields very large identified sets for causal effects. However, on the positive side, it is possible to test the linearity assumption in the data.

Nonparametric version of the identification approach

To explain the sensitivity of conventional identification via heteroskedasticity to the linearity assumption on the structural function, it is helpful to first lay out a nonparametric version of the framework before we review the linear case. Since this section is mainly concerned with giving examples of how the identification approach can fail, we specialize the dynamic set-up from (ref) to a simpler static model.

\paragraph{Nonparametric setup.}

We observe an $n$-dimensional vector $\boldsymbol{\mathbf{Y}}$ of variables that are nonlinearly related to a latent, scalar shock of interest $X$ as well as an $(m-1)$-dimensional latent vector $\boldsymbol{\mathbf{U}}$ of nuisance shocks:

where we suppress time subscripts to ease notation. The above model is a (static) nonparametric factor model, since we do not impose parametric restrictions on the unknown structural function $\boldsymbol{\mathbf{\psi}} \colon \mathbb{R}^m \to \mathbb{R}^n$.

The econometrician observes a scalar $D$ that is informative about the heteroskedasticity of the shock of interest $X$ but independent of the nuisance shocks $\boldsymbol{\mathbf{U}}$ (jointly with $X$):

This assumption implies that $D$ is a valid proxy for $X$, in the sense that $\boldsymbol{\mathbf{Y}}$ and $D$ are independent conditional on $X$. But because the variable $D$ only influences the variance and higher moments of $X$ but not its mean,

we cannot use the proxy in local projections as in (ref). For concreteness, it may be useful to think of $D$ as a binary regime indicator, which affects the conditional variance $\operatorname*{Var}(X \mid D)$ but not the conditional mean (ref), as in the original work by Rigobon2003.

If we assume that the structural function $\boldsymbol{\mathbf{\psi}}$ is linear, it is possible to achieve identification even if $D$ is unobserved, and we relax (ref) by allowing $D$ to affect the variances of nuisance shocks. See Bacchiocchi2024 and Lewis2024 for excellent reviews. However, since we are only interested in showing how the basic identification approach can fail in a nonparametric context, we maintain the stronger assumptions above. It then follows a fortiori that nonparametric identification is even more challenging under weaker assumptions.

\paragraph{Review of linear identification.} If the structural function $\boldsymbol{\mathbf{\psi}}$ in (ref) is known to be partially linear, identification of causal effects obtains under an additional relevance assumption. Thus, we temporarily assume that

where $\boldsymbol{\mathbf{\theta}}$ is the unknown vector of causal effects of $X$, while $\boldsymbol{\mathbf{\gamma}} \colon \mathbb{R}^{m-1} \to \mathbb{R}^n$ is an unknown function. Following Rigobon2004 and Lewbel2012, construct the scalar instrumental variable

where $Y_1$ is the first element of $\boldsymbol{\mathbf{Y}}$. In applications, $Y_1$ may be a policy instrument that is known to be strongly related to $X$, though it is also allowed to be correlated with the nuisance shocks. Under the linear model (ref) and the identification assumptions (ref)--(ref), $Z$ satisfies the exogeneity restriction for linear identification in Stock2018 since $E[Z \mid \boldsymbol{\mathbf{U}}]=0$. In particular, under these assumptions, a regression of $\boldsymbol{\mathbf{Y}}$ on $Y_1$ using $Z$ as an instrument identifies the (relative) causal effects of $X$:

To ensure we are not dividing by zero, we need to additionally assume the relevance conditions that (i) the shock of interest is heteroskedastic across regimes, $\operatorname*{Cov}(X^2,D) \neq 0$, and (ii) the causal effect of $X$ on $Y_1$ is nonzero, $\theta_1 \neq 0$. For completeness, we review the calculations leading to (ref) in \Cref*{app:hetero_identif_linear}.

Fragility under nonlinearity

We now argue that the simple linear identification argument fundamentally cannot be extended to nonparametric contexts.

\paragraph{Nonparametric identified set.} We first show that the nonparametric model of identification via heteroskedasticity yields a large identified set for the causal effects of $X$ on $\boldsymbol{\mathbf{Y}}$. To do this, we strengthen the independence and conditional mean assumptions (ref)--(ref) by imposing a specific model for the relationship between $X$ and $D$:

$\sigma \colon \mathbb{R} \to \mathbb{R}_+$ is a known function, and $R$ has a known distribution that is symmetric around 0. This model would, for example, be consistent with the conditionally Gaussian model $X \mid D \sim N(0,\sigma^2(D))$.

The proposition states that the identified set for $\boldsymbol{\mathbf{\psi}}$ is so large that it always contains a structural function $\boldsymbol{\mathbf{\psi}}(x, \boldsymbol{\mathbf{u}})$ that is symmetric in $x$ around 0. In particular, we can never rule out that the average marginal effect $\int \omega(x) (\partial E[\boldsymbol{\mathbf{\psi}}(x, \boldsymbol{\mathbf{U}})]/\partial x)\, dx$ is zero when the weight function $\omega(x)$ is symmetric around 0. Intuitively, the challenge is that $D$ does not affect the mean of $X$, only higher moments, so---without strong functional form restrictions on the relationship between the outcomes and the shocks---we do not have enough information to sign mean effects of shifts in the latent shock $X$. This holds even though we assume that the econometrician knows exactly how $D$ affects the dispersion of the $X$ distribution. Notice that the construction of the observationally equivalent symmetric structural function in (ref) only relies on a single (scalar) nuisance shock; hence, knowledge about the true number of shocks does not ameliorate the identification failure (see (ref) for further discussion of this point).

A careful inspection of the proof of (ref) reveals that the result is closely related to a known issue with identification via heteroskedasticity in a linear context: while the variance of the shock of interest $X$ must vary across regimes, we cannot simultaneously allow the impulse responses of $X$ to vary across regimes Lewis2024. However, in a nonparametric context this problem is even worse, since there is no fundamental distinction between “coefficients” and “shock variances” in a general nonlinear model. A priori restrictions that certain “coefficients” are independent of the regime are only meaningful once we parametrize the model, which complicates the development of an empirically useful nonparametric generalization of the identification approach.

\paragraph{Sensitivity of linear procedures.}

Because the nonparametric identified set is large, we can expect estimation procedures based on linearity of the structural function to fail to estimate causal objects in general. The next result implies that this is indeed the case for the linear instrumental variable estimator (ref) of Rigobon2004 and Lewbel2012.

(ref) shows that regressing $\boldsymbol{\mathbf{Y}}$ onto the instrument $Z$ yields a weighted average of marginal effects, but with a weight function $\check{\omega}(x)$ that cannot be guaranteed to be positive.\footnote{This is not a special case of (ref), since $Z$ does not satisfy the nonparametric proxy assumption (ref).} In fact, the weights even integrate to 0 in some cases, for example if $\Psi_1(x)=\Psi_1(-x)$ and the conditional distribution of $X$ given $D$ is symmetric around 0. In such cases, the instrument erroneously estimates a zero causal effect of $X$ on $Y_j$ for $j \geq 2$ even if $\Psi_j(x)=\theta_j x$ is a linear function with $\theta_j \neq 0$.

The weights can also be negative---and therefore cause the econometrician to get the sign of the marginal effects wrong---even in the seemingly favorable setting where (unbeknownst to the econometrician) the policy variable $Y_{1}$ simply equals the shock of interest $X$, without any nonlinearity or contamination by nuisance shocks, i.e., $\Psi_{1}(x)=x$. Assume in addition, as in Rigobon2003, that the regime indicator $D \in \lbrace 0,1\rbrace$ is binary and $E[X \mid D]=0$. Then a simple calculation shows that the weights in (ref) equal

where $f_{X|d}(x)$ is the density of $X$ conditional on regime $D=d$. Suppose the right (resp., left) tail of the $X$ distribution is fatter (resp., thinner) in regime $D=1$ than in regime $D=0$, meaning that $f_{X|1}(x) > f_{X|0}(x)$ for $x \gg 0$ and $f_{X|0}(x) > f_{X|1}(x)$ for $x \ll 0$. Then it follows from (ref) that $\check{\omega}(x)>0$ for $x \gg 0$, while $\check{\omega}(x)<0$ for $x \ll 0$. This simple example shows that the instrumental variable estimator can easily generate negative weights, even when it satisfies the exclusion and relevance conditions and the policy variable is linear in the shock. To trust that the weights are positive, we would need to have quite detailed information about the conditional shock density in the two regimes; simple moment restrictions do not suffice.

Intuitively, the problem of negative weights comes about because the Rigobon2003 and Lewbel2012 instrumental variable $Z$ defined in (ref) fails the proxy monotonicity assumption discussed earlier in connection with (ref). Because the only source of exogenous variation is the regime indicator $D$, and this indicator does not affect the mean of the latent shock $X$ but only higher moments, it is generally impossible to construct any proxy variable that is guaranteed to be monotone in $X$, unless we make strong assumptions about the structural function.

If the model is not additively separable as assumed in (ref), the instrumental variables estimator can exhibit even more pathological behavior, in that it may not equal a weighted average of marginal effects at all. As a simple example, consider a multiplicative model $\boldsymbol{\mathbf{Y}} = X\boldsymbol{\mathbf{\gamma}}(\boldsymbol{\mathbf{U}})$ with $E[\boldsymbol{\mathbf{\gamma}}(\boldsymbol{\mathbf{U}})]=\boldsymbol{\mathbf{0}}$ and impose the independence assumption (ref). In that model, $E[\boldsymbol{\mathbf{Y}} \mid X] = \boldsymbol{\mathbf{0}}$, so the marginal effect function is identically zero, but $\operatorname*{Cov}(\boldsymbol{\mathbf{Y}}, Z)=\operatorname*{Cov}(X^2,D)\operatorname*{Cov}(\gamma_1(\boldsymbol{\mathbf{U}}), \boldsymbol{\mathbf{\gamma}}(\boldsymbol{\mathbf{U}})) \neq 0$ in general, so the instrument erroneously estimates a nonzero effect.

Silver lining: testability of the linearity assumption

While the sensitivity of identification via heteroskedasticity to linearity of the structural function $\boldsymbol{\mathbf{\psi}}$ is disheartening, at least the linear model (ref) implies testable restrictions. Specifically, as noted by Rigobon2004 and Wright2012, for any $d_0,d_1$ in the support of $D$, the difference $\operatorname*{Var}(\boldsymbol{\mathbf{Y}} \mid D=d_1)-\operatorname*{Var}(\boldsymbol{\mathbf{Y}} \mid D=d_0) = [\operatorname*{Var}(X \mid D=d_1)-\operatorname*{Var}(X \mid D=d_0)]\boldsymbol{\mathbf{\theta}}\boldsymbol{\mathbf{\theta}}'$ should be a rank-1 matrix under linearity and the maintained independence assumption (ref). Other over-identification tests in more general linear models of identification via heteroskedasticity are discussed in the review article by Lewis2024. We are not aware of any thorough analysis of the power properties of these tests against nonlinear alternatives.

The ugly: identification via non-Gaussianity

A second approach to identification in linear models in the absence of direct shock measures is to assume that the structural shocks are mutually independent and non-Gaussian; see Gourieroux2017, Lanne2017, and the review article by Lewis2024. \Citet{Lewbel2024} propose a similar approach to achieve identification in cross-sectional endogenous regressions in the absence of external instruments. An earlier literature outside economics goes by the name \acf{ICA}, see Kagan1973, Comon1994, and the textbook by Hyvarinen2001.

This section delivers ugly news regarding the sensitivity of this identification approach to linearity of the structural function: once the linearity assumption is dropped, the non-Gaussianity assumption is essentially vacuous; as a consequence, estimators based on non-Gaussianity and linearity of the structural function can fail spectacularly even under mild departures from linearity. What is worse, the linearity assumption is untestable in general.

Nonparametric version of the identification approach

As in (ref), we consider the nonparametric factor model (ref). However, unlike in the case of identification via heteroskedasticity, we now do not observe any additional proxy variables that aid in identifying the latent shocks. Instead, we hope to achieve identification via restrictions on the distributions of the shocks.

\paragraph{Review of linear identification.} Assume temporarily that the number of shocks equals the number of observables, $m=n$, and that the structural function is linear:

Assume also that the $n$ shocks $(X, U_1, \dotsc, U_{n-1})$ are mutually independent, and at most one of these shocks has a Gaussian distribution. We will refer to this model as the linear \ac{ICA} model. A deep result in probability theory, the Darmois-Skitovich Theorem, says that two nontrivial linear combinations of independent variables cannot themselves be independent, unless all the underlying variables are Gaussian. In the context of the linear \ac{ICA} model, the theorem implies that any two linear combinations $\boldsymbol{\mathbf{\varsigma}}'\boldsymbol{\mathbf{Y}}$ and $\tilde{\boldsymbol{\mathbf{\varsigma}}}'\boldsymbol{\mathbf{Y}}$ of the data $\boldsymbol{\mathbf{Y}} = \boldsymbol{\mathbf{\theta}}X + \boldsymbol{\mathbf{\gamma}}\boldsymbol{\mathbf{U}}$ can be independent if and only if these linear combinations equal two different shocks in the model (up to sign and scale). Hence, the shocks in the model can be identified by searching for those linear combinations of the observed variables that are independent; once we have the shocks, we can then estimate their causal effects. See Hyvarinen2001 and Lewis2024 for reviews of estimation procedures.

The abstract identification argument above can be made less mysterious through a method of moments framework that exploits implications of shock independence for higher moments of the data Lewis2024. Nevertheless, it is clear from both the abstract argument and the more concrete moment-based approach that linearity of the structural function is being leveraged heavily.

Fragility under nonlinearity

\paragraph{Nonparametric identified set.} Unfortunately, the mere assumptions that the latent shocks are independent and non-Gaussian provide essentially no identification power in a nonparametric context. The identified set under these assumptions is so large that nearly any function of the data can be labeled a “shock”.

(ref) shows that the nonparametric factor model (ref) is very under-identified, even if we restrict the number of latent, independent shocks to equal $m=n$ and impose smoothness on the structural function $\boldsymbol{\mathbf{\psi}}$. Indeed, any element of almost any one-to-one transformation $\boldsymbol{\mathbf{\Upsilon}}(\boldsymbol{\mathbf{Y}})$ of the observables could be construed as a “shock” for some data-consistent choice of structural function $\boldsymbol{\mathbf{\psi}}$. Note that the identification issues go well beyond the familiar “labeling problem” in linear ICA analysis, where the shocks are identified up to sign and permutation, such that additional economic information is required to label each of the statistically identified shocks Lewis2024. While the challenge of identifying nonlinear factor models is well known in the broader literature (see the review by Jutten2003), it appears that the serious consequences of this fact for shock identification in macroeconometrics have not been explored previously.

The fundamental issue is that non-Gaussianity of the shocks is a vacuous assumption in the nonparametric setting: it is an innocuous normalization to assume that all shocks have uniform distributions, since we can always nonlinearly transform any shock distribution to the uniform distribution via the quantile function. In other words, we have severe identification failure as in (ref) even if the econometrician knows the exact distributions of each shock. Hence, it is no accident that the identification argument in the \ac{ICA} and structural \ac{VAR} literatures relies heavily on the linearity assumption: there is no nonlinear equivalent of the Darmois-Skitovich Theorem.

If we allow for slightly less smoothness of the structural function $\boldsymbol{\mathbf{\psi}}$, then the identification problem is even worse. As Gunsilius2023 note, any $n$-dimensional vector $\boldsymbol{\mathbf{Y}}$ can be represented as a nonlinear factor model (ref) in a single latent shock $X$ (so $m=1$ and $\boldsymbol{\mathbf{U}}=\boldsymbol{\mathbf{0}}$) using a so-called space-filling curve (e.g., Hilbert curve) construction, though the associated $\boldsymbol{\mathbf{\psi}}$ function would not be one-to-one. Hence, without restrictions on the structural function, we cannot rule out that the latent shock of interest $X$ drives all the variation in the $n$ observed variables $\boldsymbol{\mathbf{Y}}$. Indeed, given a uniform random variable $U$ on $[0,1]$, we can generate an infinite number of independent uniform random variables $\{V_{j}\}$ from the decimal expansion of $U=0.u^{1}_{1}u_{2}^{1}u_{1}^{2}u_{3}^{1}u_{2}^{2}u_{1}^{3}u_{4}^{1}u_{3}^{2}u_{2}^{3}u_{1}^{4}\dotsb$, and taking $V_{j}\equiv 0.u^{j}_{1}u^{j}_{2}\dotsb$ (and hence an infinite number of independent random variables with arbitrary distributions $F_{j}$ by taking the inverse transform $F_{j}^{-1}(V_{j})$).\footnote{This is a consequence of the fact that there exists a one-to-one function $\boldsymbol{\mathbf{\phi}}$ such that both $\boldsymbol{\mathbf{\phi}}$ and $\boldsymbol{\mathbf{\phi}}^{-1}$ are measurable between the measurable spaces $(M, \mathcal{B}_{M})$ and $([0,1], \mathcal{B}_{[0,1]})$, where $M$ is any separable complete metric space and $\mathcal{B}_{M}$ is the Borel $\sigma$-algebra on $M$ dudley02.}

\paragraph{Sensitivity of linear procedures.} Due to the nonparametric identification failure, we can expect identification approaches based on non-Gaussianity to be very sensitive to exact linearity in the structural function. The following two simple examples illustrate such sensitivity in two settings where the linearity assumption is untestable. Thus, identification via non-Gaussianity is not only fragile, it is also generally not falsifiable.

In conclusion, linearity-based \ac{ICA} identification procedures can be highly misleading under departures from a linear model, as with identification via heteroskedasticity (but unlike identification with observed shocks or proxies). In fact, the situation is arguably worse than in (ref), since even arbitrarily small structural nonlinearities can yield large biases, and the linearity assumption is not testable in general.

Identification of average marginal effects

In (ref), we reverse-engineered the weight function allowing us to interpret the linear local projection estimand as an average marginal effect. We now consider a forward-engineering problem: how can we estimate an average marginal effect with a pre-specified weight function? For simplicity, we focus on the case when the shock of interest is observed, but we drop the requirement that the shock is continuously distributed. We first consider the case without control variables before extending the analysis to allow for controls.

Identification without controls

We consider the setup from (ref), but drop time subscripts to make it clearer that our analysis applies to cross-sectional as well as time series settings. Indeed, part of our goal is to show that identification arguments from cross-section settings carry directly through to time-series contexts. Let

denote the conditional mean function from a nonparametric regression of the scalar outcome $Y$ onto the scalar variable $X$. We do not restrict the marginal distribution of $X$: it can be continuous, discrete, or mixed. Let $I\subseteq\mathbb{R}$ denote a (possibly unbounded) interval that contains the support of $X$. We are interested in summarizing $g$ by reporting its weighted average derivative, weighted by some pre-specified weight function $\omega$. With some abuse of notation, we still denote this weighted average derivative by

as in (ref), even though we don't require that $g(x)$ corresponds to some structural function. To ensure that this object is well-defined, we assume that $g$ is locally absolutely continuous on $I$. Since (ref) only defines $g$ on the support of $X$, this requires us to extend $g$ to all of $I$ in cases when there are gaps in the support of $X$, such as when $X$ is discrete. This can be done by linear interpolation: if $P(X\in (a, b))=0$ for some $(a, b)\subseteq I$, we set $g(x)=(g(b)-g(a))(x-a)/(b-a)+g(a)$ for $x\in(a, b)$. If the distribution of $X$ is discrete, this defines the derivative $g'$ as the slope between adjacent support points (and the extension will automatically be locally absolutely continuous provided that the spacing between adjacent support points is bounded away from 0).

A regression-based approach to estimating $\theta(\omega)$ first estimates the entire derivative function $g'(\cdot)$ nonparametrically (by, say, series or kernel regression), and then averages it using the weights $\omega$. The next result shows that we can alternatively estimate $\theta(\omega)$ as a weighted average of outcomes, $\theta(\omega)=E[\alpha(X)Y]$, where $\alpha$ is the Riesz representer of the linear functional $g\mapsto \theta(\omega)$.

Analogous representations for $\theta(\omega)$ are well-known in the literature if we additionally assume that $X$ is continuously distributed NeSt93. The representation is usually derived by directly applying integration by parts. Our proof instead generalizes the proof of Stein's lemma stein81, allowing us to drop the requirement that $X$ is continuously distributed and impose only very mild regularity conditions, which essentially just require that both sides of (ref) are well-defined. In particular, absolute continuity of $g$ is needed to ensure that $\theta(\omega)$ is well-defined, and in \Cref*{lemma:finite_integral_omega} in \Cref*{app:proofs}, we show that if the tails of $\omega$ or the tails of $g$ are monotone, then condition (iv) of (ref) holds provided the integral on the right-hand side of (ref) exists.

(ref) gives a recipe for constructing weighting-based estimators of $\theta(\omega)$ for particular choices of weight function $\omega$ by replacing the expectation in (ref) with a sample average and, if the function $\alpha$ is unknown, replacing $\alpha$ with an estimate.\footnote{As we discuss in (ref) below, recent results in the semiparametric literature suggest that rather than picking between this weighting-based approach and the regression-based approach to estimation of $\theta(\omega)$, it may pay off to combine them, yielding a “doubly-robust” estimator.} For instance, suppose $X$ is continuous, and let $\omega(x)=f_{X}(x)$ correspond to the density of $X$, so that $\theta(\omega)=E[g'(X)]$ is the (unweighted) average derivative. Then the required weighting is given by $\alpha(x)=-f_{X}'(x)/f_{X}(x)$, leading to the estimator of hardle1989 if one uses kernel estimators to estimate the density and its derivative. If the identification condition (ref) holds, this estimator identifies the average causal impact of increasing $X$ by an infinitesimal amount. To estimate the average impact of increasing $X$ by a fixed amount $\delta$ (i.e., the unweighted average causal effect), which corresponds to setting $\omega(x)=\frac{1}{\delta}\int_{x-\delta}^{x}f_{X}(x)\, dx$, let $\alpha(x)=-\frac{f_{X}(x)-f_{X}(x-\delta)}{\delta f_{X}(x)}$, replacing the derivative of the density by a discrete change. If we set $\omega(x)=f_{X}^{2}(x)$, so that $\alpha(x)=-2f_{X}'(x)$, and we use a leave-one-out kernel estimator for the derivative of the density, we recover the famous density-weighted average derivative estimator of PoStSt89. Finally, for $\omega(x)=E[\1{X\geq x}X]$, we get $\alpha(x)=x$, and (ref) reduces to (ref): this weighting scheme can be estimated by linear regression. As the proofs of (ref) reveal, these results are also special cases of (ref).\footnote{As a consequence, the assumption that $X_{t}$ be continuous in (ref) can be dropped.}

Identification with control variables

We now generalize the setup to allow for a vector of controls $\boldsymbol{\mathbf{W}}$. Consider the weighted average derivative

where the expectation is over the marginal distribution of $\boldsymbol{\mathbf{W}}$, and $g'$ is the derivative with respect to $x$ of the conditional mean function $g(x, \boldsymbol{\mathbf{w}}) \equiv E[Y\mid X=x, \boldsymbol{\mathbf{W}}=\boldsymbol{\mathbf{w}}]$.\footnote{Weighting by the marginal distribution of $\boldsymbol{\mathbf{W}}$ is not restrictive, since weighting schemes that use other forms of averaging across $\boldsymbol{\mathbf{w}}$ can be recovered by defining $\omega$ appropriately.} To ensure this object is well-defined, we assume that for each $\boldsymbol{\mathbf{w}}$, the weights $\omega$ are zero outside the interval $I_{\boldsymbol{\mathbf{w}}}$ containing the conditional support of $X$ given $\boldsymbol{\mathbf{W}}=\boldsymbol{\mathbf{w}}$, and that we can extend $g(\cdot, \boldsymbol{\mathbf{w}})$ to $I_{\boldsymbol{\mathbf{w}}}$ such that $g(\cdot, \boldsymbol{\mathbf{w}})$ is locally absolutely continuous on $I_{\boldsymbol{\mathbf{w}}}$, such as by linearly interpolating across any gaps in the support.

Like in the case without covariates, a regression-based estimator of $\theta(\omega)$ first estimates the derivative of the regression function $g'(x, \boldsymbol{\mathbf{w}})$, and then averages the estimated derivative function using the weights $\omega$ and the marginal distribution of the covariates. The next result shows that we can alternatively estimate $\theta(\omega)$ by taking weighted averages of the outcome.

As discussed in (ref), the representation (ref) is well-known if the distribution of $X$ is continuous conditional on $\boldsymbol{\mathbf{W}}$. The novelty of (ref) is to drop the continuity requirement and relax the regularity conditions.

If $X\in\{0,1\}$ is a binary treatment variable, and we additionally assume that $X$ is as good as randomly assigned conditional on $\boldsymbol{\mathbf{W}}$, then $g(1,\boldsymbol{\mathbf{w}})-g(0,\boldsymbol{\mathbf{w}})$ corresponds to the conditional \ac{ATE} for individuals with $\boldsymbol{\mathbf{W}}=\boldsymbol{\mathbf{w}}$. In this case, the average derivative simplifies to

which corresponds to a weighted average of conditional \acp{ATE}. By letting $\alpha(X, \boldsymbol{\mathbf{W}})=X/P(X=1\mid \boldsymbol{\mathbf{W}})-(1-X)/P(X=0\mid \boldsymbol{\mathbf{W}})$, (ref) recovers the classic result that we can estimate the (unweighted) \ac{ATE} by inverse probability weighting. If $X$ is continuous with density $f_{X}(x\mid \boldsymbol{\mathbf{W}})$ conditional on $\boldsymbol{\mathbf{W}}$, letting $\alpha(x, \boldsymbol{\mathbf{W}})=-f_{X}'(x\mid \boldsymbol{\mathbf{W}})/f_{X}(x\mid \boldsymbol{\mathbf{W}})$ recovers the average derivative $E[g'(X, \boldsymbol{\mathbf{W}})]$.

For both of these special cases, there is a wealth of papers studying how to best implement regression-based or weighting-based approaches to estimating $\theta(\omega)$, or combinations of both. Recent influential results in the cross-sectional literature ceinr22,ccddhnr18 highlight the advantages of combining both approaches using the Neyman orthogonal moment condition

where $\mu(X, \boldsymbol{\mathbf{W}}, g)=g(1,\boldsymbol{\mathbf{W}})-g(0,\boldsymbol{\mathbf{W}})$ for the \ac{ATE} and $\mu(X, \boldsymbol{\mathbf{W}}, g)=g'(X, \boldsymbol{\mathbf{W}})$ for the average derivative. This moment condition is orthogonal in the sense that it is insensitive to small perturbations in $g$, in contrast to the regression-based moment condition $\theta(\omega)= E[\mu(X, \boldsymbol{\mathbf{W}}, g)]$. As a result, an orthogonal method-of-moments estimator based on (ref) that plugs in first-stage estimates of $g$ and $\alpha$ can be viewed as a debiased version of the plug-in estimator utilizing the regression-based moment condition. Actually, the moment condition (ref) is not only orthogonal, but also doubly robust---insensitive to large perturbations in either $g$ or $\alpha$ so that the orthogonal method-of-moments estimator remains consistent so long as any one of the first-stage estimators is consistent for $\alpha$ or $g$, even if the other estimator is inconsistent. For the binary treatment case, the orthogonal method-of-moments estimator corresponds to the classic augmented inverse probability weighted estimator of rrz94.

For i.i.d.\ data, it is popular to combine the orthogonal moment condition with cross-fitting ceinr22,ccddhnr18.\footnote{Take a sample sum of the moment condition (ref) over the first half of the sample, plugging in estimates $\hat{\alpha}_{2}$ and $\hat{\gamma}_{2}$ of $\alpha$ and $\gamma$ constructed using the second half of the sample, and add to it a sample sum of the moment condition over the second half of the sample, but where we plug in estimates $\hat{\alpha}_{1}$ and $\hat{\gamma}_{1}$ from the first half of the sample.} This allows for regularity conditions that are weak enough to accommodate a variety of first-step estimators of $\alpha$ and $\gamma$, including kernels, series, as well as lasso, random forests, or other machine learning estimators, provided that these estimators converge sufficiently fast. Because of this flexibility, the approach is known as debiased machine learning. \Citet{cns22} and HiWa21 develop alternatives to this approach that bypass the need to explicitly estimate the Riesz representer $\alpha$.

These approaches all deliver estimators of $\theta(\omega)$ that converge, under appropriate regularity conditions, at the usual parametric rate (square root of sample size) even if the first-stage estimators are based on complicated nonparametric or machine learning algorithms. Recent work by goncalves24nonparametric, BaWe24, and ballarin2025 aims to adapt these semiparametric approaches to time series contexts. One potential worry is that even in the absence of covariates, given the small sample sizes typically available in macroeconomic applications, estimates of average marginal effects relying on machine-learning or nonparametric first-step estimates of the shock density and the structural function may yield estimates that are too noisy and sensitive to the choice of first-stage tuning parameters. When covariates are needed to argue that the observed variable $X$ is exogenous, the data requirements become even more severe.

The practical challenges associated with fully nonparametric estimation motivates studying what the simple \ac{OLS} local projection (ref) estimates when the true regression function is nonlinear. Extending the analysis of (ref), we now allow the researcher to control for covariates more flexibly by considering the partially linear regression

and $\Gamma$ is a linear space of control functions that contains the constant function $1$. This covers the case with a linear adjustment by letting $\Gamma=\{a+\boldsymbol{\mathbf{w}}'\boldsymbol{\mathbf{b}}\colon a\in\mathbb{R}, \boldsymbol{\mathbf{b}}\in\mathbb{R}^{\dim(\boldsymbol{\mathbf{w}})}\}$ be the class of linear functions of $\boldsymbol{\mathbf{w}}$, as well as the semiparametric partially linear model that lets $\Gamma$ be a large class of “nonparametric” functions. By the projection theorem, the estimand $\beta$ in this regression is given by

where $\pi(\boldsymbol{\mathbf{W}}) \equiv \operatorname*{argmin}_{\gamma_{0}\in\Gamma}E[(X-\gamma_{0}(\boldsymbol{\mathbf{W}}))^{2}]$ denotes the projection of $X$ onto $\Gamma$ (if $X$ and $\boldsymbol{\mathbf{W}}$ are independent, the estimand (ref) reduces to that in (ref); we assume the projection exists). We denote the true conditional expectation (the propensity score, if $X$ is binary) by $\pi^{*}(\boldsymbol{\mathbf{W}})\equiv E[X\mid \boldsymbol{\mathbf{W}}]$.

The next result uses (ref) to generalize (ref) to the case with covariates.

To interpret this result, it is useful to first consider the case where the class $\Gamma$ of control functions is rich enough so that condition (a) holds. This is trivially the case, for instance, if $\boldsymbol{\mathbf{W}}$ just consists of a single set of fixed effects. In this case, we obtain an analogue of (ref): (partially) linear regression identifies a weighted average of marginal effects, with hump-shaped weights---this holds regardless of the distribution of $X$, be it discrete, continuous or mixed. If the conditional distribution of $X$ given $\boldsymbol{\mathbf{W}}$ is continuous, then the weighting $\omega(x, \boldsymbol{\mathbf{W}})=\omega^{*}(x, \boldsymbol{\mathbf{W}})/\operatorname*{Var}(X-\pi^{*}(\boldsymbol{\mathbf{W}}))$ varies smoothly with $x$; if there are mass points, as in the case of a discrete or mixed distribution, then the weight function jumps discontinuously at the mass points. If the treatment $X$ is discrete with support $0, 1, \dotsc$, we recover the result in Angrist1999 that regression estimates a weighted average of the causal effects of increasing $X$ by one unit,

Now suppose that condition (a) is violated, because the specification for $\Gamma$ is not sufficiently flexible to model the true propensity score $\pi^{*}(\boldsymbol{\mathbf{W}})$. In general, this may lead to omitted variable bias, as the partially linear model may not be sufficiently flexible to account for all confounding due to $\boldsymbol{\mathbf{W}}$. Condition (b) prevents this scenario, ensuring that any bias due to confounding is accounted for. As a simple example of when the condition holds, consider the case where the true conditional mean function has a multiplicative form: $g(X, \boldsymbol{\mathbf{W}})=X g'(\boldsymbol{\mathbf{W}})+\gamma_{0}(\boldsymbol{\mathbf{W}})$, with $\gamma_{0}\in\Gamma$. The partially linear model is misspecified, because the marginal effect is not constant, but varies with $\boldsymbol{\mathbf{W}}$. But condition (b) holds with $\gamma_{1}(\boldsymbol{\mathbf{w}})=0$ and $\lambda(x, \boldsymbol{\mathbf{w}})=\pi^{*}(\boldsymbol{\mathbf{w}})\varphi(x)$, where $\varphi(x)$ is an arbitrary density function. Because the marginal effect varies only with $\boldsymbol{\mathbf{W}}$ but not with $X$, (ref) simplifies to

If $\pi=\pi^{*}$, we obtain a generalization of the angrist98 result for binary treatments: the weights are proportional to the conditional variance of $X$, $\operatorname*{Var}(X\mid \boldsymbol{\mathbf{W}})$. But if $\pi\neq \pi^{*}$, the weight function may be negative for some values of $\boldsymbol{\mathbf{W}}$. Consider, for instance, a panel data scenario where $\Gamma$ is linear, and $\boldsymbol{\mathbf{W}}$ consists of unit and time fixed effects. Then the assumption $\gamma_{0}\in\Gamma$ amounts to a parallel trends assumption: in the absence of the treatment, the average differences in outcomes for different units are constant and do not depend on the time period. If the treatment effects are heterogeneous, so that $g'(\boldsymbol{\mathbf{W}})$ depends on the unit and time period, then this two-way fixed effects regression still estimates a weighted average of marginal effects, but with weights that are negative if $E[X^{2}\mid \boldsymbol{\mathbf{W}}]<\pi(\boldsymbol{\mathbf{W}})\pi^{*}(\boldsymbol{\mathbf{W}})$. In the context of a binary treatment and two-way fixed effects regressions, this result has been noted in dCDH20 and gb21. (ref) generalizes this to an arbitrary treatment distribution and a general regression specification.

Another example where condition (b) of (ref) holds is when $X$ is bounded below by some baseline value (say, 0), $X\geq 0$, and the baseline outcome model is correctly specified, $g(0,\boldsymbol{\mathbf{W}})\in\Gamma$. Then condition (b) holds with $\gamma_{0}(\boldsymbol{\mathbf{W}})=g(0,\boldsymbol{\mathbf{W}})$, $\gamma_{1}=0$, and $\lambda(x, \boldsymbol{\mathbf{W}})=P(X\geq x\mid \boldsymbol{\mathbf{W}})$. In this case, the weights simplify to $\omega(x, \boldsymbol{\mathbf{W}})=E[\1{X\geq x}(X-\pi(\boldsymbol{\mathbf{W}}))\mid \boldsymbol{\mathbf{W}}]/\operatorname*{Var}(X-\pi(\boldsymbol{\mathbf{W}}))$.

Thus, if the marginal effects are constant, the partially linear model is doubly robust: the regression estimand is consistent for this constant treatment effect so long as either $\pi^{*}\in\Gamma$ or $g(0,\boldsymbol{\mathbf{W}})\in\Gamma$, as noted, for instance, in RoMaNe92. But this double robustness doesn't fully extend to the case with heterogeneous marginal effects. If the researcher gets the model for $X$ right, in the sense that $\pi=\pi^{*}$, then the partially linear regression estimates an average marginal effect. However, if the researcher gets it wrong, and only gets the outcome model under no treatment right, so that only condition (b) of (ref) holds, then weights on some of the true marginal effects $g'(X, \boldsymbol{\mathbf{W}})$ may be negative, risking a sign-reversal. This asymmetry has been noted by GoldsmithPinkham2024 for the case with a binary treatment $X$.\footnote{GoldsmithPinkham2024 also show that in regressions on multiple mutually exclusive treatment indicators, the regression estimand on a given treatment contains an additional contamination bias term corresponding to non-convex average effects of the other treatments.} (ref) shows that the result is general. The upshot of this asymmetry is that in cases where the treatment variable $X$ is only conditionally exogenous---whether the data is cross-sectional, panel, or time series---it pays off to conduct a sensitivity analysis with respect to the functional form of the control specification, in addition to the standard sensitivity analysis with respect to the set of controls.

Conclusion

We have shown that conventional linear methods for identifying causal effects in applied time series analysis based on observed shocks or proxies are robust to misspecification: they estimate a positively weighted average of the true nonlinear causal effects, irrespective of the extent of nonlinearities in the underlying \ac{DGP}. By contrast, identification approaches that exploit heteroskedasticity or non-Gaussianity of latent shocks are highly sensitive to violations of the assumed linear functional form of the structural model. Moreover, while linear identification via heteroskedasticity provides some testable restrictions, identification via non-Gaussianity is generally unfalsifiable despite the potential for severe biases.

Our results suggest that it is worthwhile for applied researchers who are interested in average marginal effects to expend the effort involved in constructing direct measures of shocks, or at least proxies that are credibly (approximately) monotonically related to the latent shock of interest. Shock measurement via narrative approaches or detailed institutional knowledge is admittedly highly work-intensive and in many applications may be practically impossible. Nevertheless, when the strategy is feasible, it affords an insurance against functional form misspecification which is not matched by the other identification approaches that we analyze. This is similar to the robustness of linear regression for estimating average treatment effects in randomized control trials relative to estimators that rely on parametric adjustment for nonlinearities or selection in observational data. Our results do not directly speak to other identification approaches like identification via long-run or sign restrictions, and we leave these for future work.\footnote{As a referee pointed out to us, it may be possible to analyze long-run restrictions by extending (ref) to apply to the instrumental variable specification in Shapiro1988. However, since the derivations in the latter paper rely on linearity of the structural function, it is not immediately clear what assumptions would guarantee convex weights.}

When reporting impulse responses from linear specifications with observed shocks, we recommend that researchers routinely report the implicit weight function, which is easy to compute via standard regression software. If the weight function associated with conventional local projections or \acp{VAR} is deemed to be unattractive, one can target other causal summaries as discussed in (ref), though further analysis is required on the best practices for doing so in macroeconomic applications.

More broadly, we hope that our paper will boost the agenda spearheaded by White2006, White2009, Angrist2011, Angrist2018, Rambachan2021, and goncalves24nonparametric,Goncalves2024 that seeks to draw lessons for macroeconometrics from the microeconometric treatment effect literature. The nonparametric framework used in the treatment effect literature contains useful lessons for empirical work in macroeconomics, even though conventional approaches to nonparametric estimation or debiased machine learning methods are impractical due to the much smaller data sets typical in macroeconomics.

\phantomsection \addcontentsline{toc}{section}{References}