Chapter 20: Multivariate Time Series Methods

Franz X. Mohr, Created: December 9, 2019, Last update: December 9, 2019

Load required packages

library(bvartools)
library(dplyr)
library(ggplot2)
library(spam)
library(tidyr)

Data

US macro data until 2007Q4

# Specify URL
url <- "https://web.ics.purdue.edu/~jltobias/second_edition/Chapter20/code_for_exercise_1/US_macrodata.csv"
# Load data
us_macro_2007 <- read.csv(url, col.names = c("date", "Dp", "u", "r"))
# Omit NA values
us_macro_2007 <- na.omit(us_macro_2007)
# Transform into time-series object
us_macro_2007 <- ts(us_macro_2007[, -1], start = c(1959, 2), frequency = 4)

Exercise 20.1

Prepare the data

temp <- gen_var(us_macro_2007, p = 2, deterministic = "const")
y <- temp$Y
x <- temp$Z

k <- NROW(y)
tt <- ncol(y)
z <- kronecker(t(x), diag(1, k))
m <- ncol(z)

Specify the priors

mu_prior <- matrix(0, m) # Prior mean
v_i_prior <- diag(1, m) # Prior variance
diag(v_i_prior)[m - k + 1:k] <- 1 / 10 # Deterministics

df_prior <- k + 3
df_post <- tt + df_prior
scale_prior = diag(1, k)

Initial values

ols <- tcrossprod(y, x) %*% solve(tcrossprod(x))
epsilon <- y - ols %*% x
sigma <- tcrossprod(epsilon) / tt
sigma_i <- solve(sigma)

Gibbs sampler

iterations <- 21000 # Iterations
burnin <- 1000 # Burnin

# Storage for coefficients
store <- iterations - burnin
draws_b <- matrix(NA, m, store)
draws_sigma <- matrix(NA, k * k, store)

# Start Gibbs sampler
for (draw in 1:iterations) {
  # Draw parameters
  b <- post_normal_sur(y = y, z = z, sigma_i = sigma_i,
                       a_prior = mu_prior, v_i_prior = v_i_prior)
  
  # Error
  epsilon <- y - matrix(b, k) %*% x
  scale_post <- scale_prior + tcrossprod(epsilon)
  sigma_i <- matrix(rWishart(1, df_post, solve(scale_post))[,,1], k)
  sigma <- solve(sigma_i)
  
  # Save post-burnin draws
  if (draw > burnin) {
    draws_b[, draw - burnin] <- b
    draws_sigma[, draw - burnin] <- sigma
  }
}

Create a bvar object for IRF

draws_a <- draws_b[1:(m - k), ] # Draws of non-deterministic coefficients
draws_det <- draws_b[m - k + 1:k, ] # Draws of deterministic coefficients

ex_20_1 <- bvar(y = y, x = x, A = draws_a, C = draws_det, Sigma = draws_sigma)

ex_20_1 <- thin(ex_20_1) # Thin the output

Impuse response functions

Dp <- irf(ex_20_1, impulse = "r", response = "Dp", n.ahead = 20, type = "oir", keep_draws = TRUE)
u <- irf(ex_20_1, impulse = "r", response = "u", n.ahead = 20, type = "oir", keep_draws = TRUE)
r <- irf(ex_20_1, impulse = "r", response = "r", n.ahead = 20, type = "oir", keep_draws = TRUE)

Dp <- colMeans(Dp)
u <- colMeans(u)
r <- colMeans(r)

temp <- cbind(Dp, u, r) %>%
  as.data.frame() %>%
  mutate(time = 1:n() - 1) %>%
  pivot_longer(cols = c("Dp", "u", "r")) %>%
  mutate(name = factor(name, levels = c("Dp", "u", "r"),
                       labels = c("Inflation", "Unemployment", "Interest rate")))

ggplot(temp, aes(x = time, y = value)) +
  geom_line() +
  facet_wrap(~ name, scales = "free")

Exercise 20.6: Time varying parameter VAR

Data: US macro data until 2007Q4

# Specify URL
url <- "https://web.ics.purdue.edu/~jltobias/second_edition/Chapter20/code_for_exercise_5/US_macrodata.csv"
# Load data
us_macro_2007 <- read.csv(url, col.names = c("date", "Dp", "u", "r"))
# Omit NA values
us_macro_2007 <- na.omit(us_macro_2007)
# Transform into time-series object
us_macro_2007 <- ts(us_macro_2007[, -1], start = c(1959, 2), frequency = 4)

