## Bayesian threshold autoregressive models

Autoregressive (AR) models are some of the most widely used models in applied economics, among other disciplines, because of their generality and simplicity. However, the dynamic characteristics of real economic and financial data can change from one time period to another, limiting the applicability of linear time-series models. For example, the change of unemployment rate is a function of the state of the economy, whether it is expanding or contracting. A variety of models have been developed that allow time-series dynamics to depend on the regime of the system they are part of. The class of regime-dependent models include Markov-switching, smooth transition, and threshold autoregressive (TAR) models.

TAR (Tong 1982) is a class of nonlinear time-series models with applications in econometrics (Hansen 2011), financial analysis (Cao and Tsay 1992), and ecology (Tong 2011). TAR models allow regime-switching to be triggered by the observed level of an outcome in the past. The **threshold** command in Stata provides frequentist estimation of some TAR models.

In Stata 17, the **bayesmh** command supports time-series operators in linear and nonlinear model specifications; see **[BAYES] bayesmh**. In this blog entry, I want to show how we can fit some Bayesian TAR models using the **bayesmh** command. The examples will also demonstrate modeling flexibility not possible with the existing **threshold** command.

**TAR model definition**

Let \(\{y_t\}\) be a series observed at discrete times \(t\). A general TAR model of order \(p\), TAR(\(p\)), with \(k\) regimes has the form

\[

y_t = a_0^{j} + a_1^{j} y_{t-1} + \dots + a_p^{j} y_{t-p} + \sigma_{j} e_t, \quad {\rm if} \quad

r_{j-1} < z_t \le r_{j}

\]

where \(z_t\) is a threshold variable, \(-\infty < r_0 < \dots < r_k =

\infty\) are regime thresholds, and \(e_t\) are independent standard normally distributed errors. The \(j\)th regime has its own set of autoregressive coefficients \(\{a_i^j\}\)’s and standard deviations \(\sigma_j\)’s. Different regimes are also allowed to have different orders \(p\). In the above equation, this can be indicated by replacing \(p\) with regime-specific \(p_j\)’s.

In a TAR model, as implied by the definition, structural breaks happen not at certain time points but are triggered by the magnitude of the threshold variable \(z_t\). It is common to have \(z_t = y_{t-d}\), where \(d\) is a positive integer the called the delay parameter. These models are known as self-exciting TAR (SETAR) and are the ones I am illustrating below.

**Real GDP dataset**

In Beaudry and Koop (1993), TAR models were used to model gross national product. The authors demonstrated asymmetric persistence in the growth rate of gross national product, with positive shocks (associated with expansion periods) being more persistent than negative shocks (recession periods). In a similar approach, I use the growth rate of real gross domestic product (GDP) of the United States as an indicator of the status of the economy to model the difference between expansion and recession periods.

Quarterly observations on real GDP, measured in billions of dollars, are obtained from the Federal Reserve Economic Data repository using the **import fred** command. I consider observations only between the first quarter of 1947 and second quarter of 2021. A quarterly date variable, **dateq**, is generated and used with **tsset** to set up the time series.

. import fred GDPC1 Summary ------------------------------------------------------------------------------- Series ID Nobs Date range Frequency ------------------------------------------------------------------------------- GDPC1 301 1947-01-01 to 2022-01-01 Quarterly ------------------------------------------------------------------------------- # of series imported: 1 highest frequency: Quarterly lowest frequency: Quarterly . keep if daten >= date("01jan1947", "DMY") & daten <= > date("01apr2021", "DMY") (3 observations deleted) . generate dateq = qofd(daten) . tsset dateq, quarterly Time variable: dateq, 1947q1 to 2021q2 Delta: 1 quarter

I am interested in the change of real GDP from one quarter to the next. For this purpose, I generate a new variable, **rgdp**, to measure this change in percentages. Positive values of **rgdp** indicate economic growth or expansion, while values close to 0 or negative are associated with stagnation or recession. Below, I use the **tsline** command to plot the time series. There, some of the recession periods, indicated by sharp drops in GDP, are visible, including the latest one in 2020.

