State-Space

1 Introduction

This article is looking at state-space models in macroeconomic forecasting, particularly a regime switching model. Regime switching models use the hypothesis that there are hidden states or characteristic equations that model the data generating process. Thus this is truly a subset of Hidden Markov Models with two states. An example would be a market rally vs normal market movement. Here, a hypothetical new equation could be used to model rally (as the behavior during these times is markedly different than during normal operation).

In the case of macroeconomic forecasting, the states one might be interested in would be expansion and contraction in the Gross Domestic Product (GDP). Knowing if we are moving into an expansion or recession hidden Markov state would be very interesting to the Federal Reserve in order to change policy in order to steer the economy. However, it must be mentioned that these types of models are hard. Why? Well we really don’t know if the two hidden states exist as we have modeled them (hence hidden). Additionally, these models have a problem with identifiability, i.e. we really don’t know if we can mathematically determine when the data are representing one hidden state versus another.

Regardless of these challenges, I will represent some of these models and the power of using Bayesian modeling for these types of problems. Here I’ll model these this process, largely borrowing from work already done by James Savage in his blog post(“Regime-Switching,” n.d.).

2 Data

First we are going to pull the GDP data using the ‘fredr’ package. We’ll get the GDP since 2010 and return the percent change quarter over quarter.

library(patchwork)
gdp_dat <- fredr(series_id = "GDPC1",
                                 observation_start = as.Date("2010-01-01"),
                                 units = "pch")

Always important when modeling data is to make sure that you visualize the data. Here we can build a simple graph to ensure that our data match our expectations and that we don’t see anything strange as shown in Figure ??.

(gdp_trend <- gdp_dat %>% 
    ggplot(aes(date, value))+
    geom_line(color = uncg_pal[[1]])+
    geom_hline(yintercept = 0, lty = "dashed", color = "red")+
    labs(
        title = "US GDP (% Change)",
        y = "%",
        caption = "Data: FRED St. Louis Federal Reserve"
    ))

Figure 2.1: US GDP Evolution Since Jan 1, 2010

Everything looks as expected, so we can move onto modeling.

library(highcharter)
hchart(gdp_dat, "line", hcaes(x = date, y = value),
       color = uncg_pal[[1]]) %>% 
    hc_title(
    text = "US GDP (% Change)",
    useHTML = TRUE) %>% 
    hc_tooltip(table = TRUE, sort = TRUE) %>% 
    hc_yAxis(tickPositions = c(-2, -1,0,1,2),
           gridLineColor = "#B71C1C",
           labels = list(format = "{value} %", useHTML = TRUE)) %>% 
    hc_credits(
    enabled = TRUE,
    text = "Source: FRED",
    href = "https://fred.stlouisfed.org/series/GDPC1") %>% 
    hc_add_theme(hc_theme_tufte2())

Figure 2.2: And an interactive looks at GDP

3 Method

In order to model our regime switching, state-space model, we will use Bayesian methods. Additionally, we will parameterize our model using the Stan language. In order to get started we will first load ‘rstan’ and set some general defaults to help us with our modeling.

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

3.1 Model

We will use a Bayesian modeling approach to model regime switching. Here there are two hypothesized regimes, a random walk state which will be represented as:

\[\text{GDP}_{t}\sim N(\alpha_{rw}, \sigma_{rw})\]

The second state will be a momentum state which will have an auto-regressive term as follows:

\[\text{GDP}_{t} \sim N(\alpha_{mo} + \rho*GDP_{t-1}, \sigma_{mo})\]

Additionally, we will have our different hidden states. As shown in Figure ??, there are four different transitions that can occur. A random walk can move to another random walk state or transition to a momentum state. This is similar for the momentum state which can transition out to random walk or stay in the momentum state.

[Random Walk]-[Momentum]
[Momentum]-[Random Walk]
[Momentum]-[Momentum]
[Random Walk]-[Random Walk]

3.1.1 Our Model

Now we can code our model using Stan.

writeLines(readLines("regime_switch.stan"))