# Genenerate data matrices
temp <- gen_var(us_macro_2007, p = 2, deterministic = "const")
y <- temp$Y
x <- temp$Z

k <- nrow(y) # Number of endogenous variables
tt <- ncol(y) # Number of obs
m <- nrow(x) # Number of regressors
n <- m * k # Number of estimated coefficients

Estimation

Priors

# Initial conditions
a0_mu_prior <- matrix(0, n)
b0_v_i_prior <- diag(1, n)
diag(b0_v_i_prior)[n - k + 1:k] <- 1 / 10 # Prior for determinsitic terms

# Priors for measurement equation errors
sigma_df_prior <- k + 3
sigma_df_post <- tt + sigma_df_prior
sigma_scale_prior <- diag(1, k)

# Priors for state equation errors
q_shape_prior <- 3
q_shape_post <- q_shape_prior + tt / 2
q_rate_prior <- rep(.01^2, n)
q_rate_prior[n - k + 1:k] <- .1^2 # Determinsitic

Initial values

y_vec <- matrix(y)

z <- matrix(0, k * tt, n * tt)
for (i in 1:tt) {
  z[(i - 1) * k + 1:k, (i - 1) * n + 1:n] <- kronecker(t(x[, i]), diag(1, k))
}
z <- as.spam(z)

beta0 <- tcrossprod(y, x) %*% solve(tcrossprod(x))
sigma_temp <- tcrossprod(y - beta0 %*% x) / tt
sigma_i_temp <- solve(sigma_temp)
sigma_i <- kronecker(diag(1, tt), sigma_i_temp)
sigma_i <- as.spam(sigma_i)

beta0 <- matrix(beta0)

q_i <- kronecker(diag(1, tt), diag(100, n))
q_i <- as.spam(q_i)

h <- diag(1, m * k * tt)
diag(h[-(1:(m * k)), 1:(m * k * (tt - 1))]) <- -1
h <- as.spam(h)

Gibbs sampler

iterations <- 21000
burnin <- 1000

# Data containers for draws
draws_beta <- NULL
for (i in 1:tt) {
  draws_beta[[i]] <- matrix(NA, n, iterations - burnin)
}
draws_sigma <- matrix(NA, k * k, iterations - burnin)

# Start Gibbs sampler
for (draw in 1:iterations) {
  # Draw parameters
  HQiH <- crossprod(h, q_i) %*% h
  XSi <- crossprod(z, sigma_i)
  K_beta <- HQiH + XSi %*% z
  mu_beta <- matrix(solve(K_beta, HQiH %*% kronecker(matrix(1, tt), beta0) + XSi %*% y_vec))
  beta <- mu_beta + solve.spam(chol(K_beta), rnorm(n * tt))
  
  # Measurement eq. variance-covariance matrix
  u <- matrix(y_vec - z %*% beta, k)
  sigma_i_temp <- rWishart(n, sigma_df_post, solve(sigma_scale_prior + tcrossprod(u)))[,, 1]
  sigma_i <- as.spam(kronecker(diag(1, tt), sigma_i_temp))
  
  # State eq. error variances
  beta <- matrix(beta, n)
  v <- apply(beta - cbind(beta0, beta[, -tt]), 1, function(x) {sum(x^2)})
  q_post <- cbind(q_shape_post, q_rate_prior + v / 2)
  q_i_temp <- apply(q_post, 1, function(x) {rgamma(1, shape = x[1], rate = x[2])})
  diag(q_i) <- rep(q_i_temp, tt)
  
  # Initial condition
  K_beta0 <- b0_v_i_prior + q_i[1:n, 1:n]
  mu_beta0 <- solve(K_beta0, b0_v_i_prior %*% a0_mu_prior + q_i[1:n, 1:n] %*% beta[, 1])
  beta0 <- mu_beta0 + solve(chol(K_beta0), rnorm(n))
  
  # Save draws
  if (draw > burnin) {
    pos_draw <- draw - burnin
    
    for (i in 1:tt) {
      draws_beta[[i]][, pos_draw] <- as.matrix(beta[, i])
    }
    
    draws_sigma[, pos_draw] <- solve(sigma_i_temp)
  }
}

Impulse responses

Create bvar obejcts for further analysis

# For 1975Q1
tt_1975 <- which(time(us_macro_2007) == 1975) - 2
bvar_1975 <- bvar(y = y, x = x,
                   A = draws_beta[[tt_1975]][1:(n - k),],
                   C = draws_beta[[tt_1975]][n - k + 1:k, ],
                   Sigma = draws_sigma)