. generate double rgdp = 100*D1.GDPC1/L1.GDPC1 (1 missing value generated) . tsline rgdp

A TAR model with two regimes estimates one threshold value \(r\) that can be visualized as a horizontal line separating, somewhat informally, expansion from recession periods. In Bayesian TAR, the threshold \(r\) is a random variable with distribution estimated from a prior and observed data.

**Bayesian TAR specification**

Before I show how to specify a Bayesian TAR model in Stata, let me first fit a simpler Bayesian AR(1) model for **rgdp** using the **bayesmh** command. It will serve as a baseline for comparison with models with structural breaks.

I use the fairly uninformative, given the range of **rgdp**, **normal(0, 100)** prior for the two coefficients in the **{rgdp:}** equation and the **igamma(0.01, 0.01)** prior for the variance parameter **{sig2}**. I also use Gibbs sampling for more efficient simulation of the model parameters.

. bayesmh rgdp L1.rgdp, likelihood(normal({sig2})) > prior({rgdp:}, normal(0, 100)) block({rgdp:}, gibbs) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}, gibbs) > rseed(17) dots Burn-in 2500 .........1000.........2000..... done Simulation 10000 .........1000.........2000.........3000.........4000......... > 5000.........6000.........7000.........8000.........9000.........10000 done Model summary ------------------------------------------------------------------------------ Likelihood: rgdp ~ normal(xb_rgdp,{sig2}) Priors: {rgdp:L.rgdp _cons} ~ normal(0,100) (1) {sig2} ~ igamma(0.01,0.01) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_rgdp. Bayesian normal regression MCMC iterations = 12,500 Gibbs sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 296 Acceptance rate = 1 Efficiency: min = .9584 avg = .9755 Log marginal-likelihood = -478.07327 max = 1 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. dev. MCSE Median [95% cred. interval] -------------+---------------------------------------------------------------- rgdp | rgdp | L1. | .1239926 .0579035 .000591 .1240375 .0096098 .237592 | _cons | .6762768 .0805277 .000805 .676687 .5186598 .8378632 -------------+---------------------------------------------------------------- sig2 | 1.344712 .1098713 .001117 1.339746 1.144444 1.571802 ------------------------------------------------------------------------------

The posterior mean estimate for the AR(1) coefficient **{L1.rgdp}** is 0.12, indicating a positive serial correlation for **rgdp**. This implies some degree of persistence in the real GDP growth. The posterior mean estimate of 1.34 for **{sig2}** suggests a volatility level above 1 percent, if the latter is measured by standard deviation. However, this simple AR(1) model does not tell us how the persistence and volatility change depending on the state of the economy.

Before continuing, I save the simulation and estimation results for later reference.

. bayesmh, saving(bar1sim, replace) note: file bar1sim.dta saved. . estimates store bar1

To describe a two-state model for the economy, I want to specify the simplest two-regime SETAR model with order of 1 and delay of 1. In the next section, I discuss the choice of order and delay parameters.

The model can be summarized by the following two equations:

\begin{align}

{\bf rgdp_t} = a_0^{1} + a_1^{1} {\bf rgdp_{t-1}} + \sigma_{1} e_t, \quad {\rm if} \quad

y_{t-1} < r \\

{\bf rgdp_t} = a_0^{2} + a_1^{2} {\bf rgdp_{t-1}} + \sigma_{2}

e_t, \quad {\rm if} \quad y_{t-1} \ge r

\end{align}

To specify the regression portion of this model with **bayesmh**, I use a substitutable expression with conditional logic,

**cond(L1.rgdp<{r}, {r1:a0}+{r1:a1}*L1.rgdp, {r2:a0}+{r2:a1}*L1.rgdp)**