//https://khakieconomics.github.io/2018/02/24/Regime-switching-models.html
data {
  int N;
  vector[N] y;
  int sticky;
}
parameters {
  vector<lower = 0, upper = 1>[2] p;
  real<lower = 0> rho;
  vector[2] alpha;
  vector<lower = 0>[2] sigma;
  real<lower = 0, upper = 1> xi1_init;
  real y_tm1_init;
}
transformed parameters {
  matrix[N, 2] eta;
  matrix[N, 2] xi;
  vector[N] f;

  // fill in etas
  for(t in 1:N) {
    eta[t,1] = exp(normal_lpdf(y[t]| alpha[1], sigma[1]));
    if(t==1) {
      eta[t,2] = exp(normal_lpdf(y[t]| alpha[2] + rho * y_tm1_init, sigma[2]));
    } else {
      eta[t,2] = exp(normal_lpdf(y[t]| alpha[2] + rho * y[t-1], sigma[2]));
    }
  }

  // work out likelihood contributions

  for(t in 1:N) {
    // for the first observation
    if(t==1) {
      f[t] = p[1]*xi1_init*eta[t,1] + // stay in state 1
             (1 - p[1])*xi1_init*eta[t,2] + // transition from 1 to 2
             p[2]*(1 - xi1_init)*eta[t,2] + // stay in state 2
             (1 - p[2])*(1 - xi1_init)*eta[t,1]; // transition from 2 to 1

      xi[t,1] = (p[1]*xi1_init*eta[t,1] +(1 - p[2])*(1 - xi1_init)*eta[t,1])/f[t];
      xi[t,2] = 1.0 - xi[t,1];

    } else {
    // and for the rest

      f[t] = p[1]*xi[t-1,1]*eta[t,1] + // stay in state 1
             (1 - p[1])*xi[t-1,1]*eta[t,2] + // transition from 1 to 2
             p[2]*xi[t-1,2]*eta[t,2] + // stay in state 2
             (1 - p[2])*xi[t-1,2]*eta[t,1]; // transition from 2 to 1

      // work out xi

      xi[t,1] = (p[1]*xi[t-1,1]*eta[t,1] +(1 - p[2])*xi[t-1,2]*eta[t,1])/f[t];

      // there are only two states so the probability of the other state is 1 - prob of the first
      xi[t,2] = 1.0 - xi[t,1];
    }
  }

}
model {
  // priors
  p ~ beta(sticky, 2);
  rho ~ normal(1, .1);
  alpha ~ normal(0, .1);
  sigma ~ cauchy(0, 1);
  xi1_init ~ beta(2, 2);
  y_tm1_init ~ normal(0, .1);


  // likelihood is really easy here!
  target += sum(log(f));
}

Now we can compile the model.

model_regime <- stan_model(file = "regime_switch.stan")

3.1.2 Our Priors

The priors set for the model should be derived from subject matter expertise. As I don’t have a ton of that regarding these types of problem, we can go with relatively uninformative priors.

One prior that I would like to talk about is the transition probability, which is parameterized as \(p\) in the model. This represents the probability at any given step that the state will transition to another hidden state. This is currently modeled as a beta distribution. Beta distributions have several nice properties, one of which is that they are bounded between 0 and 1 and thus make a nice prior for probabilities. What matters for this model is that the value of \(p\) indicates how “sticky” a given state is (i.e. how likely or probably a shift to another state is). Higher values indicate that the state is “sticky” and unlikely to change while lower values indicate that the system favors change.

Figure 3.1 shows different values of the first shape parameter (often called \(\alpha\)) for the distribution. As this value increases, the distribution becomes more and more skewed with higher probabilities favored. The other interesting thing about the Beta distribution is that if both values of the shape parameters are less than 1, the distribution becomes bathtub-shaped.

