Nonparametric Time Varying IV-SVARs: Estimation and Inference

Finance and Economics Discussion Series Federal Reserve Board, Washington, D.C. ISSN 1936-2854 (Print) ISSN 2767-3898 (Online) Nonparametric Time Varying IV-SVARs: Estimation and Inference Robin Braun, George Kapetanios, Massimiliano Marcellino 2025-004 Please cite this paper as: Braun, Robin, George Kapetanios, and Massimiliano Marcellino (2025). “Nonparametric Time Varying IV-SVARs: Estimation and Inference,” Finance and Economics Discussion Series 2025-004. Washington: Board of Governors of the Federal Reserve System, https://doi.org/10.17016/FEDS.2025.004. NOTE: Staff working papers in the Finance and Economics Discussion Series (FEDS) are preliminary materials circulated to stimulate discussion and critical comment. The analysis and conclusions set forth are those of the authors and do not indicate concurrence by other members of the research staff or the Board of Governors. References in publications to the Finance and Economics Discussion Series (other than acknowledgement) should be cleared with the author(s) to protect the tentative character of these papers.

Nonparametric Time Varying IV-SVARs: Estimation and Inference ∗ Robin Braun George Kapetanios Federal Reserve Board King’s College, London robin.a.braun@frb.gov george.kapetanios@kcl.ac.uk Massimiliano Marcellino Bocconi University, IGIER and CEPR massimiliano.marcellino@unibocconi.it December 9, 2024 Abstract Thispaperstudiestheestimationandinferenceoftime-varyingimpulseresponse functions in structural vector autoregressions (SVARs) identified with external instruments. Building on kernel estimators that allow for nonparametric time variation, we derive the asymptotic distributions of the relevant quantities. Our estimators are simple and computationally trivial and allow for potentially weak instruments. Simulations suggest satisfactory empirical coverage even in relatively small samples as long as the underlying parameter instabilities are sufficiently smooth. We illustrate the methods by studying the time-varying effects of global oil supply news shocks on US industrial production. JEL classification: C14, C32, C53, C55 Keywords: Time-varyingparameters, Nonparametricestimation, StructuralVAR, External instruments, Weak instruments, Oil supply news shocks, Impulse response analysis ∗WethankseminarparticipantsattheEuropeanCentralBank,BocconiUniversity,BankofEngland, King’s College, the Friendly Faces workshop, University of Konstanz, University of East Anglia, the PhiladelphiaFed,theIAAE2022,ESEM2022,theCFE,andaconferenceatHarvardUniversityforuseful comments. WealsothankJonasArias,RalfBrüggemann,AmbrogioCesa-Bianchi,ThorstenDrautzburg, Sophocles Mavroeidis, Aaron Metheny, Silvia Miranda-Agrippino, Pascal Paul, Mikkel Plagborg-Moller, Jim Stock, Mark Watson and Yiru Wang for insightful comments. Marcellino thanks Ministero dell Istruzione, Universita e Ricerca - PRIN Bando 2017 prot. 2017TTA7TYC for financial support. The views expressed in this paper are those of the authors and do not necessarily reflect those of the Board of Governors of the Federal Reserve System. 1

1 Introduction Instrumental variable (IV) identification of structural vector autoregressions (IV-SVARs) has become increasingly popular to study dynamic causal effects in empirical macroeconomics, including those of monetary policy (Gertler and Karadi, 2015; Caldara and Herbst, 2019; Jarociński and Karadi, 2020), oil price shocks (Känzig, 2021) or technology shocks (Miranda-Agrippino et al., 2024) among many. At the same time, important refinements to the methodology have been developed, building on the pioneering contributions of Stock (2008), Stock and Watson (2018) and Mertens and Ravn (2013). This includes the conduct of Bayesian inference (Arias et al., 2021; Giacomini et al., 2022), the establishment of the connection to local projections and the robustness to noninvertibility (Plagborg-Møller and Wolf, 2021, 2022; Forni et al., 2023), and, importantly, robust inference under weak identification (Montiel-Olea et al., 2021). Furthermore, Paul (2020) introduced the possibility of allowing for time-varying parameters in a Bayesian setting, while Inoue et al. (2024a,b) leverage the path estimator of Müller and Petalas (2010) for a similar purpose. In this paper we contribute to the literature of VARs identified by external instruments, developing estimators for IV-SVARs with slowly changing parameters aimed at capturing instabilities salient in macroeconomic relationships (Stock and Watson, 1996). While we take no stance on what causes the parameter changes, often discussed factors include institutional modifications, technological developments, economic trends such as globalization, or an evolving policy toolkit. Our paper complements and extends previous work on time-varying IV-SVARs in various ways. First, instead of assuming a Gaussian process for the model coefficients, we take a nonparametric approach that relies on persistence and smoothness assumptions on the pattern of parameter evolution. Formally, we build on classical kernel-based estimators introducedby Giraitisetal. (2014) andadapted forIVestimationin Giraitiset al.(2021). 2

Besides its nonparametric nature, our frequentist inference procedure for the model parameters and structural impulse response functions (IRFs) is computationally trivial and scales easily to larger dimensions and sample sizes. Furthermore, unlike the Bayesian alternative, inference can be robustified to account for potentially weak identification and easily handles very persistent time-series. Second, we provide results for two estimators that cater to different needs of the researchers, namely, (1) the classical IV-SVAR and (2) the internal instrument VAR estimator proposed by Plagborg-Møller and Wolf (2021). Under shock invertibility of the model, the IV-SVAR may be a powerful device, as it allows the researcher to back out the structural shocks up to a known constant. Hence, it is possible to construct time-varying IRFs that remain comparable over time in response to shocks of constant scale. Shock invertibility can be tested. If rejected, it is still possible to rely on the internal instrument estimator, which allows to estimate relative IRFs consistently in the absence of invertibility. However, given that the shock scale remains unknown in the internal instrument VAR, one cannot set a comparable shock size across time without further assumptions on the relationship between the shock of interest and the external instrument (see Paul (2020)). Our main object of interest are IRFs. In order to conduct inference for the relevant quantities, we proceed in two steps. First, we derive the asymptotic theory for the correspondingreduced-formparametersthatcharacterizethejointdynamicsoftheendogenous time series and the external instrument. We then either rely on an application of the delta method or follow Montiel-Olea et al. (2021) in constructing confidence sets via an inversion of the Anderson Rubin test statistic. The latter has the advantage of providing confidence-set robustness to a situation in which the instrument is only weakly correlated with the shock of interest, see, for example, Staiger and Stock (1997). This feature can be particularly important when a smaller bandwidth of the kernel estimator lowers the effec- 3

tive sample size considerably, even in larger sample sizes. Both methods are accompanied by closed-form solutions that allow for computationally efficient implementation. In order to understand the finite sample properties of the proposed method, we include a Monte Carlo exercise. Here, we calibrate a data-generating processes based on timevarying estimates of the global oil market VAR by Kilian (2009). We are able to obtain satisfactory empirical coverage if the evolution of parameters is sufficiently smooth, particularly for the weak-IV robust confidence sets. When selecting the bandwidth with a data-driven method that targets out-of-sample model fit, we document only a small deterioration in empirical coverage. We illustrate the methodology revisiting estimates of the transmission of oil supply news shocks on US industrial production (IP). Building on Känzig (2021), we study a timevarying IV-SVAR that includes monthly macroeconomic variables for the global crude-oil market,aswellasUSmining-andmanufacturingIP.Foridentification,themodelrelieson an external instrument that leverages futures price movements around OPEC production quota announcements. A constant parameter model suggests that these shocks transmit as a cost-push shock to the US economy, where manufacturing production declines with increasing real oil prices. However, our methodology reveals strong time-variation in the estimated impulse response functions, which seems to align with the shale-oil revolution. US mining output, which includes extraction of oil- and gas, reacts more strongly and quickly nowadays than in the past. Furthermore, US manufacturing output no longer declines, challenging that the oil-market specific shock still transmits as cost-push shock to the US. Our finding complement recent evidence presented in Bjørnland and Skretting (2024) on the time-varying impact of oil-price shocks. However, in contrast to that paper, we identify an oil-market specific shock by instrumental variables instead of exclusion restrictions, and use the proposed kernel based methods instead of Bayesian techniques. 4

Related Literature Our paper builds on the seminal work of Cogley and Sargent (2005) and Primiceri (2005) who introduced TPV into VARs by letting coefficients evolve according to a random walk in a Bayesian setting. Paul (2020) extends their framework to achieve identification of IRFs by external instruments, including the instrument as a regressor. Inoue et al. (2024a,b) also estimate time-varying IRFs identified by IV, but rely on the frequentist path estimator of Müller and Petalas (2010). While this allows for less prior dependence, the underlying implementation still relies on a Gaussian random walk assumption to obtain point estimates and standard errors. Our methodological approach is distinct from these papers along various dimensions. First, we rely on an entirely frequentist, nonparametric approach that leverages kernel estimators. This approach has several practical advantages. First, it avoids a parametric choice for the law of motion underlying the SVAR parameters and instead relies on nonparametric smoothness conditions. Second, it is computationally simple and can handleverylargedatasetsaswellaspersistenttimeseries. Insuchcircumstances, existing methods may struggle as they require to loop trough each observation and potentially need to deal with non-stationary draws implying explosive IRFs. Third, similar to a constant parameter VAR, our framework provides very simple formulations for standard errors. We leverage those to conduct robust inference valid under weak instruments (see also Inoue et al. (2024a)). Finally, unlike previous papers, we cover both the standard IV-SVAR model as well as the internal IV estimator (Plagborg-Møller and Wolf, 2021).1 The theory behind the TVP kernel estimators in IV-SVARs largely builds on earlier work of Giraitis et al. (2014), Giraitis et al. (2018) and Giraitis et al. (2021). However, our paper provides additional results that are required to accommodate identification 1The estimators by Paul (2020) and Müller and Petalas (2010) are compared with our proposal in Appendix E, where we find that with simulated data they perform similarly, while with actual data our method seems better capable of handling very persistent time series. 5

via external instruments, including the asymptotic distribution of the covariance matrix estimator, the construction of confidence sets, and the joint distribution of neighboring estimators. The joint distribution allows for inference of IRFs that are comparable over time when studying relative IRFs in the internal instrument VAR. We note that we are not the first to leverage kernel-based estimators to introduce TVP into VARs, see e.g. Kapetanios et al. (2019) and Hipp (2020). However, unlike these papers we focus on external instrument identification. Finally, we would also like to relate our paper to Amir-Ahmadi et al. (2023) who also allow for a time-varying relationship between the instrument and the structural shock of interest. Unlike their paper, however, we allow the IRFs to be time-varying. Outline The paper is organized as follows. Section 2 develops the methodology for kernel-based inference in TVP IV-SVARs. Section 3 presents results of Monte Carlo simulations. Section 4 studies the transmission of oil supply news shocks on US industrial production, and section 5 recapitulates. Proofs and additional empirical results are gathered in the online supplementary material. This also includes a comparison of our methodology to a Bayesian- and a path estimator. 2 Methodology Inthissection,westartrevisitinginstrumentalvariableidentificationofImpulseResponse Functions in a constant-parameter SVARs. We then generalize the model towards timevarying coefficients, discuss normalization of the shock size across time, and inference of reduced form quantities via kernel based methods. Finally, we show how the results can be leveraged to compute confidence sets of IRFs. 6

2.1 Identification of VAR impulse response functions via external instruments Consider the n-variate SVAR(p) model given by: y = ν +A y +A y +...+A y +u , u ∼ (0,Σ) (1) t 1 t−1 2 t−2 p t−p t t u = Bε , ε ∼ (0,I), (2) t t t where y = (y ,...,y )′ is a n×1 vector of endogenous time series, ν is a n×1 vector t 1t nt of intercepts, A ,i = 1,...,p are n × n matrices of autoregressive coefficients and the i error terms u and ε are, for simplicity, assumed to be i.i.d. white noise with covariance t t matrix Σ and I respectively. Equation (1) describes the reduced form VAR dynamics t n of y as a function of lagged realizations and a vector of n × 1 error terms u with full t t covariance matrix Σ. Equation (2) relates the prediction errors u to n × 1 structural t shocks ε whose elements are orthogonal and standardized to unit variance. The n×n t matrix B is the contemporaneous impact matrix and reflects the immediate responses of the variables y to the structural shocks ε . For the moment, we assume that the t t model is stable, which implies that the SVAR(p) has a MA(∞) representation given by y = µ + (cid:80)∞ C (A)Bε = µ + (cid:80)∞ Θ ε , where µ = E(y ) and the n × n t y j=0 j t−j y j=0 j t−j y t coefficient matrices Θ = C (A)B, are the structural impulse response functions (IRFs). j j The reduced form MA(∞) matrices C (A) can be computed recursively from C (A) = j j (cid:80)j C (A)A with C (A) = I and A = 0 for i > p. i=1 j−1 i 0 n i The main focus of this paper is the computation of impulse responses to a single shock. Without loss of generality, let this shock be ordered first in the system (ε ) and call it 1,t the target shock. Corresponding IRFs are then given by picking elements in the MA(∞) matrices: ∂Y i,t+k = λ = e′C (A)Be , (3) ∂ε k,i i k 1 1,t where e denotes the ith column of the identity matrix I . Hence, equation (3) defines i n 7