where **{r1:a0}** and **{r1:a1}** are the coefficients for the first regime and **{r2:a0}** and **{r2:a1}** are the coefficients for the second regime.

The regime-specific variance of the normal likelihood can be similarly

specified by the expression

**cond(L1.rgdp<{r},{sig1},{sig2})**

Instead of assuming a fixed threshold value for \(r\), with 0 being a natural choice, I consider \(r\) to be a hyperparameter with **uniform(-0.5, 0.5)** prior. I thus assume that the threshold is within a half percentage point of 0. Given the range of **rgdp** and that 0 separates positive from negative growth, this seems to be a reasonable assumption. Using uninformative priors for \(r\) without any restrictions on its range will make the model unstable because of the possibility of collapsing one of the regimes, that is, having a regime with 0 or only a few observed points. The priors for coefficients and variances stay the same as in the previous model.

. bayesmh rgdp = (cond(L1.rgdp<{r}, > {r1:a0}+{r1:a1}*L1.rgdp, > {r2:a0}+{r2:a1}*L1.rgdp)), > likelihood(normal(cond(L1.rgdp<{r}, {sig1}, {sig2}))) > prior({r1:}, normal(0, 100)) block({r1:}) > prior({r2:}, normal(0, 100)) block({r2:}) > prior({sig1}, igamma(0.01, 0.01)) block({sig1}) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}) > prior({r}, uniform(-0.5, 0.5)) block({r}) > rseed(17) init({sig1} {sig2} 1) dots Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done Simulation 10000 .........1000.........2000.........3000.........4000......... > 5000.........6000.........7000.........8000.........9000.........10000 done Model summary ------------------------------------------------------------------------------ Likelihood: rgdp ~ normal(<expr1>,<expr2>) Priors: {r1:a0 a1} ~ normal(0,100) {r2:a0 a1} ~ normal(0,100) {sig1 sig2} ~ igamma(0.01,0.01) {r} ~ uniform(-0.5,0.5) Expressions: expr1 : cond(L1.rgdp<{r},{r1:a0}+{r1:a1}*L1.rgdp,{r2:a0}+{r2:a1}*L1.rgdp) expr2 : cond(L1.rgdp<{r},{sig1},{sig2}) ------------------------------------------------------------------------------ Bayesian normal regression MCMC iterations = 12,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 296 Acceptance rate = .3554 Efficiency: min = .04586 avg = .1205 Log marginal-likelihood = -415.46111 max = .2235 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. dev. MCSE Median [95% cred. interval] -------------+---------------------------------------------------------------- r1 | a0 | -1.327623 .7508926 .020208 -1.284013 -2.794005 .2010012 a1 | -.8866186 .3328065 .009604 -.8828538 -1.582993 -.2458004 -------------+---------------------------------------------------------------- r2 | a0 | .5521038 .0739373 .002338 .5516799 .4093326 .6992918 a1 | .3015663 .0555794 .001702 .3020909 .1883468 .4070365 -------------+---------------------------------------------------------------- r | -.4834577 .0205573 .00096 -.4903714 -.4995296 -.416813 sig1 | 6.752686 2.201104 .066593 6.356055 3.684466 12.27644 sig2 | .678672 .0580648 .001228 .676249 .5714164 .7989442 ------------------------------------------------------------------------------

The model takes less than a minute to run. There are no obvious convergence problems reported by **bayesmh**, and the average sampling efficiency of 12% is good.

The threshold parameter is estimated to be about \(-0.48\). Although this is close to the lower limit of \(-0.5\) set by the prior, it is still strictly greater than \(-0.5\) because of its high precision—MCSE is less than 0.001.

The autoregression coefficients are negative in the first regime, **r1**, and positive in the second, **r2**. The second regime thus has much higher persistency. Also notable is the much higher variability in the first regime, about 6.75, in comparison with the second, 0.68.

I save the last estimation results and use the **bayesstats ic** command to compare the SETAR(1) and baseline AR(1) models.