p <-  seq(0,1, length=100)
plot(p, dbeta(p, .5, .5), ylab="density", type ="l", col=4, ylim=c(0,10))
lines(p, dbeta(p, 2, 2), type ="l", col=5)
lines(p, dbeta(p, 5, 2), type ="l", col=3)
lines(p, dbeta(p, 10, 2), col=2) 
lines(p, dbeta(p, 20, 2), col=1) 
legend("topleft", c("Beta(.5,.5)","Beta(2,2)","Beta(5,2)","Beta(10,2)", "Beta(20,2)"),lty=c(1,1,1,1,1),col=c(4,5,3,2,1))

Figure 3.1: Beta Distribution for Different Values of Shape Parameters

I have elected to make this a parameter in our model so that we can see if the results change for stickier state changes.

3.2 Data Prep

Stan requires that we format our data in a list format with names exactly matching those in the Stan file.

stan_dat <- list(N = nrow(gdp_dat),
                                 y = dplyr::pull(gdp_dat, value),
                                 sticky = 5)

3.3 Running Model

Now that we have our model compiled, we can sample from the posterior distribution. It is always a good practice to set the seed for reproducibility. I am turning off the updates from Stan as I know the code run. If you want to see the progress, you can remove the refresh = 0 line.

fit_regime <- sampling(model_regime, data = stan_dat, 
                                             chains = 2, refresh = 0,
                                             seed = 42)

I am going to go ahead and sample from a model with a stickier probability function as well.

fit_regime_sticky <- sampling(model_regime, data = 
                                                                list(N = nrow(gdp_dat),
                                                                         y = dplyr::pull(gdp_dat, value),
                                                                         sticky = 20), 
                                             chains = 2, refresh = 0,
                                             seed = 42)

3.4 Diagnostics

As with all (Bayesian) modeling, it is important to assess that your model actually fits your data. Here we will look at the traceplots to verify that there aren’t any clearly diverging chains which would indicate a chain mixing/ convergence issue.

traceplot(fit_regime, pars = c("alpha", "rho", "sigma"))

Figure 3.2: Traceplots of Parameters

Additionally, we can look at the pairs plots of the parameters and confirm that we have general convergence.

pairs(fit_regime, pars = c("alpha","rho","sigma"))

Figure 3.3: Pairs Plot for Parameter Estimates

These results look pretty good, though the \(\sigma\)s for the momentum state look strange which indicates that we might not be modeling our data generating process well. We could go back and look at our parameterization where we have the cauchy distribution prior on \(\sigma\).

Finally, we can print the output of the model which shows our parameter estimate but don’t look at those yet! Here we look at our effective samples (n_eff) and our Rhat (Gelman-Rubin statistic) which help us understand how many effective samples we have from our posterior (higher is better) and our chain mixing (Rhat <= 1.1 is good).

print(fit_regime, pars = c("alpha", "rho", "p", "sigma"))

## Inference for Stan model: regime_switch.
## 2 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=2000.
## 
##          mean se_mean   sd  2.5%   25%  50%  75% 97.5% n_eff Rhat
## alpha[1] 0.21    0.00 0.12 -0.03  0.13 0.21 0.29  0.42  1011    1
## alpha[2] 0.02    0.00 0.07 -0.13 -0.02 0.02 0.07  0.15  1547    1
## rho      0.92    0.00 0.09  0.74  0.85 0.92 0.98  1.10  1123    1
## p[1]     0.77    0.00 0.12  0.49  0.69 0.78 0.86  0.95  1497    1
## p[2]     0.75    0.00 0.15  0.40  0.66 0.78 0.86  0.95  1255    1
## sigma[1] 0.59    0.00 0.15  0.37  0.48 0.56 0.66  0.97   981    1
## sigma[2] 0.27    0.01 0.18  0.10  0.21 0.26 0.31  0.44   554    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jan 21 16:49:08 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

All in all the model is fit.

4 Results

Now we can look at our parameter estimates as shown in Table ??.

broom::tidyMCMC(fit_regime, 
                                pars = c("alpha", "rho", "p", "sigma"),
                                conf.int = TRUE,
                                conf.level = .9) %>% 
    mutate_if(is.numeric, round, digits = 3) %>% 
    knitr::kable(caption = "Estimates and 90% Interval Estimates") %>% 
    kableExtra::kable_styling()

