Structural Analysis of Vector Autoregressive Models

Structural Analysis of Vector Autoregressive Models

Introduction

This set of lecture notes are meant to discuss fundamental concepts on the mechanisms of structural vector autoregressive models with many economic theory driven examples which are drawn from several seminal studies presented in the existing literature. The particular perspective we follow permits to present several mathematical and theoretical concepts in an informal way while discussing important concepts by presenting how some long-standing problems in the applied macroeconometrics literature have been tackled. Moreover, our discussion of selected studies from the literature in the form of examples that present key aspects in the identification and estimation of SVAR models, main theoretical econometric concepts such as matrix representations and asymptotic theory convergence results; are motivated from a research-based pedagogical approach that encourages the reader to seek a comprehensive understanding of these concepts while fostering further discussions on current state-of-the art methodological approaches and debates in the literature. Note that a full treatment of these topics are presented in the "Time Series Analysis" book of hamilton1994time. Additional material can be found in the books of gourieroux1997time, lutkepohl2005new and kilian2017structural.

More specifically, stock2001vector describe the key job tasks of applied macro-econometricians, which shall also serve as the main learning objectives motivating the preparation of this set of lecture notes. There are summarized as below:

Vector Autoregressive models are a popular econometric tool for estimating the dynamic interaction between the variables included in the system. Moreover, several test statistics are constructed using VAR models such that: Wald tests for Granger Causality, impulse response functions (IRFs) and forecast error variance decomposition (FEVDs). Therefore, inference on these statistics is typically based on either on first-order asymptotic approximations or on bootstrap methods. However, the deviation from i.i.d innovations as in the case of conditional heteroscedasticity, invalidates a number of standard inference procedures such that the application of these methods may lead to conclusions that are not in line with the true underlying dynamics. Consequently, in many VAR applications there is need for inference methods that are valid when innovations are only serially correlated but not independent. Nevertheless, the VAR representation of the economy is suitable for estimation/testing as well as forecasting purposes but is not adequate for policy analysis which is the purpose of this course (thus, our interest in SVARs).

The Identification Problem

Roughly speaking, the identification problem of structural parameters in linear simultaneous equations models are closely related\footnote{The notion of coefficient restrictions was originally extended to show the equivalence relationship between identifiability and instrumental variables estimation, that is, the restrictions required for identification give rise to instrumental variables required for estimation.}. Seminal papers discussing the identification problem include among others sargan1983identification and dufour2003identification. According to hausman1983identification, necessary and sufficient conditions for identification with linear coefficient and covariance restrictions are developed in a limited information context. In particular, imposing covariance restrictions facilitate identification iff they imply that a set of endogenous variables is predetermined in the equation of interest - which generalizes the notion of recursiveness for structural learning and causal recovery. Under full information, covariance restrictions imply that residuals from other equations are predetermined in a particular equation, and, under certain conditions, can facilitate system identification. This implies that in the general case, FIML first order conditions show that if a system of equations is identifiable as a whole, covariance restrictions cause residuals to behave as instruments.

Furthermore, imposing exclusion restrictions and the normalization $\beta_{ii} = 1$, then the classical structural econometric specification for the simultaneous equation model $\boldsymbol{Y} \boldsymbol{B}^{\prime} + \boldsymbol{Z} \boldsymbol{\Gamma}^{\prime} = \boldsymbol{U}$, where the vector $\boldsymbol{y}_i$, for $i \in \left\{ 1,..., G \right\}$, includes $G$ jointly dependent random variables. A key result for system identification purposes is the fact that the relative triangularity of equations $(i,j)$ is precisely equivalent to a zero in the $(i,j)-$th position of $\boldsymbol{B}^{-1}$, denoted as $\boldsymbol{B}^{-1}_{ij}$. As a result, equations $(i,j)$ are relatively recursive if and only if there are no paths by which a shock to $\boldsymbol{u}_j$ can be transmitted to $\boldsymbol{y}_i$.

Next, we consider the equivalence between instrumental variables and covariance restrictions.

A well-known problem to the literature of structural VAR models is that imposing the assumption of Gaussian errors renders non-identifiable structural parameters due to the nonlinearities and complex dependence structure of the system. Additional identifying restrictions are needed to ensure robust identification and estimation of SVARs. Based on these concerns, the framework of lanne2017identification shows that in practice the Gaussian case is an exception in that a SVAR model whose error vector consists of independent non-Gaussian components is, without any additional restrictions, identified and leads to essentially unique responses. The identification scheme employs the non-Gaussianity of structural shocks while an MLE estimator for the parameters of the non-Gaussian SVAR model is consistent and asymptotically normally distributed. Moreover, the scheme allows to impose additional identifying restrictions which can be formulated into testable hypotheses. Further extensions of the non-Gaussian identification strategy includes cases such as when the system includes only strictly positive components (see, nyberg2022structural), cases when a GMM estimation approach is employed (see, lanne2021gmm and keweloh2021generalized) as well as specific cases when structural shocks are leptokurtic (e.g. see lanne2023identifying) or heavy-tailed (e.g., see anttonen2023statistically).

Furthermore, while the reduced-form VAR model can be seen as a convenient description of the joint dynamics of a number of time series that facilitates forecasting, the SVAR model is more appropriate for answering economic questions of theoretical and practical interest. The main tools in analyzing the dynamics in SVAR models are the impulse response function and the forecast error variance decomposition. Moreover, usually when fitting a SVAR model, the joint distribution of the error terms is almost always (either explicitly or implicitly) assumed to follow a multivariate Gaussian distribution which might render the system non-identifiable especially under the presence of non-linearities, higher-order dynamics and heavy-tails. The positive news is that the joint distribution of the reduced-form errors is fully determined by their covariances only but this which implies that the structural errors of the system cannot be identified unless we impose additional linear restrictions (see, lanne2017identification). However the use of suitable orthogonal transformations or even permutations (due to being exchangeable random variables), can be shown to induce observationally equivalent SVAR processes.

Furthermore, we discuss some key identification strategies in the following sections. Generally, speaking understanding all the available tools in our disposal is important when considering relevant empirical and theoretical driven questions in relation to the impact of structural shocks on economic outcomes. Therefore, in this set of lecture notes we present some key concepts and recent developments in the literature with focus on the implementation of econometric and statistical tools for the identification and estimation of structural vector autoregressive models, (SVAR). In particular, we discuss several different approaches that are commonly employed for the structural identification of these interdependent systems in relation to economic phenomena such as business cycle fluctuations and macroeconomic conditions in order to obtain policy recommendations. On the other hand, exogenous variation to macroeconomic conditions such as the adverse effect of climate change across interconnected economies requires novel econometric methods to ensure robust identification and estimation\footnote{In particular, the omission of such an important component as the climate when considering the structural analysis of time series models could lead to non-causality, leading to inaccurate model estimates and forecasts (see, lutkepohl1982non). Further properties on Noncausal vector AR processes are presented by davis2020noncausal.} (e.g., see pesaran2000structural).

The stability of cointegrated VAR processes can be assessed with respect to the long-run relations as well as with respect to the short-run parameters. Usually the long-run equilibrium relations are based on economic theory and so we expect them to be relatively more stable after a perturbation of structural shocks (such as a change in the monetary policy framework), than the short-run coefficients which are data-driven. Consider for example the study of the transmission of monetrary policy shocks to macroeconomic variables using structural cointegrated VAR models. Roughly speaking the identification of SVAR models is achieved by imposing restrictions on the covariance matrix of the residuals of a reduced-form VAR to provide an economic interpretation of shocks while a structural cointegrated VAR model is identified by imposing restrictions on the cointegrating space. Consequently, based on a structural identification scheme of the covariance matrix then metrics such as impulse-responses, variance-decompositions and historical decompositions can be constructed.

In particular, rubio2010structural verifies the identification for linear Gaussian models by checking whether two parameter points are observetionally equivalent iff they have the same reduced-form representation. Moreover, a global identification can be verified by checking whether the model is identified at different points in the parameter space prior to the estimation step. On the other hand, a well-know fact in the literature of nonstationary autoregressive models, asymptotic distribution discontinuities arise in the case when the parameter is at the boundaries of the parameter space.

Following rubio2010structural

Generally, it is well-known that identification in parametric models implies that a specific model structure is identifiable if there is no other structure which is observationally equivalent (see, also sargent1976observational). For instance, the identification scheme proposed by lanne2017identification requires the construction of observationally equivalent SVAR processes when comparing two $MA(\infty)$ representations of the model. In conjecture with conventional rank conditions that ensure the valid identification of SVAR models under non-Gaussianity the identification of structural parameters can be extended into more complex model configurations regardless of the presence of nonlinear or nonstationary dynamics.

More specifically, the main intuition of the identification methodology of rubio2010structural, is to establish general rank conditions for global identification of identified and exactly identified models. Moreover, bognanni2018class highlights that the key conceptual properties which ensure SVAR tractability include: $\textit{(i)}$ observational equivalence of points $( \boldsymbol{A}, \boldsymbol{F} )$ and $( \boldsymbol{A Q}, \boldsymbol{F Q} )$ for $\boldsymbol{Q} \in \mathcal{O}_n$ in the structural model, and $\textit{(ii)}$ the ability to reparametrize the structural model from $( \boldsymbol{A}, \boldsymbol{F} )$ to $( \boldsymbol{B}, \boldsymbol{H}, \boldsymbol{Q} )$ and separate identifiable elements $( \boldsymbol{B}, \boldsymbol{H} )$ and non-identifiable elements $\boldsymbol{Q}$. Thus, the above definitions imply the existence of a unique mapping of reduced-form parameters which is invariant to scaling and reordering.

Another issue discussed by forni2014sufficient is the presence of sufficient information to ensure that the conditions of existence and uniqueness are ensure for identification purposes. Specifically forni2014sufficient present necessary and sufficient conditions under which a VAR contains sufficient information to estimate the structural shocks. Based on this theoretical result the authors propose a test statistic for detecting information deficiency and a procedure to amend a deficient VAR model. A relevant example is the bivariate VAR with unemployment and labour productivity which is found to be informentionally deficient. The building block of this framework is a representation of the economy where $q$ mutually orthogonal structural shocks affect macro variables through square-summable impulse-response functions (see, also kocikecki2018global, kocikecki2023solution).

In practice the number of observable variables $n$ is very large, so that the econometrician needs to reduce it in order to estimate a VAR. The VAR information set is then spanned by an $s-$dimensional sub-vector of $\boldsymbol{x}_t^{*}$, or, more generally, an $s-$dimensional linear combination of $\boldsymbol{x}_t^{*}$, say $\boldsymbol{z}_t^{*} = \boldsymbol{C} \boldsymbol{x}_t^{*}$.

Linear Vector Autoregressions

Cointegration and Vector Autoregressive Processes

Integration Order

An important property of $I(1)$ variables is that there can be linear combinations of these variables that are $I(0)$. If this is so then these variables are said to be cointegrated. Notice that econometric cointegration analysis can be used to overcome difficulties associated with stochastic trends in time series and applied to test whether there exist combinations of non-stationary series that are themselfs stationary.

Cointegrated Vector Autoregressive Models

A VAR has several equivalent representations that are valuable for understanding the interactions between exogeneity, cointegration and economic policy analysis. To start, the levels form of the $s-$th order Gaussian VAR for $x$ is

where $K$ is an $N \times N_0$ matrix of coefficients of the $N_0$ deterministic variables $q_t$. Suppose that we have a vector $Y_t = \left[ y_{1t}, y_{2t},..., y_{nt} \right]^{\prime}$ that does not satisfy the conditions for stationarity. One way to achieve stationarity might be to model $\Delta y_t$, rather than $y_t$ itself. However, differencing can discard important information about the equilibrium relationships between the variables. This is because another way to achieve stationarity can be through linear combinations of the levels of the variables. Thus, if such linear combinations exist then we have cointegration and the variables are said to be cointegrated. The notion of Cointegration has some important implications: (i) It implies a set of dynamic long-run equilibria between the variables, (ii) Estimates of the cointegrating relationships are super-consistent, they converge at rate $T$ rather than $\sqrt{T}$, and (iii) Modelling cointegrated variables allows for separate short-run and long-run dynamic responses. Further details on Vector Autoregression and Cointegration can be found in the corresponding Chapter of watson1994vector.

Cointegration and Dynamic Inference from ADLM

Consider a general autoregressive distributed lag ARDL (p,q) model where a series, $y_t$, is a function of a constant term, $\boldsymbol{\alpha}_0$, past values of itself stretching back $p$periods, contemporaneous and lagged values of an independent variable, $x_t$, of lag order $q$, and indepenendent, identically distributed error term:

The contemporaneous effect of $x_t$ on $y_t$ is given by $\boldsymbol{\beta}_0$. Moreover, the magnitude of $\boldsymbol{\alpha}_1$ informs us about the memory property of $y_t$. Assuming that $0 < \alpha_1 < 1$, larger values indicate that movements in $y_t$ take longer to dissipate. The long-run effect (or long-run multiplier) is the total effect that a change in $x_t$ has on $y_t$. A simple model that incorporates such dynamic effects, is the distributed lag model.

Open Problems Various open problems remain in the literature such as the use of the cross-section as a mechanism for beliefs elicitation and identification of structural shocks. Related studies from the finance literature where the cross-section is used for improving the predictability of portfolio returns include among others chinco2019sparse, freyberger2020dissecting and kyle2023beliefs, but in the case of modelling macroeconomic fundamentals using cross-sectional information is still a growing literature. Relevant studies in this direction include forni2017dynamic, kong2019weak, yamamoto2023cross and hannadige2023forecasting among others.

Vector Autoregressions: Prediction and Granger Causality

Stability Conditions in VAR$(p)$ Model

Relation between Dynamic Structural Models and Vector Autoregressions

Consider that these interelated equations can be written in the following form:

Therefore, by pre-multiplying the above dependent variable with $\boldsymbol{B}_0$ we can obtain a VAR representation given by the following expression:

In other words, we can view the above representation as a special case of the Dynamic Structural system equation since we eliminate the interelations of the dependent variable to a reduced form of a VAR.

Main Assumptions

Consider an $n-$dimensional covariance stationary zero-mean vector stochastic process $\boldsymbol{x}_t$ of observable variables, driven by $q-$dimensional unobservable vector process $u_t$ of structural shocks.