. bayesmh, saving(bster1sim, replace) note: file bster1sim.dta not found; file saved. . estimates store bster1 . bayesstats ic bar1 bster1 Bayesian information criteria ---------------------------------------------- | DIC log(ML) log(BF) -------------+-------------------------------- bar1 | 929.311 -478.0733 . bster1 | 782.0817 -415.4611 62.61216 ---------------------------------------------- Note: Marginal likelihood (ML) is computed using Laplace–Metropolis approximation.

The SETAR(1) model has a lower DIC and higher log-marginal likelihood than the AR(1) model. Of course, we expect the more complex and flexible SETAR(1) model to provide a better fit based on the likelihood alone. Note, however, that the marginal likelihood incorporates, in addition to the likelihood, the priors on model parameters and thus, indirectly, model complexity as well.

For comparison, I also perform frequentist estimation of the same model using the **threshold** command.

. threshold rgdp, regionvars(l.rgdp) threshvar(l1.rgdp) Searching for threshold: 1 (running 237 regressions) .................................................. 50 .................................................. 100 .................................................. 150 .................................................. 200 ..................................... Threshold regression Number of obs = 296 Full sample: 1947q3 thru 2021q2 AIC = 30.1871 Number of thresholds = 1 BIC = 44.9485 Threshold variable: L.rgdp HQIC = 36.0973 --------------------------------- Order Threshold SSR --------------------------------- 1 -0.3881 319.0398 --------------------------------- ------------------------------------------------------------------------------ rgdp | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- Region1 | rgdp | L1. | -.7989782 .12472 -6.41 0.000 -1.043425 -.5545315 | _cons | -.9796166 .2473612 -3.96 0.000 -1.464436 -.4947977 -------------+---------------------------------------------------------------- Region2 | rgdp | L1. | .2910881 .0755517 3.85 0.000 .1430094 .4391667 | _cons | .5663334 .0987315 5.74 0.000 .3728231 .7598436 ------------------------------------------------------------------------------

The estimates for regression coefficients are similar to the Bayesian ones: negative in the first regime and positive in the second regime. The threshold estimate is \(-0.39\), somewhat higher than the posterior mean estimate in the Bayesian model. A limitation of the **threshold** command is the lack of error-variance estimates for the two regimes.

**Autoregression order selection**

In the previous example, I fit the simplest SETAR model of order 1 and delay of 1. In general, these parameters are unknown, and one may not have a good prior choice for them. One solution is to fit models of different orders and compare them. A better solution is to consider one Bayesian model in which the orders are included as hyperparameters and are thus estimated along with all other parameters.

The following extension of the previous Bayesian model considers as possible choices orders from 1 to 4 for each regime. Two additional discrete hyperparameters, **p1** and **p2**, indicate the regime orders. Both regimes are assumed to be at least of order 1. These hyperparameters thus take values in the set \(\{1,2,3,4\}\) according to some prior probabilities. I use the **index(0.2,0.5,0.2,0.1)** prior to set my highest expectation, 0.5, on order 2, then equal probabilities of 0.2 on orders 1 and 3, and finally probability of 0.1 on order 4. Orders 2, 3, and 4 are turned on and off using indicator variables as multipliers to the coefficients **b2**, **b3**, and **b4**, separately for each regime.