Table 4.1: Estimates and 90% Interval Estimates

term	estimate	std.error	conf.low	conf.high
alpha[1]	0.207	0.119	0.005	0.390
alpha[2]	0.020	0.071	-0.106	0.134
rho	0.915	0.092	0.765	1.068
p[1]	0.765	0.123	0.538	0.940
p[2]	0.747	0.147	0.457	0.935
sigma[1]	0.588	0.155	0.394	0.881
sigma[2]	0.273	0.185	0.137	0.403

We see that during the momentum state that there is a strong auto-regressive parameter with \(\rho\) equal to nearly one. As we saw in our convergence diagnostics, \(\sigma\) for our momentum state has a very large standard error which is concerning.

4.1 Interpretation

It’s best to visualize the output of this model.

state_summary <- as.data.frame(fit_regime, pars = "xi") %>% 
  gather(par, value) %>% 
  filter(grepl(",1", par)) %>%
  mutate(time = parse_number(str_extract(par, "[0-9]{1,3}[,]"))) %>% 
  group_by(par) %>% 
  summarise(time = first(time), 
            mean = mean(value),
            lower = quantile(value, .05),
            upper = quantile(value, .95)) %>% 
  add_column(date = gdp_dat$date) %>% 
  ggplot(aes(x = date, y = mean)) +
    geom_line(color = uncg_pal[[1]])+
    geom_hline(yintercept = .5, color = uncg_pal[[5]], lty = "dashed")+
  geom_ribbon(aes(ymin = lower, ymax = upper), 
                        alpha = 0.3,  fill = uncg_pal[[2]]) +
  labs(title ="Probability of being in random walk state\n vs momentum state",
       y = "Probability",
       subtitle = "US GDP (90% Credibility Intervals Shown)")

state_summary/gdp_trend

Figure 4.1: State Probability and Actual GDP Evolution

So here we can see that our model isn’t great. There are only a few periods in which the probability of being in a random walk is high (sometime around Q1 2012 and since the end of 2018). This could indicate that the current growth we are seeing now is in a random walk random state while previous periods were probably in this rally state. This might be a harbinger of a slow down, but this model doesn’t fit the data that well.

5 Next Steps

The fun part with these types of models is that you don’t have to stop here. Two states may not be appropriate. Additionally, there is definitely more information in the world than just GDP. This model could be enhanced by using different times series models as well as bring in additional data to use to help forecast GDP and identify the the hidden state. We also saw the issue with identifiability with our model not really being able to discern distinct states.

6 And what about the sticker model?

The somewhat sticker model indicates very similar results to the less stickier prior. No big changes other than a slightly higher probability of a random walk state 2011-2012.

as.data.frame(fit_regime_sticky, pars = "xi") %>% 
  gather(par, value) %>% 
  filter(grepl(",1", par)) %>%
  mutate(time = parse_number(str_extract(par, "[0-9]{1,3}[,]"))) %>% 
  group_by(par) %>% 
  summarise(time = first(time), 
            mean = mean(value),
            lower = quantile(value, .05),
            upper = quantile(value, .95)) %>% 
  add_column(date = gdp_dat$date) %>% 
  ggplot(aes(x = date, y = mean)) +
    geom_line(color = uncg_pal[[1]])+
    geom_hline(yintercept = .5, color = uncg_pal[[5]], lty = "dashed")+
  geom_ribbon(aes(ymin = lower, ymax = upper), 
                        alpha = 0.3,  fill = uncg_pal[[2]]) +
  labs(title ="Probability of being in random walk state\n vs momentum state",
       y = "Probability",
       subtitle = "US GDP (90% Credibility Intervals Shown)")

Figure 6.1: State Probability with Stickier Prior

References

“Regime-Switching.” n.d. https://khakieconomics.github.io/2018/02/24/Regime-switching-models.html.

---

Intermediate Macroeconomics
me.dewitt.jr@gmail.com

Code licensed under the BSD 3-clause license
Text licensed under the CC-BY-ND-NC 4.0 license

See the Repository

Copyright © 2020 Michael DeWitt