# For 2005Q1
tt_2005 <- which(time(us_macro_2007) == 2005) - 2
bvar_2005 <- bvar(y = y, x = x,
                   A = draws_beta[[tt_2005]][1:(n - k),],
                   C = draws_beta[[tt_2005]][n - k + 1:k, ],
                   Sigma = draws_sigma)

Obtain impulse responses for inflation

irf_1975 <- irf(bvar_1975, impulse = "r", response = "Dp",
                n.ahead = 20, type = "oir", keep_draws = TRUE)

irf_2005 <- irf(bvar_2005, impulse = "r", response = "Dp",
                n.ahead = 20, type = "oir", keep_draws = TRUE)

# Difference between the IRF
diff_irf_Dp <- irf_2005 - irf_1975
diff_irf_Dp <- t(apply(diff_irf_Dp, 2, quantile, probs = c(.05, .5, .95)))
diff_irf_Dp <- diff_irf_Dp %>%
  as.data.frame() %>%
  mutate(var = "Inflation",
         time = 1:n() - 1)

# Mean IRF
irf_1975 <- colMeans(irf_1975)
irf_2005 <- colMeans(irf_2005)

irf_Dp <- cbind(irf_1975, irf_2005) %>%
  as.data.frame() %>%
  mutate(var = "Inflation",
         time = 1:n() - 1)

Obtain impulse responses for unemployment

irf_1975 <- irf(bvar_1975, impulse = "r", response = "u",
                n.ahead = 20, type = "feir", keep_draws = TRUE)

irf_2005 <- irf(bvar_2005, impulse = "r", response = "u",
                n.ahead = 20, type = "feir", keep_draws = TRUE)

# Difference between the IRF
diff_irf_u <- irf_2005 - irf_1975
diff_irf_u <- t(apply(diff_irf_u, 2, quantile, probs = c(.05, .5, .95)))
diff_irf_u <- diff_irf_u %>%
  as.data.frame() %>%
  mutate(var = "Unemployment",
         time = 1:n() - 1)

# Mean IRF
irf_1975 <- colMeans(irf_1975)
irf_2005 <- colMeans(irf_2005)

irf_u <- cbind(irf_1975, irf_2005) %>%
  as.data.frame() %>%
  mutate(var = "Unemployment",
         time = 1:n() - 1)

Obtain impulse responses for interest rate

irf_1975 <- irf(bvar_1975, impulse = "r", response = "r",
                n.ahead = 20, type = "feir", keep_draws = TRUE)

irf_2005 <- irf(bvar_2005, impulse = "r", response = "r",
                n.ahead = 20, type = "feir", keep_draws = TRUE)

# Difference between the IRF
diff_irf_r <- irf_2005 - irf_1975
diff_irf_r <- t(apply(diff_irf_r, 2, quantile, probs = c(.05, .5, .95)))
diff_irf_r <- diff_irf_r %>%
  as.data.frame() %>%
  mutate(var = "Interest rate",
         time = 1:n() - 1)

# Mean IRF
irf_1975 <- colMeans(irf_1975)
irf_2005 <- colMeans(irf_2005)

irf_r <- cbind(irf_1975, irf_2005) %>%
  as.data.frame() %>%
  mutate(var = "Interest rate",
         time = 1:n() - 1)

Figure 20.2: Mean IRF

temp <- bind_rows(irf_Dp, irf_u, irf_r) %>%
  # Bring it into long form
  pivot_longer(cols = starts_with("irf_"),
               names_to = "Period") %>%
  # For better labelling in the plot
  mutate(var = factor(var, levels = c("Inflation", "Unemployment", "Interest rate")),
         Period = factor(Period, levels = c("irf_1975", "irf_2005"),
                         labels = c("1975", "2005")))

# Plot
ggplot(temp, aes(x = time, y = value, colour = Period)) +
  facet_wrap(~ var, scales = "free") +
  geom_line() +
  theme(legend.position = "bottom", legend.title = element_blank())

Figure 20.3: Plot differences between IRF

temp <- bind_rows(diff_irf_Dp, diff_irf_u, diff_irf_r)
names(temp) <- c("ymin", "y", "ymax", "var", "time")

temp <- temp %>%
  mutate(var = factor(var, levels = c("Inflation", "Unemployment", "Interest rate")))

# Plot
ggplot(temp, aes(x = time)) +
  geom_ribbon(aes(ymin = ymin, ymax = ymax, fill = "a"), alpha = .3, show.legend = FALSE) +
  geom_line(aes(y = y, colour = "a"), show.legend = FALSE) +
  facet_wrap(~ var, scales = "free")

Work in progress