. bayesmh rgdp = (cond(L1.rgdp<{r}, > {r1:a0} + {r1:a1}*L1.rgdp + ({p1}>1)*{r1:b2}*L2.rgdp + > ({p1}>2)*{r1:b3}*L3.rgdp + ({p1}>3)*{r1:b4}*L4.rgdp, > {r2:a0} + {r2:a1}*L1.rgdp + ({p2}>1)*{r2:b2}*L2.rgdp + > ({p2}>2)*{r2:b3}*L3.rgdp + ({p2}>3)*{r2:b4}*L4.rgdp)), > likelihood(normal(cond(L1.rgdp<{r}, {sig1}, {sig2}))) > prior({p1}, index(0.2,0.5,0.2,0.1)) block({p1}) > prior({p2}, index(0.2,0.5,0.2,0.1)) block({p2}) > prior({r1:}, normal(0, 100)) block({r1:}) > prior({r2:}, normal(0, 100)) block({r2:}) > prior({sig1}, igamma(0.01, 0.01)) block({sig1}) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}) > prior({r}, uniform(-0.5, 0.5)) block({r}) > rseed(17) init({sig1} {sig2} 1 {p1} {p2} 2) dots Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done Simulation 10000 .........1000.........2000.........3000.........4000......... > 5000.........6000.........7000.........8000.........9000.........10000 done Model summary ------------------------------------------------------------------------------ Likelihood: rgdp ~ normal(<expr1>,<expr2>) Priors: {p1 p2} ~ index(0.2,0.5,0.2,0.1) {r1:a0 a1 b2 b3 b4} ~ normal(0,100) {r2:a0 a1 b2 b3 b4} ~ normal(0,100) {sig1 sig2} ~ igamma(0.01,0.01) {r} ~ uniform(-0.5,0.5) Expressions: expr1 : cond(L1.rgdp<{r},{r1:a0} + {r1:a1}*L1.rgdp + ({p1}>1)*{r1:b2}*L2.rgd p + ({p1}>2)*{r1:b3}*L3.rgdp + ({p1}>3)*{r1:b4}*L4.rgdp,{r2:a0} + {r2:a1}*L1.rgdp + ({p2}>1)*{r2:b2}*L2.rgdp + ({p2}>2)*{r2:b3}*L3.rgd p + ({p2}>3)*{r2:b4}*L4.rgdp) expr2 : cond(L1.rgdp<{r},{sig1},{sig2}) ------------------------------------------------------------------------------ Bayesian normal regression MCMC iterations = 12,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 293 Acceptance rate = .3534 Efficiency: min = .0167 avg = .04996 Log marginal-likelihood = -415.62492 max = .2163 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. dev. MCSE Median [95% cred. interval] -------------+---------------------------------------------------------------- r1 | a0 | -1.382746 .7689649 .043474 -1.395605 -2.953167 .0933557 a1 | -.8994067 .3164027 .021135 -.8966563 -1.53448 -.2834772 b2 | -.6850185 8.58136 .561206 .0350022 -18.77976 16.92638 b3 | -1.115146 9.72345 .595084 -.5968546 -20.72936 17.06076 b4 | .1783556 10.30925 .572088 -.0286035 -18.84217 21.22304 -------------+---------------------------------------------------------------- r2 | a0 | .4381789 .0809618 .003256 .4369241 .2837359 .6029897 a1 | .3064629 .0549879 .002723 .3078107 .1969199 .409214 b2 | .1311621 .0531822 .004115 .1343153 .0269019 .2253627 b3 | -.2968566 9.603545 .515804 -.1204644 -19.08613 18.83395 b4 | -.9700926 9.811401 .462162 -1.123427 -20.09342 18.72596 -------------+---------------------------------------------------------------- p1 | 1.1602 .3812486 .018329 1 1 2 p2 | 1.9858 .1902682 .012683 2 1 2 r | -.4845632 .0183435 .000932 -.4903276 -.4994973 -.4211905 sig1 | 6.850306 2.403665 .078689 6.427896 3.629047 12.61115 sig2 | .6557395 .0570403 .001226 .6538532 .5533153 .7733018 ------------------------------------------------------------------------------

The model takes about 2 minutes to run and has a good average sampling efficiency of 5%. Posterior median estimates for the order parameters are 1 for the first regime and 2 for the second. We saw that the first regime is more volatile. During recessions, having shorter order is consistent with having higher volatility.

Note that the parameters **b2**, **b3**, and **b4** are not the actual autocorrelation coefficients for the series. To summarize the autoregression coefficients for the first regime, **r1**, we need to include the order indicators for **p1** from the model specification.