where $\boldsymbol{C} (L) = \displaystyle \sum_{j=0}^{\infty} C_j L^j$ is an one-sided polynomial in the lag operator $L$ in infinite order. These shocks are orthogonal white noises such that $u_t \sim ( 0, \boldsymbol{\Sigma}_u )$, where $\boldsymbol{\Sigma}_u$ is diagonal. Furthermore, notice that a VAR model based on a subset of the variables implies a different $MA(\infty)$ representation than its counterpart based on the joint distribution. We begin by considering various examples regarding the identification and estimation of SVAR models using properties of structural shocks.

If the roots of $\mathsf{det} [ \boldsymbol{C}(z) ]$ are outside the unit disc, we have invertability in the past, that is, the inverse representation depends only on non-negative powers of $L$, and we have fundamentalness. If at least one of the roots of $\mathsf{det} [ \boldsymbol{C}(z) ]$ is inside the unit disc, then invertability and non-fundementalness holds.

The identification of structural vector autoregressions is a crucial step before proceeding with constructing test statistics and forecasts (see Section (ref)). A suitable approach found in the literature, is the method of recursive identification\footnote{The theoretical analysis of recursive identification methods is given by soderstrom1978theoretical and chen2010new } (see, chen2014recursive) which implies exact identification of the structural parameter $B_0^{-1}$ (e.g., see kim2000exchange). Moreover, the presence of stochastic singularity are crucial on whether a system is identified or unidentifiable. Based on these considerations komunjer2011dynamic employ restrictions implied by observational equivalence to establish conditions for avoiding the presence of stochastic singularity. To ensure that these structural econometric models can be still partially identified (see, also phillips1989partially) regardless of the possible presence of stochastic singularity, then the econometrician can exploit the existence and uniqueness of local identification of DSGE models based on the corresponding linearized solution. Lastly, given the prominent role of expectations in news aspects such as the information flow and anticipated shocks can affect the identification of structural shocks (e.g., optimal inter-temporal decisions discount future tax obligations, see leeper2013fiscal and mertens2013dynamic).

Weak Exogeneity in $I(2)$ VAR Systems

The notion of weak exogeneity is important when considering the structural analysis of cointegrating regression models. Under the assumption of weak exogeneity estimation of the cointegration parameters in conditional models is accordingly influenced. In particular, for the VAR model allowing for $I(1)$ variables under the assumption of Gaussian errors implies the use of a cointegrated VAR specification which is suitable for modelling long-run equilibrium dynamics for multivariate time series. Within this stream of literature one is interested to analyze the conditions under which a subset of equations is weakly exogenous with respect to the cointegration parameters (see, paruolo1997asymptotic, paruolo2000asymptotic, paruolo1999weak, dolado1992note and pesaran2000structural). Moreover, tchatoka2013finite propose unified exogeneity test statistics and examine the pivotality property under strict exogeneity (see, also white2014granger) and characterize the finite-sample distributions of the statistics under $H_0$, including when identification is weak and the errors are possibly non-Gaussian. Assume that the vectors $\boldsymbol{U}_t = [ \boldsymbol{u}_t, \boldsymbol{V}_t ]^{\prime}$ for $t = 1,..., T$, have the same nonsingular covariance matrix defined below

where $\boldsymbol{\Sigma}_V$ has dimension $G$ (see, also khalaf2014identification). Then, the covariance matrix of the reduced-form disturbances $\boldsymbol{W}_t = [ \boldsymbol{v}_t, \boldsymbol{V}_t^{\prime} ]^{\prime} = [ \boldsymbol{u} + \boldsymbol{V} \boldsymbol{\beta}, \boldsymbol{V} ]$ takes the following form

where $\boldsymbol{\Omega}$ is positive definite matrix. The exogeneity hypothesis can be expressed as $\mathcal{H}_0: \boldsymbol{\delta} = \boldsymbol{0}$.

Local Projections

Application: Measuring the Impact of Fiscal Policy

This example is based on the study of mertens2010measuring who consider the consequences of anticipation effects for VAR-based estimates of the impact of government spending shocks. We consider a bivariate time series representation of a vector consisting of control variables $z_t$, and government spending $\mathsf{g}_t$. We assume transitory government spending shocks, such that $\mu_{ \mathsf{g} } (L)$ is a stable polynomial. The identification is facilitated by assuming that information can arrive with any anticipation horizon between 1 and $q$ periods. Thus, the vector of observables $\boldsymbol{v}_t = [ \mathsf{g}_t, z_t ]^{\prime}$ is formulated as below

Substituting to the steady-state equilibrium solution of the system, the MA representation is

Two popular methods to estimate impulse responses when a shock has already been identified include: local projections and distributed lag models. In particular, one can show that the dynamic response from a VAR with the shock embedded as an endogenous variable is equivalent to that of a VAR with the shock included as an exogenous variable only when that shock has no serial correlation (e.g., see plagborg2021local, olea2021inference). The result follows because a VAR with a shock as an exogenous variable (VAR-X) can be seen as a multivariate generalization of a DLM. A comparison of these cases can be found in alloza2020dynamic. Moreover, lusompa2023local propose an efficient estimation method for Local projections when residuals are autocorrelated. Lastly, the simultaneous impact of monetary and fiscal policy is discussed by bruneau2003monetary while the relation of fiscal policy and persistent inflation is studied by bianchi2023fiscal.

Estimation Method

Recently, it has been shown that the method of local projections has advantages in drawing inference of impulse responses, especially with respect to uniform validity (e.g., see inoue2020uniform). We demonstrate an example of the estimation method from the framework of xu2023local. The VAR$(p)$ model is of the form $y_t = A_1 y_{t-1} + ... + A_p y_{t-p} + u_t$, $t = 1,...,n$, where $u_t$ is serially uncorrelated shock.

To measure the responses of the first endogenous variable $y_{1t}$ to shocks after $h$ propagation periods, where $h \geq 1$, the LP method runs the regression below (see, xu2023local).

Alternative, we can estimate the following system of equations

Notice that by definition it holds that $\theta_{1,p} (h) = 0$. Recall that for a matrix $x$, we denote its Frobenius norm, such that $|x| = [ \mathsf{trace} ( x^{\prime} x ) ]^{1/2}$. Suppose that $\left\{ \widehat{u}_t(h): t = p^{*},..., n - h \right\}$ are OLS residuals of the $\mathsf{VAR} ( p^{*} - 1 )$ regression, using the data $\left\{ y_t: t = p^{*},..., n - h \right\}$. In addition, related distributional assumptions for the error term of the model can be imposed.

The OLS estimator is given by

Notice that these local projection estimates correspond to an in-sample estimates for the VAR$(p)$ model based on $p$ lags. On the other hand, if we are interested to obtain out-of-sample forecast sequences then we need to implement a forecasting scheme such as a moving window which is used to obtain pseudo-out-of-sample forecasts. The window size should be larger than the number of lags to have good statistical properties for the forecast sequences. Relevant applications of local projections include the construction of impulse response functions.

Impulse Response Functions

A major advantage of modelling interdependent systems such as Structural Vector Autoregressions is that these system representations can be used to construct impulse response functions and forecast error variance decompositions, which are useful for investigating the dynamics within the system as well as the statistical properties of econometric specifications for forecasting purposes (see, baillie2013estimation). Thus, reduced-form impulse responses allows us to find the response of one variable to an impulse in another variable. Usually impulse response functions can be constructed based on the stable VAR$(p)$ process. An alternative method for estimating impulse response functions has been proposed in the literature, using the method of local projections as in jorda2005estimation, barnichon2019impulse, gorodnichenko2020forecast, montiel2021local and plagborg2021local. In addition, the use of exogenous instrumentation as an identification strategy of forecast errors is examined by olea2021inference and plagborg2022instrumental. The construction of robust confidence intervals for IV-based local projections is studied by noh2017impulse, gafarov2018delta, ganics2021confidence and koo2022impulse among others.

Asymptotic Results for VAR Processes with Known Order

Suppose that $\boldsymbol{\beta}$ is an $( n \times 1)$ vector of parameters and $\hat{\boldsymbol{\beta}}$ is an estimator such that

Moreover, let $\big( g_1 ( \boldsymbol{\beta}),..., g_m ( \boldsymbol{\beta}) \big)^{\prime}$ be a continuously differentiable function with values in $m-$dimensional Euclidean space and $\displaystyle \frac{ \partial g_i }{ \partial \boldsymbol{\beta}^{\prime} } = \frac{ \partial g_i }{ \partial \beta_j }$, for $i \in \left\{ 1,..., m \right\}$. Then, it holds that (see, lutkepohl1990asymptotic)

Forecast Error Variance Decomposition

Generally, the Forecast Error Variance Decomposition gives the proportional contribution of each of the structural shocks to the forecast error variance of each variable at different horizons. Recall that the $h-$period forecast error of $y_t$ is given by

Hence, the $h-$period ahead forecast error of the $j-$th component of $y_t$ is given by

Moreover, the variance of the $h-$period forecast error of the $j-$th component of $y_t$ is given by

which holds since $\Sigma_w = I_K$. Notice that the sum of the terms $\left\{ \theta_{jk,0}^2 + ... + \theta_{jk,h-1}^2 \right\}$, is thus the $k-$th shock to the forecast error variance of the $h-$period forecast of the $j-$th variable. Then, we have that

As an example of application of the above metrics see kilian2009impact (and kilian2009not) who consider the effects of three structural shocks (oil supply shock, aggregate demand shock, and oil-specific demand shock), in the US stock market. Another useful metric is the construction of historical decompositions which facilitate the quantification of the contribution of a given structural shock to the historically observed fluctuations in the variables included in the SVAR model. In particular, one might be interested to find to what extent the fiscal policy shock explains recessionary conditions (e.g., see cevik2023s) during uncertain times (e.g., see bloom2009impact). This is useful because we can concentrate on explaining a historical event which we know a priori might have induced a structural change to the macroeconomic system, instead of using for example the average contribution of the fiscal policy shock to the business cycle fluctuations as given by the FEVD (see, kilian2014quantifying).

Suppose that $y_t$ is weakly stationary, and we have observations from 1 to $t$, such that for any $t \in \mathbb{N}$

The MA coefficients $\Theta_s$ converge to zero the further to the past, so $\hat{y}_t = \sum_{s=0}^{t-1} \Theta_s w_{t-s}$ approximates $y_t$.

Identification of Structural Vector Autoregression Models

While theory-based identification through sign restrictions guarantees economic interpretation of the identified shocks (e.g., see antolin2018narrative and mavroeidis2021identification), data-based identification offers diagnostic tools for testing the otherwise just identifying structural model specifications, while remain silent about the economic interpretation of the resulting structural shocks.

Identification using Short-run and Long-run Restrictions

\paragraph{Recursive Identification} Usually recursive identification with $\boldsymbol{B}_0^{-1}$ being a lower diagonal matrix was particularly popular in the earlier SVAR literature which implies exact identification. In particular, when $B_0^{-1}$ is lower triangular, so it is the matrix $B_0$ and thus the SVAR is recursive.

\paragraph{Non-Recursive Models}

A non-recursively identified SVAR model resembles a traditional system of simultaneous equation models. The associated identifying restrictions practically generate moment conditions that can be used to estimate the SVAR model. Furthermore, using such moment conditions one can employ suitably imposed rank conditions as a mechanism for identification (e.g., see blanchard2003empirical). Although exact identification is not directly a testable hypothesis, we can employ the $J-$test for over identifying restrictions.

The GMM estimator is obtained by

Dynamic Structural Factor Models with Overidentifying Restrictions

As an example of formal statistical testing for overidentifying restrictions we present the framework proposed by han2018estimation who develops a new estimator for the impulse response functions (IRFs) in structural factor models with a fixed number of over-identifying restrictions (see, also Chapter 12 in dhrymes2013mathematics). The proposed identification scheme nests the conventional just-identified recursive scheme as a special case. The first step before formulating the test statistic is to present the identification of structural shocks in the structural factor models setting (see, also koistinen2022estimation).

Consider the following structural factor model for $t = 1,...,T$ such that $\boldsymbol{X}_t = \boldsymbol{\Lambda} \boldsymbol{F}_t + \boldsymbol{e}_t$,

where $\boldsymbol{X}_t = \big[ X_{1t},..., X_{nt} \big]^{\top}$ is an $n-$dimensional vector, $\boldsymbol{F}_t$ is an $r-$dimensional static factor, $\boldsymbol{\Lambda}$ is an $( n \times r )$ factor loading matrix and $e_t = \big[ e_{1t},...., e_{nt} \big]^{\top}$ is an $n-$dimensional idiosyncratic error term. Moreover, $\boldsymbol{G}$ is an $( r \times q )$ matrix of rank $q$, $\boldsymbol{\eta}_t$ is the $q-$dimensional reduced-form shock and $\boldsymbol{\zeta}_t$ is the $q-$dimensional structural shock. We also have that $\mathbb{E} ( \boldsymbol{\eta}_t \boldsymbol{\eta}_t^{\top} ) = \boldsymbol{I}_q$ and $\boldsymbol{A}$ is a $q \times q$ nonsingular matrix. Setting $q \leq r$, then the proposed framework allows to to obtain dynamic factors. Furthermore, by assuming the stationarity of $\boldsymbol{F}_t$, we have that

where $\boldsymbol{\Psi}_j$ is the coefficient matrix in the vector moving average representation. Moreover, denote with $\boldsymbol{\mathcal{F}}_t = \big[ \boldsymbol{F}_{t-1}^{\top},..., \boldsymbol{F}_{t-p}^{\top} \big]^{\top}$ and $\boldsymbol{\Phi} = \big[ \boldsymbol{\Phi}_1,..., \boldsymbol{\Phi}_p \big]$. Rearranging equations gives

where $\boldsymbol{\Pi} = \boldsymbol{\Lambda} \boldsymbol{\Phi}$, $\boldsymbol{\Theta} = \boldsymbol{\Lambda} \boldsymbol{G}$ and $\boldsymbol{\Gamma} = \boldsymbol{\Theta} \boldsymbol{A}$. In other words, for the $i-$th cross-section, the factor model can be represented as $\underline{X}_i = F \lambda_i + \underline{e}_i$, where $\underline{X}_i = [ X_{i1},..., X_{iT} ]^{\top}$ and $F = [ F_1,..., F_T ]^{\top}$ and $\underline{e}_i = [ e_{i1},..., e_{iT} ]^{\top}$. Furthermore, let $\hat{F}_t$ denote the standard PCA estimator for $F_t$, such that $\frac{1}{T} \sum_{t=1}^T \hat{F}_t \hat{F}^{\top}_t = \boldsymbol{I}_r$.