the IRFs λ as the dynamic effect a unit standard deviation shock in ε on variable i, k,i 1t k periods ahead. It is important to note that, without further assumptions, IRFs are not identified. The reason is that the same reduced form dynamics of the VAR forecast errors u = Bε are t t ˜ obtained for any alternative structural model B = BQ where Q is an orthogonal rotation matrix ({Q : Q′Q = I ,Q′ = Q}). To see this, note that both models imply the same n reduced form covariance matrix Σ = BB′ = BQQ′B = B ˜ B ˜′. In this paper, we rely on u an identification strategy that involves an instrumental variable z for the target shock t (Stock and Watson, 2012; Mertens and Ravn, 2013). Assumption 1 (External Instrument). Let z be an instrument for the first shock. The t stochastic process {(ε ,z )}∞ satisfies: t t t=1 1. E(z ε ) = α ̸= 0, t 1,t 2. E(z ε ) = 0 for j ̸= 1. t j,t Assumption 1 allows to identify b = Be up to scale and sign normalization, since: 1 1     E[ε z ] (cid:20) (cid:21) α 1,t t Γ = E(u t z t ) =B    = b 1 b 2     = αb 1 , E[ ε z ] 0 2:n,t t where ε = [ε ,...,ε ]′. In words, the correlation between the external instrument 2:n,t 2t nt and reduced form prediction error is proportional to the first column of the impact matrix b . 1 In order to obtain interpretable magnitudes, there are two popular approaches to normalize IRFs in IV-SVARs. The first is known as the unit shock standardization (Stock and Watson, 2016) or relative IRFs. Here, the shock variance is re-normalized to yield IRFs ˜ that increase the first variable by unity on impact (say b = 1). In that case, it holds 11 that Γ = E(z u ) = α, implying that ˜ b = Γ/e′Γ is the first column of the rescaled 11 t 1,t 1 1 impact matrix measuring the response to a target shock with unidentified standard deviation Var(ε˜ ) = b2 . Corresponding IRFs as function of reduced form parameter are 1t 11 8