. bayesstats summary {r1:a0} {r1:a1} (a2:({p1}>1)*{r1:b2}) > (a3:({p1}>2)*{r1:b3}) (a4:({p1}>3)*{r1:b4}) Posterior summary statistics MCMC sample size = 10,000 a2 : ({p1}>1)*{r1:b2} a3 : ({p1}>2)*{r1:b3} a4 : ({p1}>3)*{r1:b4} ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. dev. MCSE Median [95% cred. interval] -------------+---------------------------------------------------------------- r1 | a0 | -1.382746 .7689649 .043474 -1.395605 -2.953167 .0933557 a1 | -.8994067 .3164027 .021135 -.8966563 -1.53448 -.2834772 -------------+---------------------------------------------------------------- a2 | .0517708 .2540154 .013244 0 -.2375046 .9287818 a3 | -.0014845 .038323 .000658 0 0 0 a4 | 0 0 0 0 0 0 ------------------------------------------------------------------------------

The autocorrelation estimates for orders 2 through 4 are very close to 0, as we expect given that the estimate for **p1** is 1.

Similarly, the autoregression coefficients for the second regime have essentially estimates of 0 for orders 3 and 4.

. bayesstats summary {r2:a0} {r2:a1} (a2:({p2}>1)*{r2:b2}) > (a3:({p2}>2)*{r2:b3}) (a4:({p2}>3)*{r2:b4}) Posterior summary statistics MCMC sample size = 10,000 a2 : ({p2}>1)*{r2:b2} a3 : ({p2}>2)*{r2:b3} a4 : ({p2}>3)*{r2:b4} ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. dev. MCSE Median [95% cred. interval] -------------+---------------------------------------------------------------- r2 | a0 | .4381789 .0809618 .003256 .4369241 .2837359 .6029897 a1 | .3064629 .0549879 .002723 .3078107 .1969199 .409214 -------------+---------------------------------------------------------------- a2 | .1306131 .0480967 .003019 .1336667 0 .2179723 a3 | -.0008258 .0089415 .000611 0 0 0 a4 | 0 0 0 0 0 0 ------------------------------------------------------------------------------

The delay \(d\) is another important parameter in SETAR models. Thus far, we considered a delay of one quarter period, which may not be optimal. Although it is possible to incorporate \(d\) as a hyperparameter in a single Bayesian model similarly to what I have done with the order parameters, to avoid an overly complicated specification, I run three additional models with \(d=2\), \(d=3\), and \(d=4\) by using **L2.rgdp**, **L3.rgdp**, and **L4.rgdp**, respectively, as threshold variables and compare them with the model with \(d=1\).