Identification via Conditional Heteroscedasticity

The presence of structural change in a sample can be exploited as an identification mechanism within a structural system context, especially when the heteroscedasticity of structural variances across the two regimes is exploited (see, hodoshima1988estimation, choi2002structural, lanne2008identifying as well as bacchiocchi2011new). In other words, the presence of heteroscedasticity in the error terms $\varepsilon_t$ of the system can be employed for identifiability purposes. Although, the conditional heteroscedastiticy approach often requires to use a pre-testing procedure before estimating the model as in meitz2021testing, lutkepohl2021testing, bertsche2022identification and guay2021identification; we follow the study of bruggemann2016inference who proposed a framework for the identification of SVARs with conditional heteroscedasticity of unknown form.

Let $( u_t, t \in \mathbb{Z} )$ be a $K-$dimensional white noise sequence defined on a probability space $( \Omega, \mathcal{F}, \mathbb{P} )$, such that each $u_t = \left( u_{1t},..., u_{2t} \right)^{\prime}$ is assumed to be measurable with respect to $\mathcal{F}_t$, where $( \mathcal{F}_t )$ is a sequence of increasing $\sigma-$fields of $\mathcal{F}$. Suppose that we observe a data sample $( y_{-p+1},..., y_0, y_1,..., y_T )$ of sample size $T$ plus $p$ pre-sample values from the following DGP for the $K-$dimensional time series $y_t = \left( y_{1t},..., y_{Kt} \right)^{\prime}$ such that

where $A(L) = I_k - A_1 L - A_2 L^2 - ... - A_p L^p, \ \ \ A_p \neq 0$.

Consider a $K-$dimensional time series such that $y_t = \left( y_{1t},..., y_{Kt} \right)$ where

or $A(L) y_t = \mu + u_t$, in compact representation. Denote with $\boldsymbol{y} = \mathsf{vec} \left( y_1,..., y_T \right)$ to be $(KT \times 1)$ vector. Moreover, the parameter $\boldsymbol{\beta}$ is estimated by $ \widehat{\boldsymbol{\beta}} = \mathsf{vec} \left( \widehat{A}_1,..., \widehat{A}_p \right)$, via the multivariate OLS estimator

Assume that the process $\boldsymbol{y}_t$ is stable, then it has a vector moving-average representation (VMA) s.t.

where $\Phi_j, j \in \mathbb{N}$, is a sequence of (exponentially fast decaying) $( K \times K)$ coefficient matrices with

Notice that the standard estimator of $\boldsymbol{\Sigma}_u$ is given by $ \boldsymbol{\Sigma}_u = \frac{1}{T} \sum_{t=1}^T \widehat{u}_t \widehat{u}_t^{\prime}$, where $\hat{u}_t = y_t - \widehat{A}_1 y_{t-1} - ... - \widehat{A}_p y_{t-p}$ are the residuals obtained from the estimated VAR(p) model.

Moreover, we set $\boldsymbol{\sigma} = \mathsf{vech} ( \boldsymbol{\Sigma}_u )$ and $\widehat{\boldsymbol{\sigma}} = \mathsf{vech} ( \widehat{\boldsymbol{\Sigma}}_u )$. The $\mathsf{vech}$-operator is defined to stack column wise the elements on and below the main diagonal of the matrix. A particular useful way to consider the development of the asymptotic theory is to obtain the joint limiting distribution using an unconditional central limit theorem as below

such that the covariance matrix is partitioned as below

with the following analytical forms

Next, we focus on the residual-based moving block bootstrap resampling method. Various studies in the literature have demonstrated that block bootstrap methods are suitable for capturing dependencies in time series data. Specifically, we are interested in applying the moving block bootstrap technique for the residuals obtained from a fitter VAR$(p)$ model to approximate the limiting distribution of

Identification via Non-Gaussianity of Structural Shocks

Generally, exploiting the presence of non-Gaussianity especially in a multivariate system which is more likely to include non-linearities and complex dependencies of unknown forms between variables is an identification method which uniquely determines the structural parameter. However, identification up to signed-permutations reduces the identified set and thus to achieve point identification it is necessary to have a mechanism for selecting a specific permutation as in lanne2017identification (see, also hallin2015r). Recently, chan2023large propose a fast algorithm for identifying the exact permutation scheme based on multiple sign and ranking restrictions.

Assume that the matrix $\boldsymbol{B} \equiv \boldsymbol{S} \boldsymbol{C}$ where $\boldsymbol{C} \in \mathbb{R}^{d \times d}$ is orthogonal and $\boldsymbol{S} \in \mathbb{R}^{d \times d}$ is such that $\boldsymbol{\Sigma}_u = \boldsymbol{S} \boldsymbol{S}^{\prime}$. Thus, given a large number of time series observations $\left\{ \boldsymbol{u}_t \right\}_{t=1,...,n}$, then the statistical problem is to consistently estimate the matrix product $\boldsymbol{S} \boldsymbol{C}$ and determine the distribution of the structural shocks $\boldsymbol{\varepsilon}_t$. Since, the matrix $\boldsymbol{S}$ can be directly estimated by matching from the covariance matrix $\boldsymbol{\Sigma}_u$, using the Cholesky decomposition, then the identification problem of the structural parameter matrix, boils down to the consistent estimation of the orthogonal matrix $\boldsymbol{C} \in \mathbb{R}^{d \times d}$. When structural shocks $\boldsymbol{\epsilon}_t$ are assumed to be Gaussian, one can only identify the space of the mixed signals without separating the individual components. The only way to ensure consistent identification is by assuming that structural shocks are non-Gaussian and then proceeding by consistently estimating the orthogonal matrix $\boldsymbol{C}$. To consistently estimate the individual components of the structural shocks we assume joint non-Gaussianity and that the orthogonal matrix $\boldsymbol{C}$ is consistently estimated. The main intuition of exploiting the non-Gaussianity of structural shocks as an identification mechanism as in lanne2017identification can be seen via the concept of time-reversibility (see, findley1986uniqueness, tong2005time, chan2006note and lawrance1991directionality\footnote{Professor Tony Lawrance is an Emeritus Professor at the Department of Statistics, University of Warwick, UK.}) in time series analysis. The notion of time-series reversibility was proposed by economists and financial engineers to address the asymmetry property of business cycles. The presence of asymmetric signals can induce enough skewness to reject the null hypothesis of time-reversibility. Business cycle asymmetry can be tested with respect to whether macroeconomic fluctuations are time irreversible and can also be exploited for identification purposes as in virolainen2020structural (see, nyberg2023risk).

\paragraph{Proof of Proposition 1 in lanne2017identification}

Suppose that the following two conditions hold:

