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.
67,076 characters · 14 sections · 68 citation commands
Smooth transition vector autoregressive (STVAR) models that allow for gradual shifts in parameter values have become popular due to their ability to capture nonlinear dynamics in time series data Hubrich+Terasvirta:2013. They extend the conventional linear vector autoregressive (VAR) model by allowing for smooth transitions between regimes, each of which is characterized by a linear vector autoregression with different coefficients or error covariance matrices. Each observation is the sum of a weighted average of the conditional means of the regimes and a random error, whose covariance matrix is a weighted average of the covariance matrices of the regimes. Like other nonlinear structural VAR models, the structural counterparts of STVAR models enable tracing out the causal effects of economic shocks, which may depend on the initial state of the economy or on the sign or size of the shock.
Different STVAR models are obtained by specifying the transition weights or the error distribution in various ways. In this paper, we introduce the Gaussian smooth transition vector autoregressive (GSTVAR) model, where the transition weights are similar to the mixing weights of the Gaussian mixture vector autoregressive (GMVAR) model of Kalliovirta+Meitz+Saikkonen:2016. For a $p$th order model, the transition weights in both models depend on the full distribution of the preceding $p$ observations, which enables capturing complicated switching dynamics. Specifically, the greater the relative weighted likelihood of a regime is, the greater its transition weight is, which facilitates associating statistical characteristics and economic interpretations to the regimes. However, in contrast to the GMVAR model, which incorporates unobserved discrete regime switches, our GSTVAR model has the advantage that it allows for capturing gradual shifts in the dynamics of the data. Moreover, the identified shocks of the structural GSTVAR model can be recovered from the fitted model and used in impulse response analysis to further reflect the properties of the data, which is not the case with the structural GMVAR model.
Probably the most popular among the STVAR models put forward in the previous literature is the logistic STVAR (LSTVAR) model Anderson+Vahid:1998 characterized by logistic transition weights. In the LSTVAR model, the interpretation of the regimes is clear, and it can accommodate exogenous as well as endogenous switching variables. These properties are shared by the often employed threshold VAR (TVAR) model Tsay:1998, albeit the latter features discrete regimes instead of allowing for smooth transitions between them. Hence, these models are well-suited, when the transition weights can be expected to depend on the level of some specific switching variables. In contrast, our GSTVAR model facilitates capturing more complicated switching dynamics and results in more clearly data-driven regimes, but it may not always be straightforward to interpret the regimes economically.
Besides STVAR and TVAR models, there are a number of alternative approaches to modeling nonlinearities in multivariate time series, including so-called Markov-switching VAR (MS-VAR) models Krolzig:1997 and time-varying parameter VAR (TVP-VAR) models (Cogley+Sargent:2001; Cogley+Sargent:2005; Primiceri:2005; Koop+Leon-Gonzalez+Strachan:2009; among others). The former incorporate discrete regime switches that are random and unobserved, and the switching probabilities depend on only the preceding regime. While they can flexibly capture nonlinearities in the data, unlike our GSTVAR model, they do not accommodate gradual shifts in the dynamics nor complicated switching dynamics. As to TVP-VAR models, compared to our GSTVAR model, they have the major limitation that the stochastic processes governing the changes in parameter values are exogenous to the variables included, which may lead to overlooking important endogenous dynamics.
We apply the structural GSTVAR model to find out about the macroeconomic effects of severe weather. This issue was recently studied by Kim+Matthes+Phan:2022 by a STVAR model with transition weights defined by a linear time trend, which only allows a single switch from one regime to the other. Following them, we include an indicator of severe weather in addition to a number of macroeconomic variables in the model and estimate it on monthly U.S. data from 1961:1 to 2022:3. In line with Kim+Matthes+Phan:2022, we find evidence in favor of nonlinearity, and a GSTVAR model with two regimes is deemed sufficient. However, our regimes are characterized by features quite different from theirs. One regime dominates the earlier part of the sample period, particularly the volatile periods of the 1970s and 1980s, and also prevails later during the Financial crisis and the COVID-19 crisis, while the latter part of the sample period is predominantly characterized by the other regime.
In both regimes, a recursively identified positive severe weather shock decreases GDP, consumer prices, and the interest rate, but the effects are clearly stronger in the regime prevailing in the earlier part of the sample period and in certain crisis periods. In contrast to Kim+Matthes+Phan:2022, who found the effects of severe weather stronger in the latter part of the sample period, we thus conclude that the U.S. economy has adapted to the changing distribution of weather related shocks over time. However, since strong effects are found also in the crisis periods, it seems that the apparent adaptation does not provide sufficient resilience when the economy is in a vulnerable state.
The rest of this paper is organized as follows. Section (ref) first presents our framework for STVAR models and then discusses their stationarity and ergodicity. We also introduce the structural STVAR model and some tools useful in empirical analysis based on it. Finally, the new Gaussian STVAR (GSTVAR) model is introduced by specifying the transition weights and the error distribution. In Section (ref), model selection and estimation of the parameters by the method of maximum likelihood are discussed. Section (ref) contains the empirical application to the macroeconomic effects of severe weather shocks. Section (ref) concludes the paper. Further details can be found the in the appendices, and the introduced methods are implemented to the accompanying R package sstvars sstvars, available in the CRAN repository.
The STVAR model with $M$ regimes and autoregressive order $p$ considered in this paper can be written as
where $\phi_{1,0},...,\phi_{M,0}\in\mathbb{R}^{d}$, $m=1,...,M$, are the intercept parameters, $A_{1,i},...,A_{M,i}\in\mathbb{R}^{d\times d}$, $i=1,...,p$, are the autoregressive matrices, and $\Omega_1,...,\Omega_M$ are the positive definite $(d\times d)$ covariance matrices of the regimes. The serially uncorrelated reduced form innovations $u_t$ follow some distribution with zero mean and conditional covariance matrix $\Omega_{y,t}$.
The transition weights $\alpha_{m,t}$ are assumed to be $\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$. These weights convey the relative proportions of the regimes at each point in time and determine how the process shifts between them. Specifically, when the process is completely in one of the regimes, the transition weight $\alpha_{m,t}$ of that regime equals unity, while the weights of the rest of the regimes are equal to zero. As the process begins a shift towards another regime, the transition weight of the emerging regime increases due to the dynamics captured in the transition weight function through the preceding observations. Concurrently, the weight of the previously dominant regime decreases.
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}$, and its conditional covariance matrix is $\Omega_{y,t} = \text{Cov}(y_t|\mathcal{F}_{t-1}) = \sum_{m=1}^M \alpha_{m,t}\Omega_m$. That is, the conditional mean is a weighted sum the regime-specific means $\mu_{m,t}$ with the weights given by the transition weights $\alpha_{m,t}$, whereas the conditional covariance matrix is a weighted sum of the regime-specific conditional covariance matrices $\Omega_m$.
For standard asymptotic inference, the STVAR process must be stationary and ergodic. To verify that this is indeed the case, we rely on the sufficient condition for stationarity and ergodicity derived by Kheifets+Saikkonen:2020 Saikkonen:2008. However, their parametrization is not quite the same as ours, so their result must be slightly modified for our purposes (see Appendix (ref) for details). In terms of the companion form AR matrices of the regimes defined as
Kheifets+Saikkonen:2020 show that if the following condition holds, the STVAR process is ergodic stationary (both strictly and second-order).
Here
denotes the joint spectral radius (JSR) of a finite set of square matrices $\mathcal{A}$ with $\mathcal{A}^j=\lbrace A_1A_2...A_j:A_i\in\mathcal{A}\rbrace$ and $\rho(A)$ is the spectral radius of the square matrix $A$. As Kheifets+Saikkonen:2020 note, Condition (ref) is not necessary for ergodic stationarity of the process, meaning that if it does not hold, we just cannot use the result to state whether the process is ergodic stationary or not.
It is worth mentioning that a necessary condition for Condition (ref) is that the usual stability condition is satisfied for each of the regimes, as $\max(\rho(\boldsymbol{A}_1),...,\rho(\boldsymbol{A}_M))\leq \rho(\lbrace \boldsymbol{A}_1,...,\boldsymbol{A}_M \rbrace)$ Kheifets+Saikkonen:2020. Therefore, the following condition, which is analogous to Corollary 1 of Kheifets+Saikkonen:2020, is necessary for Condition (ref).
Since the sufficient Condition (ref) is computationally costly to verify with reasonable accuracy Chang+Blondel:2013, it is useful to employ the more easily verified necessary Condition (ref) in numerical estimation (see Section (ref)). As our estimation procedure produces a set of alternative local solutions, Condition (ref) can be checked for each of them after the estimation, and the best local solution is selected among those for which stationarity can be verified.
A number of methods for bounding the JSR have been proposed in the literature, many of which are discussed by Chang+Blondel:2013. The accompanying R package sstvars sstvars implements the branch-and-bound method of Gripenberg:1996. The JSR toolbox in MATLAB Jungers:2023, in turn, automatically combines various methods in the estimation of the JSR to enhance computational efficiency.
To conduct structural analysis, we need to find the structural counterpart of the STVAR model that involves orthogonal, serially uncorrelated structural errors, or shocks, $e_t$. This amounts to finding a non-singular ($d\times d)$ impact matrix $B_t$ such that the conditional covariance matrix of $e_t = B_t^{-1}u_t$, conditional on $\mathcal{F}_{t-1}$, is a diagonal matrix (typically normalized to the identity matrix). In other words, $B_t$ must be such that $B_t^{-1}\Omega_{y,t}B_t'^{-1}=I_d$. However, without additional identifying restrictions $B_t$ satisfying this equation is not unique, and as discussed by, e.g., Kilian+Lutkepohl:2017 (Kilian+Lutkepohl:2017, Chapter 18), imposing such restrictions in nonlinear structural VAR models may not be straightforward. In our empirical application in Section (ref), we consider a recursive model, where $B_t$ is obtained by a lower-triangular Cholesky decomposition of the conditional covariance matrix $\Omega_{y,t}$.
Conventional impulse responses can be computed separately for each regime in the structural STVAR model, as has been done in much of the nonlinear structural VAR literature. However, by thus precluding future regime switches, this approach changes the structure of the model and makes the impulse response analysis subject to Lucas critique. Moreover, the expected effects of the structural shocks generally depend on the initial values of the variables as well as on the sign and size of the shock, which is not accounted for by the conventional impulse responses. An appropriate alternative is the generalized impulse response function (GIRF) of Koop+Pesaran+Potter:1996, defined as
where $h$ is the horizon and $\mathcal{F}_{t-1}=\sigma\lbrace y_{t-j},j>0\rbrace$ as before. The first term on the right side of ((ref)) is the expected realization of the process at time $t+h$ conditional on a structural shock of sign and size $\delta_j \in\mathbb{R}$ in the $j$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 conditional on the previous observations only. The GIRF thus expresses the expected difference in the future outcomes when the structural shock of sign and size $\delta_j$ in the $j$th element hits the system at time $t$ as opposed to all shocks being random. An interesting feature of the structural STVAR model is that besides the generalized impulse response functions of the observable variables, it produces the GIRFs of the transition weights $\alpha_{m,t}$, $m=1,...,M$, which yield information on the dynamic effects of the shocks on regime switches. They are obtained by replacing $y_{t+h}$ with $\alpha_{m,t}$ on the right side of Equation ((ref)).
Typically STVAR models have a $p$-step Markov property, and this is also the case with our GSTVAR model to be introduced in Section (ref). Thus conditioning on (the $\sigma$-algebra generated by) the $p$ previous observations $\boldsymbol{y}_{t-1}=(y_{t-1},...,y_{t-p})$ is effectively the same as conditioning on $\mathcal{F}_{t-1}$ at time $t$ and later. To estimate the GIRFs conditional on the economy being in a specific regime, say $\tilde{m}$, when the shock arrives, we select the length $p$ histories $\boldsymbol{y}_{t-1}$ from the data that indicate the dominance of this regime in the time period $t$. That is, we select the histories $\boldsymbol{y}_{t-1}$ for which the transition weight $\alpha_{\tilde{m},t}$ is large, say, larger than $0.75$, indicating that the process is largely in Regime $\tilde{m}$. For each of these histories $\boldsymbol{y}_{t-1}$, we compute the GIRFs to the structural shock $e_{it}$ recovered from the data (the $i$th element of the estimate of $B_t^{-1}(y_t - \mu_{y,t})$), so they reflect the properties of the data. To then make GIRFs to shocks with different signs and sizes comparable, they are all scaled to correspond to some fixed (say, unit) instantaneous increase in one of the variables. Finally, the distribution of GIRFs obtained for all histories $\boldsymbol{y}_{t-1}$ can be superimposed in a so-called “shotgun plot" Inoue+Kilian:2016. When a low level of opacity is used, the darkness of each region in the figure indicates the concentration of GIRFs in it (for an example, see Figure (ref)). A detailed description of our Monte Carlo algorithm for estimating the GIRF is presented in Appendix (ref).\footnote{Another way to estimate GIRFs conditional on the economy being in a specific regime when the shock arrives is to generate the initial values $\boldsymbol{y}_{t-1}=(y_{t-1},...,y_{t-p})$ from the stationary distribution of that regime. Comparison of regime-dependent GIRFs for different signs and sizes of the shock then reveals asymmetries with respect to these features within the regime.}
Similarly to the conventional impulse response functions being unsuitable for impulse response analysis in the structural STVAR model (due to their inability to capture asymmetries in the effects of the shocks or to take future switches in the regime into account), the conventional forecast error variance decomposition is unsuitable for tracking the relative contribution of each shock to the variance of the forecast errors. Instead, the generalized forecast error variance decomposition (GFEVD) of Lanne+Nyberg:2016 can be used. It is defined for variable $i$, shock $j$, and horizon $h$ as
where $h$ is the chosen horizon and $\text{GIRF}(l,\delta_j,\mathcal{F}_{t-1})_i$ is the $i$th element of the related GIRF. The GFEVD is otherwise similar to the conventional forecast error variance decomposition but with GIRFs in the place of conventional impulse response functions, and it can be interpreted in a similar manner to the conventional forecast error variance decomposition.
Specifying a particular STVAR (and the related structural STVAR) model from the definition in Section (ref) amounts to specifying the distribution of the reduced form innovations $u_t$ (or the structural shocks $e_t$) and the transition weights $\alpha_{m,t}$. In this paper, we propose a Gaussian STVAR model with standard normal distributions for the structural errors $e_t$. Hence, the conditional distribution of $y_t$, conditional on $\mathcal{F}_{t-1}$, is Gaussian and characterized by the density function
That is, the conditional distribution is the $d$-dimensional Gaussian distribution with mean $\mu_{y,t}$ and covariance matrix $\Omega_{y,t}$.
The GSTVAR model has the advantage that, as the conditional distribution is Gaussian, the stationary distributions of the regimes corresponding to $p$ consecutive observations are known, and, hence, the weighted relative likelihoods of the regimes can be used as transition weights. In this specification, the transition weights depend on the full distribution of the preceding $p$ observations, and they are defined identically to the mixing weights of the GMVAR model Kalliovirta+Meitz+Saikkonen:2016. Denoting $\boldsymbol{y}_{t-1}=(y_{t-1},...,y_{t-p})$, the transition weights are defined as
where $\alpha_1,...,\alpha_M$ are transition weight parameters that satisfy $\sum_{m=1}^M \alpha_m=1$ and $n_{dp}(\cdot;\boldsymbol{1}_p\otimes \mu_m, \boldsymbol{\Sigma}_{m,p})$ is the density function of the $dp$-dimensional normal distribution with mean $\boldsymbol{1}_p\otimes \mu_m$ and covariance matrix $\boldsymbol{\Sigma}_{m,p}$. The symbol $\boldsymbol{1}_p$ denotes a $p$-dimensional vector of ones, $\otimes$ is Kronecker product, $\mu_m=(I_d - \sum_{i=1}^pA_{m,i})^{-1}\phi_{m,0}$, and the covariance matrix $\boldsymbol{\Sigma}_{m,p}$ is given in Lutkepohl:2005, but using the parameters of the $m$th regime. In other words, $n_{dp}(\cdot;\boldsymbol{1}_p\otimes \mu_m, \boldsymbol{\Sigma}_{m,p})$ corresponds to the density function of the stationary distribution of the $m$th regime.
The transition weights are thus weighted ratios of the stationary densities of the regimes corresponding to the preceding $p$ observations. The specification of the transition weights is appealing, as it states that the greater the weighted relative likelihood of a regime is, the greater the weight of this regime is. The regimes are, hence, formed based on the statistical properties of the data and are not affected by the choice of switching variables. This is a convenient feature for forecasting, and it also facilitates associating statistical characteristics and economic interpretations to the regimes.
Our GSTVAR model has a number of desirable features compared to popular alternative nonlinear VAR models. As already discussed, in line with the GMVAR model of Kalliovirta+Meitz+Saikkonen:2016, the transition between the regimes depends on the full distribution of multiple past observations. By contrast, in the Markov-switching VAR (MS-VAR) model,the transition probabilities depend only on the most recent regime, whereas in the logistic STVAR (LSTVAR) and threshold VAR (TVAR) models, the regime is determined by the lagged value(s) of some observable variable(s) or exogenous switching variable(s). On the other hand, like the LSTVAR model, the GSTVAR model facilitates gradual shifts between the regimes, whereas in the TVAR, GMVAR and MS-VAR models each observation is generated from a single regime. Hence, the GSTVAR model combines several advantages of the GMVAR, TVAR, LSTVAR and MS-VAR models. A more detailed comparison of the models is given in Appendix (ref).
The parameters of the GSTVAR model can be estimated by the method of maximum likelihood (ML). We collect the parameters to the vector $ \boldsymbol{\theta}=(\phi_{1,0},...,\phi_{m,0},\varphi_1,...,\varphi_M,\sigma,\alpha) $, where $\varphi_m=(\text{vec}(A_{m,1}),....,\text{vec}(A_{m,p}))$, $m=1,...,M$, $\sigma=(\text{vech}(\Omega_1),...,\text{vech}(\Omega_M))$, and $\alpha=(\alpha_1,...,\alpha_{M-1})$ contains the transition weight parameters (notice that $\alpha_M=1-\sum_{m=1}^{M-1} \alpha_m$).
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 $n_d(y_t;\mu_{y,t},\Omega_{y,t})$ is the $d$-dimensional conditional density of the process, conditional on $\mathcal{F}_{t-1}$, at time $t$, given in Equation ((ref)). Hence, the ML estimator of $\boldsymbol{\theta}$ maximizes the log-likelihood function $L_t(\boldsymbol{\theta})$ over the parameter space specified below.
To ensure ergodic stationarity of the process, we assume that Condition (ref) in Section (ref) holds. Moreover, it is assumed that the true parameter value is an interior point of a compact subset of the parameter space, which is a standard condition for asymptotic normality of the ML estimator. Thus, given ergodic stationarity of the process, there is no particular reason to expect that the standard asymptotic results of consistency and asymptotic normality would not apply to the ML estimator. Finally, to achieve identification of the parameters such that the regimes cannot 'relabelled' to obtain the same model with different parameter vectors,we order the transition weight parameters $\alpha_m$, $m=1,...,M$ in a decreasing order. That is, we assume
The constraints imposed on the parameter space are summarized in the following assumption.
Finding the ML estimate amounts to maximizing the log-likelihood function in ((ref)) over a high dimensional parameter space satisfying the constraints in Assumption (ref). Due to the complexity of the log-likelihood function, numerical optimization methods are required. The maximization problem can, however, be challenging in practice due to the complicated dependence of the transition weights on the preceding observations, which induces a large number of modes of the log-likelihood function and large areas of the parameter space, where it is flat in multiple directions.
We follow Virolainen:2022, Virolainen:2025 (and others) and employ a two-phase estimation procedure that is run for a large number of times.\footnote{For example, in our empirical application presented in Section (ref), we use $6000$ estimation rounds to estimate the four-variate GSTVAR with $p=4$, $M=2$.} In the first phase, a genetic algorithm is used to find parameter values that lie close to the local maximum points of the log-likelihood function. Since genetic algorithms tend to converge slowly near local solutions, a gradient based variable metric algorithm is run for each of the starting values, resulting in a number of alternative local solutions. Some of the estimation rounds may end up in saddle points or near-the-boundary points that are not local solutions, and some of the local solutions may be inappropriate for statistical inference (for instance, there might be only few observations from some of the regimes). Because Condition (ref) included in Assumption (ref) is computationally costly to verify, we recommend using the necessary Condition (ref) to restrict the parameter space.
After the estimation rounds have been run, the researcher can choose among the appropriate local solutions the one that maximizes the log-likelihood. Then, it can be verified that the selected local solution satisfies the sufficient Condition (ref) for ergodic stationarity if the more easily checked necessary Condition (ref) is used to restrict the parameter space in estimation. The accompanying R package sstvars sstvars employs a modified genetic algorithm that works similarly to that described in Virolainen:2022.
To study the finite-sample properties of the ML estimator, we conduct a small scale Monte Carlo study, which is discussed in detail in Appendix (ref). Because estimation is computationally demanding, we consider a simple bivariate GSTVAR model with $M=2$ and $p=1$. We consider two specifications that differ with respect to the values of the autoregressive parameters and intercepts (see Table (ref) in Appendix (ref)). According to the results (see Table (ref) in Appendix (ref)) the estimator is slightly biased in small samples, but the bias vanishes and estimation accuracy increases with the sample size. Also, estimation accuracy does not seem to depend on the specification.
To select the number of regimes and the autoregressive order of the GSTVAR model, a suitable strategy is to start with a relatively simple specification and then build up to more complicated models if necessary. In particular, it may be useful to begin with linear Gaussian VAR models (i.e., GSTVAR models with only one regime) to evaluate to what extent they can capture the relevant characteristics of the data. Then, the order of the model or the number of regimes can be increased, if needed.
It is well known that testing for linearity against a model with multiple regimes, in general, poses a nonstandard testing problem because the model is identified only under the alternative Davies:1977. This is the case also in our setup, and, therefore, we follow the previous literature (e.g., Kalliovirta+Meitz+Saikkonen:2016 and Virolainen:2022, Virolainen:2025) and recommend using information criteria and residual-based diagnostic checks to compare the fit of models with different orders and numbers of regimes, as well as to study their adequacy in capturing the autocorrelation structure, conditional heteroskedasticity, and distribution of the data. The derivation of formal diagnostic tests is beyond the scope of this paper, but graphical devices, including quantile-quantile plots as well as autocorrelation and cross-correlation functions of the residuals and their squares, can be used to this end.
For a number of reasons, GSTVAR models of a relatively low autoregressive order are preferable. Firstly, the estimation of GSTVAR models with a large $p$ can be tedious in practice. Secondly, if the order $p$ is large, the number of parameters increases vastly when the number of regimes is increased, which is likely to substantially decrease estimation accuracy. Thirdly, unlike linear VARs, decreasing the autoregressive order may actually improve the fit because the transition weights are calculated using the whole joint distribution of the preceding $p$ observations. Consequently, the sensitivity of the transition weights to individual observations and, hence, responsiveness to changes in the dynamics of the data decrease with $p$.
It is also advisable to be conservative with the number of regimes $M$ because if the number of regimes is chosen too large, some of the parameters in the model are not identified Kalliovirta+Meitz+Saikkonen:2016. Having too many regimes may also lead to overfitting, as the number of parameters may become overly large compared to the number of observations in some of the regimes. Moreover, increasing the number of regimes substantially increases the complexity of the surface of the log-likelihood function, making estimation of the parameters challenging in practice.
In this section, we apply our GSTVAR model to study the effects of severe weather on the U.S. macroeconomy. This issue was recently addressed by Kim+Matthes+Phan:2022, who employ a two-regime recursive STVAR model with transition weights defined by a linear time trend ($\alpha_{2,t}=t/T$ and $\alpha_{1,t}=1-\alpha_{2,t}$ in (ref)). They argue that allowing for time-varying coefficients in the model is important because the strength of the economic impact of severe weather may have diminished due to actions potentially taken to make the economy more adaptable to climate change. However, they find little evidence in favor of adaptation to severe weather in the U.S, but on the contrary, their results suggest that the severe weather shock has a significant effect only towards the end of the sample (spanning from 1963 to 2019).
Using transition weights defined by a linear time trend in a two-regime STVAR model only allows one smooth transition from one regime to the other during the sample period. This approach may be too simple to describe the adaptation of the economy, as the effects of severe weather may vary across different states of the economy due to variation in factors such as fiscal flexibility and sector-specific vulnerabilities. However, it is not obvious which exogenous or lagged endogenous variables would be capable of capturing such variation in the transition weights in a nonlinear SVAR model. Therefore, we address this issue with our GSTVAR model in which the regimes and shifts between them are formed based on the statistical properties of the data through the full distribution of the preceding $p$ observations. Due its more data-driven nature, this approach plausibly better describes the evolution of the joint dynamics of severe weather and U.S. macroeconomy through time than an STVAR model with transition weights defined by a linear time trend or depending only on the level of exogenous or lagged endogenous variables.
We consider a four-variable monthly U.S. data set that comprises indicators of severe weather and economic activity as well as consumer price inflation and an interest rate variable. While our sample period (from 1961:1 to 2022:3, $T=735$ observations) is somewhat longer than than in Kim+Matthes+Phan:2022, the main conclusions based on their shorter sample period (from 1963:4 to 2019:5) remain the same (see Appendix (ref) for the subsample results). The end point of our sample period is determined by the availability of the monthly GDP growth data, but as argued by Kim+Matthes+Phan:2022, it is important to use data of higher frequency than quarterly because some weather effects can be short-lived.
As an indicator of the frequency of severe weather and the extent of sea level rise, we use the Actuaries Climate Index (ACI), developed by actuarial organizations in the United States and Canada ACI:2023. Following Kim+Matthes+Phan:2022, we seasonally adjust the ACI series with the standard Census Bureau X-13 seasonal adjustment algorithm. For a more detailed description of ACI, see Kim+Matthes+Phan:2022.
As an aggregate measure of real economic activity, we use the monthly GDP growth rate constructed by the Federal Reserve Bank of Chicago from a collapsed dynamic factor analysis of a panel of $500$ monthly measures of real economic activity and quarterly real GDP growth MGDP:2023. While Kim+Matthes+Phan:2022 use both the industrial production index and unemployment rate to measure real economic activity at the monthly frequency, we opt for the monthly GDP growth rate instead. It has the advantage of being a more comprehensive measure of real economic activity, including agricultural output, tourism, and services, which can be significantly affected by the severe weather shocks. Moreover, using only a single variable to measure real economic activity is likely to facilitate estimation by reducing the dimension, and thereby also the number of parameters, of the model.
Finally, we include the monthly growth rate of the consumer price index (CPI) and the effective federal funds rate (RATE) that is replaced by the Wu+Xia:2016 shadow rate for the zero lower bound periods. The CPI and federal funds rate series were downloaded from the Federal Reserve Bank of S.t. Louis database and Wu+Xia:2016 shadow rate from the Federal Reserve Bank of Atlanta database. The time series are depicted in the top four panels of Figure (ref), where the shaded areas indicate the U.S. recessions defined by the NBER.
To select the order of the GSTVAR model, we start by examining the partial autocorrelation functions of the U.S. time series (shown in Figure (ref) of Appendix (ref)). They suggest that the autoregressive order $p=4$ might be adequate. Therefore, we fit two-regime GSTVAR models with autoregressive orders $p=1,....,5$ and, to compare them, compute the values of three information criteria (AIC, BIC and HQIC). As seen in Table (ref), the order $p=4$ minimizes the AIC, whereas the order $p=2$ minimizes the HQIC and BIC. The table also contains similar results for linear VAR (i.e., GSTVAR with $M=1$) models, which are clearly inferior to the two-regime GSTVAR models in terms of the information criteria. Moreover, the constancy of the AR matrices and intercepts as well as the constancy of AR matrices only are clearly rejected by the Wald test (see Appendix (ref)). While adding a third regime to the model might improve the fit, the number of parameters compared the number of observations in each regime would increase substantially, possibly leading to overfitting. Moreover, since incorporating too many regimes in the model could also result in identification issues (see Section (ref)), we confine ourselves to two-regime models.
Based on graphical residual diagnostics, presented in Appendix (ref), the two-regime GSTVAR model with $p=4$ lags captures the autocorrelation structure of the data reasonably well. Conditional heteroskedasticity and marginal distribution of the series are not completely captured, but the inadequacies are not particularly severe. Finally, the sufficient Condition (ref) for ergodic stationarity holds for the selected model, as the upper bound of the joint spectral radius of the matrices $\boldsymbol{A}_m$ ((ref)), $m=1,2$ (see Section (ref)), is found to be strictly less than one ($0.995$). Hence, we proceed with the two-regime fourth-order model.
The estimated transition weights indicate the relative importance of each regime. The time series of the weights of Regime 1 are presented in the bottom panel of Figure (ref), and the corresponding weights of Regime 2 are, of course, obtained by subtracting these weights from unity. The shifts between the regimes are relatively fast, with the switch from one regime to the other typically taking from one to three months to be completed. Regime 1 mainly dominates in the 1960s, 1970s, and 1980s, but obtains large weights also during the Financial crisis and from the beginning of the COVID-19 crisis onward (excluding a short period in 2021). Overall, Regime 1 thereby mostly prevails in the earlier sample and Regime 2 in the later sample. The stationary standard deviations of ACI, GDP growth rate, CPI growth rate, and interest rate in Regime 1 are $0.59, 1.17, 0.47$, and $9.57$, respectively, whereas they are $0.52, 0.29, 0.26$, and $3.46$ in Regime 2. These differences and the dominance of Regime 1 in the volatile periods of 1970s and 1980s as well as during the Financial crisis and the COVID-19 crisis, suggest that it could represent more turbulent times compared to Regime 2.
We are predominantly interested in the economic impact of severe weather, and, therefore, a severe weather shock (ACI shock) must be identified. This amounts to imposing a unique structure of the impact matrix $B_t$ in ((ref)) governing the contemporaneous relationships of the shocks, so that one of the shocks can be labelled the ACI shock. Following Kim+Matthes+Phan:2022, our key identification restriction is that the ACI shock is the only shock that can instantaneously affect the ACI. This identification restriction seems reasonable, as it states that the other (macroeconomic) shocks do not affect severe weather or sea level within a month, but allows them to affect the ACI in the long run. We establish the identification by placing the ACI first in the vector of variables (with the rest of the variables in the order GDP, CPI, RATE) and imposing a recursive lower-triangular structure on the impact matrix $B_t$, obtained by the Cholesky decomposition of the estimated conditional error covariance matrix $\Omega_{y,t}$.
To study the effects of the ACI shock, we compute the generalized impulse response functions of the variables to it. As pointed out in Section (ref), since the transition weights are endogenously determined and the regime can shift as a result of a shock, the impulse responses generally depend on the initial values as well as on the sign and size of the shock. We are particularly interested in finding out about the potential state-dependence of the effects of the ACI shock. To capture them, we take all such histories of length $p$ from the data that indicate the dominance of a given regime, and compute the GIRFs for each regime separately. Specifically, for Regime $\tilde{m}$, we take the histories $\boldsymbol{y}_{t-1}$ for which the corresponding transition weight $\alpha_{\tilde{m},t}$ is greater than $0.75$. With the threshold $0.75$, the given regime is clearly dominant but a reasonably large amount of histories are still included in both regimes (Regimes 1 and 2 have $115$ and $603$ such histories, respectively). To closely match the properties of the data, for each history, the GIRF to the corresponding recovered structural shock is computed. Finally, the GIRFs are scaled so that the instantaneous response of the ACI is $0.3$, making the responses comparable.
Figure (ref) presents the sample distribution (over the histories) of the GIRFs of the variables, including the transition weights, to the identified ACI shock $h=0,1,...,36$ quarters ahead. Each light gray curve with low opacity depicts the GIRF corresponding to one history $\boldsymbol{y}_{t-1}$, so the darkness of a region in the figure indicates a greater concentration of GIRFs in that area. The left and right columns depict the responses in Regime 1 and Regime 2, respectively. In Regime 1, the vast majority of the GIRFs follow virtually the same sample paths. This is the case because the ACI shock has almost no effect on the transition weights of Regime 1, as shown in the bottom left panel. Although the majority of the GIRFs follow similar sample paths also in Regime 2, they exhibit more variation depending on the history $\boldsymbol{y}_{t-1}$ (and the related sign and size of the shock).
In both regimes, a positive ACI shock mainly decreases GDP, consumer prices, and the interest rate. However, there is a major difference in the GIRFs between the regimes: the effects of the ACI shock on all macroeconomic variables are clearly stronger in Regime 1. This may be the case because Regime 1 appears to be attributed to more volatile times, including recessions and crisis periods. Thus, it might represent a state of the economy that is more vulnerable to weather related shocks, and therefore, a similar shock is likely to result in greater damage in Regime 1 than in Regime 2. Other differences between the regimes are minor.
As Regime 1 dominates the earlier part and Regime 2 is predominant in the latter part of the sample period, it seems reasonable to conclude that the effect of severe weather shocks has diminished over time. This suggests that the U.S. economy has adapted to the effects of climate change. In contrast, Kim+Matthes+Phan:2022 found the effects of the ACI shock stronger at the end than at the beginning of their sample period from 1963:4 to 2019:5, which they interpreted as insufficient adaptation of the U.S. economy to the changing distribution of weather related shocks. While our results are thus, in general, contradictory to theirs, there are some short periods in the latter part of our sample where Regime 1 with stronger impact of the ACI shock prevails. These particularly overlap with the major crises, which suggests that the potential adaptation may not provide sufficient resilience against weather related shocks in a turbulent state of the economy.
To assess the relative importance (or economic significance) of the ACI shock, we compute the generalized forecast error variance decompositions using all the histories $\boldsymbol{y}_{t-1}$, $t=1,...,T$, in the data as the initial values with the corresponding sign and size of the recovered structural shock. The rows labeled “All data" in Table (ref) present the relative contribution of the ACI shock to the forecast error variance of the variable in question at each horizon, while the rows labeled “Regime 1" and “Regime 2" contain the corresponding figures for each regime (i.e., for the histories with a transition weight of the regime in question higher than $0.75$ in time period $t$). The ACI shock dominates the ACI index, but the other shocks become relatively more important at longer horizons. For the remaining variables, the ACI shock plays a lesser, yet clearly notable, role, contributing approximately 5--17%, 6--8% and 2--15% to the forecast error variance of GDP, CPI and the interest rate in the entire sample period, respectively. The ACI shock also seems to have some importance for regime shifts, with its relative contribution to the forecast error variance of the transition weight $\alpha_{1,t}$ hovering between 13 and 15%. The ACI shock is economically significant also in both regimes. Compared to Regime 1, its relative importance is, in general, greater for GDP and the transition weight in Regime 2, whereas the opposite holds for CPI and the interest rate, especially at short horizons.
To check the robustness of our findings with respect to the sample period, we also consider the subsample period from 1963:4 to 2019:5 considered by Kim+Matthes+Phan:2022. The results general conclusion remain intact in that the regimes of the fitted GSTVAR model quite similar to our main specification. Also, the estimated GIRFs and are similar, with the exception that the interest rate variable increases in response to a positive ACI shock in Regime 1. The details, including the evolution of the transition weights and the GIRFs can be found in Appendix (ref).
For comparison, we also report the results of the GMVAR model Kalliovirta+Meitz+Saikkonen:2016 as well as the recursively identified linear Gaussian SVAR model (obtained as special case of GMVAR and GSTVAR models with $M=1$) for the entire sample period. The mixing weights of the fitted GMVAR model, presented in Figure (ref) in Appendix (ref) are rather similar to the transition weights of our GSTVAR model, as expected due to the similar function forms of the weight functions. Also the GIRFs, presented in Figure (ref) in Appendix (ref), are quite similar to those in our GSTVAR model. In particular, in line with our GSTVAR model, the GMVAR model displays clearly stronger responses of GDP, inflation and the interest rate to the ACI shock in the more turbulent Regime 1 than in Regime 2. Finally, the finding that a positive ACI shock causes a decrease in all GDP, inflation and the interest is quite robust, as it emerges in all models considered, including the recursive Gaussian SVAR model, whose impulse responses are also depicted in Figure (ref) in Appendix (ref).
We have introduced a Gaussian smooth transition vector autoregressive (GSTVAR) model with transition weights adopted from the Gaussian mixture VAR (GMVAR) model of Kalliovirta+Meitz+Saikkonen:2016. However, compared to the GMVAR model, our model has the advantage of great flexibility in that it enables capturing gradual shifts in the dynamics of the data, whereas the GMVAR model involves discrete regimes. On the other hand, switching between the regimes in the $pt$h-order GSTVAR model depends on the full distribution of the preceding $p$ observations, in contrast to previously introduced other smooth transition VAR models, in which it is governed by specific switching variables. Hence, our model facilitates associating the regimes with statistical information in the data in a versatile manner. The structural counterpart of the GSTVAR model has the desirable feature that it produces the dynamic effects of the structural shocks on the transition weights governing regime switching.
We have discussed estimation, model selection and diagnostic checking as well as conducting structural analysis in the GSTVAR model. Moreover, by making use of the results of Saikkonen:2008 and Kheifets+Saikkonen:2020, we have introduced a sufficient condition for stationarity and ergodicity of the model. Due to the complexity of the log-likelihood function, estimation by the method of maximum likelihood can be quite tedious in practice, and following Virolainen:2022, Virolainen:2025 (and others), we propose a two-phase procedure in which a genetic algorithm is used to find starting values for a gradient based estimation method. The introduced methods, including a modified version of a genetic algorithm, are implemented in the CRAN distributed R package sstvars sstvars. A potential area of future research is an extension of our model that accommodates conditional heteroskedasticity by incorporating autoregressive conditional heteroskesticity or stochastic volatility to the shocks.
We have applied our method to study the macroeconomic effects of severe weather shocks in the U.S. Our monthly data set covers the period from 1961:1 to 2022:3 and contains an indicator of the frequency of severe weather and the extent of sea level rise in addition to a number of macroeconomic variables. Our structural GSTVAR model has two regimes, one of which (Regime 2) prevails in the latter part of the sample period, while the other regime (Regime 1) mainly dominates its earlier part, particularly the turbulent times of 1970s and 1980s, but also prevails in the later sample during the Financial crisis and the COVID-19 crisis. Kim+Matthes+Phan:2022 have recently also studied the economic impact of severe weather in the U.S., and following their lead, we have identified the severe weather shock recursively, placing the severe weather indicator first in the vector of variables. A positive weather shock was found to decrease GDP, consumer prices, and the interest rate in both regimes, but the effects are stronger in Regime 1. Hence, our results suggest that the U.S. economy, with the exception of certain crisis periods, has adapted to the changing distribution of weather related shocks. This finding contradicts the main result of Kim+Matthes+Phan:2022, who found the impact of the weather shock to get stronger over time, interpreting this as lack of adaptation to changing weather.