. bayesmh rgdp = (cond(L2.rgdp<{r}, > {r1:a0} + {r1:a1}*L1.rgdp + ({p1}>1)*{r1:b2}*L2.rgdp + > ({p1}>2)*{r1:b3}*L3.rgdp + ({p1}>3)*{r1:b4}*L4.rgdp, > {r2:a0} + {r2:a1}*L1.rgdp + ({p2}>1)*{r2:b2}*L2.rgdp + > ({p2}>2)*{r2:b3}*L3.rgdp + ({p2}>3)*{r2:b4}*L4.rgdp)), > likelihood(normal(cond(L2.rgdp<{r}, {sig1}, {sig2}))) > prior({p1}, index(0.2,0.5,0.2,0.1)) block({p1}) > prior({p2}, index(0.2,0.5,0.2,0.1)) block({p2}) > prior({r1:}, normal(0, 100)) block({r1:}) > prior({r2:}, normal(0, 100)) block({r2:}) > prior({sig1}, igamma(0.01, 0.01)) block({sig1}) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}) > prior({r}, uniform(0, 1)) block({r}) > rseed(17) init({sig1} {sig2} 1 {p1} {p2} 2) > burnin(5000) nomodelsummary notable Burn-in ... Simulation ... Bayesian normal regression MCMC iterations = 15,000 Random-walk Metropolis–Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 293 Acceptance rate = .3544 Efficiency: min = .01077 avg = .05785 Log marginal-likelihood = -453.03074 max = .1904 . bayesmh rgdp = (cond(L3.rgdp<{r}, > {r1:a0} + {r1:a1}*L1.rgdp + ({p1}>1)*{r1:b2}*L2.rgdp + > ({p1}>2)*{r1:b3}*L3.rgdp + ({p1}>3)*{r1:b4}*L4.rgdp, > {r2:a0} + {r2:a1}*L1.rgdp + ({p2}>1)*{r2:b2}*L2.rgdp + > ({p2}>2)*{r2:b3}*L3.rgdp + ({p2}>3)*{r2:b4}*L4.rgdp)), > likelihood(normal(cond(L3.rgdp<{r}, {sig1}, {sig2}))) > prior({p1}, index(0.2,0.5,0.2,0.1)) block({p1}) > prior({p2}, index(0.2,0.5,0.2,0.1)) block({p2}) > prior({r1:}, normal(0, 100)) block({r1:}) > prior({r2:}, normal(0, 100)) block({r2:}) > prior({sig1}, igamma(0.01, 0.01)) block({sig1}) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}) > prior({r}, uniform(0, 1)) block({r}) > rseed(17) init({sig1} {sig2} 1 {p1} {p2} 2) > burnin(5000) nomodelsummary notable Burn-in ... Simulation ... Bayesian normal regression MCMC iterations = 15,000 Random-walk Metropolis–Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 293 Acceptance rate = .338 Efficiency: min = .006822 avg = .0667 Log marginal-likelihood = -472.66834 max = .2068 . bayesmh rgdp = (cond(L4.rgdp<{r}, > {r1:a0} + {r1:a1}*L1.rgdp + ({p1}>1)*{r1:b2}*L2.rgdp + > ({p1}>2)*{r1:b3}*L3.rgdp + ({p1}>3)*{r1:b4}*L4.rgdp, > {r2:a0} + {r2:a1}*L1.rgdp + ({p2}>1)*{r2:b2}*L2.rgdp + > ({p2}>2)*{r2:b3}*L3.rgdp + ({p2}>3)*{r2:b4}*L4.rgdp)), > likelihood(normal(cond(L4.rgdp<{r}, {sig1}, {sig2}))) > prior({p1}, index(0.2,0.5,0.2,0.1)) block({p1}) > prior({p2}, index(0.2,0.5,0.2,0.1)) block({p2}) > prior({r1:}, normal(0, 100)) block({r1:}) > prior({r2:}, normal(0, 100)) block({r2:}) > prior({sig1}, igamma(0.01, 0.01)) block({sig1}) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}) > prior({r}, uniform(0, 1)) block({r}) > rseed(17) init({sig1} {sig2} 1 {p1} {p2} 2) > burnin(5000) nomodelsummary notable Burn-in ... Simulation ... Bayesian normal regression MCMC iterations = 15,000 Random-walk Metropolis–Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 293 Acceptance rate = .3749 Efficiency: min = .003091 avg = .03948 Log marginal-likelihood = -484.88072 max = .1626

To save space, I show only the estimated log-marginal likelihoods of the models,

\(d = 1\) | \(d = 2 \) | \(d = 3 \) | \(d = 4 \) |

\(-416\) | \( -453 \) | \(-473 \) | \(-485 \) |

A delay of 1 gives us the highest log-marginal likelihood, thus validating our initial choice.

**Final model**

Here is our final model, which seems to provide the best analysis of the dynamics of **rgdp**.

