Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,647 characters · 15 sections · 81 citation commands
Linear structural vector autoregressive (SVAR) models have become a standard tool in empirical macroeconomics, as they can effectively capture constant dynamic relationships among the variables and facilitate tracing out the causal effects of economic shocks. Macroeconomic systems may, however, exhibit variation in their dynamics, induced by crises, policy shifts, or business cycle fluctuations, for example. Hence, also the effects of the shocks may depend on the state of the economy. Such features cannot be accommodated by linear SVAR models, and therefore nonlinear SVAR models that are able to capture such features are often employed Kilian+Lutkepohl:2017. In nonlinear SVAR models, the structural shocks are typically identified with conventional methods that impose economically interpretable but restrictive assumptions, such as zero contemporaneous interactions among some of the variables, which are not always plausible in practice. To address this issue, statistical identification methods relying on the statistical properties of the data can be used Kilian+Lutkepohl:2017.
There are two main branches in the statistical identification literature: identification by heteroskedasticity (Rigobon:2003, Lanne+Lutkepohl:2010, Bacchiocchi+Fanelli:2015, Lutkepohl+Netsunajev:2017, Lewis:2021, Virolainen:2025, and others) and identification by non-Gaussianity (Lanne+Meitz+Saikkonen:2017, Lanne+Luoto:2021, Lanne+Liu+Luoto:2023, and others). To the best of our knowledge, this paper is the first to examine identification by non-Gaussianity in nonlinear SVAR models. Under certain statistical conditions, both types of statistical information typically suffice to identify the shocks in a linear SVAR model. As pointed out by Lutkepohl+Netsunajev:2017, among others, identification by heteroskedasticity without additional restrictions, however, has the major drawback in nonlinear SVAR models that it imposes time-invariant (relative) impact effects of the shocks. We show that identification by non-Gaussianity, in turn, facilitates recovering the shocks without such undesirable restrictions.
This paper contributes to the literature on identification by non-Gaussianity by extending the framework of Lanne+Meitz+Saikkonen:2017 to structural smooth transition vector autoregressive (STVAR) models, which is a major class of nonlinear SVAR models Hubrich+Terasvirta:2013. The STVAR model can flexibly capture nonlinear data generating dynamics by accommodating multiple regimes and gradual as well as abrupt shifts between them, governed by the transition weights. In contrast to its linear counterpart, the impact matrix of the structural STVAR model should generally allow for time-variation in the impact responses of the variables to the shocks, which complicates identification. Nevertheless, similarly to Lanne+Meitz+Saikkonen:2017, it turns out that identification is achieved when the shocks are mutually independent and at most one of them is Gaussian. We show that under this condition, the shocks are readily identified up to ordering and signs when the impact matrix of the STVAR model is defined as a weighted sum of the impact matrices of the regimes. While we focus on exogenous, logistic, and threshold transition weights, our results can be extended to other suitable weight functions as well.
In line with the statistical identification literature, external information is required to label the identified structural shocks as economic shocks. Our nonlinear setup has the additional complication that the same shock should be assigned to the same column of the impact matrix across all regimes. Moreover, our experience shows that identification is often weak in the sense that there are often multiple local solutions with a fit close to the global solution. We find such local solutions to be often largely quite similar to each other, with the differences frequently but not exclusively related to different ordering and signs of the columns of the regime-specific impact matrices. Nonetheless, since different local solutions may produce different results in structural analysis, we recommend addressing weak identification by adopting a blended identification strategy that combines identification by non-Gaussianity with supplementary identifying information Carriero+Marcellino+Tornese:2024.
Our methods are demonstrated in an empirical application to the macroeconomic effects of climate policy uncertainty (shocks) that considers monthly U.S. data from 1987:4 to 2024:12. Following Khalil+Strobel:2023 and Huang+Punzi:2024, we measure climate policy uncertainty (CPU) with the CPU index of Gavriilidis:2021, which is constructed based on the amount of newspaper coverage on topics related to CPU. We are interested in studying how the effects of the CPU shock vary depending on the level of economic policy uncertainty (EPU). Therefore, we fit a two-regime structural logistic STVAR model using the first lag of the EPU index Baker+Bloom+Davis:2016 as the switching variable. We find that a positive CPU shock reduces production and raises inflation in times of both low and high EPU, but its effects, particularly on inflation, are stronger in the periods of high EPU. Our results are, hence, in line with the previous literature suggesting that a positive CPU shock reduces production and raises inflation (Khalil+Strobel:2023 and Huang+Punzi:2024), while Fried+Novan+Peterman:2022 found that it decreases production.
The rest of this paper is organized as follows. Section (ref) presents the framework of reduced form STVAR models. Section (ref) discusses identification of the shocks in structural STVAR models, presents our identification results, and discusses the problem of labeling the shocks as well as how to address weak identification by adopting a blended identification strategy. Section (ref) discusses stationarity of the model and proposes estimating its parameters with a penalized likelihood-based estimator using a three-step procedure. Section (ref) presents the empirical application and Section (ref) concludes. Further details can be found in the appendices, and the introduced methods have been implemented to the accompanying R package sstvars sstvars, which is available via the CRAN repository.
Let $y_t$, $t=1,2,...$, be the $d$-dimensional time series of interest and $\mathcal{F}_{t-1}$ denote the $\sigma$-algebra generated by the random vectors $\lbrace y_{t-j}, j>0 \rbrace$ (i.e., $\mathcal{F}_{t-1}$ contains the information about the history of the process). We consider STVAR models with $M$ regimes and autoregressive order $p$ assumed to satisfy
where $\phi_{1},...,\phi_{M}\in\mathbb{R}^{d}$ are the intercept parameters; $A_{1,i},...,A_{M,i}\in\mathbb{R}^{d\times d}$, $i=1,...,p$, are the autoregression matrices; $u_t$ is a martingale difference sequence of reduced form innovations; $\Omega_{y,t}$ is the positive definite conditional covariance matrix of $u_t$ (conditional on $\mathcal{F}_{t-1}$), which may depend on $\alpha_{1,t},...,\alpha_{M,t}$ and $y_{t-1},...,y_{t-p}$; and $\nu$ collects any other parameters that the distribution of $u_t$ depends on.
The transition weights $\alpha_{m,t}$ are assumed to be either exogenous (nonrandom) or $\mathcal{F}_{t-1}$-measurable functions of $\lbrace y_{t-j}, j=1,...,p \rbrace$, and to satisfy $\sum_{m=1}^{M}\alpha_{m,t}=1$ at all $t$. Thus, we accommodate the types of transition weights typically used in macroeconomic applications, including exogenous, logistic, and threshold weights discussed below. The transition weights express the proportions of the regimes the process is in at each point of time and determine how the process shifts between them.
It is easy to see that, conditional on $\mathcal{F}_{t-1}$, the conditional mean of the above-described process is $\mu_{y,t} \equiv E[y_t|\mathcal{F}_{t-1}] = \sum_{m=1}^M \alpha_{m,t}\mu_{m,t}$, a weighted sum the regime-specific means $\mu_{m,t}$ with the weights given by the transition weights $\alpha_{m,t}$. The discussion on the specification of the conditional covariance matrix $\Omega_{y,t}$ is postponed to Section (ref). The linear vector autoregressive (VAR) model is obtained as a special case by assuming constant autoregressive dynamics across the regimes, i.e., $A_{m,i}=A_i$ and $\phi_{m}=\phi$ for all $i=1,...,p$ and $m=1,...,M$.\footnote{See Hubrich+Terasvirta:2013 for a survey on STVAR literature, including also specifications more general than ours.}
A popular specification is obtained by assuming two regimes ($M=2$) and logistic transition weights. In the logistic STVAR model Anderson+Vahid:1998, the transition weights vary according to a logistic function as
where $c\in\mathbb{R}$ is the location parameter, $\gamma > 0$ is the scale parameter, and $z_t$ is the switching variable, which we assume to be a lagged endogenous variable up to the lag $p$ (that is, $z_t\in\{y_{it-j}, i=1,...,d, j=1,...,p \}$). The location parameter $c$ determines the midpoint of the transition function, i.e., the value of the switching variable when the weights are equal. The scale parameter $\gamma$, in turn, determines the smoothness of the transitions (the smaller $\gamma$ is, the smoother the transition is), and it is assumed strictly positive so that $\alpha_{2,t}$ is increasing in $z_t$. In the special case when $\gamma\rightarrow \infty$, regime-switches are discrete and the logistic weights reduce to the threshold weights of Tsay:1998 (see Appendix (ref) for details).
A structural STVAR model is obtained from the reduced form model defined in Section (ref) by identifying the serially and mutually uncorrelated structural shocks $e_{t}=(e_{1t},...,e_{dt})$ $(d\times 1)$ from the reduced form innovations $u_t$. Specifically, the structural shocks are recovered from the reduced form innovations with the transformation
where $B_{y,t}$ is an invertible ($d\times d$) impact matrix that governs the contemporaneous relationships of the shocks and may depend on $\alpha_{1,t},...,\alpha_{M,t}$ and $y_{t-1},...,y_{t-p}$. In other words, assuming a unit variance normalization for the structural shocks, the identification problem amounts to finding an impact matrix $B_{y,t}$ such that the conditional covariance matrix $\Omega_{y,t}=B_{y,t}B_{y,t}'$. However, $B_{y,t}$ is not generally uniquely identified without further restrictions, as many such decompositions exist. Various identification methods have been proposed to resolve this issue Kilian+Lutkepohl:2017.
Lanne+Meitz+Saikkonen:2017, among others, have shown that if the shocks are mutually independent and at most one of them is Gaussian, a static impact matrix is readily identified (up to ordering and signs of its columns) without additional restrictions. We extend this result to a time-varying impact matrix, and therefore, it is useful to parametrize the structural model directly with $B_{y,t}$. Hence, we make the identity $u_t=B_{y,t}e_t$ explicit in Equation ((ref)) as
where $e_t$ are independent and identically distributed structural errors with identity covariance matrix and a distribution that may depend on the parameter $\nu$.
To incorporate time-variation in the impact matrix, we specify its functional form. A natural specification for structural STVAR models is to assume a constant impact matrix for each of the regimes and define the impact matrix of the process to be their weighted sum as
where $B_1,...,B_M$ are invertible $(d\times d)$ impact matrices of the regimes. Also the impact matrix $B_{y,t}$ needs to be invertible to ensure positive definiteness of the conditional covariance matrix, which does not automatically follow from the invertibility of $B_1,...,B_M$. Nevertheless, it turns out that $B_{y,t}$ defined in ((ref)) is invertible for all $t$ almost everywhere in $[B_1:...:B_M] \in\mathbb{R}^{d\times dM}$, as is stated in the following lemma (which is proven in Appendix (ref)).
The result of Lemma (ref) holds almost everywhere in $[B_1:...:B_M] \in\mathbb{R}^{d\times dM}$, which means that the set where the matrices $B_1,...,B_M$ are such that $B_{y,t}$ is singular for some $t$ has Lebesgue measure zero. In other words, invertible matrices $B_1,...,B_M$ generally lead to an invertible impact matrix $B_{y,t}$, excluding some special cases of $B_1,...,B_M$ for which this result does not hold (e.g., if $B_2=-B_1$, then $B_{y,t}=0$ when $\alpha_{1,t}=\alpha_{2,t}=0.5$). To support intuition, note that the result of Lemma (ref) implies that if the matrices $B_1,...,B_M$ are drawn by random from a continuous distribution, the probability of obtaining matrices such that $B_{y,t}$ is singular for some $t$ is zero.
Assuming an impact matrix of the form ((ref)), the conditional covariance matrix of $y_t$, conditional on $\mathcal{F}_{t-1}$, is obtained as
where $\Omega_m\equiv B_mB_m'$ and $\Omega_{m,n}\equiv B_mB_n'$. The conditional covariance matrix of $y_t$ can thus be described as a weighted sum of $M^2$ matrices with the weights varying in time according to the transition weights $\alpha_{m,t}$, $m=1,...,M$. If the process is completely in one of the regimes, i.e., $\alpha_{m,t}=1$ for some $m\in\lbrace 1,...,M\rbrace$, then $\Omega_{y,t}$ reduces to $\Omega_{m}$, implying that this constitutes the conditional covariance matrix of Regime $m$. When the process is not completely in any of the regimes, $\Omega_{y,t}$ depends on both the regime-specific covariance matrices and the cross terms. Hence, this specification of the structural model is different to the conventional reduced form STVAR specification in which the conditional covariance matrix of $y_t$ is a weighted sum of the covariance matrices of the regimes. Nevertheless, our model specification facilitates exploiting non-Gaussianity of the shocks in their identification as is shown in the next section.
The key identification assumption is that the structural shocks $e_t=(e_{1t},...,e_{dt})$ are mutually independent and at most one of them is Gaussian. We also assume that $e_t$ is an IID sequence with zero mean and identity covariance matrix as is formally stated in the following assumption.
The assumption of a unit variance is merely a normalization since the variance of the shocks is captured by the impact matrix $B_{y,t}$. Also, due to the time-variation of the impact matrix, Assumption (ref) does not imply that the reduced form innovations $u_t$ are independent nor that they are identically distributed. The assumption of strictly positive density on sets of positive (Lebesgue) measure mainly rules out bounded shock distributions.
It is helpful to distinguish two layers of identification. First, for any fixed $t$ and conditional on $\mathcal{F}_{t-1}$, the impact matrix $B_{y,t}$ in (ref) is identified up to column permutations and sign changes by independent component analysis (ICA) under Assumption (ref) Lanne+Meitz+Saikkonen:2017. Second, under the structural parametrization $B_{y,t}=\sum_{m=1}^M \alpha_{m,t} B_m$ in (ref), identification of the regime-specific matrices $B_1,...,B_M$ (up to column ordering and signs) additionally requires identification of the transition weight parameters, or in the case of exogenous weights, sufficient variation in the known weights. The next lemma, partly similar to Proposition 1 in Lanne+Meitz+Saikkonen:2017, therefore formalizes the per-$t$ ICA result, and, in the case of logistic weights ((ref)) establishes identification of the weight function parameters (and, for completeness, identification of the AR parameters for logistic and exogenous weights). Specific weight functions are assumed here for simplicity, but the result can be extended to other suitable weight functions as well. Identification of $B_1,...,B_M$ then follows from these ingredients.
Lemma (ref) is proven in Appendix (ref). Condition (ref) implies distinguishable regimes, whereas Condition (ref) guarantees a certain small amount of variation in exogenous weights. Lemma (ref) concludes that at each $t$, conditionally on $\mathcal{F}_{t-1}$, the impact matrix $B_{y,t}$ is unique up to ordering and signs of its columns under Assumption (ref). That is, at each $t$‚ changing the ordering or signs of the columns of $B_{y,t}$ would lead to an observationally equivalent model, but changing $B_{y,t}$ in any other way would lead to an observationally distinct model. It follows that if the impact matrix is time-invariant as in Lanne+Meitz+Saikkonen:2017, i.e., $B_{y,t}=B$ for some constant matrix $B$, the structural shocks are identified up to ordering and signs. However, when the impact matrix varies over time, two complications arise.
First, because $B_{y,t}$ ((ref)) is identified only up to column ordering and signs at each $t$, its unique identification requires constraints on the regime-specific impact matrices $B_1,...,B_M$ such that any reordering or sign changes in their columns lead to observationally distinct $B_{y,t}$ at some $t$. Second, since $B_{y,t}$ is not a matrix of constant parameters but a function of parameters, it needs to be shown that the parameters in its functional form are identified. Also, in addition to $B_1,...,B_M$, the impact matrix $B_{y,t}$ depends on the transition weights $\alpha_{m,t}$, $m=1,...,M$, which must therefore be identified as well. To that end, we assume either the logistic or exogenous transition weights considered in Lemma (ref).
The following proposition, proven in Appendix (ref), establishes unique identification of $B_1,...,B_M$.
Proposition (ref) essentially states that if the structural shocks are independent and at most one of them is Gaussian (Assumption (ref)), the impact matrices $B_1,...,B_M$, and hence, the structural shocks are identified. The result holds for almost every $[B_1:...:B_M] \in\mathbb{R}^{d\times dM}$ because the identification may fail in some special cases, but this set has Lebesgue measure zero.\footnote{Beyond possible singularity of the impact matrix $B_{y,t}$ (see Lemma (ref)), the failure of the identification in a measure zero set only concerns the possibility that fixing the ordering and signs of the columns of $B_1$ does not fix the ordering and signs of the columns of $B_2,...,B_M$. }
Condition (ref) of Proposition (ref) states that the ordering and signs of the columns of $B_1$ should be fixed (e.g., by assuming that the first nonzero entry in each column is positive and the first nonzero entries are in a decreasing order, given that none of them are equal), which fixes the ordering and signs of the columns of $B_{y,t}$ for all $t$. If the shocks follow, for example, the skewed $t$-distributions defined in Equation ((ref)) (in Section (ref) discussing estimation), any fixed ordering and signs of the columns of $B_1$ can be assumed without loss of generality. This is because reordering the columns would just reorder the shocks and swapping a sign of a column corresponds to swapping the sign of the related shock and skewness parameter. Condition (ref) states that exogenous transition weights should, for some $t$, be strictly positive for all the regimes, which is used to establish that Condition (ref) fixes the ordering and signs of the columns of $B_{y,t}$ for all $t$. For logistic weights, this is achieved via the assumption $\gamma<\infty$, guaranteeing that the regime-switches are not discrete.\footnote{To accommodate discrete regime-switches, which are obtained as a special case of the logistic weights ((ref)) with the smoothness parameter tending to infinity, we establish the identification of the shocks in the threshold VAR model of Tsay:1998 in Appendix (ref).}
Somewhat surprisingly, our experience shows that identification appears to be often weak in the sense that there are multiple local maxima of the (penalized) log-likelihood function (discussed in Section (ref)) with (penalized) log-likelihoods close to each other and to the global maximum. The parameter values related to each such local maximum seem to be often largely quite similar to each other, but have some differences, frequently but not exclusively related to different ordering and signs of the columns of $B_2,...,B_M$. Since such different local solutions may produce different results in structural analysis, this weak identification should be appropriately addressed, for instance, by combining Proposition (ref) with supplementary identifying information as discussed in the next section.
As in the linear SVAR model of Lanne+Meitz+Saikkonen:2017, the statistically identified structural shocks do not necessarily have economic interpretations, and labeling them as economic shocks requires external information. However, labeling the shocks based on the estimates of the impact matrices $B_1,...,B_M$ might not always be straightforward, as the same shock must be associated with the same column of the impact matrix $B_m$ in all regimes. For example, if a positive supply shock should increase output and decrease prices on impact in all regimes, labeling the $i$th shock as the supply shock requires the $i$th column of all $B_1,...,B_M$ to satisfy such signs. Moreover, as discussed in Section (ref), our experience shows that identification is often weak in the sense that there are multiple local solutions with fit close to the global optimum. While we find such local solutions to be often largely quite similar to each other, they may produce different impulse response functions even when the shock of interest is associated with the same column of all $B_1,...,B_M$.
To formally address weak identification and the problem of labeling the shocks, we recommend employing a "blended identification" strategy that combines identification by non-Gaussianity with supplementary identifying information Carriero+Marcellino+Tornese:2024. Specifically, we propose imposing overidentifying restrictions that are sufficient to yield a unique local solution (among the ones with fit close to the global solution) and facilitate labeling the shocks of interest. Since different local solutions may yield different results, the restrictions should be economically reasonable and serve to exclude less plausible alternative solutions.
As a simple example, in a bivariate system of output and prices, it may be reasonable to assume that for each variable, one of the shocks has greater impact effect on that variable than the other shock, and that this effect has the same sign across regimes. If one of the shocks is to be interpreted as a demand shock, it might be useful to further assume that the shock with the largest effect on output also moves prices in the same direction in all regimes. Similarly, if the other shock is intended to represent a supply shock, it could be useful to assume that it moves output and prices in the opposite directions in all regimes.\footnote{ Also other forms of information can be incorporated into the blended identification strategy, such as narrative restrictions and zero impact effect restrictions (the latter of which are testable due to statistical identification), for instance. } See our empirical application in Section (ref) for an example that involves more variables.
As a practical consideration, note that our estimation method, described in Section (ref) and implemented to the accompanying R package sstvars sstvars, produces a set of alternative local solutions. By filtering these local solutions based on whether they satisfy the imposed overidentifying restrictions, it is straightforward to determine which restrictions are sufficient in the considered specification to exclude all but one of the local maximums (among the ones with a fit close to the presumed global optimum). Finally, to label the shocks of interest, the imposed restrictions can often be used to uniquely associate each of them to the same column of the impact matrix across all regimes.
The parameters of the structural STVAR model discussed in Section (ref) can be estimated by the method of maximum likelihood (ML). To obtain a well-behaving estimator with desirable asymptotic properties such as consistency, the parameter space is often restricted to the region where the model is ergodic stationary. Building on the results of Saikkonen:2008 Kheifets+Saikkonen:2020, it can be shown that when the transition weights $\alpha_{m,t}$ are logistic (ref) or of the threshold form, a sufficient condition for ergodic stationarity is that the joint spectral radius of the companion form AR matrices of the regimes is strictly less than one (see Theorem (ref) in Appendix (ref)).\footnote{ For logistic models (with $M=2$ assumed), we additionally require that the matrix $B_1^{-1}B_2$ has no negative real eigenvalues to rule out singular convex combinations of $B_1$ and $B_2$ (see Appendices (ref) and (ref) for details). See Appendix (ref) for the definition of the threshold weights.} Because this condition is computationally demanding to verify Chang+Blondel:2013, it is not particularly useful for restricting the parameter space during estimation. Therefore, we instead impose the usual stability condition for each regime (Condition (ref) in Appendix (ref)), as it is a necessary condition for the sufficient one. The sufficient condition can then be checked after estimation for the solutions of interest.
Maximizing the log-likelihood function can be challenging in practice due to its high multimodality, induced by the nonlinear dynamics.\footnote{The log-likelihood function is presented in Section (ref) for shocks that follow skewed $t$-distributions.} Imposing the stability condition can make estimation particularly difficult when the data is persistent, as in such cases numerical optimization algorithms frequently gravitate to the boundary of the parameter space, where they perform poorly. To overcome this issue, we propose to allow for unstable estimates and to maximize the penalized log-likelihood function, obtained by adding a penalty term to the log-likelihood function that penalizes parameter values falling outside or close to the boundary of the stability region. In this way, the optimization algorithm can explore the parameter space also outside the stability region, facilitating improved performance, while sufficient penalization eventually steers the algorithm back to the stability region.\footnote{Penalized likelihood-based estimation has been previously applied to time series models by Nielsen+Rahbek:2024. They impose nonnegativity constraints in the autoregressive conditional heteroskedasticity model, allowing the estimation algorithm to explore the boundary of the parameter space. Nielsen+Rahbek:2024 use penalization also to determine which parameters should be set to zero.}
Before introducing the penalized log-likelihood function, the log-likelihood function is presented. We assume that each shock $e_{it}$, $i=1,...,d$, follows the skewed $t$-distribution introduced by Hansen:1994 with zero mean, unit variance, $\nu_i>2$ degrees of freedom, and skewness controlled by the parameter $\lambda_i\in (-1, 1)$. The skewed $t$-distribution can flexibly capture fat tails and skewness, and as the Gaussian distribution is obtained as a special case, with $\nu_i=\infty$ and $\lambda_i=0$, the plausibility of the identifying Assumption (ref) can be (informally) assessed based on the estimates of these parameters.
The density function of the skewed $t$-distribution, $st(\cdot; \nu_i, \lambda_i)$, is given as Hansen:1994:
where $\mathbbm{1} \lbrace e_{it}<-a_i/b_i \rbrace$ is an indicator function that takes the value one if $e_{it}<-a_i/b_i$ and zero otherwise, and $\mathbbm{1}\lbrace e_{it}\geq -a_i/b_i\rbrace$ is an indicator function that takes the value one if $e_{it}\geq -a_i/b_i$ and zero otherwise. The constants $a_i, b_i$, and $c_i$ are defined as $a_i = 4\lambda_i c_i\left(\frac{\nu_i - 2}{\nu_i - 1}\right)$, $b_i = (1 + 3\lambda_i^2 - a_i^2)^{1/2}$, and $c_i = \frac{\Gamma \left(\frac{\nu_i + 1}{2}\right)}{(\pi(\nu_i - 2))^{1/2}\Gamma\left(\frac{\nu_i}{2}\right)}$, where $\Gamma(\cdot)$ is the Gamma function.
Indexing the observed data as $y_{-p+1},...,y_0,y_1,...,y_T$, the conditional log-likelihood function, conditional on the initial values $\boldsymbol{y}_0=(y_0,...,y_{-p+1})$, is given as
where $I_{d,i}$ is the $i$th column of the $d$-dimensional identity matrix, $\mu_{y,t}$ is the conditional mean defined in Section (ref) and $B_{y,t}$ is the impact matrix defined in Equation ((ref)). The parameters of the model are collected to the vector $\boldsymbol{\theta}=(\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu)$, where $\varphi_m=(\text{vec}(A_{m,1}),...,$ $\text{vec}(A_{m,p}))$, $m=1,...,M$, $\sigma=(\text{vec}(B_1),...,\text{vec}(B_M))$, $\alpha$ contains the parameters governing the transition weights, and $\nu=(\nu_1,...,\nu_d,\lambda_1,...,\lambda_d)$ contains the degrees-of-freedom and skewness parameters. With logistic weights ((ref)), $\alpha=(c,\gamma)$, whereas with exogenous weights, $\alpha$ is omitted from the parameter vector.
The penalized log-likelihood function is then defined as
where $L(\boldsymbol{\theta})$ is defined in ((ref)) and $P(\boldsymbol{\theta})\geq 0$ is a penalization term that penalizes the log-likelihood function from parameter values that are close to entering or are in an uninteresting region of the parameter space. To focus on avoiding unstable estimates, we define the penalization term as
where $|\rho(\boldsymbol{A}_m(\boldsymbol{\theta}))_i|$ is the modulus of the $i$th eigenvalue of the companion form AR matrix of Regime $m$, the tuning parameter $\eta\in (0, 1)$ determines how close to the boundary of the stability region the penalization starts, and $\kappa>0$ determines the strength of the penalization.
Whenever the companion form AR matrix of a regime has eigenvalues greater than $1 - \eta$ in modulus, the penalization term ((ref)) is greater than zero, and it increases in the modulus of these eigenvalues. Our penalization term incorporates the multiplicative coefficient $Td$ in an attempt to standardize the strength of penalization with respect to the number of observations $T$ and variables $d$, as the magnitude of the log-likelihood typically increases with these quantities. In our empirical application (Section (ref)), we use the tuning parameter values $\eta = 0.05$ and $\kappa = 0.2$, thus, penalizing parameter values outside the stability region significantly, while maintaining flexibility near the boundary.
Lanne+Meitz+Saikkonen:2017 propose a three-step estimation procedure for their linear SVAR model identified by non-Gaussianity. In the first step, least squares (LS) estimation produces preliminary estimates of certain parameters, which serve as initial values in numerical optimization algorithms used in the subsequent steps to maximize the log-likelihood function. The nonlinear dynamics of our structural STVAR model, however, add substantial challenges to the estimation problem by inducing a large number of modes to the (penalized) log-likelihood function. To address these challenges, we propose a modified version of the three-step estimation procedure in which multimodality is explicitly taken into account.
Our three-step estimation procedure proceeds with the following steps:
Step (ref) obtains initial estimates for the AR and weight function parameters, and it is comparable to Step 1 of Lanne+Meitz+Saikkonen:2017 but with the nonlinearity of the model explicitly accounted for. Notably, the NLS estimation often produces estimates that do not satisfy the usual stability condition for each of the regimes, and thus allowing for instability (but penalizing it) facilitates utilization of the NLS estimates in our three-step procedure. Step (ref) is comparable to Step 2 of Lanne+Meitz+Saikkonen:2017, but the multimodality of the (penalized) log-likelihood function is taken into account by making use of a robust estimation algorithm that is able to escape from local maxima. Step (ref) finalizes the estimation by initializing a standard gradient based optimization algorithm from the initial estimates obtained from the previous steps.
Due to the high multimodality of the (penalized) log-likelihood function, we recommend running Steps (ref) and (ref) a large number of times to improve the reliability of the results. Since a large number of estimation rounds are run, the estimation procedure produces a set of estimates, some of which presumably corresponding to different local maximums of the (penalized) log-likelihood function. If there are multiple solutions with (penalized) log-likelihoods close to the greatest found (penalized) log-likelihood, the blended identification strategy discussed in Section (ref) can be employed to address weak identification.
To assess the performance of the penalized maximum likelihood (PML) estimator as well as our three-step estimation procedure, we conduct a small-scale Monte Carlo study, which is discussed in detail in Appendix (ref). Since estimation is computationally demanding, we consider the simple bivariate (structural) LSTVAR model with $p=1$ and $M=2$. To evaluate how the PML estimator performs when the AR matrices are in the penalization region, we consider two specifications of both of LSTVAR models: one where the AR matrices are well inside the stability region and one where they lie clearly in the region where the penalization term is strictly positive. According to the results (see Table (ref) in Appendix (ref) for details), the estimation accuracy is slightly better when the AR matrices are inside the stability region. Both specifications exhibit some bias in small samples. However, since both the bias and the standard deviations of the estimates diminish as the sample size increases, our results align with (possible) consistency of the PML estimator.
Our empirical application studies the macroeconomic effects of climate policy uncertainty (CPU). In the related literature, Fried+Novan+Peterman:2022 develop a dynamic general equilibrium model incorporating beliefs about future climate policy, and find that an increased climate policy transition risk shifts investments toward cleaner capital while reducing total investments and output. Khalil+Strobel:2023 reach similar conclusions in their dynamic stochastic general equilibrium (DSGE) model, attributing the effects to financial institutions' aversion to uncertainty. Huang+Punzi:2024 also use a DSGE model and find that CPU shocks cause firms to delay investments, which reduces output and raises inflation. These findings are supported by their recursively identified Bayesian SVAR model.
We are interested in studying how the effects of the CPU shock vary depending on the level of general economic policy uncertainty (EPU), as the level of uncertainty may affect the behavior of economic agents. Such variation can be accommodated in our structural STVAR model by specifying a logistic transition weight function with the (lagged) level of EPU as the switching variable. This approach allows us to capture potential state-dependent effects, providing more detailed results than linear SVAR analysis.
We consider a monthly U.S. dataset consisting of five variables and covering the time period from 1987:4 to 2024:12. Following Khalil+Strobel:2023 and Huang+Punzi:2024, we measure CPU with the climate policy uncertainty index (CPUI) of Gavriilidis:2021, which is constructed based on the amount of newspaper coverage on climate policy uncertainty topics. In order to control for general economic policy uncertainty, we also include the economic policy uncertainty index (EPUI) of Baker+Bloom+Davis:2016, which is based on newspaper coverage on topics related to economic policy uncertainty as well as on components quantifying the present value of future scheduled tax code expirations and disagreement among professional forecasters over future government purchases and consumer prices.
As a measure of real economic activity, we include the log of industrial production index (IPI), which is detrended by taking first differences. For measuring the price level, we use the log of consumer price index (CPI), likewise detrended by taking first differences. Finally, the monetary policy stance is measured with the effective Federal funds rate that is replaced by the Wu+Xia:2016 shadow rate for the zero-lower-bound periods (RATE). The series of the included variables are presented in Figure (ref).\footnote{The CPUI and EPUI data are retrieved from \url{https://www.policyuncertainty.com}, whereas IPI, CPI, and the Federal funds rate are retrieved from the Federal Reserve Bank of St. Louis database and the Wu+Xia:2016 shadow rate from the Federal Reserve Bank of Atlanta's website. The availability of the CPUI data determines the beginning of our sample period.}
The structural STVAR model is specified from ((ref))-((ref)) by defining the transition weights and distributions of the shocks. We employ the logistic transition weights with two regimes given in Equation ((ref)). They facilitate smooth transitions between the regimes and enable analyzing how the effects of the CPU shock vary depending on the level of EPU, whose first lag is specified as the switching variable. We estimate the structural LSTVAR model by PML with the three-step procedure proposed in Section (ref), and the shocks are assumed to follow skewed $t$-distributions. The autoregressive order $p=2$ is selected based on the Hannan-Quinn information criterion. The relatively low autoregressive order avoids overfitting, particularly in the regime accommodating fewer observations.
To deal with the problem of weak identification discussed in Section (ref) and (ref), we implement a blended identification strategy that imposes the following overidentifying restrictions on the impact matrices $B_1$ and $B_2$. First, for each variable, one of the shocks has a greater impact effect on that variable in both regimes than the other shocks, this shock is unique to that variable, and the impact effect has the same sign in both regimes. Second, the shock that has the greatest impact effect on output moves inflation in the same direction in both regimes (consistent with a demand shock). Third, the shock that has the greatest impact effect on inflation moves output in the opposite direction in both regimes (consistent with a supply shock). Fourth, the shock that has the greatest impact effect on interest rate moves output in the opposite direction in both regimes (consistent with a monetary policy shock). The shock that has the greatest impact effect on CPUI is deemed as the CPU shock, and to support our identification, we check that the recovered CPU shock aligns with several major historical events affecting CPU. In particular, we find that President Trump's election is accompanied by large positive CPU shocks in November 2016 and 2024. Moreover, there is a large negative CPU shock in January 2021, when President Biden signed executive orders to re-enter the Paris Agreement and launch a broad federal climate strategy.
The location parameter estimate $\hat{c}=16.16$, roughly implying that Regime 1 dominates when the level of EPU is not high, while Regime 2 prevails in times of high EPU. Hence, we label Regime 1 as the Low EPU Regime and Regime 2 as the High EPU Regime. The scale parameter estimate $\hat{\gamma}=64.23$, indicating that transitions between the regimes are quite fast, which can be observed also from the bottom panel of Figure (ref), where the evolution of the fitted transition weight of the High EPU Regime is depicted.
The estimated vectors of the degrees-of-freedom and skewness parameters are $(2.77, 2.71, 3.00, 3.15, 2.03)$ and $(0.38, 0.41, -0.16, 0.15, 0.11)$, respectively, in the order of the variables CPUI, EPUI, IPI, CPI, and RATE. Clearly, none of the shocks is close to Gaussian, suggesting that Assumption (ref) is plausible. Then, we check whether our model satisfies the stationarity condition by computing an upper bound for the joint spectral radius of the companion form AR matrices of the regimes (see Theorem (ref) in Appendix (ref)). The obtained upper bound is strictly less than one ($0.98$), implying that our model is ergodic stationary.\footnote{In addition, we checked that the estimated $B_1^{-1}B_2$ does not have real negative eigenvalues, so the additional condition mentioned in Footnote (ref) for logistic weights is satisfied.} Finally, based on graphical residual diagnostics, the overall adequacy of our model seems reasonable (for details, see Appendix (ref)).
The effects of the shocks in our structural LSTVAR model may depend on the initial state of the economy as well as on the sign and size of the shock, making the conventional way of calculating impulse responses unsuitable. Therefore, we consider the generalized impulse response function (GIRF) Koop+Pesaran+Potter:1996 that accommodates such features, defined as
where $h$ is the horizon. The first term on the right side of ((ref)) is the expected realization of the process at time $t+h$ conditionally on a structural shock of sign and size $\delta_i \in\mathbb{R}$ in the $i$th element of $e_t$ at time $t$ and the previous observations. The latter term on the right side is the expected realization of the process conditionally on the previous observations only. The GIRF thus expresses the expected difference in the future outcomes when the $i$th structural shock of sign and size $\delta_i$ arrives at time $t$ as opposed to all shocks being random. Since our model has the $p$-step Markov property, the conditioning set $\mathcal{F}_{t-1}$ can be replaced by the lag vector $\boldsymbol{y}_{t-1}=(y_{t-1},...,y_{t-p})$. The GIRF also facilitates tracing out the effects of the shocks on the transition weights $\alpha_{m,t}$, $m=1,...,M$, by replacing $y_{t+h}$ with $\alpha_{m,t+h}$ on the right side of Equation ((ref)).
To study the state-dependent effects of the CPU shock, we follow a procedure similar to Lanne+Virolainen:2025 and calculate the GIRFs conditional on a given regime dominating when the shock arrives. Specifically, for each Regime $m\in \lbrace 1, 2\rbrace$, we take all the length $p$ histories $\boldsymbol{y}_{t-1}$ from the data for which the corresponding transition weight $\alpha_{m,t}$ is greater than $0.75$, indicating that Regime $m$ is clearly dominant. For each such history, we take the corresponding CPU shock recovered from the data, and compute the GIRF using the Monte Carlo algorithm described in Lanne+Virolainen:2025, Appendix B. Finally, GIRFs corresponding to shocks of different sign and size are made comparable by scaling them to correspond to a $5$ point instantaneous increase of CPUI.
Figure (ref) illustrates the distribution of the GIRFs in each regime in a so-called "shotgun plot", which depicts each GIRF using a level of opacity such that the darkness of a region displays the concentration of GIRFs in it Lanne+Virolainen:2025. The response of CPUI to a positive CPU shock follows a similar pattern in both regimes in the vast majority of the GIRFs. EPUI increases in both regimes, but response is stronger in the High EPU Regime.
Consistent with the previous literature (Fried+Novan+Peterman:2022, Khalil+Strobel:2023, and Huang+Punzi:2024), a positive CPU shock decreases production in both regimes. Huang+Punzi:2024 attribute the decline in output to firms postponing investment, whereas Fried+Novan+Peterman:2022 and Khalil+Strobel:2023 emphasize reallocation of capital away from carbon-intensive sectors. We find that the contraction is slightly stronger in the High EPU Regime than in the Low EPU Regime, possibly reflecting greater risk aversion among firms and financial intermediaries during periods of elevated economic policy uncertainty.
In both regimes, a positive CPU shock increases prices, in line with Huang+Punzi:2024 and Khalil+Strobel:2023. Huang+Punzi:2024 find that a positive CPU shock increases inflation by acting as a negative supply shock, whereas according to Khalil+Strobel:2023, a positive CPU shock may decrease aggregate demand due to risk-averse households raising savings, while firms raise their prices preemptively to lower the probability of being stuck with a negative markup. However, we find that the inflationary effects are stronger in the High EPU Regime than in the Low EPU Regime. One possible explanation is that the firms might be more risk averse in the periods of high EPU and thereby react more sensitively.
The response of the interest rate variable seems consistent with the inflationary effects of the CPU shock: it increases in both regimes, but clearly more so in the High EPU Regime. Finally, consistent with the GIRFs of EPUI, a positive CPU shock increases the transition weights of the High EPU Regime, more so in the High EPU Regime, where there is first a short-term decline in the weights. In other words, a positive CPU shock generally drives the economy towards the High EPU Regime. Overall, our results suggest that a positive CPU shock reduces production and raises inflation in both regimes, but its effects, particularly on inflation, are stronger under high economic policy uncertainty.
Linear structural vector autoregressive models can be identified statistically by non-Gaussianity, provided that the shocks are mutually independent and at most one of them is Gaussian Lanne+Meitz+Saikkonen:2017. This paper is the first to extend identification by non-Gaussianity to nonlinear SVAR models. Specifically, we have shown that structural smooth transition vector autoregressive models, featuring a time-varying impact matrix defined as a weighted sum of regime-specific impact matrices, are likewise identified under these assumptions. While we have focused on exogenous, logistic, and threshold transition weights typically used in macroeconomic applications, our results can be extended to other suitable weight functions as well.
In addition to establishing the identification of the shocks, we have discussed the problem of labeling them. We have also employed blended identification Carriero+Marcellino+Tornese:2024 to address weak identification, which in our experience is often but not exclusively related to the ordering and signs of the columns of the regime-specific impact matrices. We have proposed estimating the model parameters by the method of penalized maximum likelihood, and we introduced a three-step estimation procedure by adapting the approach of Lanne+Meitz+Saikkonen:2017 to nonlinear SVAR models. Building on the results of Saikkonen:2008, we have also provided a sufficient condition for ergodic stationarity of our structural STVAR model. The introduced methods are implemented in the accompanying R package sstvars sstvars, which is available via the CRAN repository.
In an empirical application to U.S. data from 1987:4 to 2024:12, we have studied the macroeconomic effects of the climate policy uncertainty shock using a two-regime logistic STVAR model. As a measure of climate policy uncertainty, we employ the CPU index Gavriilidis:2021 constructed based on the amount of newspaper coverage on climate policy uncertainty topics. To distinguish between periods of low and high economic policy uncertainty, we use the first lag of the EPU index Baker+Bloom+Davis:2016 as the switching variable. Our results are in line with the previous literature suggesting that a positive CPU shock reduces production and raises inflation (Khalil+Strobel:2023 and Huang+Punzi:2024), while Fried+Novan+Peterman:2022 found that it decreases production. However, we found that while a positive CPU shock reduces production and raises inflation in times of both low and high EPU, its effects, particularly on inflation, are stronger during the periods of high EPU.