A key step is to demonstrate the existence and uniqueness of an equivalent SVAR representation. In particular, one needs to show that the structural shocks of the original model can be determined as a function of the structural shocks of the observetionally equivalent SVAR process multiplied by an unknown permutation matrix $\boldsymbol{M}$ such that $ \boldsymbol{\varepsilon}_t = \boldsymbol{M} \boldsymbol{\varepsilon}_t^{*}, \ \ \ \boldsymbol{M} = \boldsymbol{B}^{-1} \boldsymbol{B}^{*}$. To ensure a valid identification scheme, we need to examine the properties of the elements of this permutation matrix. Thus, we assume that at most one column of $\boldsymbol{M}$ may contain more than one nonzero element (see, Lemma A.1 in lanne2017identification and kagan1973characterization\footnote{For the interested reader, the property is due to the theory of locally compact abelian groups which implies that any finitely generated group can be expressed as a direct product of cyclic groups. As a result, the duality result allows to construct a distributional equivalent topological space using properties from non-Gaussian topological spaces. Specifically, any compact Hausdorff abelian group can be completely described by the purely algebraic properties of its dual group. Corollary 3 in morris1979duality defines the observetionally equivalent class of topologically isomorphic groups. These theorems remain agnostic regarding the exact properties model estimators are assumed to have ex-ante to the imposed specification, which implies that using a non-Gaussian distributional assumption is sufficient.} Suppose say that the $k-$th column of $\boldsymbol{M}$ has at least two nonzero elements, $m_{i k}$ and $m_{jk}$ such that $i \neq j$. Then, it holds that

where the random variable $\varepsilon_{k,t}^{*}$ is assumed to be Gaussian with positive variance. Moreover, for all $\ell \in \left\{ 1,..., n; \ell \neq k \right\}$, it must hold that $m_{ik} m_{j \ell} = 0$ because only the $k-$th element of $\boldsymbol{M}$ could have more than one nonzero element. Notice that the purpose of the identification scheme is to distinguish or identify the structural shocks. Due to the independence assumption of the error terms and by the fact that each row of $M$ has exactly one non-zero element, then we can conclude that there exist a permutation matrix $\boldsymbol{P}$ and a diagonal matrix $\boldsymbol{D} = \mathsf{diag} \left\{ d_1,..., d_n \right\}$ with non-zero diagonal elements such that $\boldsymbol{M} = \boldsymbol{D} \boldsymbol{P}$. Rearranging terms implies that $ \boldsymbol{\varepsilon}_t = \boldsymbol{B}^{-1} \boldsymbol{B}^{*} \boldsymbol{\varepsilon}_t^{*}$ and to have observational equivalence it should hold that

since $\boldsymbol{D}$ is a diagonal matrix then $\boldsymbol{D}^{-1}$ is also a diagonal matrix and $\boldsymbol{P} \boldsymbol{P}^{\prime} = \boldsymbol{I}$ where $\boldsymbol{P}$ represents a permutation matrix within the same equivalence class.

Proposition 1 in lanne2017identification characterizes a class of observationally equivalent SVAR processes that differ only with respect to the ordering and scaling of the structural shocks in the vector $\varepsilon_t$. A (rescaled) error vector of these set of permuted structural shocks, consists of exactly the same elements, whose ordering varies and thus these SVAR processes induces asymptotically equivalently impulse responses functions. Thus, the invariance transformation property allows to identify structural shocks using an observationally equivalent class of SVAR processes based on mutual independence.

Testing for Independence of Structural Shocks

The independence property of structural shocks $\varepsilon_t$ can serve to uniquely identify the structural parameter $B$ it at most one component of $\varepsilon_t$ is Gaussian. Moreover, in non-Gaussian frameworks, one has to consider the entire dependence structure to implement a robust structural analysis. Moreover, although the original identification scheme of lanne2017identification is based on the assumption of at least one Gaussian shock for full identification, the particular result is extended by maxand2020identification who show that the same identification scheme is valid in the case of multiple Gaussian structural shocks while the remaining structural shocks are non-Gaussian. In addition some further studies as in crucil2023monetary and bruns2023testing consider the case of non-Gaussian proxy SVARs. In particular, crucil2023monetary extending the approach of schlaak2023monetary employs an instrumental regression to estimate the structural shock of interest while staying agnostic for the remaining structural shocks.

Testing for the independence of the identified structural shocks obtained from orderings (e.g., a set of permutations of all shocks) of the VAR residuals is an important validation step for the underline assumptions. When a suitable test statistic of independence for the identified shocks doesn't find evidence against the hypothesis of independence then this is an indication that using a non-Gaussian MLE approach should render valid model estimates. Motivated from the presence of such non-Gaussian and heavy-tailed error terms in SVAR models davis2023time study the role of uncertainty in economic fluctuations. Specifically, disaster shocks are shown to have short term economic impact arising mostly from feedback dynamics. In particular, the authors show that the coefficient estimate on the infinite variance regressor would be consistent but not asymptotically normal. As a result, the coefficients of impulse response functions whether there are computed from the VAR or by local projections are likely not to be asymptotically normally distributed. A possible drawback of the particular identification scheme is that exploiting the non-Gaussianity property precludes cases of non-identification. Recently, there is a growing interest in testing not only the distributional assumptions of structural shocks but also indirectly testing the absence of non-identification. In particular, drautzburg2023refining mentiones that criteria such as an inverted GMM test statistic based on co-dependence of higher moments, can be employed to rule out incorrect identification schemes when the data are non-Gaussian. When the data are Gaussian, there is no co-dependence of higher moments and the estimator becomes close to the plug-in estimator in moon2012bayesian (see, also granziera2018inference). Thus, the approach of drautzburg2023refining exploits independence for identification, but does so allowing for weak identification, which indeed remains valid in the Gaussian limit where there is complete lack of identification. Further issues related to identification via non-Gaussianity is the modelling of climate-economic systems which is relevant to the structural analysis of climate shocks to macroeconomic variables.

Identification via External Instruments

Overall, SVARs are used in the applied macroeconometrics literature to estimate the dynamic causal effects of structural shocks. Parameters in these models have traditionally been point-identified using zero or long-run restrictions. A commonly used approach in empirical macroeconometrics research is to construct time series that capture exogenous changes affecting the macroeconomy. The dynamic effects of these time series can be often estimated using distributed lag regressions. However, these variables are correlated with a particular structural shock and uncorrelated with other structural shocks (e.g., target versus non-target structural shocks). In particular, external instruments impose second moment restrictions that identify SVAR shocks and associated impulse response coefficients.

A major known issue of SVAR-IV models is the weak identification problem which appears in the case that instruments are only weakly correlated with the variable of interest. Specifically, building on the IV regression literature, stock2021inference proposed weak-instrument robust inference methods for impulse response coefficients. Their framework demonstrates how an external instrument can be used to identify the structural shock of interest, its impulse response coefficients, historical effect on the variables in the VAR and contribution to forecast error variances. A relevant identification scheme combines sign restrictions with external instruments is proposed by braun2023identification.

Identification via Sign-Restrictions and Set-Identified SVARs

Several studies consider the set-identification of SVAR models based on a combination of economic theory and data-driven approaches such as sign restrictions. The key idea behind sign restrictions is to characterise a shock through imposing restrictions on the responses of some variables, while remaining agnostic about others. In particular, inoue2013inference, inoue2016joint study the construction of confidence sets for set-identified parameters. Moreover, giacomini2021robust and giacomini2022robust consider robust inference in Proxy-SVARs with external instruments. Sign restrictions usually imply partial identification (see, dufour2005projection). Additionally, it has several drawbacks including the fact that different structural shocks in the model might be characterized by the same set of sign restrictions. Nevertheless, to illustrate the implementation of the particular identification scheme we follow the framework proposed by gafarov2018delta who contribute to the literature of set-identified SVARs by proposing a novel delta-method interval for the coefficients of the IRF.

Consider the $n-$dimensional SVAR with $p-$lags such that $\boldsymbol{Y}_t = \boldsymbol{A}_1 \boldsymbol{Y}_{t-1} + ... + \boldsymbol{A}_p \boldsymbol{Y}_{t-p} + \boldsymbol{B} \boldsymbol{\varepsilon}_t$. Then, the object of interest is the $k-$th period ahead structural impulse response function of variable $i$ to a particular shock $j$ such that $ \lambda_{k,i,j} (A,B) \equiv e_i^{\prime} \boldsymbol{C}_k (A) \boldsymbol{B}_j$, where $\boldsymbol{B}_j \equiv \boldsymbol{B} \boldsymbol{e}_j$ and $e_i$ and $e_j$ denote the $i-$th and $j-$th column of $\boldsymbol{I}_n$. A second object of interest is the vector of reduced-form VAR parameters

Overall, a common practice in empirical macroeconomics literature is to set equality and inequality restrictions to set-identify the structural IRFs. Specifically, an example of an equality restriction in a monetary VAR is that prices do not react contemporaneously to monetary policy shocks. An example of an inequality restriction is that contractionary monetary policy shock cannot increase prices. Let $\mathcal{R} (\mu) \subset \mathbb{R}^n$ be the set of values of $\boldsymbol{B}_j$ that satisfy the inequality and equality restrictions such that

where $Z(\mu)$ is a matrix of dimension $( n \times m_z )$ that collects the equality restrictions specified by the researcher and $S(\mu)$ is a matrix of dimension $( n \times m_s )$ that collects the inequality restrictions.

Therefore, gafarov2018delta suggested to exploit the specific form of the directional derivative in the SVAR framework via the following expression

\paragraph{Evaluating the posterior of sign-identified VAR Models}

Consider the $d-$variate reduced-form VAR$(p)$ model such that (see, inoue2013inference)

for $t = 1,...,T$, where $e_t \overset{ \textit{i.i.d} }{ \sim } \mathcal{N} ( 0_{n \times 1}, \Sigma )$ and $\Sigma$ is positive definite. We write the above econometric model such that

where $Y = \big[ y_1 \ y_2 \ ... \ y_T \big]^{\prime}_{ T \times k }$, $X = \left[ X_1 \ X_2 \ ... \ X_T \right]$, $X_t = [ 1 \ y_{t-1}^{\prime} \ ... \ y_{t-1}^{\prime} ]$ and $B = [ c \ B_1 \ B_2 \ ... \ B_p ]^{\prime}$. Moreover, we follow the conventional approach of specifying a normal-inverse Wishart prior distribution for the reduced-form VAR parameters.

where $N_0$ is an $np \times np$ positive definite matrix, $S_0$ is an $n \times n$ covariance matrix and $\nu_0 > 0$.

Identification in Proxy SVARs under Possible Nonstationarity

Suppose that $\left\{ Y_t: t = -p+1,...,T \right\}$ be a $( d \times 1 )$ vector of observed variables that follows a SVAR process $ \boldsymbol{Y}_t = \boldsymbol{\mu} + \sum_{j=1}^p \boldsymbol{A}_j \boldsymbol{Y}_{t-j} + \boldsymbol{u}_t, \ \ \ \boldsymbol{u}_t = \boldsymbol{B} \boldsymbol{\varepsilon}_t$, where $\boldsymbol{u}_t$ is the reduced-form error terms and $\boldsymbol{A}_j$ are $( d \times d)$ coefficient matrices for $j \in \left\{ 1,..., p \right\}$. Therefore, to study the IRFs with respect to the structural shocks $\varepsilon_{it}$, we study the estimation of the $i-$th column of the structural parameter $B$ denoted by $\boldsymbol{B} = ( \boldsymbol{b}_1,..., \boldsymbol{b}_d )^{\prime}$ with external instruments $\big\{ \boldsymbol{Z}_{t} = ( z_{1t},..., z_{kt} )^{\prime} \in \mathbb{R}^k : t = 1,..., T \big\}$.

Long-run restrictions are a popular identification method of SVAR models since the seminal contribution of blanchard1988dynamic, which are commonly used to investigate the hours debate\footnote{In particular, the hours debate is concerned with the short-run effect of a technology shock on hours initiated by gali1999technology and christiano2003happens. Using a bivariate SVAR framework $Y_{1t}$ corresponds to log productivity and $Y_{2t}$ corresponds to log hours. Further studies which investigate the impact of technology shocks to business cycles include francis2005technology, feve2009response and chaudourne2014understanding. For example, feve2009response impose a long-run identifying restriction which implies that the composite non-technology shock has no long-run effect on labor productivity. This means that the upper triangular element of $\boldsymbol{A}(L)$ in the long-run must be zero, that is $A_{12} = 0$, so to uncover this restriction from the estimated VAR$(p)$ model, the matrix $\boldsymbol{A}(1)$ is obtained by the Choleski decomposition of $\boldsymbol{C}(1) \boldsymbol{\Sigma} \boldsymbol{C}(1)$.} of gali1999technology as well as the impact of technology shocks on productivity and output (see, also chari2008structural).

Then, it can be shown that the DGP can be reformulated as below

Denote with $\boldsymbol{X}_t = \big[ \left( \boldsymbol{Y}^{*}_{1,t-1} - \boldsymbol{\mu}_1 \right)^{\prime}, \boldsymbol{1}, \boldsymbol{D}_t^{\prime} \big]^{\prime}$, then the original SVAR process (which includes a model intercept) and the above reformulation based on the demeaned original series has an one-to-one transformation. Moreover, this equivalent representation implies that the least-squares residual $\hat{\boldsymbol{u}}_t$ from the VAR process is numerically equivalent to that obtained from regressing $\Delta \boldsymbol{Y}_t$ on $\boldsymbol{X}_t$. Therefore, we shall use the conventional SVAR representation for practical estimation purposes such as to obtain the OLS residuals of the model whereas we use the above reformulation for the asymptotic theory analysis of the corresponding model estimator under the presence of both near-unit and cointegrated regressors.

Next, for the near unit root regressors we denote with $\boldsymbol{J}_c (s)$ an $( k_1 \times 1 )$ vector OU process such that it satisfies the following stochastic differential equation $d \boldsymbol{J}_c (s) = \boldsymbol{C} \boldsymbol{J}_c (s) + d \boldsymbol{B}_u (s)$. Denote with $\boldsymbol{\mathcal{D}}_n$ the diagonal matrix $ \boldsymbol{\mathcal{D}}_n =

$. Moreover, it holds that

which converges into the following stochastic matrix

Therefore, the following converges results holds

\paragraph{Proof of VEC Representation}

Our goal is to demonstrate that the data generating process with three equations representing the case of both near unit roots and cointegrated regressors, is equivalent to the original SVAR formulation.

where the above regressors correspond to their demeaned counterparts. In other words we start by rewritting the regressors specific represenation to the equivalent general SVAR process.

Recall that we have that $\boldsymbol{u}_t := \boldsymbol{P} \boldsymbol{e}_t$ which implies by the change of basis rule

Thus in matrix notation we have that $\Delta \boldsymbol{Y}^{*}_t = \boldsymbol{M} \left( \boldsymbol{Y}^{*}_{t-1} - \boldsymbol{\mu} \right) + \boldsymbol{P} \boldsymbol{\Psi} (L)^{-1} \boldsymbol{e}_t$. Multiplying both sides by $\boldsymbol{P} \boldsymbol{\Psi} (L) \boldsymbol{P}^{-1}$, we obtain

Identification with Occasionally Binding Constraints

At the Zero Lower Bound (e.g., see eggertsson2014can and aruoba2022svars) the ability of policy makers to react to macroeconomic shocks is limited and therefore the response of the economy to shocks will change when the policy instrument hits the ZLB. Thus, according to mavroeidis2021identification since the difference in the behaviour of the macroeconomic variables across the ZLB and non-ZLB regimes is only due to the impact of monetary policy, the switch across regimes provides information about the causal effects of policy. Therefore, the ZLB identifies the causal effect of policy during the unconstrained regime. However, if monetary policy remains partially effective during the ZLB regime, for example, through the use of unconventional policy instruments, then the difference in the behaviour of the economy across regimes will depend on the difference in the effectiveness of conventional and unconventional policies. In this case, we obtain only partial identification of the causal effects of monetary policy, but we can still get information bounds on the relative efficiency of unconventional policy.

From the econometric perspective a major challenge when the restriction mechanism involves OBCs is that it generates censoring of one of the dependent variables. Specifically, once the censoring mechanism is triggered then one can allows for some of the coefficients for the remaining variables to change. According to aruoba2022svars the specification of a dynamic multivariate model with censoring and regime-dependent coefficients faces two challenges: parsimony and the existence of a unique reduced form. These models distinguish between a shadow rate, $y_{1t}^{*}$, and the actual interest rate $y_{1,t} = \mathsf{max} \left\{ y_{1t}^{*}, c \right\}$ by defining the endogenous variable $s_t = \boldsymbol{1} \left\{ y_{1,t} > c \right\}$. This implies that the regime-dependency of the coefficients is equivalent to capturing nonlinearities in decision rules that arise in DSGE models with occasionally-binding constraints. Moreover, duffy2022cointegration mentions that in an autoregressive model, an occasionally binding constraint naturally give rise to nonlinearity in the form of links between unit roots and stochastic trends. A prototypical model from the VECM literature is

On Causal Discovery and Identification

Causal ordering provides a mechanism for causal effect justification in identification and estimation of SVARs. In particular, moneta2011causal studies statistical methodologies for causal inference in SVAR models while other authors propose a graph-theoretic method for causal analysis in vector autoregressions which allows to select the causal (contemporaneous) order of a SVAR model. Thus, a key role to the theory of causal discovery plays the framework for recursive causal models proposed by kiiveri1984recursive. The recursive causal modelling approach utilizes the notion of conditional independence (pioneered by wright1934method and wold1960generalization), for system identification purposes. Moreover, swanson1997impulse discuss the estimation of impulse response functions based on a causal approach to residual orthogonalization in vector autoregressions and propose a search algorithm as statistical mechanism for non-recursive structural models while hausman1983identification discuss the use of IV regression as an identification mechanism for systems of simultaneous equations

In this direction, silva2017learning examine structural learning using instrumental variables under the non-Gaussianity assumption. Their IV causal discovery method provides a set of potential causal effects, from which a unique causal discovery is attainable if at least two instrumental variables (under the same conditioning set) are present in the true underline graph. Moreover, to facilitate computational flexibility this methodology focuses on completeness, that is, on the characterization of relation as causals when instrumental variables satisfy the given Graphical criteria, rather than the identification of all causal effects induced from faithfulness assumptions. A graphical approach to estimating high-dimensional VAR models is presented by fragetta2011effects and bertsche2023directed.

Identifying Aggregate Fluctuations

A particular stream of literature focuses on the identifying aggregated fluctuations through sectoral and firm-specific dynamics. Specifically, gabaix2011granular proposed a framework for investigating the granular origins of aggregate fluctuations. In addition, gabaix2023granular illustrate when the market is sufficiently concentrated, then one can use the collection of idiosyncratic shocks to individual micro units, at each time period, as an instrument for endogenous aggregate variables. Moreover, a full-information approach to Granular Instrumental variables is proposed by baumeister2023full while a framework for heterogeneity-robust granular instruments is presented by qian2023heterogeneity. Lastly, banafti2022inferential establish inferential theory for GIVs in high dimensional settings.

Granular Instrumental Variables

The global market clearing condition is given by $y_{St} = d_t$, where $y_{St} := \boldsymbol{S}^{\prime} \boldsymbol{y}_{\cdot t} = \sum_{i=1}^N S_i y_{it}$, where $\boldsymbol{S}$ is the $( N \times 1)$ vector of shares that are normalized such that $\sum_{i=1}^N S_i = 1$.

Feasible Granular Instrumental Variables

Consider the following panel simultaneous equations model with factor error structure

Asymptotic distribution for structural parameters in factor augmented regressions in time series and panel models is already well-developed in the literature. When both $N$ and $T$ are large, then it can be shown that the sampling error from estimating the high dimensional precision matrix, the factors, as well as the instrument is negligible in the asymptotic distribution of the structural parameters. Moreover, in the case of supply elasticity, the estimator will additionally depend on the estimated (potentially high dimensional) precision matrix. Therefore, the contribution of banafti2022inferential to the GIV literature focuses around relaxing the homogenous loadings assumption which is overly restrictive and propose constructing an estimate of the instrument using PCA or iterative OLS-PCA methods.

Extensions

Extension 1: VARMA Representation

Suppose that the $k-$dimensional process $\left\{ \boldsymbol{Y}_t: t \in \mathbb{Z}\right\}$ has the following VARMA representation

where $\left\{ \varepsilon_t \right\}$ is the innovation process, a sequence of independent and identically distributed random variables with zero-mean and an invertible covariance matrix $\Sigma$ (see, also gourieroux2020identification).

The following representation also holds

(Sketch Proof)

\paragraph{Forecasting in VARMA models}

Modelling and forecasting a large set of macroeconomic variables can be done using VARMA representations. However, to ensure a robust estimation, the identification of uniquely parametrized VARMA models requires imposing restrictions on the parameter space to ensure that the model is uniquely identified (see, dias2018estimation).

Define the lag polynomials such that

We shall say that the model is unique if there is only one pair of stable and invertible polynomials $A(L)$ and $B(L)$, respectively which satisfies the canonical MA representation (see, also poskitt2006identification)

Suppose that our aim is to construct a statistic for testing the null hypothesis $\boldsymbol{\theta} = \boldsymbol{\theta}_0$ against the alternative $\boldsymbol{\theta} \neq \boldsymbol{\theta}_0$. Choosing $p_0 < p_1$ and $q_0 < q_1$ allows one to test the adequency of the specified VARMA coefficients in $\boldsymbol{\theta}_0$, while contemplating the possibility of possibly higher-order VARMA models.

Therefore, this null hypothesis is thus invariant under the group of affine transformation $\boldsymbol{\varepsilon}_t \mapsto \boldsymbol{M} \boldsymbol{\varepsilon}_t$ if and only if $\boldsymbol{M} \boldsymbol{A}_i \boldsymbol{M}^{-1} = \boldsymbol{A}_i$ for all $i = 1,...p_0$ and $\boldsymbol{M} \boldsymbol{B}_j \boldsymbol{M}^{-1} = \boldsymbol{B}_j$ for all $j = 1,..., q_0$, that is, iff each $\boldsymbol{A}_i$ and each $\boldsymbol{B}_j$ commutes with any invertible matrix $\boldsymbol{M}$, which holds true iff they are proportional to the $( k \times k )$ identify matrix (see, hallin2004rank).

Write $\mathcal{H}^{(n)} ( \boldsymbol{\theta}_0, \boldsymbol{\Sigma}, f )$ for the hypothesis under which an observation $\boldsymbol{X}^{(n)} := \big( \boldsymbol{X}_1^{(n)},..., \boldsymbol{X}_n^{(n)} \big)^{\prime}$ is generated by the $\mathsf{VARMA}(p_0, q_0)$ model with parameter value $\boldsymbol{\theta}_0$ satisfying Assumption A and innovation process satisfying Assumption B. Our objective is to test $\mathcal{H}^{(n)} ( \boldsymbol{\theta}_0 ) := \bigcup_{\boldsymbol{\Sigma} } \bigcup_{ f } \mathcal{H}^{(n)} ( \boldsymbol{\theta}_0, \boldsymbol{\Sigma}, f )$ against $\bigcup_{ \boldsymbol{\theta} \neq \boldsymbol{\theta}_0 } \mathcal{H}^{(n)} ( \boldsymbol{\theta} )$. Consequently, $\boldsymbol{\Sigma}$ and $f$ play the role of nuisance parameters; note that the unions, in the definition of $\mathcal{H}^{(n)} ( \boldsymbol{\theta}_0 )$, are taken over all possible values of $\boldsymbol{\Sigma}$ and $f$.

Let $\boldsymbol{A}(L)$ and $\boldsymbol{B}(L)$ be such that $\boldsymbol{A}_i = \boldsymbol{0}$ for $i = p_0 + 1,..., p_1$, and $\boldsymbol{B}_i = \boldsymbol{0}$ for $i = q_0 + 1,..., q_1$, and consider the sequences of linear difference operators

where the vector

is such that $\mathsf{sup}_n ( \boldsymbol{\uptau}^{(n)} )^{\prime} \boldsymbol{\uptau}^{(n)} < \infty$. These operators define a sequence of VARMA models such that

hence, the sequence of local alternatives $\mathcal{H}^{(n)} \left( \boldsymbol{\theta}_0 + n^{-1/2} \boldsymbol{\uptau}^{(n)}, \boldsymbol{\Sigma}, f \right)$.

\paragraph{Testing for Invertibility}

A key assumption for valid identification, estimation and forecasting with VARMA models (see, gourieroux2020identification, sims2012news and lutkepohl2002forecasting), is the invertability of underline polynomial representations. In particular, chen2017testing propose a test for invertability or fundamentalness of structural vector autoregressive moving average models generated by non-Gaussian independent and identically distributed structural shocks. Moreover, it can be proved that in these models under certain regularity conditions the Wold innovations are a martingale difference sequence if and only if the structural shocks are fundamental. Thus, the particular representation provides a mechanism for testing the presence of invertability. In other words, chen2017testing convert the statistical problem of polynomial invertability of VARMA processes into a statistical testing problem for the martingale difference sequence property of the Wold innovations. An MLE approach for non-invertible ARMA models is studied by meitz2013maximum. Further issues on identification and estimation of non-invertible SVARMA models are presented by funovits2020identifiability.

Extension 2: Non-Linear Dynamic System

A more challenging task is the econometric identification in nonlinear systems. In addition, for a non-linear dynamic system excited by a non-Gaussian disturbances, the system equation governing the probability distribution of the system response contains an infinite number of terms which makes identification and estimation considerable complex. Therefore, a pertrubation scheme is proposed to obtain an approximate solution to the system, thereby the effects of excitation non-Gaussianity can be evaluated (see, cai1992response, grigoriu1995linear and sengupta1989efficient).

Consider the unit-root VAR model, $\boldsymbol{A} (z) \boldsymbol{y}_t = \boldsymbol{\varepsilon}_t$, which crucially depends on $\boldsymbol{A} (z)^{-1} (L)$, thus we address the issue of of the inversion of $\boldsymbol{A} (z)$ around a pole $z = z_0$. Denote with

be a matrix polynomial of order $n$ and degree $K$ and $z_0$ such that $\mathsf{deg} \big[ \boldsymbol{A} (z) \big] = 0$. Therefore,

Then, an analytical expression for the matrix function $\boldsymbol{A} (z)^{-1}$ is given by

where $m$ is the order of the pole and the first term above represents the principal part while $\boldsymbol{M}(z) = \sum_{ j = 0 }^{ + \infty } \boldsymbol{N}^{j} ( z - z_0 )^j$, represents the regular part. Using the fact that $\boldsymbol{A} (z) \boldsymbol{A} (z)^{-1} = \boldsymbol{I}_n$

Notice that the Unit-Root VAR model is a linear non-homogeneous difference-equation system in matrix form whose solution can be formally written as below

where the first term is a particular solution of the non-homogeneous equation and the second is the so-called complementary solution. However, both depend on the operator $\boldsymbol{A}^{-1} (L)$, and thus on the matrix $\boldsymbol{A}^{-1} (z)$ since the algebras of the polynomial functions of the lag operator $L$ and of the complex variable $z$ are isomorphic. These expressions are useful when deriving statistical tests for invertability.

Extension 3: Automatic Inference for Finite Order VARs

Selecting lag length is useful when $p$ is unknown, which can directly impact the estimation of the effects of structural shocks. Notice that we only discuss lag selection procedures for stationary settings\footnote{A relevant framework for model selection in explosive time series regressions is presented by tao2020model. }. Roughly speaking the selection of the optimal lag length is achieved via the use of information criteria such as the AIC, SIC and HQC (e.g., see ivanov2005practitioner and choi2012model).

In this section, we focus on the framework proposed by kuersteiner2005automatic who consider automatic inference for infinite order vector autoregressions. In particular, infinite order dynamics in cointegration models has been examined by several studies. Specifically, $y_t$ has an infinite order VAR representation

where $\mu = C(1)^{-1} \mu_y$ and $C(L)^{-1} = \xi (L)$ such that $\xi (L) = \displaystyle \left( I - \sum_{j=1}^{\infty} \xi_j L^j \right)$. Therefore, the approximate model with VAR coefficient matrices $\left\{ \xi_{1,p},..., \xi_{p,p} \right\}$ has the following formulation

where $\Sigma_{v,p} = \mathbb{E} \big[ v_{t,p} v_{t,p}^{\prime} \big]$ the mean squared prediction error of the approximating model. Moreover, it can be shown that the parameters $( \xi_{1,p},..., \xi_{p,p} )$ are root$-n$ consistent and asymptotically normal for $\xi (p) = ( \xi_1,..., \xi_p )$ under the condition that the lag order $p$ does not increase too quickly, that is, $p$ is chosen such that $p^3 / n \to 0$.

Consider the sequential testing procedure proposed by ng2001lag, then a Wald test for the null hypothesis that the coefficients of the last lag $p$ are jointly zero, is expressed by

Extension 4: Locally Trimmed OLS Procedure for VARs

A particular methodology which has been proposed in the literature to deal with the problem of outliers when estimating vector autoregression models, is the least-trimmed squares procedure. In particular, agullo2008multivariate proposed a multivariate least-trimmed squares estimator and prove Fisher-consistency under the assumption of elliptically symmetric error distributions. Moreover, cavaliere2009robust consider robust methods for estimation and unit root testing in autoregressions with infrequent outliers whose number, size, and location can be random and unknown\footnote{In particular, the authors show that standard inferene based on ordinary least squares estimation of an augmented Dickey-Fuller (ADF) regression may not be reliable, because (a) clusters of outliers may lead to inconsistent estimation of the autoregressive parameters, and (b) large outliers induce a jump component in the asymptotic distribution of UR test statistics.}.

Let $\left\{ y_t | t \in \mathbb{Z} \right\}$ be a $d-$dimensional stationary time series. The vector autoregressive model of order $p$, denoted by VAR$(p)$, is given by the following expression

with $y_t$ a $d-$dimensional vector, the intercept parameter $\mathcal{B}_0$ a vector in $\mathbb{R}^d$ and the slope parameters $\left\{ \mathcal{B}_1, ...., \mathcal{B}_p \right\}$ are matrices in $\mathbb{R}^{d \times d}$. Throughout the paper, $B^{\prime}$ will stand for the transpose of the matrix $B$. Moreover, the $d-$dimensional error terms $\epsilon_t$ are assumed to be independently and identically distributed with a density function of the form

with $\Sigma$ a positive definite matrix, called the scatter matrix, and $g$ a positive function. If the second moment of $\epsilon_t$ exists, $\Sigma$ will be proportional to the covariance matrix of the error terms. Existence of a second moment, will not be required for the robust estimator. We focus on the unrestricted VAR$(p)$ model, where no restrictions are put on the parameters $\left\{ \mathcal{B}_1, ...., \mathcal{B}_p \right\}$. Suppose that the multivariate time series $y_t$ is observed for $t=1,...,T$. Then, the vector autoregression model has the following multivariate regression representation

for $t \in \left\{ p+1,..., T \right\}$ with $x_t = \left( 1, y_{t-1}^{\prime}, ..., y_{t-p}^{\prime} \right)^{\prime} \in \mathbb{R}^q$, where $q = pd + 1$. Moreover, the matrix $\mathcal{B} = \left( \mathcal{B}_1, ...., \mathcal{B}_p \right)^{\prime} \in \mathbb{R}^{q \times d}$ contains all the unknown matrices of coefficients. In matrix notation, we have that $X = \left( x_{p+1},...,x_{T} \right)^{\prime} \in \mathbb{R}^{ n \times q}$ is the matrix of regressors and $Y = \left( y_{p+1},..., y_{T} \right)^{\prime}$ is the matrix of predictants, where $n = T - p$ is the sample size after excluding the estimation for the lag components of the model. The classical least squares estimator for the regression parameter $\mathcal{B}$ is estimated by

and the scatter matrix $\Sigma$ is estimated by the following expression

\paragraph{The Multivariate Least Trimmed Squares Estimator}

The unknown parameters of the VAR$(p)$ will be estimated via the multivariate regression model. To allow for the robust estimation of the model under the presence of outliers the Multivariate Least Trimmed Squares (MLTS) estimator is employed, which is based on the idea of the Minimum Covariance Determinant estimator. The way that the MLTS estimator works is that it selects the subset of $h$ observations which minimize the determinant of the covariance matrix of residuals based on a least squares fit from the particular subsample.

Consider the dataset $Z = \left\{ \left( x_t, y_t \right), t = p+1,..., T \right\} \subset \mathbb{R}^{d+q}$. We denote with $\mathcal{H}$ the collection of all subsets of size $h$, such that:

Then, for any subset $H \in \mathcal{H}$, we denote with $\widehat{\mathcal{B}}_{\text{OLS}}$ the classical least squares estimator based on the observations of the particular subset given by the following expression

where $X_H$ and $Y_H$ are submatrices of $X$ and $Y$ respectively, consisting of the rows of $X$ and $Y$ that corresponds to the subsample $H$.

Then, the corresponding scatter matrix estimator computed from this subset is constructed as below

The MLTS estimator is now defined as below

Moreover, the associated estimator of the scatter matrix of the error terms is given by

where $c_{\alpha}$ is the correction factor to obtain consistent estimation of $\Sigma$ for the error distribution of the model, and $\alpha$ the trimming proportion for the MLTS estimator, that is, $\alpha \approx \left( 1 - \frac{h}{n} \right)$. In the case of multivariate error terms it has been shown that $ c_{\alpha} = (1 - \alpha ) \big/ F_{ \chi^2_{d+2} \left( q_{\alpha} \right)}$, where $F_{\chi^2_{d+2}}$ is the cumulative distribution function of a $\chi^2$ distribution with $q$ degrees of freedom, and $q_{\alpha} = \chi^2_{q, 1 - \alpha}$ is the upper $\alpha$ quantile of this distribution.

Equivalent characterizations of the MLTS estimator are given by agullo2008multivariate. The particular study shows that any $\widetilde{\mathcal{B}} \in \mathbb{R}^{d \times q}$ minimizing the sum of the $h$ smallest squared Mahlanobis distance of its residuals, subject to the condition that $\text{det} \left[ \Sigma \right] = 1$ is a solution of (ref). More formally,

where $d_{1:n} \left( \mathcal{B}, \Sigma \right),...,d_{n:n} \left( \mathcal{B}, \Sigma \right)$ is the ordered sequence of the residual Mahalanobis distances is

for some $\mathcal{B} \in \mathbb{R}^{d \times q}$. Notice that the MLTS estimator minimizes the sum of the $h$ smallest squared distances of its residuals, and is therefore the multivariate extension of the Least Trimmed Squares (LTS) estimator of rousseeuw1984least. The particular MLTS estimator has been found to have low efficiency so the literature has proposed the reweighted version of the estimator to improve the performance of the MLTS estimator. Therefore, since the efficiency of the MLTS estimator is quite low, the re-weighted version of the particular estimator is employed to improve the performance of the MLTS.

The Reweighted Multivariate Least Trimmed Squares (RMLTS) estimates are defined as

and $q_{\delta} = \chi^2_{q, 1 - \delta}$.

Since outliers tend to have large residuals with respect to the initial robust MLTS estimator, then it implies a large residual Mahalanobis distance $d_j^2 \left( \mathcal{ \widehat{ \mathcal{B}} }_{\text{MLTS}}, \widehat{\Sigma}_{\text{MLTS}} \right)$. If the later is above the critical value $q_{\delta}$, then the observation is flagged as an outlier. The final RMLTS is then based on those observations not having been detected as outliers. For instance, one can set $\delta = 0.01$ and take as trimming proportion for the initial MLTS estimator $\alpha = 25 \%$.

\paragraph{Determining the autoregressive order}

To select the order $p$ of a vector autoregressive model, information criteria are computed for several values of $p$ and an optimal order is selected by minimizing the criterion. Most information criteria are in terms of the value of the log likelihood $\ell_k$ of the VAR$(p)$ model. Using the model assumption for the distribution of the error terms, we obtain that

with $n = T - k$. When error terms are multivariate normal the above leads to the following expression

The log likelihood will depend on the autoregressive order via the estimate of the covariance matrix of the residuals. In particular, for the classical least squares estimator we have that

where $\widehat{\epsilon}^{\prime}_t (k)$ are the residuals corresponding to the VAR$(p)$ model. Using trace properties, the last term of the log-likelihood function corresponds to $(n - p) p / 2$ for the OLS estimator. In order to prevent that outliers might affect the optimal selection of the information criteria, we estimate $\Sigma$ by the RMLTS estimator given by the following expression:

with $J(k)$ as above and $m(k)$ the number of elements in $J(k)$. In this case, the last term of the log-likelihood function equals now $\left( m(k) - p \right) p \big/ \left( 2 c_{\delta} \right)$.

\paragraph{Open Problems} Some possible extensions of the particular approach is to propose a framework which employes the LTS method that accounts for outliers in multivariate time series models, specifically to the setting of structural vector autoregressions.

Extension 5: Time-Varying VAR Models

The identification and estimation of time-varying VAR models is useful when modelling dynamically macroeconomic events such as monetary policy regimes (see, primiceri2005time, koop2013large) and chan2016large as well as time-varying instrumental variable estimation as in kapetanios2019large, petrova2019quasi and giraitis2021time. A relevant setting is the framework of factor augmented SVAR models (e.g., see eickmeier2015classical).

Financial Connectedness

The study of network topology and financial connectedness is used for modelling systemic risk, spillover effects and financial contagion using the VAR framework (see, discussion in katsouris2023limit). Estimating systemic risk measures implies identifying shocks which are transmitted via system-wide connectedness. In particular, diebold2014network, diebold2012better proposed a framework for modelling the network topology via the forecast error variance decomposition (FEVD) of an estimated VAR process, as a transformation method that summarizes connectedness between a set of institutions. Thus this approach allows to construct measures of directional connectedness which are analogous to bilateral imports and exports for each of a set of $N$ countries. Recent implementations are presented by barigozzi2017network, barunik2018measuring, demirer2018estimating, korobilis2018measuring, geraci2018measuring, rehman2018precious, yoon2019network, barigozzi2021time, laborda2021volatility, mensi2023time and barunik2024persistence.

Consider, an $N$ dimensional covariance-stationary series $\boldsymbol{Y_t}=(y_{1,t},...,y_{N,t})$, which is generated from a VAR$(p)$ process for $t=1,...,T$. Then, based on the Wold decomposition theorem we have that

such that $|\Phi(z)|$ lie outside the unit-circle, then the VAR process has the MA($\infty$) representation expressed as $\boldsymbol{Y_t}= \boldsymbol{\Theta(L)}\boldsymbol{\epsilon_t}$, where $\boldsymbol{\Theta(L)}$ is an $(N x N)$ matrix of infinite lag polynomials and $\boldsymbol{\epsilon}_t$ is considered to be a white-noise generated by the covariance matrix $\boldsymbol{\Sigma}$. Therefore, the $H-$step generalized variance decomposition\footnote{See also, lanne2016generalized who employ a GFEV decomposition for linear and nonlinear multivariate models. Moreover, katsouris2023estimating develops methods for the identification and estimation of financial networks.} matrix D$^{gH}$ = $[d_{ij}^{gH}]$ is given by

where $e_j$ are selection matrices e.g., $e_j=(0,...1,...0)^{'}$, $\boldsymbol{\Theta}_h$ is the coefficient matrix evaluated at the $h$-lagged shock vector and $\sigma_{ij}$ is the $j$th diagonal element of $\boldsymbol{\Sigma}$. Modelling the asymmetric responses of the economy to monetary policy shocks over the business cycle due to the impact of the monetary transmission mechanism (e.g., see eickmeier2015classical, morley2012asymmetric and veldkamp2005slow), requires to employ specifications which allow for time-varying effects (see, paul2020time).

Dynamic Causal Effects

The two main types of identification presented so far in the literature consists of point-identification and set-identification. Specifically, point-identified SVARs are interpreted as achieving identification through internal instruments. Consequently, in these models, structural shocks are the interventions of interest, and thus the goal is to estimate the dynamic causal effect of these shocks on macroeconomic outcomes. In contrast, modern microeconometric identification stategies rely heavily on external sources of variation that provide quasi-experiemnts to identify causal effects. Such external variation might be found, for example, in firm-specific characteristics that introduce as-if randomness in the variable of interest (the treatment). The use of such external instruments in microeconometrics has proven highly productive and has yielded compelling estimates of causal effects. Therefore, the use of external instruments has opened a new and rapidly growing research avenue in macroeconometrics, in which credible identification is obtained using as-if random variation in the shock of interest that is distinct from - external to - the macroeconomic shocks hitting the economy (stock2018identification).

Identification and Estimation of Causal Effects

Causal Effects and IV Regression

Following the paper of stock2018identification, a starting point for formulating the theory on identification and estimation of dynamic causal effects in macroeconomics, is that the expected difference in outcomes between the treatment and control groups in a randomized experiment with a binary treatment is the average treatment effect. Roughly speaking, if a binary treatment $X$ is randomly assigned, then all other determinants of $Y$ are independent of $X$, which implies that the (average) treatment effect is

In the linear model $Y = \alpha + \beta X + u$, where $\beta$ is the treatment effect, random assignment implies that $\mathbb{E} \left( u | X \right) = 0$ so that the population regression coefficient is the treatment effect. When randomization is conditional on covariates $W$, then the treatment effect for an individual with covariates $W = w$ is estimated by the outcome of a random experiment on a group of subjects with the same value of $W$,

Therefore the identification and estimation of dynamic treatment effects allows to investigate the VAT behaviour of firms in relation to imposed by the government tax rates. The intuition goes as follows. Governmental taxation is the pillar of fiscal policy and decisions on the optimal tax policy are usually taken in relation to their effectiveness on various macroeconomic outcomes, such as output fluctuations as well as structural changes in agent's behaviour over the business cycle (aggregate fluctuations). From the econometrics perspective to facilitate inference requires to derive the asymptotic behaviour of the dynamic treatment effect estimator across a panel of economics which have undergone a period of fiscal adaption to new tax policies. In other words, one would be interested to obtain the asymptotic behaviour of these estimators based on the trajectory of the control against the experimental groups in relation to a common taxation policy as well as to obtain counterfactual outcomes when a control group has adapted an updated fiscal policy across the business cycle. In this direction, we are interested to incorporate tax policy variables such as taxation revenues which are often considered as policy outcomes. Furthermore, due to the fact that these variables might be positive correlated during the business cycle, suitable instrumental variables that can ensure consistent model estimation and inference procedures can be employed. In particular, vegh2015tax use the VAT rates across economic during the business cycle as a valid instrumental variable. Moreover, the authors use a novel tax database which combines both corporate and personal income tax rates as well as data on value-added taxes. A relevant study which discusses the effects of global tax reform is presented by gomez2023measuring.

Counterfactual Analysis

Causal Effect Inference with Time Series Data

In particular, the parameter of interest is the ATE $:= \mathbb{E} \left[ Y_1 - Y_0 \right]$. Moreover, the Conditional Average Treatment Effect (CATE) defined by

captures how the average treatment effect changes across sub-populations with covariate value $X = x$. Furthermore, under the unconfoundness assumption $\tau(x)$ can be nonparametrically identified by

Consider the $k-$step ahead forecast of $Y_t$ given the information available up to time $t^{*}$ which is given by the following expression

Notice that by definition, $Y_{ t^{*} + k } (0)$ is the potential outcome expected at time $t^{*} + k$ in case the intervention does not occur. Therefore, an estimator of the point causal effect can be defined as

Potential Outcome Framework

Let $D \in \left\{ 0, 1 \right\}$ be a binary treatment indicator and $Y_d$ be the potential outcomes for $d = 0,1$ in the status quo environment (see, also hsu2022counterfactual). In particular, $Y_1$ is the outcome if an individual is exogenously assigned to the treatement $(D = 1)$ and $Y_0$ is the outcome in the absence of treatement $( D = 0)$. Then, the actual observed outcome is

We observe a $d-$dimensional vector of pretreatment covariates $X = \left( X_1,..., X_k \right)$ in the status quo environment and $X^{\star} = \left( X^{\star}_1,..., X^{\star}_k \right)$ in the counterfactual environment which is of the same dimension as the set of covariates $X$.

The corresponding treatment indicator, outcome, and potential outcomes in the counterfactual environment are denoted as $D^{*}, Y^{*}$ and $Y_d^{*}$ for $d = 0,1$, respectively with

However, notice that since the treatment has not been implemented in the counterfactual environment yet, neither $D^{*}$ nor $Y^{*}$ is observed in our model (see, hsu2022counterfactual). We adopt the former approach such that a closed-form solution exists

where $\widehat{ \mathbb{E}} \left( Y_d | X = X_j^{*} \right)$ is the Nadaraya-Watson estimator such that

The asymptotic properties for $\widehat{\delta}^{*}$ can be derived under weaker regularity conditions stated below. Furthermore based on regularity conditions it can be proved that

Therefore, compared to the semiparametric efficiency bound of the ATE estimator we have that

Time Series Experiments and Causal Effects

Estimating Contemporaneous Effect of Treatment

Time series causal effects are defined as a comparison between the potential outcomes at a fixed point in time. Thus, the primary object of interest is the temporal average of these causal effects. Therefore, the only source of randomness in our formulation is the randomization of the treatment. Once we begin the experiment, the randomization reveals a particular outcome path; however, this does not change the set of potential outcomes, and therefore, does not alter what we would have observed had we administered a different treatment. Thus, in this sense, we are treating the potential outcomes as predetermined, and the only source of randomness comes from which path we observe.

Then, the temporal average treatment effect of the paths $w_{1:T}$ and $w_{1:T}^{\prime}$ is

Moreover, we define the potential outcomes just intervening on treatment the last $j$ periods as

In other words, this "marginal" potential outcome represents the potential or counterfactual level of regime A in country $i$, if we let welfare spending run its natural course up to $t-j-1$ and just set the last $j$ lags of spending to $x_{t-j:t}$. Thus, we can define an important quantity of interest, that is, the contemporaneous effect of treatment (CET) of $X_{it}$ on $Y_{it}$ such that:

Furthermore, following bojinov2019time specifically in longitudinal studies and survey studies with a follow-up wave of cross-sectional datasets, we are particularly interested to investigate the treatment effect of survey participants as well as their membership in certain endogenously generated groups. This implies that we are modelling the presence of time-varying exposure to dynamic causal estimands. Therefore, to construct a feasible structural econometric framework we consider the treatment and potential outcome path with time-varying effects. At each time step, $t = 1,...,T$, we assume that each individual (sampling unit), is exposed to either treatment, $W_t = 1$, or control, $W_t = 0$, and then we measure the corresponding outcome variable of interest.

Although we focus on binary treatments, our results can be generalized to multiple treatments as well. Thus, the random "treatment path" is given by $W_{1:t} = \big( W_1,..., W_t \big)$. As a result, the potential path for the treatment path $w_{1:t}$ can corresponds to the following sequence

Null Hypothesis of Temporal Causal Effects

Following the framework proposed by bojinov2019time, to access whether the treatment has a statistically significant effect, we consider the sharp null of no temporal causal effects

In particular, this will null hypothesis will be tested against a portmanteu alternative. Thus, invoking the sharp null hypothesis implies that $Y_t ( w_{1:t}^{\text{obs}} ) = Y_t ( w_{1:t}^{\prime} )$ for all $w_{1:t}^{\prime}$. Moreover, the particular formulation of the null hypothesis means "the null of no temporal causality at lag $p \geq 0$ of the treatment on the outcome" holds, such that $H_{0,p}: \tau_{t,p} = 0$, for all $t = 1,2,...,T$.

In the case where $p = 0$, then $\widehat{\tau}_{t,0}$ estimation error can be written as below

and under the sharp null hypothesis we obtain that

In particular, whenever $p_t(1) = p_t(0)$ the estimator of the upper bound to the variance for the contemporaneous causal effect, is equal to the variance under the sharp null of no treatment effect. Then, the $p$ lagged error variance $Var^R \big( u_{t-p,p} | \mathcal{F}_{T, t-p-1} \big)$ equals to

We begin by rewriting the within-unit matching estimator as below

Consider the conditional independence property which implies that $\big\{ Y_{it}(1), Y_{it}(0) \big\} \perp \boldsymbol{X}_i | \boldsymbol{U}_i$. By the law of iterated expectations this implies that

for $x = 0,1$.

Bias-Reduced Doubly Robust Estimation

According to vermeulen2015bias the proposed approach may more generally lend itself better to small-sample inference. For instance, suppose that interest lies in the marginal causal effect $\tau = \mathbb{E} [ Y(1) - Y(0) ]$. Because the proposed estimation strategy does not require acknowledging the uncertainty of the estimated nuisance parameters (up to first order), we foresee that it may potentially lend itself better to randomization inference. How such randomization inference could be accomplished and how it performs in small to large samples is indeed relevant questions for further investigation.

Marginal Treatment Effects

Consider $\textit{i.i.d}$ data $\left\{ \boldsymbol{Z}_i = \left( Y_i, A_i, \boldsymbol{X}_i \right), i = 1,..., n \right\}$, where $Y_i$ is the outcome of interest, $A_i$ is a dichotomous treatment taking values zero and one and $\boldsymbol{X}_i$ is a sufficient set of covariates to control for confounding of the treatment effect, in the sence that $Y( a ) \perp A| \boldsymbol{X}$ for $a \in \left\{ 0, 1\right\}$. Specifically, $Y(a)$ denotes the counterfactual outcome for treatment level $a \in \left\{ 0, 1\right\}$, which is linked to the observed data through the consistency assumption (i.e., $Y(a) = Y$ iff $A=a$).

Notice that the selection on observables assumption says that policies are independent of potential outcomes after appropriate conditioning. Moreover, the conditioning set includes the unknown error vector of the model, $\varepsilon_t$, such that $Y_{t,1}(d), Y_{t,2}(d),..., \perp \varepsilon_t | z_t$. This is because, conditional on the set of instruments $z_t$, randomness in $D_t$ is due exclusively due to randomness induced from the unknown $\varepsilon_t$ vector. Therefore, in order to convert the above assumption into a testable hypothesis which has an identifiable parameter space, we consider the sharp null hypothesis such that $Y_{t,j} ( d^{\prime} ) = Y_{t,j} (d) = Y_{t+j}$.

Then, if we substitute the observed with the potential outcomes we obtain the following testable conditional independence assumption:

Notice that $D_t$ is not necessarily a binary treatment status but can be a set of different taxation policies which implies that we have a continuum of treatment status.

The following theorem of Sims and Granger causality holds.

The Augmented Synthetic Control Method

The synthetic control method (SCM) is a popular approach for estimating the impact of a treatment on a single unit in panel data settings (see, ben2021augmented). In other words is a weighted average of control units that balances the treated units pre-treatment outcomes and other covariates as closely as possible. We define the treated potential outcome as $Y_{it} = Y_{it} (0) + \tau_{it}$, where the treatment effects $\tau_{it}$ are fixed parameters. Since the first unit is treated, the key estimand of interest is

Then, the observed outcomes are as below

Furthermore, we assume that the assignmemnt process is sequentially probabilistic, meaning that any assignment vector may be realized with positive probability at time $t$ given the history of the observable stochastic process up to time $t-1$. Let $\left\{ \mathcal{F}_t \right\}_{ t \geq 1}$ denote the natural filtration generated by the realized $\left\{ w_t, y_t \right\}_{ t \geq 1}$.

Panel Data Regression for Dynamic Causal Effects

Panel Experiments and Dynamic Causal Effects

Following the framework proposed by bojinov2021panel we consider a design-based framework which provides a generalization of the finite population literature in cross-sectional causal inference and time series experiments to panel experiments.

Identifying Dynamic Causal Effects with Panel Data: Applications

Advanced Topics

Testing for Unstable Root in Structural ECMs

We follow the framework proposed by boswijk1994testing. Consider the single-equation error correction model of a time series $\left\{ y_t \right\}$ conditional upon the $( k \times 1 )$ vector time series $\left\{ z_t \right\}$, for $t \in \left\{ 1,... T \right\}$ as

where $\left\{ v_t \right\}$ is an innovation process relative to $\left\{ z_t, y_{t-j}, z_{t-j}, j = 1,2,... \right\}$ with positive variance $\omega^2$. Notice that $\theta$ and $\beta_j$ are $( k \times 1)$ parameter vectors such that $j \in \left\{ 1,..., p-1 \right\}$, such that $\theta$ defines the long-run equilibrium relation $y = \theta^{\prime} z$, the deviations from which lead to a correction of $y_t$ by a proportion of $\lambda$, the adjustment or error correction coefficient.

The conditional model is said to be stable if all roots of the characteristic equation

are outside the unit circle. In other words, stability of the model implies that the disequilibrium error $\left( y_t - \theta^{\prime} z_t \right)$ is a stationary process, even though $z_t$ and $y_t$ are integrated of order one and hence nonstationary. Specifically, if the model is stable, then $x_t = ( y_t, z_t^{\prime} )^{\prime}$ is cointegrated with cointegrating vector $( 1, - \theta )^{\prime}$. Thus, the purpose of this exercise is to develop a class of tests for the null hypothesis that the characteristic equation has a unit root, so that the model is unstable, against the alternative hypothesis of stability. Furthermore, the single-equation conditional model can be seen as a special case if a structural error correction model. This is a system of Error Correction Equations for a $( g \times 1 )$ vector of time series $\left\{ y_t \right\}$ conditional upon $\left\{ z_t \right\}$ such that (see, boswijk1994testing):

where $\left\{ \boldsymbol{v}_t \right\}$ is an innovation process with a positive-definite covariance matrix $\boldsymbol{\Omega}$ such that the above expression corresponds to a parametrization of a conditional model of $y_t$ given $z_t$ (if $\boldsymbol{\Gamma}_0 \neq \boldsymbol{I}_g)$. Next we consider the identification of these matrix parameters. In particular, the above model implies $g$ cointegrating relationships $\boldsymbol{\Gamma} \boldsymbol{y} + \boldsymbol{B} \boldsymbol{z} = \boldsymbol{0}$, provided that it is stable, such that the characteristic equation is expressed as below

has all roots outside the unit circle.

In other words, unless the parameters are restricted in some way, the structural model is not identified. Specifically, we identify the long-run relations by imposing restrictions of the usual form such that

where $\boldsymbol{\Gamma}_i$ and $\boldsymbol{B}_i$ denote the $i-$th row of $\boldsymbol{\Gamma}$ and $\boldsymbol{B}$, respectively, and where $\boldsymbol{R}_i$ is a known matrix of approximate order. Thus, the rank condition for identification of the $i-$th long-run relation

Then, the remaining parameters are identified by the normalization $\boldsymbol{\Gamma}_{0,ii} = 1$, and the restriction that $\boldsymbol{\Lambda} = \mathsf{diag} \big( \lambda_1,..., \lambda_g \big)$, that is, where $\boldsymbol{\Lambda}$ is a diagonal matrix. In practice, this means that only the disequilibrium error of the $i-$th long-run relation appears in the $i-$th structural error correction equation which means that the $i-$th equation is (over)-identified if the $i-$th long-run equation is.

Consequently, under this null hypothesis, there is no error correction in the $i-$th equation, which suggests that the $i-$th row of the system $\boldsymbol{\Gamma} y + \boldsymbol{B} \boldsymbol{z} = 0$ is not a cointegrating relationship. Specifically, let $z_t$ be generated by the following expression:

Assumption 1, can be interpreted as a form of exogeneity condition, because it states that the cointegration properties of the conditional model carry over to the full VAR system. Understanding well these concepts are essential especially if one is interested to apply advanced identification techniques in SVARs such as the Markov switching regime (see, lanne2002threshold and lanne2010structural) as well as identification in Gaussian Mixture Vector Autoregressive Models (see, kalliovirta2016gaussian and meitz2023mixture). In particular, kalliovirta2016gaussian develop a framework for Gaussian mixture vector autoregression which is a nonlinear vector autoregression and is designed for analyzing time series that exhibit regime-switching dynamics\footnote{The property of long memory and regime switching is examined in the study of diebold2001long while a regime switching panel data regression model is proposed by cheng2019regime. For a forecasting application see nyberg2018forecasting.} (see, kalliovirta2015gaussian and virolainen2020structural). Next, we consider another example motivated from the structural break literature specifically in the context of cointegrated models (see, also hansen2002testing).

Testing Cointegration Rank when some Cointegrating Directions are Changing

Consider the following VECM of order $p$ with a structural break such that

Then, the date of the break is characterized by the fraction $\lambda_0$ of the sample size $n$ such that $k_0 = \floor{\lambda_0}$, where $\lambda_0 \in [ \underline{\lambda}, \bar{\lambda} ]$ and $0 < \underline{\lambda} < \bar{\lambda} < 1$. Moreover, the VECM is assumed to satisfy the stability condition over the two regimes as given by the Assumption below.

Empirical Likelihood Estimation Approach

Based on the aforementioned assumptions the estimation of the unknown parameter vector $\boldsymbol{\beta}_o$ based on the OLS estimator $\hat{\boldsymbol{\beta}}$ is given by

where $\boldsymbol{\Lambda} = \boldsymbol{\Omega}_1^{-1} \boldsymbol{\Omega}_2 \boldsymbol{\Omega}_1^{-1}$, are defined as $(p+1) \times (p+1)$ matrices.

To construct an empirical likelihood function, the estimation equations are defined as below:

for a generic parameter $\boldsymbol{b} \in \mathbb{R}^{p+1}$. By using the Lagrange multipliers method, we have that $\hat{\boldsymbol{\lambda}} = \hat{\boldsymbol{\lambda}} ( \boldsymbol{b} ) \in \mathbb{R}^{p+1}$ is the solution of the following set of equations:

Then, the corresponding empirical log-likelihood ratio is given by

and it holds that $\ell ( \boldsymbol{\beta}_0 ) \overset{d}{\to} \chi_p^2$, as $T \to \infty$.

High Dimensional Generalized Empirical Likelihood Estimation

In this section, we discuss relevant aspects to the empirical likelihood estimation for high dimensional dependent data, which is applicable to time series regression models (see, chang2015high). Further frameworks related to identification and estimation of high-dimensional time series with SVAR models are proposed by krampe2023structural, zhang2023statistical and adamek2022local. In particular, denote with $\boldsymbol{\theta} = ( \theta_1,..., \theta_p )^{\prime}$ be a $p-$dimensional parameter taking values in a parameter space $\Theta$. Consider a sequence of $r-$dimensional estimating equation such that

for some $r \geq p$, Then, the model information regarding the data and the data parameter is summarized by moment restrictions below:

where $\theta_0 \in \Theta$ is the true parameter. Furthermore, in order to preserve the dependence structure among the underlying data, we employ the blocking technique. Let $M$ and $L$ be two integers denoting the block length and separation between adjacent blocks, respectively. Then, the total number of blocks is

Then, the EL estimator $\boldsymbol{\theta}_o$ is $\widehat{\boldsymbol{\theta}}_{EL} = \mathsf{arg max}_{ \theta in \Theta } \ \mathsf{log} \ \mathcal{L} ( \boldsymbol{\theta} )$. Consequently, the maximization problem can be carried out more efficiently by solving the corresponding dual problem, which implies that $\widehat{\boldsymbol{\theta}}_{EL}$ can be obtained as below:

for any $\boldsymbol{\theta} \in Q$ and $\mathcal{N}$ an open interval containing zero.

High Dimensional VARs with Common factors

We follow the framework proposed by miao2023high who study high-dimensional vector autoregressions (VARs) augmented with common factors that allow for strong cross-sectional dependence. This approach allows to incorporate in a unified framework a convinient mechanism for accommodating the interconnectedness and temporal co-variability that are often present in large dimensional systems.

Consider the $N-$dimensional vector-valued time series $\left\{ Y_t \right\} = \left\{ \left( y_{1t},..., y_{NT} \right)^{\prime} \right\}$, the high-dimensional VAR model of order $p$ with CFs given by

where $\boldsymbol{A}_1^0,..., \boldsymbol{A}_p^0$ are the $( N \times N )$ transition matrices and $\boldsymbol{u}_t$ is an $N-$dimensional vector of unobserved idiosyncratic errors. Moreover, the analytical framework allows for both the number of cross-sectional units $N$ and the number of time periods $T$ to pass to infinity. The lag length is also allowed to (slowly) grow to infinity with $(N,T)$. Estimation then is a natural high-dimensional problem.

Then, the $N-$dimensional VAR$(p)$ process $\left\{ Y_t \right\}$ can be rewritten in a companion form as an $Np-$dimensional VAR$(1)$ process with common factors such that:

As a result, the reverse characteristic polynomial of $Y_t$ can be written as below:

In the low-dimensional framework, the process is stationary if $\mathcal{A}(z)$ has no roots in and on the complex unit circle, or equivalently the largest modules of the eigenvalues of $\boldsymbol{\Phi}$ is less than 1. Therefore, to achieve identification, we shall study the Gram or signal matrix $\boldsymbol{S}_X = \boldsymbol{X}^{\prime} \boldsymbol{X}/ T$ and its population counterpart $\Sigma_X = \mathbb{E} \left( \boldsymbol{X}^{\prime}_t \boldsymbol{X}_t \right)$. In other words, one can study the deviation bounds for the Gram matrix, under the Gaussianity assumption and boundedness of the spectral density function.

In order to ensure that the matrix $\Sigma_X$ is well-behaved, we write $\boldsymbol{X}_{t+1}$ as a moving average process of infinite order MA$( \infty )$ such that

Eigenvalue Analysis:

Fast Algorithm for Detection of Breaks in Large VAR models

Currently in the literature there are two approaches in terms of the identification of break-points. Firstly, there is a so-called "top down" procedure, in the sense that one tests all the data to determine if there is at least one change-point and iterates the procedure in the intervals immediately to the "left" and "right" of the most recently detected change-point. However, the particular methodology has certain weaknesses compared to a number of other thresholding procedures. A second approach is the so-called "bottom up" procedure motivated by the observation in the presence of multiple change-points or to mitigate the effects of inadequately controlled drift in the "baseline" mean value, it may be useful to compare a candidate change-point at $j$ to an appropriate "local" background $(j,k)$, where $i < j < k$. Similar approaches are the Wild Binary Segmentation (WBS) methodology which uses a random set of possible backgrounds and an apparently determined threshold (see, fang2020segmentation).

Illustrative Examples

The formulation of the test statistic presented in the above example, implies that is the (smallest) time index $k \geq n_0 := \floor{nr}$, for some small $r > 0$, at which the maximum in the definition of the test statistic $C_n$ is attained. On the other hand, in order to robustify these CUSUM-type test statistics for detecting changes not only in the middle of the sample, the following statistic is proposed

which considers the maximal deviation from the average of the cumulated sums, the maximum being taken over all possible subsamples $\boldsymbol{Y}_{ n \ell}, i \leq \ell \leq j$, for $1 \leq i < j \leq n$. To be more precise, the particular algorithm implies that each subsample induces a corresponding CUSUM statistic and then the maximal test statistic is taken over all these subsamples. Therefore, it can be proved

Further studies related to testing for structural breaks in covariance structures include kao2018testing. In particular, $\widehat{\boldsymbol{\Sigma}}$ be the sample covariance matrix such that $\widehat{\boldsymbol{\Sigma}} = \frac{1}{T} \sum_{t=1}^T y_t y_t^{\prime}$. Then, for a given $\uptau \in [0,1]$, we define a point in time $\floor{ T \uptau }$, and we use the subscripts $\uptau$ and $1 - \uptau$ to denote quantities calculated using the subsamples $t = 1,..., \floor{ T \uptau }$ and $t = \floor{ T \uptau } + 1,...,T$, respectively. In particular, we consider the sequence of partial sample estimators

We denote with $w_t = \mathsf{vec} \big( y_t y_t^{\prime} \big)$ and $\bar{w}_t = \mathsf{vec} \big( y_t y_t^{\prime} - \boldsymbol{\Sigma} \big)$.

Main Asymptotic theory results

In particular, if no serial dependence is present, a possible choice is the full sample estimator

Alternatively, one could use the sequence of partial sample estimators such that

Furthermore, to accommodate for the case in which $\Psi_{\ell} \equiv \mathbb{E} \big( \bar{w}_t \bar{w}_{t- \ell} \big) \neq 0$ for some $\ell$, we propose a weighted sum-of-covariance estimator with bandwidth $m$, such that

Alternatively, we can employ the following matrix

A similar expression holds for $\widehat{\Psi}_{\ell, 1 - \uptau}$. Furthermore, denote with

where $\Omega_T$ is positive definite uniformly in $T$, and $\Omega_T \to \Omega$ with $\left\lVert\Omega\right\rVert < \infty$.

Furthermore, a growing literature develops methodologies for testing structural breaks in factor models as in breitung2011testing, baltagi2021estimating and bai2024likelihood among others (see, also bai1998estimating). Regarding the literature relevant to structural break testing around the Vector Autoregressive framework of interest there are two main approaches. On the one hand, the approach of safikhani2022fast implies a fast and scalable algorithm for detection of structural breaks in large VAR models. The particular modelling environment corresponds to a piecewise linear Vector Autoregression model. Thus, testing for structural breaks takes the form of a time-varying VAR model and therefore in practise the statistical problem of interest is to test for the presence of time-varying vector autoregressive matrix coefficients across the blocks. Moreover, since the model is parametrized in a form of a large VAR model with a corresponding companion matrix, then the focus is on the fast and accurate break detection algorithm with respect to the dimensionality of the model. On the other hand, gourieroux2020identification proposed an econometric framework for identification and estimation of the structural vector autoregressive model where the main focus is how the distributional assumptions on the shocks affect the identification procedure. In practice, based on the non-Gaussianity assumption of the structural shocks the model is identifiable and the estimation of model parameters is achieved via the use of the likelihood function. Therefore, testing for structural breaks in such a setting could be implemented using an appropriate specification which captures time-variation.

Large VAR Models: Model Formulation

For the remaining of this section, we follow the framework proposed by safikhani2022fast. Their procedure is found to be scalable to big high-dimensional time series datasets with a computational complexity that can achieve $O ( \sqrt{n} )$, where $n$ is the length of the time series (sample size), compared to an exhaustive procedure that requires $O(n)$ steps.

Suppose that there exist $m_0$ break points $0 < t_1 < ... t_{m_0} < T$, with $t_0 = 0$ and $t_{ m_0 + 1 } = T$ where $j \in \left\{ 1,..., m_0 + 1\right\}$. Notice that if there are no break points, then $m_0 = 0$ and so $( t_0 = 1, t_1 = T )$ which corresponds to the full sample period.

Thus, we can assume that $m_0 \geq 1$, which implies that $( t_0 = 1, t_2 = T )$ and so $j = \left\{ 1, 2 \right\}$. In other words, for $t_{j-1} \leq t < t_j$ (e.g., $t_0 \leq t < t_1$, i.e., $t = 1,..., T$) we have the following model:

Since in the case that $m_0 = 1$ then $j \in \left\{ 1,2\right\}$ then we have a two-regime model estimation:

There is a single break-point and we fit a different VAR model before the break, i.e. $M_0$ pre-break based on the observations up to the unknown break-point, and a we fit $M_1$ post-break based on the observations from the unknown break-point and onward until the end of the full sample (see, safikhani2022fast).

Based on the aforementioned setting proposed by safikhani2022fast, the unknown parameters in each segment are: the number of break-points $m_0$, their locations $t_j, j = 1,..., m_0$, the VAR parameters $\boldsymbol{\Phi}_{(q, j)}$ together with the covariance matrix. Therefore, according to , the statistical problem is then to detect the break points $t_j$, in a computational efficient manner that is also scalable for very large values of $T$. Additionally, we are also interested to estimate accurately the VAR parameters $\boldsymbol{\Phi}_{( \ell, j)}$, under a high-dimensional environment $( p >> T )$.

A Block Segmentation Scheme Based Algorithm

The main idea of the BSS is to partition the time points into blocks of size $b_n$ and fix the VAR parameters within each block. To this end, define a sequence of time points $q = r_0 < r_1 < ... < r_{ k_n } = T + 1$ which play the role of the end points for the blocks. In other words, the total number of blocks is calculated as:

where $k_n$ is the total number of blocks based on the (effective) sample size and $b_n$ is the length size of each block which is kept fixed across all blocks. Notice that throughout, the dimensions of the vector $\boldsymbol{y}_t$ remains fixed having $p$ elements as well as the (effective) sample size being $n$.

The theoretical results presented in safikhani2022fast deal with the high-dimensional case such that the sparsity levels $d_{ k_j }$ increase with the sample size, $T$. Specifically, we define with $p \equiv p(n)$ and $m_0 \equiv m_0 (n)$ and $d_{ k_j } \equiv d_{ k_j } (n)$, where $n = T - q + 1$. In other words, $n$ is the effective sample size or the length of the time series for estimation of model parameters purposes, since the number of lags used for constructing additional regressors affects the sample size (see, safikhani2022fast).

The minimum distance between two consecutive break points is denoted by

Let $\boldsymbol{\mathcal{Y}} \in \mathbb{R}^{ n \times p}$ and $\boldsymbol{\mathcal{X}} \in \mathbb{R}^{ n \times k_n p q }$ and $\boldsymbol{\theta} \in \mathbb{R}^{k_n p q \times p }$ and $\boldsymbol{E} \in \mathbb{R}^{ n \times p }$. In other words, based on this parametrization, $\theta_i \neq 0$ for $i \geq 2$ implies a change in the VAR coefficients. Therefore, for $j \in \left\{ 1,..., m_0 \right\}$, the structural break points $t_j$ can be estimated as block-end time point $r_{i-1}$, where $i \geq 2$ and $\theta_i \neq 0$.

Then, the model is formulated as below

where $\pi_b = k_n p^2 q$ and $\boldsymbol{\Theta} = \mathsf{vec} ( \boldsymbol{\theta} ) \in \mathbb{R}^{ \pi_b \times 1 }$. The initial estimate of the parameter $\boldsymbol{\Theta}$ is obtained:

Notice that the above high dimensional problem (convexity is preserved) uses fused lasso penalty with two $\ell_1$ penalties controlling the number of break points and the sparsity of the VAR model.

Denote the set of indices of blocks with nonzero jumps and corresponding estimated change points obtained from solving the high-dimensional problem as below:

The total number of estimated change points in this step corresponds to the cardinality of the set $\mathcal{A}_n$, where $\hat{m} = \left| \mathcal{A}_n \right|$. Moreover, the block size $b_n$ acts as a tuning parameter that regulates the number of model parameters to be estimated. In addition to this, the algorithm applies a local screening step. Since the set $\mathcal{A}_n$ of candidate change points tends to overestimates their numbers, a screening step to eliminate redundant break points is required. According to safikhani2022fast, the main idea goes is to estimate the VAR parameters locally on the left and right side of each selected break point in the first step and compare them to one VAR parameter estimated from combining the left and right of the selected break point as one large stationary segment. Now, if the selected break point is close to a true break point, the sum of squared errors calculated assuming stationarity around the true break point will be much larger compared to the sum of squared errors calculated from two separate VAR parameter estimates on the right and left of the selected break point. Therefore, we can get consistent estimates of the number of break points by minimizing a localized information criterion (LIC) comprising of the sum of squared errors and a penalty term on the number of break points.

Next, we describe the localized screening step in more details. Recall that $\mathcal{A} = \left\{ \hat{t}_1, \hat{t}_2,..., \hat{t}_{ \hat{m} } \right\}$ is the set of candidate break points selected in the first step. Then, for each subset $A \subset \mathcal{A}_n$, we define the following local VAR parameter estimates for some $\widehat{t}_i \in A$ such that

where $a_n$ corresponds to the neighbourhood size in which the VAR parameters are estimated. The procedure of safikhani2022fast, can be directly implemented for identifying and dating the presence of multiple break-points in multivariate time series. In contrast to the literature of SVAR models, the assumption is that the model parameters are identifiable but a possible parameter instability could exist.

Consistency of the BSS Estimator

To establish the consistency properties of the BSS-based estimator relevant regularity conditions are required which we briefly present below.

Sieve Bootstrap for Functional Vector Autoregressions

According to paparoditis2018sieve, bootstrap procedures for Hilbert space-valued time series proposed so far in the literature, are mainly attempts to adapt, to the infinite dimensional functional framework, of bootstrap methods that have been developed for the finite dimensional time series case. Specifically, paparoditis2018sieve shows, that under quite general assumptions, the stochastic process of Fourier coefficients obeys a so-called vector autoregressive representation and this representation plays a key role in developing a bootstrap procedure for the functional time series at hand. To capture the essential driving functional parts of the underlying infinite dimensional process, the first $m$ functional principal components are used and the corresponding $m-$dimensional time series of Fourier coefficients is bootstrapped using a $p$th order vector autoregression fitted to the vector time series of sample Fourier coefficients. In this way, a $m-$dimensional pseudo-time series of Fourier coefficients is generated which imitates the temporal dependence structure of the vector time series of sample Fourier coefficients. Thus, using the truncated Karhunen-Loeve expansion, these pseudo-coefficients are then transformed to functional bootstrap replicates of the main driving, principal components, of the observed functional time series.

The bootstrap procedure

The main building block of the procedure proposed by paparoditis2018sieve is to generate pseudo-replicates $X_1^{*},...,X_n^{*}$ of the functional time series at hand by first bootstrapping the $m-$dimensional time series of Fourier coefficients $\xi_t = ( \xi_{1,t}, \xi_{2,t},..., \xi_{m,t} )$, $t =1,2,...n$, corresponding to the first $m$ principal components. This $m-$dimensional time series of Fourier coefficients is bootstrapped using the autoregressive representation $\xi_t$. Then, the generated $m-$dimensional pseudo-time series of Fourier coefficients is then transformed to functional principal pseudo-components by means of truncated Karhunen-Loeve expansion $\sum_{j=1}^n \xi_{j,t} v_j$. Adding to this, an appropriately resampled functional noise leads to the functional pseudo-time series $X_1^{*}, X_2^{*} ,..., X_n^{*}$. However, since the $\xi_t'$s are not observed, we work with the time series of estimates scores. The functional sieve bootstrap algorithm of paparoditis2018sieve is:

Furthermore, an important issue for bootstrap-based inference methods is to investigate the validity of the functional sieve bootstrap applied in order to approximate the distribution of some statistic $T_n = T \left( X_1,...,X_n \right)$ of interest, when the bootstrap analogue $T_n^{*} = T \left( X_1^{*},...,X_n^{*} \right)$ is used. Notice that establishing validity of a bootstrap procedure for time series heavily depends on two issues. On the dependence structure of the underlying process which affects the distribution of the statistic of interest and on the capability of the bootstrap procedure used to mimic appropriately this dependence structure.

Recall the definition given by paparoditis2018sieve such that

and observe that $\xi^{*}_t = \sum_{l=0}^{\infty} \widehat{\Psi}_{l,p}(m) e^{*}_{t-l}$, where $\widehat{\Psi}_{0,p}(m) = I_m$ and the power series

converges for $|z| \leq 1$. Thus, following paparoditis2018sieve we can decompose each sequence as below

where for each $t \in \mathbb{Z}, \left\{ e_{s,t}^{*} , s \in \mathbb{Z} \right\}$ is an independent copy of $\left\{ e_{s}^{*} , s \in \mathbb{Z} \right\}$. Notice that,

By Minskowski's inequality, we have that

Thus, evaluating the first expectation term, we obtain that using $\left\lVertA\right\rVert_F^2 = tr \left( A A^{\top} \right)$ and the submultiplicative property of the Frobenius matrix norm, that

where $\widehat{\Sigma}_{e,p}(m) = \widehat{\Sigma}_{e,p}^{1 / 2}(m) \widehat{\Sigma}_{e,p}^{1 / 2}(m)$. An identical expression appears for the second expectation term on the right-hand side of (ref). Applying Minkowski's inequality again, we get by the exponential decay of $\left\lVert \widehat{\Psi}_{l,p} (m) \right\rVert_F$, such that the following probability bound holds (see, paparoditis2018sieve)

Conclusion

There are various statistical identification procedures for SVAR models such as imposing sign-restrictions, short-run restrictions, long-run restrictions, identification by higher-order moments as well as identification by conditional heteroscedasticity and changes in volatility (e.g., see, lanne2008identifying, lanne2008statistical. More precisely, by exploiting the heteroscedasticity of the errors $\varepsilon_t$, then the identification of the structural parameters of the system can be achieved. In particular, lanne2010structural assume Markov switching and a smooth transition in the covariance matrix of the error term $\varepsilon_t$ of the model. Specifically, the structural analysis of VARs in this case is based on the Markov regime switching property to identify shocks, which holds due to the fact that the reduced form error covariance matrix varies across states. Furthermore, a different stream of literature considers identification in SVAR models by assuming that the error terms are non-Gaussian and mutually independent (e.g., see lanne2017identification and gourieroux2020identification). For example, lanne2010structural assume that the errors of the model are independent over time with a distribution that is a mixture of two Gaussian distributions with zero means and diagonal covariance matrices, one of which is an identify matrix and the other one has positive diagonal elements, which for identifiability have to be distinct. Roughly speaking there are specific systems where the use of the non-Gaussianity property is a more suitable identification strategy of structural shocks of the SVAR system which we leave for discussion in a follow-up paper.

In this note we have reviewed key statistical properties of VAR processes, cointegrated VAR processes and SVAR processes as well as suitable identification strategies for the structural parameters of these structural econometric specifications (see, also mumtaz2018vars). Moreover, we have discusses several empirical-driven applications from the applied macroeconometrics literature and we studied the proposed identification and estimation methods employed to those cases. Then, motivated from the literature of dynamic causal effects (see, mertens2013dynamic, stock2018identification and budnik2023identifying), we have motivated the main statistical and econometric tools for identification of structural econometric models, commonly used for counterfactual analysis purposes, and then presented some examples of time series regression models, both important for understanding suitable techniques for policy-making evaluation in macroeconometrics. We leave any further presentation of other relevant topics for further research discussions. Lastly, some advanced topics were also presented such as break testing for unstable roots in structural error correction models, estimation techniques for the optimal lag order, the estimation of a high dimensional VAR model as well as a fast algorithmic procedure presented in the literature for detecting structural breaks in large VAR models.

In the Appendix of this set of lecture notes we present an example of the maximum likelihood estimation for $\alpha-$stable autoregressive processes (see, andrews2009maximum) as well as the main asymptotic theory results for inference for the VEC(1) model from guo2023inference.