then given by: λ ˜ = e′C (A)Γ/e′Γ. (4) k,i i k 1 The second approach is to normalize the standard deviation of the shock to unit variance (Var(ε ) = 1), yielding absolute IRFs. Here, one is required to incorporate additional 1,t informationofthereducedformcovariancematrixΣtorecoverαandhenceb . Exploiting 1 invertibility of the model Σ = BB′ yields the following quadratic form: Γ′Σ−1Γ = (αBe )′(BB′−1(αBe ) = α2. 1 1 √ Normalizing α > 0, one can back out b = Γ/α = Γ/ Γ′Σ−1Γ and define the absolute 1 IRFs as the following function of reduced form parameters: √ λ = e′C (A)Γ/ Γ′Σ−1Γ. (5) k,i i k At this point, it is worth discussing a key difference between the two definitions of im- ˜ pulse response functions: λ does not rely on invertibility of the model, which is the k,i assumption that structural shocks can be recovered as a function of the VAR prediction errors ε = B−1u . As shown in Plagborg-Møller and Wolf (2021), augmenting the VAR t t with the external instrument z allows for consistent estimation of relative impulse ret sponse functions λ ˜ , even if invertibility does not hold. Specifically, for y˜ = [z ,y′]′, the k,i t t t resulting internal instrument VAR model reads: ˜ ˜ ˜ ˜ y˜ = A y˜ +A y˜ +...+A y˜ +u˜ , u˜ ∼ (0,Σ), (6) t 1 t−1 2 t−2 p t−p t t (cid:16) (cid:17) and robust relative IRFs are obtained by λ ˜ = e′ C A ˜ Γ ˜ /(e′Γ ˜ ) for Γ ˜ = e′chol(Σ ˜ ). k,i 1+i k 2 t 1 On the other hand, without invertibility it is no longer possible to identify absolute IRFs (λ ).2 As we will discuss in the next subsection, the ability to recover the shock up to a k,i known constant will be an important advantage when it comes to studying time-varying 2See Plagborg-Møller and Wolf (2022) for detailed analysis on how the shock can be set-identified, however. 9

impulse response functions that remain comparable over time. 2.2 Introducing time-varying coefficients Introducing time-varying coefficients into the SVAR reads: y = A y +A y +...+A y +B ε , ε ∼ (0,I ) (7) t 1t t−1 2t t−2 pt t−p t t t n where E (ε ) = 0, E (ε ε′) = I and E (u u′) = Σ = B B′. Also, let Γ = E (z u ) and t t t t t n t t t t t t t t t t update the IV assumption to the time-varying case: Assumption 2 (External Instrument). Let z be an instrument for the first shock. The t stochastic process {(ε ,z )}∞ satisfies: t t t=1 1. E (z ε ) = α ̸= 0, t t 1,t t 2. E (z ε ) = 0 for j ̸= 1. t t j,t At this point, one approach would be to impose a specific parametric assumption about how time variation is generated, e.g. via a random walk, allowing for likelihood based inferenceusingtheKalmanfilter(Primiceri,2005;Paul,2020)orquasilikelihoodmethods as Müller and Petalas (2010). Instead, in this paper we follow a nonparametric approach along the lines of Giraitis et al. (2014, 2018), which assumes a bound on the degree of time variation that can be allowed for in order to conduct valid asymptotic inference via kernel-based estimators: Assumption 3. Let β = vec(A ) for A = [A ,...,A ], σ = vech(Σ ), and θ = t t t 1t pt t t t [β′,Γ′,σ′]′. Then: t t t (cid:16)s(cid:17) sup ||θ −θ ||2 = O ,||θ || < ∞, for all t. t t+j t T j≤s Assumption 3 states that the model parameters are bounded and that changes to those parameters are restricted to be small. The rate is assumed to be of the order T−1 but in previous work (see, e.g. Giraitis et al. (2018)), a relaxation to an order given by T−γ, 0 < γ ≤ 2, has been shown to be feasible. Such an order is equivalent to a mild Lipschitz condition on the smoothness of the parameters and is much milder then existing 10

conditions in the time-varying literature. Note that, unlike most other existing work, it is not assumed that parameters are smooth deterministic functions of time but, instead, we place a restriction on their differences. For simplicity we assume that parameters are a sequence of deterministic constants, though allowing for smooth stochastic processes is also feasible. The theory in all existing work (such as, e.g., Giraitis et al. (2018)) has been developed with a single rate of change for all parameter processes. As a result, we alignwiththissetting. Extendingtheanalysistodifferentratesfordifferentparametersis possible but is outside the scope of the paper. Assumption 3 enables consistency and rate results as presented in Theorem 1-3. The results are similar in nature to those presented in, e.g., Giraitis et al. (2018). The assumption requires that the deviation in parameters is shrinking with T, so the idea behind consistency is close to an infill asymptotic setup. Under Assumption 3, Giraitis et al. (2018) show that the MA(∞) representation can be expressed as: ∞ (cid:88) y = C (A )B ε +o(1). (8) t k t t t−k k=0 Equation (8) states that, under assumption 3, the MA(∞) representation of the TVP- SVAR is asymptotically given by that of a fixed-coefficient model, but replacing A and B with their time-varying counterparts. Under the instrumental variables assumption 2, time-varying IRFs to a shock of size one standard deviation are given by:3 (cid:113) λ = e′C (A )b = e′C (A )Γ / Γ′Σ−1Γ , (9) k,i,t i k t 1t i k t t t t t where b = B e is the first column of B . For the unit shock normalization, the corre- 1t t 1 t sponding time-varying IRFs are: λ ˜ = e′C (A )Γ /e′Γ , (10) k,i,t i k t t 1 t effectively measuring IRFs to a re-normalized target shock ε˜ with variance b2 . 1,t 11,t 3Alternatively, one might pursue a simulation based approach to obtain a more accurate picture as advocated in Koop et al. (1996), which is based on the exact MA(∞) representation. 11

In a time-varying parameter VAR, the choice of IRF normalization is not just a matter of preference. To unpack the differences, consider a toy model of supply and demand for quantities and prices y = [q ,p ]′: t t t   supply: q = α p +σ ε εs t t t 1t 1t t   ∼ (0,I ),   2 demand: q = β p +σ ε εd t t t 2t 2t t where α and β are supply- and demand elasticities, respectively, while [σ ,σ ] are the t t 1t 2t   1 −α t volatilities of the structural shocks. In this case, A =  , and t   1 −β t   −σ β /(α −β ) σ α /(α −β ) 1t t t t 2t t t t B = A−1diag([σ ,σ ]′) =  . t t 1t 2t   −σ 1/(α −β ) σ /(α −β ) 1t t t 2t t t Assume the availability of an instrument for the supply shocks. Then, the absolute IRFs λ leverage estimates of Γ and Σ to identify b , the first column of B . This IRF is k,i,t t t 1t t a function of both, time-varying elasticities and volatilities facilitating interpretation by holding constant the scale of the shock to unity throughout time. For example, in the case of a monetary policy shock, this means that estimates of λ will not only reflect k,i,t variation in the elasticities but also shock volatility, which summarizes how successful a central bank is in steering interest rates. This may be time-varying for many reasons, e.g. the design of new policies or temporary policy constraints such as the zero lower bound. (cid:18) (cid:19) Relative IRFs, on the other hand, yield Γ t /e′ 1 Γ t = 1 1/β t in our toy-model, and hence only depend on the structural elasticities. In this case, the IRF corresponds to a shock with a time-varying scale ε˜ ∼ (0,σ2 (β /(α −β ))2). Depending on the applica- 1t 1t t t t tion at hand, this might still be a very useful quantity to study. It is also worth noting that unlike absolute IRFs, relative IRF does not depend on the shock volatilities. Hence, smoothness as stated in Assumption 3 is not necessary for σ , and our inference procet dures for relative IRFs could be adjusted to allow for other forms of heteroskedasticity in 12

structural shocks, e.g. as in Gonçalves and Kilian (2004), but we do not consider this in the current paper. Following the discussion in section 2.1, the computation of λ relies on shock invertibilk,i,t ityandhence,maynotalwaysbeanoptionfortheresearcher. However,wheninvertibility is a concern, it is still possible to obtain relative IRFs to a fix shock size under stronger assumption about the relationship between the external instrument and the target shock (Paul, 2020). Specifically, assuming that E(z ε ) = α is constant as in Assumption 1, all t 1t the time-variation observed in Γ = αb can be attributed to differences in the scale 1,t 11,t of the shock (b ). In that case, normalizing the IRFs to increase the first variable by 11,t unity at a fixed time point t is sufficient to obtain responses to a constant shock size b over time: λ ˜ = e′C (A )Γ /e′Γ . (11) k,i,t,t b i k t t 1 t b Specifically, Γ /e′Γ = αb /(αb ) = b /b measures the impact effect of the target t 1 t b 1t 11,t b 1t 11,t b shock normalized to have variance b2 . Hence, although the shock volatility is still 11,t b unidentified, it is constant throughout the sample and hence remains comparable across time. Note that this is generally not the case if α itself was subject to time-variation. t Summing up our discussion, once time-variation is introduced into the model, a tradeoff arises. Under invertibility, an IV-SVAR allows to identify the scale of the shock throughoutthesampleandstudyabsoluteIRFsofaconstantshocksize(λ ). Whenever k,i,t invertibility does not hold, unit shock IRFs based on a internal instrument VAR may be a useful alternative. While relative IRFs require less assumption on the smoothness of shock volatilities, they requires stronger assumptions on α to obtain (relative) IRFs t ˜ that remain comparable across time (λ ). In practice, we therefore recommend to k,i,t,t b pre-test for shock-invertibility. If there is no evidence against invertibility of the target shock in a given application, proceeding with the IV-SVAR estimator may be preferable as it requires minimal assumptions on α . However, if shock invertibility is rejected, t 13

informative results under the stronger assumption α = α may still be obtained based on t the internal instrument VAR. A pre-test for shock invertibility that can be readily applied in a time-varying set-up is described in Plagborg-Møller and Wolf (2022). Specifically, the testable prediction is that under invertibility of the shock, z should not Granger t cause y in an instrument augmented VAR. t 2.3 Joint inference for the reduced form parameters ˜ In order to conduct inference for λ and λ we proceed in two steps. We start k,i,t k,i,t,t b deriving the joint asymptotic distribution of kernel-based estimators of the reduced form parameters A , Γ , Σ and Σ , both for the IV-SVAR and internal instrument VAR. In t t t t b a second step, we construct confidence sets for the impulse responses either by the Delta method or an inversion of the Anderson and Rubin test statistic. StartingwiththeIV-SVARgiveninequation(7),letβ = vec(A )andx = [y′ ,y′ ,...,y′ ]. t t t t−1 t−2 t−p Then, the kernel estimator is given by: (cid:34) (cid:35)−1(cid:34) (cid:35) T T (cid:88) (cid:88) β ˆ = I ⊗ w (H)x x′ w (H)vec(x y′) , (12) t n t,j j j t,j j j j=1 j=1 T 1 (cid:88) ˆ Γ = w (H)uˆ z , (13) t t,j j j H j=1 T 1 (cid:88) Σ ˆ = w (H)uˆ uˆ′, (14) t H t,j j j j=1 whereuˆ = y −(I ⊗x′)β ˆ andw (H) = K(|t−j|/H)isaKernelfunctiontoensuremore j j n j t t,j distant observations get discounted when forming the estimate at time t. To establish theoretical properties of the estimator, we make the following two assumptions on the error term and kernel: Assumption 4. ε = (ε ,··· ,ε )′ is an iid process such that E[ε4 ] < ∞. z is a t 1t nt i1 t stationary, α-mixing process with exponentially declining mixing coefficients, such that E[z4] < ∞. Further, E[y4] < ∞ for i = 1,··· ,n. 1 i0 Assumption 5. K is a non-negative bounded function with a piecewise bounded deriva- ˙ (cid:82) tive K(x) such that K(x)dx = 1. If K has unbounded support, we assume in addition 14

that K(x) ≤ Cexp(−cx2), |K ˙ (x)| ≤ C(1+x2)−1, x ≥ 0, for some C > 0, c > 0. Concerning Assumption 4, two remarks are in order. First, as in previous work on timevarying regressions using kernel estimators, we assume the errors to be iid processes. Contemporaneous work by Giraitis et al. (2024a) shows that this can be relaxed to allow for martingale difference processes. However, a rigorous treatment is technically demanding and we do not think that, for the purposes of our analysis, much would be gained. Second, Assumption 4 precludes the presence of serial correlation in the error term, and the associated need for potential heteroskedasticity and autocorrelation robust (HAR) corrections. We note three things in this context. First, it is reasonable to expect that the presence of enough lags in the VAR model will soak up serial correlation. Second, one can test for remaining residual correlation, although recent work by Giraitis et al. (2024b) notes several problems with such tests, typically associated with over-rejection. Finally, we are not aware of any rigorous work on HAR procedures for time varying models. Such work would want to account for time variation in the correction and, while this is a very interesting topic of research, we consider it beyond the scope of the current paper. Concerning Assumption 5, one suitable choice that we adapt is the Gaussian kernel (cid:104) (cid:105) K (H) ∝ exp −1 (cid:0)j−t(cid:1)2 , further normalized such that (cid:80) w = H. In Appendix A j,t 2 H j t,j we show that: Theorem 1. [joint asymptotic normality of reduced form parameters in the TVP-IV- SVAR] Under Assumption 3-5 and H = o(T1) it holds that: 2  ˆ  β −β √ t t H Γ ˆ t −Γ t  → d N(0,V θt ), ˆ vech(Σ )−σ t t for V = S Π S′ and θt t ww,t t  I ⊗Π−1 0 0  n x,t (cid:0) (cid:1) S t = − I n ⊗Π xz,t Π− x, 1 t I 0 , 0 0 S σ 15

for Π = plim 1 (cid:80)T w x x′, Π = plim 1 (cid:80)T w z x , x,t T→∞H j=1 j,t j j xz,t T→∞H j=1 j,t j j Π = plim 1 (cid:80)T w2 ξ ξ′, ξ = [vec(x u )′,(z u −Γ)′,vec (cid:0) u′u −Σ (cid:1)′ ]′, and S ww,t T→∞H j=1 j,t j j j j j j j j j t σ such that vech(Σ ) = S vec(Σ ). t σ t It is worth mentioning that the bandwidth H assumed in Theorem 1 is strictly related to the amount of time-variation permitted in Assumption 3. If a different rate is assumed in Assumption 3, say T−γ rather than T−1, a different H should be used (the more time variation, the smaller the optimal H). While our theoretical results imply certain conditions for the bandwidth (H = o(T1)), it is reasonable to allow for larger bandwidths 2 in practice, for smaller sample sizes that are common in fields such as macroeconomics. From a practical point of view, since γ is not known, a cross validation approach can be used to select H. Towards the end of this section, we describe a simple procedure that targets out-of-sample model fit for IV identified impulse response functions. A second important comment is that, as mentioned, the way we model nonparametric time variation follows GKY and it is different from the more common approach of assuming that the parameters change as a function (with at least a bounded derivative) of t/T e.g. as in Dahlhaus (1997). Hence, the proof of the results and the conditions on the optimal bandwidth differ from the usual ones. Equivalent results can be obtained for the reduced form parameters of the time-varying internal instrument VAR. For y˜ = [z ,y′]′, the underlying model reads: t t t ˜ ˜ ˜ ˜ y˜ = A y˜ +A y˜ +...+A y˜ +u˜ , u˜ ∼ (0,Σ ), (15) t 1t t−1 2t t−2 pt t−p t t t (cid:16)(cid:104) (cid:105)(cid:17) ˜ ˜ ˜ ˜ ˜ where β = vec A ,...,A and P = chol(Σ ) is the Cholesky decomposition such t 1t pt t t that P ˜ P ˜′ = Σ ˜ . As discussed above, the main object of interest based on the internal t t t (cid:16) (cid:17) instrument VAR is λ ˜ = e′ C A ˜ P ˜ /(e′P ˜ ), that is the relative IRF stank,i,t,t b 1+i k t •1,t 2 •1,t b dardized to increase the first variable by unit on date t . We start with results for the b 16

following estimator of time t reduced form coefficients: (cid:34) (cid:35)−1(cid:34) (cid:35) T T β ˆ˜ = I ⊗ (cid:88) w (H)x˜ x˜′ (cid:88) w (H)vec(x˜ y˜′) (16) t n+1 t,j j j t,j j j j=1 j=1 T Σ ˆ˜ = H−1 (cid:88) w (H)uˆ˜ uˆ˜′, (17) t t,j j j j=1 where uˆ˜ = y˜ − (I ⊗ x′)β ˆ˜ . Joint asymptotic normality between the reduced form j j n+1 j t parameters are then given as follows: Theorem 2. [joint asymptotic normality of reduced form parameters in the TVP internal instrumentVAR.]UnderAssumption3-5andH = o(T1): defineΠ ˜ = plim 1 (cid:80)T w x˜ x˜′, 2 x,t T→∞H j=1 j,t j j Π ˜ = plim 1 (cid:80)T w2 x˜ x˜′,Π = plim 1 (cid:80)T w vec(u˜ u˜′)vec(u˜ u˜′)′ , ww,t T→∞H j=1 j,t j j uu,uu,t T→∞H j=1 j,t j j j j σ˜ = vech(Σ ˜ ) and L be the n(n + 1)/2 × n2 elimination matrix such that vech(A) = t t n L vec(A). Then, the estimators β ˆ˜ and σˆ˜ are jointly asymptotically normal and asympn t t totically independent of each other. Their respective distributions are given by √ (cid:16)ˆ˜ ˜ (cid:17) d (cid:18) ˜ (cid:16) ˜ (cid:17)−1 ˜ (cid:16) ˜ (cid:17)−1 (cid:19) H β −β → N 0,Σ ⊗ Π Π Π , t t t x,t ww,t x,t √ (cid:16) (cid:17) H σˆ˜ −σ˜ → d N (cid:0) 0,L Π L′ −σ˜ σ˜′ (cid:1) . t t n+1 uu,uu,t n+1 t t Under an additional normality assumption for the errors, the asymptotic variance of (cid:16) (cid:17) σˆ˜ further reduces to 2D+ Σ ˜ ⊗Π ˜ D+′ , where Π ˜ = plim 1 (cid:80)T w2 u˜ u˜′ t n+1 t uu,t n+1 uu,t T→∞H j=1 j,t j j and D+ = (D′ D )−1D′ for D the duplication matrix such that vec(Σ ˜ ) = n+1 n+1 n+1 n+1 n+1 t ˜ D vech(Σ ). n+1 t ˜ Given that estimates of λ are based on reduced form parameters at time t and t , k,i,t,t b b the construction of corresponding confidence sets requires an expression for their joint distribution. This is particularly relevant when t and t are close, and estimates are b highly correlated by construction. Given that VAR slope and covariance parameters are asymptotically uncorrelated, and given that only covariance estimates of t are used b to construct λ ˜ , it’s sufficient to focus on the joint distribution of σˆ˜ and σˆ˜ . The k,i,t,t t t b b following Corollary gives their asymptotic joint distribution. Corollary 3. LetAssumptions3-5holdandH = o(T1). Definew (H) = [w (H),...,w (H)]′ 2 t t,1 t,T 17

and σ˜ = vech(Σ ˜ ). Let Π = plim 1 (cid:80)T vec(ξ ˜ ξ ˜′ )vec(ξ ˜ ξ ˜′ )′ for t,t b t,t b uu,uu,t,t b T→∞H j=1 wj wj wj wj (cid:104) (cid:16) (cid:17) (cid:16) (cid:17)(cid:105) ξ ˜ = w1/2(H) y˜ −x˜ Θ ˜ ,w1/2(H) y˜ −x˜ Θ ˜ . wj t,j j j t t b ,j j j t b Under these definitions, it follows that: √ (cid:16) (cid:17) H σˆ˜ −σ˜ → d N (cid:0) 0,L Π L′ −σ˜ σ˜′ (cid:1) t,t b t,t b 2(n+1) uu,uu,t,t b 2(n+1) t,t b t,t b As above, under normality assumption of the errors, the asymptotic variance of σˆ˜ further t (cid:16) (cid:17) simplifies to 2D+ Σ ˜ ⊗Π ˜ D+′ , where Π ˜ = plim 1 (cid:80)T U ˜ U ˜′ for 2(n+1) t,t b uu,t,t b 2(n+1) uu,t,t b T→∞H j=1 2 2 (cid:104) (cid:105) ˜ ˜ ˜ ˜ ˜ ˜ ˜ U = w (H)⊗(Y −XΘ ),w (H)⊗(Y −XΘ ) 2 t t t t b b and D+ = (D′ D )−1D′ for D the duplication matrix such that 2(n+1) 2(n+1) 2(n+1) 2(n+1) 2(n+1) ˜ ˜ vec(Σ ) = D vech(Σ ). t,t b 2(n+1) t,t b 2.4 Inference for impulse response functions Based on asymptotic results for the estimators of the reduced form parameters, we can rely on standard methods to construct confidence sets for the object of interest, that are the estimates of time-varying structural impulse response functions: (cid:113) λ ˆ = e′C (A ˆ )Γ ˆ / Γ ˆ′Σ ˆ−1Γ ˆ , k,i,t i k t t t t t λ ˆ˜ = e′ C (cid:16) A ˆ˜ (cid:17) P ˆ˜ /(e′P ˆ˜ ), k,i,t,t b 1+i k t •1,t 2 •1,t b ˆ ˆ ˆ ˆ˜ ˆ˜ ˆ˜ where A , Γ and Σ are based on the IV-SVAR, while A , P and P are based on t t t •1,t •1,t b the internal instrument VAR. In this paper, we discuss two approaches to construct appropriate confidence sets, either via the classical Delta method or an inversion of the Anderson Rubin (AR) test statistic as in Montiel-Olea et al. (2021). The latter fixes α under the null hypothesis and hence remains valid even under asymptotically weak instruments, that is if α → 0. √ (cid:16) (cid:17) ˆ d StartingwiththeDeltaMethod,itsapplicationyieldsthat H λ −λ → N (0,Ω ), k,t k,t k,t where Ω = J (β ,Γ ,σ )V J (β ,Γ ,σ )′. Here, V denotes the joint distribution k,t k t t t θt k t t t θt 18

of the IV-SVAR reduced form parameters which we give in Theorem 1. Furthermore, J (β ,Γ ,σ ) denotes the derivative of λ with respect to the reduced form parameters. k t t t k,t Similarly, for relative IRFs it is √ (cid:16)ˆ˜ ˜ (cid:17) d (cid:16) ˜ (cid:17) H λ −λ → N 0,Ω , k,t,t k,t,t k,t,t b b b (cid:16) (cid:17) (cid:16) (cid:17)′ ˜ ˜ ˜ ˜ ˜ ˜ ˜ where Ω = J β ,σ˜ V J β ,σ˜ for V being elements of the joint asympk,t,t b k t t,t b θt,tb k t t,t b θt,tb ˜ totic covariance matrix stated in Theorem 2 and Corollary 3, and J () is the gradient k ˜ of Ω with respect to the reduced form parameters. Analytical formulas for both k,t,t b gradients are derived in Appendix B. As documented in Montiel-Olea et al. (2021) for relative IRFs, empirical coverage rates of the Delta method can deteriorate quickly when the instrument is only weakly correlated with the target shock. This may be particularly true in TVP models where the effective sample size is fairly small. Hence, we also cover weak identification robust confidence ˜ sets. Consider the (n+1)×1 vectors L and the (n+2)×1 vector L : k,t     (cid:16) (cid:17) ˜ ˜ C (A )Γ C A P L =  k t t , L ˜ =  k t •1,t , k,t  (cid:112)  k,t,t b   Γ′Σ−1Γ e′P ˜ t t t 2 •1,t b for which it holds that λ = (e′L )/(e′ L ) and λ ˜ = (e′ L ˜ )/(e′ L ˜ ). k,i,t i k,t n+1 k,t k,t,t b 1+i k,t,t b n+2 k,t,t b To derive the AR confidence set, first note that an application of the Delta Method √ (cid:16) (cid:17) implies that H L ˆ −L → d N(0,ΩL ) where ΩL depends on the covariance matrix k,t k,t k,t k,t given in Theorem 1 and the gradient of L with respect to the reduced form parameters. k,t Similarly, a statement can be obtained for √ H (cid:16) L ˜˜ −L ˜ (cid:17) → d N(0,ΩL˜ ) based on k,t,t b k,t,t b k,t,t b the reduced form results for the internal VAR estimator. Without loss of generality, let us focus on the confidence set for λ . The null hypothesis k,i,t λ = λ implies e′L − λ e′ L = 0, a linear restriction on L (see also Fieller k,i,t 0 i k,t 0 n+1 k,t k,t (1944)). Following Montiel-Olea et al. (2021), a Wald Test statistic can be set up as q(λ ) = H(e′ i Lˆ k,t −λ0e′ n+1 Lˆ k,t )2 where ωˆ is the ijth element of Ω ˆL . Further inversion 0 ωˆii−2λ0ωˆi,n+1+λ2 0 ωˆn+1,n+1 ij k,t 19

yields the AR confidence set of coverage 1 − a, given by CSAR{λ |q(λ ) ≤ χ2 }. k,i,t k,i,t 1,1−a The inequality q(λ ) ≤ χ2 is quadratic in λ and can be solved in closed form. k,i,t 1,1−a k,i,t For details, including the gradients necessary to obtain ΩL , we refer to Appendix B. k,t A few properties are worth mentioning at this point. First, even in a weak instrument √ casewhereα = a/ H forsomefixeda, theARCSremainsvalid. Thereasonisthatthe H Wald statistic fixes λ under the null hypothesis and hence does not require consistent k,i,t estimates thereof for its validity. Second, in the strong instrument case, Montiel-Olea et al. (2021) prove that the AR confidence set converges to Delta Method implied confidence intervals. For those reasons, we generally recommend to use the weak-IV robust confidence intervals. ˜ Third, we note that for relative IRFs λ , both the Delta method and AR CS depend k,t,t b on the choice of t . Montiel-Olea et al. (2021) show that the 100%(1−a) AR confidence b set is finite only if the Wald test statistic for e′P ˜ is above its corresponding critical 2 •1,t value (χ2 ). For this reason, we recommend setting t to a point in the sample where 1,1−a b the Wald-test statistic is very high, which can be done fairly automatic way. Finally, we note that this paper covers the simplest case of a single instrument identifying one target shock. However, the reduced form results hold for a more general case of k instruments. Hence, it is possible to extend the results in this subsection to the case of r targetshocksandr instruments,whichrequiresadditionalrestrictions,e.g.asemployedin Mertens and Montiel Olea (2018). We point the interested reader to the supplementary appendix of Montiel-Olea et al. (2021), discussing how robust confidence sets can be constructed in this case. 2.5 Bandwidth selection To control for the amount of time-variation in any given application, the user is required to set a bandwidth H prior to estimation. We acknowledge that there are several ways this can be done. First, similar to a Bayesian approach, one might take on off-model 20

information about how much time-variation is reasonable to see in Impulse Response Functions over a certain time span, and select a bandwidth accordingly. In this paper, however, we pursue a purely data-driven approach selecting a bandwidth H that provides the best out-of-sample model fit. Acknowledging that IRFs are conditional forecasts, we propose to evaluate the models out-of-sample predictive performance comparingdatarealizationsy toconditionalforecastsyˆ (H) = E [y |y˜ ,...,y˜,z ,...,z ]. t+h t+h t t+h 1 t t+1 t+h ˜ Just like the IRFs identified by IV, this is a function of both the VAR slope parameters A t and covariance matrix Σ ˜ in the instrument-augmented VAR (Waggoner and Zha, 2003).4 t Based on those forecasts, we propose the following objective function to choose the bandwidth: (cid:88) n T (cid:88) −h (cid:88) hm min w (y −yˆ (H))2, i i,t+h i,t+h H i=1 t=Ts h=1 where T is the time in the sample where the pseudo out-of-sample exercise starts, h s m is the maximum forecast horizon to be included in the evaluation, and w is a variable i specific weight to account for differences in scale. For the latter, we simply use the inverse variance of AR(1) residuals in each variable. 3 Monte Carlo Simulations Inthefollowing, westudythefinitesamplepropertiesoftheproposedinferenceprocedure for time-varying impulse response function. As we expect, the performance of confidence sets will depend on the effective sample size, the speed of time-variation underlying the SVAR coefficients and the instrument strength. Overall, our findings suggest that the asymptotic theory provides a reasonable approximation in finite samples. 4Note that to evaluate the conditional expectations based on the information set up to time t, we relyonmodifiedestimatorsA ˆ˜ (H)andΣ ˆ˜ (H)definedasinequation(16)-(17),butbasedontruncated t|t t|t kernels for which we set w (H)=0,j >t. t,j 21

3.1 Data Generating Process In order to simulate from a practically relevant Data Generating Process (DGP), we follow Montiel-Olea et al. (2021) and calibrate a time-varying parameter VAR model based on actual macroeconomic data. Building on the oil market literature (Kilian, 2008, 2009), we fit the TVP VAR kernel estimators to a monthly trivariate dataset of size T = 377. The dataset includes the change in (log) global crude oil production, an index for real economic activity and the log of real oil price. A total of three lags are considered, yielding the following estimated structural VAR model: (cid:16) (cid:17) y = A ˆ (H)x +chol Σ ˆ1/2(H) Qε , ε ∼ N(0,I), t t t t t t where the bandwidth is set to H = 100 and Q is a rotation matrix set in a way that b ∝ [1,1,−1]′ resembles a supply shock in a fixed parameter model (H → ∞). Finally, 1 an instrument z is generated by the following measurement error equation: t z = ϕ ε +σ η , η ∼ N(0,1), t z 1t z t t where we consider θstrong = {ϕ = 0.86,σ = 0.06} and θweak = {ϕ = 0.48,σ = 0.71} z z z z following two parameter constellations proposed in Montiel-Olea et al. (2021) that yield a ˜ strong- and weak instrument for relative IRFs λ in the fixed parameter case. However, k,i as we will document, these parameter constellations do not necessarily translate to strong andweakinstrumentdynamicsforabsoluteIRFsλ , sinceitleveragesalltheinformation k,i coming from Γ and Σ exploiting the underlying shock invertibility assumption. The true impulse response functions for the resulting DGP are given in Figure 1 for three horizons: h = 0,10,20. As visible from the chart, they display substantial time-variation over the sample period. For example, the impact effect (h = 0) of the supply shock on the first variable (dprod ) halves over the sample, while at h = 10 and h = 20 the sign t changes at around t = 300. Similar patterns for magnitudes and signs are present in the 22

IRFs of the other two variables, although with an increased persistence. 4.5 4.5 2 4 4 1 3.5 3.5 3 0 3 2.5 F F F R 2 R2.5 R-1 I I I 1.5 2 -2 1 1.5 0.5 -3 1 0 -0.5 0.5 -4 100 200 300 100 200 300 100 200 300 Time Time Time Figure 1: True impulse response functions λ for horizons h = 0, h = 10 and h = 20. h,i,t 3.2 Empirical coverage Weproceedsimulatingatotalof5000datasetsfromtheDGP,eachofsamplesizeT = 377. Empirical coverage at 95% confidence level is then computed for TVP kernel estimates ˜ of λ (IV-SVAR) and λ (internal instrument VAR). For the latter, IRFs are reh,i,t h,i,t,t b standardized to increase the first variable by one at the fixed time point t = T/2, and b for ease of readability we drop the subscript in the remainder of this section. While we assume the lag length to be known during the Monte Carlo exercise, we explore empirical coverage under both the true bandwidth and the simple data-driven selection method described in section 2.4. Finally, to keep the discussion simple, we focus on empirical coverage at two points of time, t = T/2 and t = 3/4T. 23

Figure 2: Estimated empirical coverage at 95% confidence level obtained for λ (red) h,i,t ˜ and λ (blue) at t = 1/2T = 189 (first row) and t = 3/4T = 283 (second row). h,i,t θstrong = {ϕ = 0.86,σ = 0.06} and H is known. Confidence Sets (CS) based on the z z Delta Method (DM) are highlighted by diamonds, while Anderson Rubin confidence sets (AR) by stars. Figure 2 shows simulation results under the strong instrument parameters setting and known bandwidth. Regarding absolute IRFs λ (red), we document that the Delta h,i,t Method(DM)andAndersonRubin(AR)confidencesets(CS)largelycoincide, suggesting that the instrument remains strong. Empirical coverage for λ is very close to the h,1,t nominal confidence level (95%), while that of λ and λ is still reasonable with h,2,t h,3,t values of more than 90% at most horizons. The worst performance we document is for λ at t = T/2 = 189 for horizons h = 5 and h = 10, where coverage is about 85%. h,3,t ˜ With respect to estimates of relative IRFs λ (blue), we document that the DM- and h,i,t AR CS provide somewhat different empirical coverage. This suggests that the instrument may be weak in the time-varying case. This is not surprising, given that relative to a constant parameter setup, the kernel estimator is subject to a lower effective sample size. Given the resulting weak instrument problem, the Delta Method generally provides ˜ worse coverage than the AR CS. For the IRFs of the first variable, λ , the DM provides h,1,t 24

considerable over-coverage throughout horizons, while the AR confidence sets are close to ˜ ˜ the nominal level. On the other hand, for λ and λ the DM provides coverage that h,2,t h,3,t is generally too low, whereas coverage by the weak-IV robust confidence sets performs better and remains at or above 90%. Figure 3: Estimated empirical coverage at 95% confidence level obtained for λ (red) h,i,t ˜ and λ (blue) at t = 1/2T = 189 (first row) and t = 3/4T = 283 (second row). h,i,t θstrong = {ϕ = 0.86,σ = 0.06} and H is estimated. Confidence Sets (CS) based on the z z Delta Method (DM) are highlighted by diamonds, while Anderson Rubin confidence sets (AR) by stars. Figure 3 shows equivalent simulations when H is chosen by the data-driven method we describe in 2.5. We document very similar coverage rates with exception of λ where h,2,289 we see some deterioration to levels of about 80%. 25

Figure 4: Estimated empirical coverage at 95% confidence level obtained for λ (red) h,i,t ˜ and λ (blue) at t = 1/2T = 189 (first row) and t = 3/4T = 283 (second row). h,i,t θweak = {ϕ = 0.48,σ = 0.71}. Confidence Sets (CS) based on the Delta Method (DM) z z are highlighted by diamonds, while Anderson Rubin confidence sets (AR) by stars. For the parameter constellation θweak, results obtained under the known bandwidth are reported in Figure 4. Starting with λ (red), we find very similar results than reported h,i,t previously in the strong instrument case. Generally, the coverage remains satisfactory nearthenominalvalueof95%. Theworstcoverageisobtainedforλ withratesslightly h,3,t above 85%. Interestingly, it is still the case that only marginal differences arise between the DM- and AR confidence sets, suggesting that the weak instrument problem created byMontiel-Oleaetal.(2021)forrelativeIRFsdoesnottranslatetoabsoluteIRFs, despite the lower effective sample size in the time-varying case.5 ˜ With respect to λ (blue), the performance of both confidence sets deteriorates as h,i,t one would expect when the instruments becomes weaker. Still, the AR confidence sets perform better than the DM, remaining closer to the nominal 95% level. However, one ˜ starts to observe some over-coverage, particularly for λ . h,1,t 5Simulation results available upon request find that for absolute IRFs, a much weaker instrument is needed to note a difference between the AR and DM confidence sets, e.g. θ ={ϕ =0.48/4,σ =0.71}. z z 26

Figure 5: Estimated empirical coverage at 95% confidence level obtained for λ (red) h,i,t ˜ and λ (blue) at t = 1/2T = 189 (first row) and t = 3/4T = 283 (second row). h,i,t θweak = {ϕ = 0.48,σ = 0.71}. and H is estimated. Confidence Sets (CS) based on the z z Delta Method (DM) are highlighted by diamonds, while Anderson Rubin confidence sets (AR) by stars. Similar to the first parameter constellation for the instrument, choosing H by a datadriven method yields broadly similar results, with only minor deterioration in coverage rates for some of the IRFs (see Figure 5). In Appendix C, we provide supplementary Monte Carlo results obtained for larger effective sample sizes. Here, we interpolate coefficients linearly to obtain an equivalent shape in the time-varying coefficients, but spread out over a larger sample and hence much smoother. We choose T = 30×377 = 11310, and let the kernel bandwidth for estimation √ increase by H = 30×100 = 547. Our findings suggest that estimated empirical coverage rates get very close to the nominal size as one would expect in large effective sample sizes. 27

4 The time-varying effects of oil supply news on US industrial production Inthefollowing,weillustratetheuseofourmethodologyrevisitingtheeffectsofoilsupply news shocks on US industrial production. Our analysis builds on the work of Känzig (2021) who studies the effects of exogenous changes in oil price expectations caused by OPEC communications, an intergovernmental organization of major oil-producing nations. To capture changes in oil prices orthogonal to the business cycle, Känzig (2021) constructs an external instrument based on quotes of WTI oil price futures in a narrow window around OPEC production quota announcements. A constant parameter IV- SVAR estimated over the period from 1974 to 2017 suggest consequences for the US economy that mimic a typical supply shock; activity falls, as measured by US industrial production, while both consumer prices and inflation expectations rise. Figure 6: US petroleum consumption, production, and net imports (1950-2023). Source: US Energy Information Administration. Basedonourmethodology, weextendtheanalysistostudyinstabilitiesovertime. Alarge body of literature has found that the relationship between oil prices and US macroeco- 28

nomic conditions has changed, see for example Baumeister and Peersman (2013), Ramey and Vine (2011). Kilian (2009) notes that large part of the instability can be explained by the varying importance of supply- and demand shocks. However, even conditioning on oil-market specific supply shocks a large degree of time-variation remains (Baumeister and Peersman, 2013). A variety of potential drivers have been put forward to explain the variation, including thetime-varyingoilintensityofeconomicactivity, improvedmonetarypolicy, orchanging importance of certain sectors in the economy (Ramey and Vine, 2011). In recent history, a plausible explanation may also be the shale oil revolution. A combination of hydraulic fracturing and horizontal drilling allowed the US to sharply increase its production of crude oil and natural gas. Over the period 2005-2023, total US petroleum production more than doubled, from an average of 7.9 to 21.7 million barrels per day, as shown in Figure 6 (black line).6 This allowed to United States to transition from a large petroleum net-importer to a petroleum net-exporter in 2020 (blue line, Figure 6). Indeed, Bjørnland and Skretting (2024) document evidence that the shale oil revolution aligns well with changes in the transmission of oil-market specific shocks to the US economy. Within a Bayesian time-varying factor model, the authors identify an oil-market specific shock by exclusion restrictions, finding that US industrial production and investment reacts more positive to oil price increases since the the shale-oil revolution, boosted by activity in oil-intensive regions and industries. Our empirical findings complement those results, instead relying on kernel estimators and on the identification of oil-supply news shocks by instrumental variables. 6Petroleum production includes field production of crude oil and natural gas, as well as products produced from refining crude oil and from processing natural gas plant liquids. 29

4.1 Data and identification strategy As in Känzig (2021), our VAR model includes a measure of real oil prices, world crude oil production, a proxy of world crude oil stocks, and world industrial production.7 We augment the model further by two sub-aggregates of the index of US industrial production, that is manufacturing and mining output.8 All variables are included in log levels. To identify shocks to oil supply expectations, we rely on an updated surprise series as external instrument, made available on the homepage of Diego Kanzig. The estimation sample includes data from January 1974 to December 2023. To mitigate the effect of Covid outliers on our estimates, we use a series of dummies, thereby discounting any signal from the data between February 2020 and June 2022. We confirm that this is equivalent to setting the kernel weighting function to zero over that time period. Following the original paper, we use p = 13 lags. ForthehyperparameterH governingoveralltime-variation,thecrossvalidationprocedure discussed in section 2.5 suggests a bandwidth of H ≈ 190 when applied to the pre-covid period.9 We note that the objective function is relatively flat between 110 and 250, yielding at most 5% deterioration in the MSE relative to H = 190. For this reason, we choose a slightly lower bandwith of H = 150 allow for more time-variation in the IRFs, particularly since the shale-oil revolution takes place relatively late in the sample. A sequence of Granger causality tests presented in Table 1 provide no clear evidence at the 5% significance level that the instrument is Granger causing the endogenous variables in the model, with the exception of the very end of the sample. A constant parameter model yieldssimilarconclusionswithp-valuesof0.8and0.85fortheWaldandF-test. Therefore, 7The real price of oil is defined as the WTI price deflated by US CPI, and the proxy of world crude oil stocks is included in seasonally adjusted log levels. World industrial production is downloaded from Christiane Baumeisters homepage, see Baumeister and Hamilton (2019). 8The US industrial production index measures the combined real output of the manufacturing, mining, and electric and gas utilities industries. The variability in the latter is mostly driven by weather, and hence excluded from our analysis. 9Here,theobjectivefunctionisbasedonoutofsampleone-stepaheadconditionalforecastscomputed over the second half of the sample. 30

weproceedassumingshockinvertibility, andstudyabsoluteimpulseresponsefunctionsto ashockofunitstandarddeviationthroughouttime(λ ). Toassessthesensitivityofthe k,i,t results to our model choices, Appendix D displays estimates using different bandwidths, and relative IRFs estimated relying on the internal instrument VAR.10 Table 1: Granger causality test results computed at different points of time for the null hypothesis that z does not predict y in a VAR for y˜ = [z ,y′]′. t t t t t date t July 77 May 86 Feb 95 Dec 03 Sep 12 Jun 21 Wald Statistic 96.32 78.49 74.70 81.32 100.89 137.11 p-value 0.08 0.46 0.59 0.38 0.04 0 F Statistic 1.23 1.01 0.96 1.04 1.29 1.76 p-value 0.20 0.49 0.57 0.44 0.15 0.01 The F-test is based on (n−1)p = 78 nominator degrees of freedom (dof), and H −np−1 = 406 denominator dofs (Lütkepohl, 2005). 10We also test locally for residual autocorrelation and find no evidence thereof, see Appendix D. 31

4.2 Results Figure7: Time-varyingimpulseresponsefunctionstoanoil-supplyshockofunitvariance. Figure 7 displays estimates of time-varying IRFs to an oil supply news shock at various points in time over the sample (blue line) and compares it to the constant parameter estimates (red line). Shaded areas denote 90% confidence intervals. As expected, the constant parameter results replicate those of Känzig (2021). A supply 32

news shock raises the oil price for up to 5 years. Crude oil production is declining gradually, while crude stocks are rising reflecting precaution by market participants. Global economic activity declines, measured by world industrial production, and so does US manufacturing output. US mining output, which in large part reflects the extraction of oil and gas, increases slightly but with some lag.11 There is a striking amount of time-variation in the transmission of the shock that aligns well with the US shale-oil revolution. An oil supply news shock of constant size is estimated to have 2/3 of the price effect towards the end of the sample, compared to estimates from 1977-2003. Furthermore, the price effect is notably less persistent. Despite the overall more muted price signal, US mining output is estimated to increase by larger amounts towards the end of the sample, and react more quickly. Such a quicker reaction of US mining output aligns well with micro-evidence on larger price elasticity of shale-oil producers (Aastveit et al., 2022). Since 2012, world crude oil production is no longer estimated to decline significantly in reaction to an oil-supply news shock, but instead increases somewhat in the shortrun. This may reflect, in part, that increasing output by non-OPEC oil producers is able to offset OPEC production declines. The IRF of the world industrial production (IP) index is no longer estimated to decline but instead increases temporarily. However, since the world IP index includes mining output it is difficult to disentangle how much of the response reflects increased oil and gas output by non-OPEC states. Indeed, for US manufacturing output, the time-varying effects are less pronounced. Still, there is striking evidence that the oil price shock no longer triggers a significant decline in US manufacturing output, questioning if it still resembles a cost-push shock for the United States. To shed more light on the drivers of the aggregate time-variation documented for US 11Relative importance weights for the US industrial production index suggest that extraction of oil and natural gas (NAICS 2111) reflects currently about 70% of US mining output. 33

minining and manufacturing output, we further study responses for three digit (manufacturing) and four digit (mining) industries according to the North American Industry Classification System (NAICS). The collection of responses are obtained by augmenting the model with one industry at the time, while netting out the respective industry variation of the total manufacturing or mining index included in the baseline VAR. Figure 8 provides an overview of IRFs for selected industries that display large degrees of time-variation between 1995 and 2021. For a complete picture of each industry, we refer to Appendix D. Figure 8: Estimated impulse response functions to an oil-supply shock of unit variance for selected industries at three-(manufacturing) and four- (mining) digit NAICS level. For comparison, the red line shows point estimates from a constant parameter IV-SVAR. Within US mining output, oil and gas extraction (NAICS 2111) shows the strongest pattern. For the mid 90’s, our methodology points to a very slow and imprecisely estimated response. However, in recent times, we find a strong and rapid response which is significant at the 90% confidence interval. Within manufacturing, the response of Petroleum and Coal Products (NAICS 324) shows astrongpatternoftime-variation. Here,thedominantprocessispetroleumrefiningwhich is downstream to Oil and Gas extraction. Our estimates for 1995 show a response close to a constant parameter model, where the industry output declines significantly up to 20 months as the cost of crude increases. Nowadays, our estimates point towards no significant decline, but instead a slight increase in output after 2-3 years. 34

Estimatesfortwootherdownstreamindustries, thatisChemicalsandPlasticandRubber Products, also have varied significantly, moving from contractionary to insignificant territory. Four other industry responses stand out. For Transportation Equipment (NAICS 336), which reflects mostly the production of motor vehicles and aircrafts, oil supply news shocks tended to be strongly contractionary in line with larger operating costs of the produced goods. There is no statistical evidence this is still the case, as the industry output has become insensitive to the oil price increase. Finally, and somewhat surprisingly, the responses of Primary Metal (NAICS 331) and Fabricated Metal products (NAICS 332) have switched in sign over the short run. As they currently represent a combined 12% of manufacturing output in the US, the time-variation in those sectors is likely contributing strongly to the total manufacturing response. Summing up, our findings suggest that OPEC oil-supply news shocks transmit differently inrecenttimesthanaconstantparametermodelwouldsuggest. Theoilpriceeffectsseem to have declined and are less persistent, while there is no evidence that global crude oil production still declines. Within the US, mining output responds more positive and at a faster pace, while there is no evidence that the shock is still contractionary for the US manufacturing sector. 5 Conclusion In this paper, we develop kernel based estimators for time varying impulse response functions of structural VAR models identified by external instruments. Compared to prominent Bayesian approaches, our frequentist estimators are particularly simple to implement, computationally efficient and require no choice for the law of motion and corresponding priors. The amount of time-variation in a given dataset can be set in an automatic fashion, e.g. by optimizing out-of-sample model fit. Importantly, inference can be reliably conducted even if identification is only weak. 35

We illustrate the methodology revisiting the influential paper by Känzig (2021) on the transmission of oil-supply news shocks. We find strong patterns of time-variation, mostly aligning in time with the shale-oil revolution. While for much of the sample an oil price increase due to oil-supply news resulted in clear headwinds to US manufacturing output, this is no longer the case in more recent history. 36

References Aastveit, K. A., H. C. Bjørnland, and T. S. Gundersen (2022): “The price responsiveness of shale producers: Evidence from micro data,” Available at SSRN 4273926. Amir-Ahmadi, P., C. Matthes, and M.-C. Wang (2023): “Understanding InstrumentsinMacroeconomics-AStudyofHigh-FrequencyIdentification,” Tech.rep., Working paper. Arias, J. E., J. F. Rubio-Ramírez, and D. F. Waggoner (2021): “Inference in Bayesian Proxy-SVARs,” Journal of Econometrics, 225, 88–106. Baumeister, C. and J. D. Hamilton (2019): “Structural interpretation of vector autoregressions with incomplete identification: Revisiting the role of oil supply and demand shocks,” American Economic Review, 109, 1873–1910. Baumeister, C. and G. Peersman (2013): “Time-varying effects of oil supply shocks on the US economy,” American Economic Journal: Macroeconomics, 5, 1–28. Bjørnland, H. C. and J. Skretting (2024): “The shale oil boom and the US economy: Spillovers and time-varying effects,” Journal of Applied Econometrics. Caldara, D. and E. Herbst (2019): “Monetary policy, real activity, and credit spreads: Evidence from Bayesian Proxy SVARs,” American Economic Journal: Macroeconomics, 11, 157–92. Cogley, T. and T. J. Sargent (2005): “Drifts and volatilities: monetary policies and outcomes in the post WWII US,” Review of Economic dynamics, 8, 262–302. Dahlhaus, R. (1997): “Fitting time series models to nonstationary processes,” The annals of Statistics, 25, 1–37. 37

Fieller, E. C. (1944): “A fundamental formula in the statistics of biological assay, and some applications,” Quart. J. Pharm, 17, 117–123. Forni, M., L. Gambetti, G. Ricco, et al. (2023): “External Instrument SVAR Analysis for Noninvertible Shocks,” Working Papers 2023-03, Center for Research in Economics and Statistics. Gertler, M. and P. Karadi (2015): “Monetary policy surprises, credit costs, and economic activity,” American Economic Journal-Macroeconomics, 7, 44–76. Giacomini, R., T. Kitagawa, and M. Read (2022): “Robust Bayesian inference in proxy SVARs,” Journal of Econometrics, 228, 107–126. Giraitis, L., G. Kapetanios, and Y. Li (2024a): “Regression Modelling under General Heterogeneity,” Preprint, No. 983, Queen Mary University of London. Giraitis, L., G. Kapetanios, and M. Marcellino (2021): “Time-varying instrumental variable estimation,” Journal of Econometrics, 224, 394–415. Giraitis, L., G. Kapetanios, and T. Yates (2014): “Inference on stochastic timevarying coefficient models,” Journal of Econometrics, 179, 46–65. ——— (2018): “Inference on multivariate heteroscedastic time varying random coefficient models,” Journal of Time Series Analysis, 39, 129–149. Giraitis, L., Y. Li, and P. C. B. Phillips (2024b): “Robust inference on correlation under general heterogeneity,” Journal of Econometrics, 240, 105–120. Gonçalves, S. and L. Kilian (2004): “Bootstrapping autoregressions with conditional heteroskedasticity of unknown form,” Journal of Econometrics, 123, 89–120. Hipp, R. (2020): “On causal networks of financial firms: Structural identification via non-parametric heteroskedasticity,” Tech. rep., Bank of Canada. 38

Inoue, A., B. Rossi, and Y. Wang (2024a): “Has the Phillips Curve Flattened?” CEPR Discussion Paper No. 18846. ——— (2024b): “Local projections in unstable environments,” Journal of Econometrics, 105726. Jarociński, M. and P. Karadi (2020): “Deconstructing monetary policy surprises—theroleofinformationshocks,” American Economic Journal: Macroeconomics, 12, 1–43. Känzig, D. R. (2021): “The macroeconomic effects of oil supply news: Evidence from OPEC announcements,” American Economic Review, 111, 1092–1125. Kapetanios, G., M. Marcellino, and F. Venditti (2019): “Large time-varying parameter VARs: A nonparametric approach,” Journal of Applied Econometrics, 34, 1027–1049. Kilian, L. (2008): “Exogenous oil supply shocks: How big are they and how much do they matter for the U.S. economy?” The Review of Economics and Statistics, 90, 216–240. ——— (2009): “Not all oil price shocks are alike: Disentangling demand and supply shocks in the crude oil market,” American Economic Review, 99, 1053–1069. Koop, G., M. H. Pesaran, and S. M. Potter (1996): “Impulse response analysis in nonlinear multivariate models,” Journal of econometrics, 74, 119–147. Lütkepohl, H. (1993): “Testing for causation between two variables in higher dimensional VAR models,” in Studies in Applied Econometrics, ed. by H. Schneeweiß and K. F. Zimmermann, Springer-Verlag, Heidelberg, 75–91. ——— (2005): “New Introduction to Multiple Time Series Analysis,” Springer Books. 39

Mertens, K. and J. L. Montiel Olea (2018): “Marginal tax rates and income: New time series evidence,” The Quarterly Journal of Economics, 133, 1803–1884. Mertens, K. and M. O. Ravn (2013): “The dynamic effects of personal and corporate incometaxchangesintheUnitedStates,” American Economic Review, 103, 1212–1247. Miranda-Agrippino, S., S. H. Hoke, K. Bluwstein, et al. (2024): “Patents, News, and Business Cycles,” Staff Working Paper No. 788, Bank of England. Montiel-Olea, J. L., J. H. Stock, and M. W. Watson (2021): “Inference in structural vector autoregressions identified with an external instrument,” Journal of Econometrics, 225, 74–87. Müller, U. K. and P.-E. Petalas (2010): “Efficient estimation of the parameter path in unstable time series models,” The Review of Economic Studies, 77, 1508–1539. Paul, P. (2020): “The time-varying effect of monetary policy on asset prices,” Review of Economics and Statistics, 102, 690–704. Plagborg-Møller, M. and C. K. Wolf (2021): “Local projections and VARs estimate the same impulse responses,” Econometrica, 89, 955–980. ———(2022): “Instrumentalvariableidentificationofdynamicvariancedecompositions,” Journal of Political Economy, 130, 2164–2202. Primiceri, G. E. (2005): “Time varying structural vector autoregressions and monetary policy,” Review of Economic Studies, 72, 821–852. Ramey, V. A. and D. J. Vine (2011): “Oil, automobiles, and the US economy: How much have things really changed?” NBER macroeconomics annual, 25, 333–368. Staiger, D. and J. H. Stock (1997): “Instrumental Variables Regression with Weak Instruments,” Econometrica, 65, 557–586. 40

Stock, J. (2008): “What Is New in Econometrics: Time Series,” Tech. rep., Lecture 7. In: Short Course Lectures, NBER Summer Institute. Stock, J. H. and M. W. Watson (1996): “Evidence on structural instability in macroeconomic time series relations,” Journal of Business & Economic Statistics, 14, 11–30. ——— (2012): “Disentangling the channels of the 2007-09 recession,” Brookings Papers on Economic Activity, 43, 81–156. ——— (2016): “Dynamic factor models, factor-augmented vector autoregressions, and structural vector autoregressions in macroeconomics,” in Handbook of macroeconomics, Elsevier, vol. 2, 415–525. ———(2018): “Identificationandestimationofdynamiccausaleffectsinmacroeconomics using external instruments,” The Economic Journal, 128, 917–948. Waggoner, D. F. and T. Zha (2003): “A Gibbs sampler for structural vector autoregressions,” Journal of Economic Dynamics and Control, 28, 349–366. 41

Appendix A Proofs A.1 Proof of Theorem 1 Theorem 1 states that under Assumption 3-5 and H = o(T1), it holds that: 2   ˆ β −β t t   √   H  Γ ˆ t −Γ t   → d N(0,V θt ),     ˆ vech(Σ )−σ t t for V = S Π S′ and θt t ww,t t   I ⊗Π−1 0 0  n x,t    S t =   − (cid:0) I n ⊗Π xz,t Π− x, 1 t (cid:1) I 0   ,     0 0 S σ for Π = plim 1 (cid:80)T w x x′, Π = plim 1 (cid:80)T w z x , x,t T→∞H j=1 j,t j j xz,t T→∞H j=1 j,t j j Π = plim 1 (cid:80)T w2 ξ ξ′, ξ = [vec(x u )′,(z u −Γ)′,vec (cid:0) u′u −Σ (cid:1)′ ]′, and S ww,t T→∞H j=1 j,t j j j j j j j j j t σ such that vech(Σ ) = S vec(Σ ). t σ t For x′ = [y′ ,y′ ,...,y′ ,1]′ a 1×k vector, the model reads t t−1 t−2 t−p y′ = x′ Θ +u′ t t t t (cid:124)(cid:123)(cid:122)(cid:125) (cid:124)(cid:123)(cid:122)(cid:125)(cid:124)(cid:123)(cid:122)(cid:125) 1×n 1×k k×n y = (I ⊗x′) β + u t n t t t (cid:124)(cid:123)(cid:122)(cid:125) (cid:124) (cid:123)(cid:122) (cid:125)(cid:124)(cid:123)(cid:122)(cid:125) (cid:124)(cid:123)(cid:122)(cid:125) n×1 n×nk nk×1 n×1 y = x˜ β +u t t t t where x˜ = (I ⊗ x′) and β = vec(Θ ). Let z be a m × 1 random vector that is t n t t t t correlated with u . We wish to consider estimating Γ = E(u z′) allowing for this t t t t (cid:124)(cid:123)(cid:122)(cid:125) n×m quantity to vary over time. To do so we wish to derive the asymptotic distribution of √ 1 H (cid:80) j w t,j (H)(uˆ t,j z j ′ −E (cid:0) u j z j ′ (cid:1) ) where w t,j (H 2 ) = (cid:80) H j w˜ w˜ t, t j ,j (H (H ) ) , uˆ j = y j −x˜ j β ˆ t and 1

(cid:34) (cid:35)−1(cid:34) (cid:35) T T (cid:88) (cid:88) β ˆ = I ⊗ w (H )x x′ w (H )vec(x y′) t n t,j 1 j j t,j 1 j j j=1 j=1 where w (H) = K(|t−j|/H), (18) t,j where H → ∞, H = o(T). K(x), x ∈ (0, a) is a non-negative continuous function with finite or infinite support, such that for some C > 0 and ν > 3, K(x) ≤ C(1+xν)−1, |(d/dx)K(x)| ≤ C(1+xν)−1, x ∈ (0, a). (19) √ (cid:16) (cid:17) ˆ First, consider T β −β : t t   (cid:32) (cid:33)−1 √ (cid:16) (cid:17) 1 (cid:88) T 1 (cid:88) T H β ˆ t −β t = I n ⊗ H w t,j (H)x j x′ j  √ H w t,j (H)vec(x j u′ j ) j=1 j=1 (cid:124) (cid:123)(cid:122) (cid:125) St,xx(H) T 1 (cid:88) = S (H)√ w (H)vec(x u′) t,xx t,j j j H j=1 T 1 (cid:88) = S (H)√ w (H)(I ⊗x )u t,xx t,j n j j H j=1 Next,considerΓ ˆ = 1 (cid:80)T w (H)uˆ z′ andγˆ = vec(Γ ˆ ) = 1 (cid:80)T w (H)(z ⊗I )uˆ . t H j=1 t,j t,j j t t H j=1 t,j j n t,j (cid:16) (cid:17) Denote by γ = E [vec(u z′)] and use that uˆ = u −x˜ β ˆ −β : t t t j t t t t t (cid:32) (cid:33) √ 1 (cid:88) T H(γˆ −γ ) = √ w (H)uˆ z′ −Γ t t t,j t,j j t H j=1 (cid:32) (cid:33) T T 1 (cid:88) 1 (cid:88) (cid:16) (cid:17) = √ w (H)(z ⊗I )u −γ − w (H)(z ⊗I )x˜ H β ˆ −β t,j j n j t t,j j n j t t H H j=1 j=1 (cid:32) (cid:33) (cid:32) (cid:33) 1 (cid:88) T 1 (cid:88) T √ (cid:16) (cid:17) = √ w (H)(z ⊗I )u −γ − w (H)(z ⊗I )x˜ H β ˆ −β . t,j j n j t t,j j n j t t H H j=1 j=1 (cid:124) (cid:123)(cid:122) (cid:125) St,zx(H) 2

  S (H) 0 0 t,xx     Define S t (H) =   −S t,zx (H)S t,xx (H) I 0   , then it is:     0 0 S σ       β ˆ −β S (H) 0 0 vec(x u′)  t t   t,xx   j j  √     1 (cid:88) T   H   γˆ t −γ t    =    −S t,zx (H)S t,xx (H) I 0    √ H w t,j (H)   vec(u j z j ′ −Γ t )    j=1       σˆ −σ 0 0 S vec(u u −Σ ) t t σ j j t (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) St ξj and therefore the asymptotic covariance is given by V = S Π S′ for √ 1 (cid:80)T w ξ → θt t ww,t t H j=1 t,j j N(0,Π ). ww,t The results of the Theorem follow directly from Theorem 2.2 of Giraitis et al. (2018) (GKY18) once we account for the presence of the exogenous variable, z (Extension t 1 (E1)) and the introduction of a lag order greater than 1 (Extension 2 (E2)). The only other difference between the analysis of GKY18 and ours is that GKY18 allow for stochasticparameterprocesses. Wechoosetorestrictourselvestodeterministicsequences for the parameter processes, to simplify the presentation of our asymptotic results. We consider each extension in turn, starting with E1. There are two matters relating to proving E1. The first relates to extending Theorem 2.1 of GKY18 to this case (Result E11, (RE11)), and the second is to establish asymptotic normality as in (2.15) of GKY18 (Result E12, (RE12)). RE11 follows immediately by 3 and (6.2)-(6.3) of GKY18. RE12 relates to showing normality of term T (the first term of T ) in page 41 of the n,t;1 n,t online appendix of GKY18. Normality follows immediately by Lemma 6.2 (ii) of GKY18 using Assumption 4. Next, consider E2. The result here follows immediately by considering the companion form given by ˜ y˜ = A y˜ +ν , (20) t t t−1 t 3

  A A ... A  1t 1t pt       I 0 ... 0  where y˜ = (y′,y′ ,...,y′ )′, A ˜ =  , ν = (cid:0) (B ε )′,0,...,0 (cid:1)′ and t t t−1 t−p+1 t   t t t   0 ... ... ...       ... ... I 0 applying Theorem 2.2 of GKY18. ˆ The only result that needs to be proven is the asymptotic independence of β and σˆ . t t We revisit the proof of Theorem 2.2 of GKY18. The asymptotically relevant terms of √ (cid:16) (cid:17) √ ˆ H β −β and H(σˆ −σ )aregivenbyT andq whicharebothdefinedinpage t t t t n,t;1 n,t 41 of the online appendix of GKY18. The expectation of their cross product involves the thirdmomentsofε whicharezerobythesymmetryassumptionofTheorem2provingthe t result. The proof for independence between γˆ and σˆ can be established in an equivalent t t way. A.2 Proof of Corollary 2 We discuss the covariance term of the two differently dated estimators in the statement of the Corollary. To do this we need to extend slightly the work of GKY18. To do so we will revert to the notation of the proof of their Lemma 6.2. Recall ξ := K−1/2b′ε y′ V−1/2a tj 2,t j j−1 ψ,t0 where K = (cid:80)T w , K = (cid:80)T w2 V is defined in (2.17) of GKY18 and b, and t j=1 tj 2,t j=1 tj ψ,t0 a are vectors of constants. Then, following the proof of Lemma 6.2, following (6.28) of (cid:80) GKY18, it suffices to determine the probability limit of w w E[ξ ξ |F ]. |t−j|<h t−1j tj tj t−1j j−1 We note that E[ξ ξ |F ] = K−1/2K−1/2E(b′ε )2a′V−1/2y y′ V−1/2a tj t−1j j−1 2,t 2,t−1 j ψ,t0 j−1 j−1 ψ,t0 where E(b′ε )2 = ||b||2. Setting V ˜ := K−1 (cid:80) w w y y′ , for K = j yyc,t (1),t |t−j|<h t−1j tj j−1 j−1 (q),t (cid:80)n w w , we obtain j=1 t−qj tj j := K ˜ ||b||2a′V−1/2V ˜ V−1/2a = K ˜ ||b||2 +r , tn (1)t ψ,t0 yyc,t ψ,t0 (1)t tn 4

where K ˜ = K−1/2K−1/2K , and r = K ˜ ||b||2a′V−1/2(V ˜ − V )V−1/2. It (q)t 2,t 2,t−q (q),t cn (1)t ψ,t0 yyc,t ψ,t0 ψ,t0 remains to show that r → 0 which involves checking that tn p ˜ ||V −V || = o (1). yyc,t ψ,t0 sp p We need to consider ||K−1 (cid:80)n w w y y′ − V || which is o (1) by Lemma (1),t j=1 t−1j tj j−1 j−1 ψ,t sp p (cid:80) 6.1(i) of GKY18. This of course easily generalises to w w E[ξ ξ |F ] for |t−j|<h t−qj tj tj t−qj j−1 all finite q. A.3 Proof of Theorem 2 All the results of this Theorem follow directly from the proof of Theorem 1. Appendix B Inference for structural impulse response functions This part of the Appendix gives detailed formulas in order to compute closed form Delta Method and Anderson Rubin confidence sets for the (time-varying) IV-SVAR estimator and the internal IV-VAR estimator. B.1 Absolute Impulse Response Functions (IV-SVAR estimator) In this paper, we use the IV-SVAR estimator to recover absolute impulse response functions λ , that is the ith element of the n×1 vector λ . The corresponding function is h,i,t k,t given by: (cid:113) λ ˆ = C (A ˆ )Γ ˆ / Γ ˆ′Σ ˆ−1Γ ˆ k,t k t t t t t Building on the reduced form results given in theorem 1, we get   ˆ β −β t t   √   H  Γ ˆ t −Γ t   → d N(0,V θt ).     ˆ vech(Σ )−σ t t 5

√ (cid:16) (cid:17) ˆ d Starting with the Delta Method, as described in section 2.4, we have H λ −λ → k,t k,t N (0,Ω ), where Ω = J (β ,Γ ,σ )V J (β ,Γ ,σ )′ for β = vec(A ), and k,t k,t k t t t θt k t t t t t (cid:20) (cid:21) ∂λ ∂λ ∂λ k,t k,t k,t J (β ,Γ ,σ ) = : : k t t t ∂β ∂Γ ∂σ t t t is the n × (n2p+n+n(n+1)/2) dimensional gradient. The corresponding derivatives are stated in the following. First, note that C (A ) = J AkJ′ where J = [I ,0,...,0] k t s t s s n and   A A ... A A 1t 2t p−1,t pt     I 0 ... 0 0   n      A t =   0 I n 0 0   .     . . . . . . ... . . . 0       0 0 ... I 0 n ∂λ Hence, it is k,t = 0 for k = 0 while for k > 1: ∂βt ∂λ k,t = (cid:0) (Γ /α )′ ⊗I (cid:1) G , ∂β′ t t n k t where α = (cid:112) Γ′Σ−1Γ and G = ∂vec(C k (At)) = k (cid:80) −1 (cid:2) J(A′)k−1−m (cid:3) ⊗C (A ) (Lütkepohl, t t t t k ∂β′ t m t t m=0 (cid:104) (cid:105) 1993). Next, define ∂[Γt,αt] = ∂Γt : ∂αt where it holds that: ∂[Γ′,σ′]′ ∂[Γ′,σ′]′ ∂[Γ′,σ′]′ t t t t t t ∂Γ t = [I : 0], ∂[Γ′,σ′]′ n t t ∂α 1 t = (cid:0) Γ′Σ−1Γ (cid:1)−1/2(cid:2) 2Γ′Σ−1,− (cid:0) Γ Σ−1 ⊗Γ′Σ−1 (cid:1) D (cid:3) , ∂[Γ′,σ′]′ 2 t t t t t t t t t t t for D is the duplication matrix such that vec(Σ ) = Dvech(Σ ). Also, it holds that: t t ∂λ k,t = C (A ) (cid:2) I /α : Γ /α2 (cid:3) . ∂[Γ′,α′]′ k t n t t t t t Combining both results via the Chain rule yields the missing parts of J (): k (cid:20) (cid:21) ∂λ ∂λ ∂λ ∂[Γ ,α ] k,t k,t k,t t t : = × . ∂Γ ∂σ ∂[Γ′,α′]′ ∂[Γ′,σ′]′ t t t t t t 6

With respect to the AR confidence set, the first step is to obtain the asymptotic distribution of the (n+1)×1 vector:   C (A )Γ k t t L =   k,t   (cid:112) Γ′Σ−1Γ t t t √ (cid:16) (cid:17) forwhichitholdsthatλ = (e′L )/(e′ L ). ViatheDeltaMethodweget: H L ˆ −L → d i,t i k,t n+1 k,t k,t k,t N(0,ΩL ) for ΩL = J(2)(β ,Γ ,σ )V J(2)(β ,Γ ,σ )′ where k,t k,t k t t t θt k t t t (cid:20) (cid:21) ∂L ∂L ∂L J(2)(β ,Γ ,σ ) = k,t : k,t : k,t . k t t t ∂β ∂Γ ∂σ t t t ∂L Similar to above, it holds that k,t = 0 for k = 0 while for k > 1: ∂βt   ∂L k,t =  (Γ′ t ⊗I n )G k . ∂β′   t 0 Finally, the last step is: (cid:20) (cid:21) ∂L ∂L ∂L ∂[Γ ,α ] k,t k,t k,t t t : = × , ∂Γ ∂σ ∂[Γ′,α′]′ ∂[Γ′,σ′]′ t t t t t t where ∂[Γt,αt] is as defined above and ∂[Γ′,σ′]′ t t   ∂L k,t =  C k (A t ) 0 . ∂[Γ′,α′]′   t t 0 1 Next, consider the linear test e′L ˆ −λ e′ L ˆ = 0 with the corresponding Wald test i k,t 0 n+1 k,t statistic q(λ ) = H(e′ i Lˆ k,t −λ0e′ n+1 Lˆ k,t )2 where ωˆ is the ijth element of of Ω ˆL . The AR 0 ωˆii−2λ0ωˆi,n+1+λ2 0 ωˆn+1,n+1 ij k,t confidence set of coverage 1 − a is then given by inverting the test statistic, yielding CSAR{λ |q(λ ) ≤ χ2 }. The inversion can be solved in closed form following, e.g. k,i,t k,i,t 1,1−a footnote 14 in Montiel-Olea et al. (2021). 7

B.2 Relative Impulse Response Functions (internal IV VAR estimator) For relative impulse response functions, the corresponding function of reduced form parameters is given by: λ ˆ˜ = C (cid:16) A ˆ˜ (cid:17) P ˆ˜ /(e′P ˆ˜ ), k,t,t b k t •1,t 2 •1,t b where λ ˆ˜ = e′ λ ˆ˜ . Here, A ˆ˜ , P ˆ˜ = e′chol(Σ ˆ˜ ) and P ˆ˜ = e′chol(Σ ˆ˜ ) are based k,i,t,t b 1+i k,t,t b t •1,t 1 t •1,t b 1 t b on kernel estimates of the TVP internal instrument VAR. Starting from the reduced form results of Theorem 2 and Corollary 3, we have:   √ (cid:16)ˆ˜ ˜ (cid:17) d  ˜ (cid:16) ˜ (cid:17)−1 ˜ (cid:16) ˜ (cid:17)−1 H β −β → N 0,Σ ⊗ Π Π Π , t t t x,t ww,t x,t b   (cid:124) (cid:123)(cid:122) (cid:125) V1   √ (cid:16) (cid:17) H σˆ˜ −σ˜ → d N 0,L Π L′ −σ˜ σ˜′ , t,t b t,t b  2(n+1) uu,uu,t,t b 2(n+1) t,t b t,t b  (cid:124) (cid:123)(cid:122) (cid:125) V2 ˜ ˜ √ (cid:16)ˆ˜ ˜ (cid:17) d (cid:16) ˜ (cid:17) for σ˜ = vech(Σ ). To obtain Ω in H λ −λ → N 0,Ω , an t,t t,t k,t,t k,t,t k,t,t k,t,t b b b b b b (cid:16) (cid:17) (cid:16) (cid:17)′ ˜ ˜ ˜ ˜ application of the Delta method yields Ω = J β ,σ˜ diag(V ,V )J β ,σ˜ for k,t,t k t t,t 1 2 k t t,t b b b (cid:34) (cid:35) (cid:16) (cid:17) ∂λ ˜ ∂λ ˜ ˜ ˜ k,t k,t J β ,σ˜ = : . k t t,t b ∂β ˜′ ∂σ˜′ t t,t b The first part of J ˜ () is given by ∂λ˜ k,t = 0 for k = 0 while for k > 1: K ∂β˜′ t ∂λ ˜ (cid:18) (cid:16) (cid:16) (cid:17)(cid:17)′ (cid:19) k,t = P ˜ / e′P ˜ ⊗I G , ∂β ˜′ •1,t 2 •1,t b n k t for G = ∂vec(C k (A˜ t)) = k (cid:80) −1 (cid:2) J(A′)k−1−m (cid:3) ⊗C (cid:16) A ˜ (cid:17) . To obtain ∂λ˜ k,t we make use of the k ∂β˜ t ′ m=0 t m t ∂σ˜ t ′ ,tb Chain rule. First, consider the Gradient ∂[P˜ • ′ 1,t ,e′ 2 P˜ •1,tb ]′ . To this end, let S be selection ∂σ˜t,tb σt matrix such that σ˜ = S σ˜ and S a selection matrix such that σ˜ = S σ˜ . t σt t,t b σt t b σtb t,t b Define S a matrix of 0 and 1’s such that [P ˜′ ,e′P ˜ ]′ = S [vech(P ˜ )′,vech(P ˜ )]′. P •1,t 2 •1,t b P t t b ˜ ˜ Furthermore, let L be the elimination matrix such that vech(Σ ) = L vec(Σ ), and n+1 t n+1 t K be the mn×mn commutation matrix such that vec(A′) = vec(A) for A any m×n mn 8

matrix. Then:   ∂[P ˜ • ′ 1,t ,e′ 2 P ˜ •1,t b ]′ = S  (cid:16) L n+1 (cid:0) I (n+1)2 +K n+1,n+1 (cid:1) (cid:16) P ˜ t ⊗I n+1 (cid:17) L′ n+1 (cid:17)−1 S σt . ∂σ˜ t,t b P (cid:16) L (cid:0) I +K (cid:1) (cid:16) P ˜ ⊗I (cid:17) L′ (cid:17)−1 S  n+1 (n+1)2 n+1,n+1 t b n+1 n+1 σtb Finally: ∂λ ˜ (cid:16) (cid:17) (cid:20) (cid:16) (cid:17) (cid:16) (cid:17)2 (cid:21) k,t = C A ˜ I / e′P ˜ : P ˜ / e′P ˜ , ∂[P ˜′ ,e′P ˜ ] k t n+1 2 •1,t b •1,t 2 •1,t b •1,t 2 •1,t b ˜ and hence the second part of J () is given by: K ∂λ ˜ ∂λ ˜ ∂[P ˜′ ,e′P ˜ ]′ k,t = k,t × •1,t 2 •1,t b . ∂σ˜′ ∂[P ˜′ ,e′P ˜ ] ∂σ˜′ t,t b •1,t 2 •1,t b t,t b With respect to the AR confidence set, the first step is to obtain the asymptotic distribution of the (n+2)×1 vector:   ˜ C (A )P L ˜ =  k t •1,t  k,t,t b   e′P ˜ 2 •1,t b for which it holds that λ ˜ = (e′ L ˜ )/(e′ L ˜ ). Via the Delta Method we get: i,t 1+i k,t,t b n+2 k,t,t b √ H (cid:16) L ˆ˜ −L ˜ (cid:17) → d N(0,ΩL˜ )forΩL˜ = J ˜(2) (cid:16) β ˜ ,σ˜ (cid:17) diag(V ,V )J ˜(2) (cid:16) β ˜ ,σ˜ (cid:17)′ k,t,t b k,t,t b k,t,t b k,t,t b k t t,t b 1 2 k t t,t b where (cid:34) (cid:35) (cid:16) (cid:17) ∂L ˜ ∂L ˜ J ˜(2) β ˜ ,σ˜ = k,t : k,t . k t t,t b ∂β ˜ ∂σ˜ t t,t b ∂L˜ Similar to above, it holds that k,t,tb = 0 for k = 0 while for k > 1: ∂βt   (cid:16) (cid:17) ∂L ˜ k,t =  P ˜ • ′ 1,t ⊗I n G k . ∂β′   t 0 Finally: ∂L ˜ ∂L ˜ ∂[P ˜′ ,e′P ˜ ]′ k,t,t b = k,t,t b × •1,t 2 •1,t b , ∂σ˜ ∂[P ˜′ ,e′P ˜ ] ∂σ˜′ t,t b •1,t 2 •1,t b t,t b 9

where ∂[P˜ • ′ 1,t ,e′ 2 P˜ •1,tb ]′ is as above and ∂σ˜′ t,tb   ∂L ˜ k,t =  C k (A t ) 0 . ∂[P ˜′ ,e′P ˜ ]   •1,t 2 •1,t b 0 1 Appendix C Supplementary Monte Carlo Results This part of the Appendix illustrates the performance of the kernel based confidence sets in large samples. Specifically, we increase sample size and kernel bandwidth by a factor of √ 30 and 30 respectively, interpolating coefficients from the same data-generating process described in section 3. Figure C.9 shows the true underlying impulse response functions, which display the same dynamics just over a larger time frame. Figure C.9: True impulse response functions λ . h,i,t As reported in Figure C.10, estimated empirical coverage gets very close to the nominal 95% confidence level, in the DGP that considers a strong instrument and H to be known. 10

Figure C.10: Estimated empirical coverage at 95% confidence level obtained for λ h,i,t ˜ (red) and λ (blue) at t = 1/2T = 5655 (first row) and t = 3/4T = 8483 (second row). h,i,t θstrong = {ϕ = 0.86,σ = 0.06} and H known. Confidence Sets (CS) based on the Delta z z Method (DM) are highlighted by diamonds, while Anderson Rubin confidence sets (AR) by stars. Appendix D Supplementary empirical results D.1 Test of residual autocorrelation In this part of the supplementary material, we test if the iid assumption on the residuals may be violated during some periods of the sample. To do so, we follow the textbook treatmentinLütkepohl(2005)leveragingsimplemultivariatePortmanteauTestsadjusted fortheTVPcase. Atthispoint, wewarnthatwehavenotverifiedtheasymptoticvalidity of the test in the time-varying case, nor studied its finite sample properties. A thorough analysis thereof is beyond the scope of our paper. Let the time t sample auto-covariances be C ˆ = 1 (cid:80)T w (H)uˆ uˆ ,i = 1,...,h < it H j=i+1 t j j−i q T and the corresponding sample autocorrelation be R ˆ = D ˆ−1C ˆ D ˆ−1 where D is a it t it t t ˆ diagonal matrix with elements of C on the diagonal. Then, we aim to test the null 0t hypothesis that H : R = (R ,...,R ) = 0 vs H : R ̸= 0. The corresponding 0 hq,t 1,t hq,t 1 hq,t multivariate Portmanteau test statistic and the approximate distribution is given by 11

Q = H (cid:80)hq tr(C ˆ C ˆ−1C ˆ C ˆ−1) ≈ χ2(n2(h −p)) for large sample sizes H and h . hq,t i=1 it 0t it 0t q q Table 2 summarizes the results from a sequence of Portmanteau Tests for residual autocorrelationuptolagh = 26inourempiricalapplication. Thetestpointstonostatistical q evidence of remaining residual autocorrelation throughout the sample. Table 2: Portmanteau test results computed at different points of time for the null hypothesis that R = (R ,...,R ) = 0. hq,t 1,t hq,t date t July 77 May 86 Feb 95 Dec 03 Sep 12 Jun 21 Portmanteau Statistic 332.47 228.76 164.28 141.41 161.57 212.82 p-value 1 1 1 1 1 1 The Portmanteau test is based on n(h −p) = 468 degrees of freedom (Lütkepohl, q 2005). D.2 Robustness of the main results In this part of the supplementary material, we briefly study the robustness of our results with respect to TVP estimates obtained by relative IRFs, and the bandwidth. Figure D.11: Estimate of α obtained in the IV-SVAR. t First, Figure D.11 displays point estimates alongside 90% confidence intervals of α obt tained in the IV-SVAR model. Assuming shock invertibility, there is no strong evidence that α varies over time. This allows to fix the shock size across time and study the relat ˜ tive IRF λ . We choose t such that we have a reasonable local instrument strength. h,i,t,t b b In our case this is December 2003, which is when the Wald test statistic for the null hypothesis of e′P ˜ = 0 is the largest, also guaranteeing finite length of the AR confidence 2 •1,t sets. ˜ The top two rows of Figure D.12 show the estimates obtained for relative IRFs λ h,i,t,t b throughout time, standardized to increase the real oil price by 6.2% in December 2003 12

(H = 150), hence aligning the shock size with λ at that month. While uncertainty is h,i,t ˜ fairly large for our estimates of λ , the baseline point estimates are mostly included h,i,t,t b in the 90% confidence sets. The results are qualitatively similar, in that oil-supply news shocks are no longer clearly contractionary for the US manufacturing sector, and more recently lead to a large boost in the US mining output. The bottom two rows of Figure D.12 show estimates obtained for relative IRFs λ h,i,t obtained under different bandwidths. As expected, for a larger bandwidth (H = 190) the response of manufacturing IP is less time-varying and remains negative towards the end of the sample. Estimates obtained with a lower bandwidth H = 110, instead, are more erratic and turn positive towards the beginning and end of the sample. In any case, the alternative point estimates remain within the 68% confidence intervals obtained in our baseline model. Interestingly, IRF estimates of US mining output are less sensitive to the choice of the bandwidth. 13

Figure D.12: Time-varying impulse response functions of US manufacturing- and mining output to an oil-supply shock. The top two rows contrast baseline estimates obtained ˜ by the IV-SVAR (blue) to estimates of relative IRFs λ obtained by the internal h,i,t,t b instrument VAR (green), standardized to increase the real oil price by 6.2% in December 2003(H = 150). Thebottomtworowscontrastthebaselineresults(blue)obtainedunder H = 150 to estimates under different bandwidths (H = 110 and H = 190). Shaded areas indicate 68% and 90% confidence intervals. ˜ Finally, Figure D.13 shows the entire set of estimates of λ for all six variables in the h,i,t,t b model, and contrasts it to results obtained under the constant parameter case (red). Two findings are worth highlighting. First, the time-varying IRF estimates of the oil-market variablesarequalitativelysimilartothoseobtainedusinganIV-SVARmodel. Second,the ˜ confidence sets for λ are broadly comparable to those obtained in the fix parameter h,i,t,t b model, and at certain times even narrower. Given the smaller effective sample size for the TVP estimator, this means that estimates of the asymptotic covariance are substantially lower. 14

˜ Figure D.13: Time-varying relative impulse response functions λ to an oil-supply h,i,t,t b shock. Estimates are obtained using the internal instrument VAR, standardized to increase the real oil price by 6.2% in December 2003 (H = 150). For comparison, the red thick line indicates relative IRF estimates while dashed lines and the shaded area show corresponding 90% confidence sets. D.3 Complementary results: TVP IRF estimates by industry Thispartoftheappendixincludescomplementaryestimatesoftime-varyingIRFsforvariousindustriesinthemanufacturingandminingsector, sortedbydurablegoodsproducing industries (Figures D.14 and D.15), nondurable goods producing industries (Figures D.16 15

and D.17) and mining (Figure D.18) Figure D.14: Time-varying impulse response functions to an oil-supply shock of unit variance. Shaded areas indicate 65% and 90% confidence sets, while the red line indicates point estimates obtained under fix parameters. 16

Figure D.15: Time-varying impulse response functions to an oil-supply shock of unit variance. Shaded areas indicate 65% and 90% confidence sets, while the red line indicates point estimates obtained under fix parameters. Figure D.16: Time-varying impulse response functions to an oil-supply shock of unit variance. Shaded areas indicate 65% and 90% confidence sets, while the red line indicates point estimates obtained under fix parameters. 17

Figure D.17: Time-varying impulse response functions to an oil-supply shock of unit variance. Shaded areas indicate 65% and 90% confidence sets, while the red line indicates point estimates obtained under fix parameters. Figure D.18: Time-varying impulse response functions to an oil-supply shock of unit variance. Shaded areas indicate 65% and 90% confidence sets, while the red line indicates point estimates obtained under fix parameters. 18

Appendix E A comparison to alternative time-varying IRF estimators A major benefit of the kernel estimator is its simplicity, computational efficiency and the ability to easily include very persistent time-series into the VAR model. Existing methods, including the Bayesian VAR-X model of Paul (2020) and path estimators based on the framework in Müller and Petalas (2010) (see e.g. Inoue et al. (2024b,a)) require to loop through each time t estimate, and hence are computationally more demanding. Furthermore, the Bayesian estimator of IRFs may struggle with explosive posterior draws when persistent time-series are included in the VAR without further transformation. Finally, the methodology in Müller and Petalas (2010) does not offer a joint distribution ofparametersattimetandt ,hencemakingitdifficulttocomputeIRFsofthesameshock b size based on an internal instrument VAR. We also encountered that estimates based on the methodology of Müller and Petalas (2010) tend to favor very unstable parameters if we used the proposed default equal-probability mixture for the random walk weighting functions. Since our empirical applications include long, very persistent time-series we weren’t able to compare our method to alternatives in the empirical application.12 However, in the following, we offer a simple comparison of IRF estimates in a controlled environment, generating the time-series and instrumental variable from a bivariate VAR. We find that all estimator yield fairly similar results. Our experiment is based on time-series of size T = 500, simulated from a Model where y = A y + B ε where ε ∼ N(0,I ), B = B for t = 1,..., T and B = B for t 1 t−1 t t t 2 t 1 2 t 2       0.8 −0.05 0.1 0.2 0.5 0.4 t = T +1,...,T. We set A =  , B =  , B =  . The 2 1   1   2   0.2 0.7 1 1 −1 1 instrument is simulated from z = 0.2ε +0.1η where η = N(0,1). t 1t t t 12Another challenge we face are the large outliers based on the pandemic, which we simply dummy out in the kernel-based methods. It’s less clear how to best treat those observations in the alternative methods. 19

We compute estimates of impulse responses at t = [1/8,2/8,3/8,4/8,5/8,6/8,7/8]T, comparing the kernel based estimator of an IV-SVAR to two alternative methods: the path estimator of an IV-SVAR as proposed in Müller and Petalas (2010), and the VAR- X estimator of Paul (2020). The amount of time-variation is obtained as follows. Our kernel-based estimator relies on a bandwidth of H = T0.5, while the path estimator uses the default equal probability mixture proposed in Müller and Petalas (2010). The Bayesian VAR-X is based on a series of independent random walks for the intercepts, the regression coefficients of the instrument, and each of the autoregressive coefficients. Its variances are treated as random and given a conjugate inverse gamma prior with a mean of 0.012 and five degrees of freedom. Note that unlike the IV-SVAR estimators, estimates based on the VAR-X are unable to identify the shock variance. For that reason, we standardize the IRFs of the VAR-X estimator to increase the first variable by one at t = 1/4T, hence matching the true effects of a unit variance shock. Figures E.19, E.20 and E.21 show the results of the exercise for the kernel, path and Bayesian estimator respectively. All estimators yield fairly similar point estimates and are able to correctly detect the break at the middle of the sample. Naturally, since all methods assume that parameters are smooth, they fail around the break-point at T = 1/2T. However, as they move away from the break estimates get fairly accurate. The width of the confidence intervals seem close for the kernel and VAR-X estimator. Considerable wider intervals are obtained for the path estimator, which puts a lot of weight on very unstable parameters. However, setting up a path estimator with a tighter specification that favors more stable parameters yields a width of the confidence interval that is comparable to the kernel and VAR-X estimator.13 13The path estimator suggested in Müller and Petalas (2010) is based on a multivariate Gaussian random walk with a variance proportional to the inverse Hessian of the likelihood. Their default method choice is based on minimizing weighted average risk (WAR) relative to an equal-probability mixture of 11 values for the constant of proportionality c2/T2, with c ∈ {0,5,...,50}. When we set c ∈ {0,5/10,...,50/10}, we obtain confidence widths that are similar to that of the other two estimators. 20

Figure E.19: Kernel based estimator: Point estimates of IRFs alongside 90% confidence intervals at t = [1/8,2/8,3/8,4/8,5/8,6/8,7/8]T for y (first column) and y (second 1 2 column). The true IRFs are plotted as red lines. 21

Figure E.20: Path estimator: Point estimates of IRFs alongside 90% confidence intervals at t = [1/8,2/8,3/8,4/8,5/8,6/8,7/8]T for y (first column) and y (second column). 1 2 The true IRFs are plotted as red lines. 22

Figure E.21: Bayesian VAR-X: Posterior median estimates of IRFs alongside 90% posterior credible intervals at t = [1/8,2/8,3/8,4/8,5/8,6/8,7/8]T for y (first column) and 1 y (second column). The true IRFs are plotted as red lines. 2 23