. bayesmh rgdp = (cond(L1.rgdp<{r}, > {r1:a0}+{r1:a1}*L1.rgdp, > {r2:a0}+{r2:a1}*L1.rgdp+{r2:a2}*L2.rgdp)), > likelihood(normal(cond(L1.rgdp<{r}, {sig1}, {sig2}))) > prior({r1:}, normal(0, 100)) block({r1:}) > prior({r2:}, normal(0, 100)) block({r2:}) > prior({sig1}, igamma(0.01, 0.01)) block({sig1}) > prior({sig2}, igamma(0.01, 0.01)) block({sig2}) > prior({r}, uniform(-0.5, 0.5)) block({r}) > rseed(17) init({sig1} {sig2} 1) dots Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done Simulation 10000 .........1000.........2000.........3000.........4000......... > 5000.........6000.........7000.........8000.........9000.........10000 done Model summary ------------------------------------------------------------------------------ Likelihood: rgdp ~ normal(<expr1>,<expr2>) Priors: {r1:a0 a1} ~ normal(0,100) {r2:a0 a1 a2} ~ normal(0,100) {sig1 sig2} ~ igamma(0.01,0.01) {r} ~ uniform(-0.5,0.5) Expressions: expr1 : cond(L1.rgdp<{r},{r1:a0}+{r1:a1}*L1.rgdp,{r2:a0}+{r2:a1}* L1.rgdp+{r2:a2}*L2.rgdp) expr2 : cond(L1.rgdp<{r},{sig1},{sig2}) ------------------------------------------------------------------------------ Bayesian normal regression MCMC iterations = 12,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 295 Acceptance rate = .3497 Efficiency: min = .04804 avg = .09848 Log marginal-likelihood = -414.93784 max = .1997 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. dev. MCSE Median [95% cred. interval] -------------+---------------------------------------------------------------- r1 | a0 | -1.269802 .7325285 .024046 -1.268012 -2.746194 .2098139 a1 | -.858765 .3224316 .009838 -.8566081 -1.48599 -.1966072 -------------+---------------------------------------------------------------- r2 | a0 | .4496805 .0830125 .002535 .4501146 .2868064 .6120813 a1 | .2980119 .0562405 .001979 .2955394 .1915367 .412661 a2 | .1302317 .0417601 .001504 .1285035 .0465193 .2122086 -------------+---------------------------------------------------------------- r | -.4831157 .0213086 .000972 -.4905123 -.4996558 -.4155153 sig1 | 6.747231 2.396798 .087641 6.255731 3.666464 12.7075 sig2 | .6563023 .055926 .001251 .6510408 .5575998 .7745131 ------------------------------------------------------------------------------

In conclusion, the expansion state, **r2**, is characterized by positive trend and autocorrelation, comparatively higher persistency, and lower volatility. The recession state, **r1**, on the other hand, experiences negative trend and autocorrelation, and higher volatility.

Although SETAR(1) provides a much more detailed analysis than a simple AR(1) model, it still does not capture all the changes in the dynamics of GDP growth. For example, the expansion periods before 1985 seem to have much higher volatility than those after 1985. Alternative regime-switching models may need to be considered to address this and other aspects of the time evolution of economic growth.

**References**

Beaudry, P., and G. Koop. 1993. Do recessions permanently change output? *Journal of Monetary Economics* 31: 149–163. https://doi.org/10.1016/0304-3932(93)90042-E.

Cao, C. Q., and R. S. Tsay. 1992. Nonlinear time-series analysis of stock volatilities. *Journal of Applied Econometrics* 7: S165–S185. https://doi.org/10.1002/jae.3950070512.

Hansen, B. E. 2011. Threshold autoregression in economics. *Statistics and Its Inference* 4: 123–127. https://doi.org/10.4310/SII.2011.v4.n2.a4.

Tong, H. 1982. Discontinuous decision processes and threshold autoregressive time series modelling. *Biometrica* 69: 274–276. https://doi.org/10.2307/2335885.

——. 2011. Threshold models in time series analysis—30 years on. *Statistics and Its Inference* 4: 107–118. https://dx.doi.org/10.4310/SII.2011.v4.n2.a1.