Step one: Check data fits an ARCH model

The point of this test is to check whether the data is serially correlated and if variance is non-constant. We present the p-values from the test in the table below. Since they are all below 0.05, that tells us that the data exhibits the effects we’re looking for.

Step Two: Iterate model specifications to find good model

Iterate through different combinations of lags on the return and volatility data to find the model that produces the lowest value to two statistics commonly used to measure prediction error: Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). The tutorial goes into much greater detail than we do. We show the two best models below. We chose to use a GARCH(1,1) in our simulation for purposes of simplicity. We’re not sure we could explain an ARMA (6,6) + GARCH (1,1) model in a short post!

Step three: Simulate

We graphed the simulations above. We show a sample of the results in the table below.

Using the simulations, one can create a cumulative probability distribution of exceeding the former peak by a certain date (another neat chart from the tutorial). We find that even after five years there’s still around a 20% chance the S&P won’t have surpassed it’s former peak. We’re not sure that’s the right probability, but it’s reasonable given past history.

Of course, if we had used daily instead of monthly returns and volatility, the probabilities would have gone up. This is an artifact of the time series—namely, the annualized volatility of the monthly returns is less than that of daily ones. While it is beyond the scope of this post to explain why that is the case, the simple explanation is that there is less noise in the monthly data. Nonetheless, the point is that with a lower volatility, expected moves are lower, thus probability of exceeding some specific level will be lower too.

We have tried to reproduce the results from the tutorial faithfully, albeit with different data. Any errors are ours.

Here’s the code.

## We imported many of the functions found on the github page for the tutorial. 
## Please reference that page for the full code. These functions are preceded 
## by source('garch_fcts.R')

## NB: Often one will see the dplyr package called directly to get the filter()
## function. This is because one of other packages masks the dplyr function.


### Load packages
suppressPackageStartupMessages({
  library(tidyquant)
  library(tidyverse)
  library(FinTS)
  library(fGarch)
})

options("getSymbols.warning4.0"=FALSE)

### Load data
sp <- getSymbols("^GSPC", src = "yahoo", from = "1950-01-01", auto.assign = FALSE) %>% 
  Ad() %>% 
  `colnames<-`("price")

sp <- data.frame(date = index(sp), price = as.numeric(sp$price)) %>% 
  mutate(log_ret = log(price/dplyr::lag(price)),
         disc_ret = price/dplyr::lag(price)-1)

sp_mon <- sp %>% 
  mutate(year_mon = as.yearmon(date)) %>% 
  group_by(year_mon) %>% 
  dplyr::filter(date == max(date)) %>% 
  ungroup() %>% 
  mutate(log_ret = log(price/dplyr::lag(price)),
         disc_ret = price/dplyr::lag(price)-1) %>% 
  select(-year_mon)

## Plot annualized volatility
sp %>% 
  mutate(vol = rollapply(log_ret, width = 22, sd, align = "right", fill = NA)*sqrt(252)) %>% 
  ggplot(aes(date, vol*100)) +
  geom_line(color = "blue") +
  labs(x = "",
       y = "Volatility (%)",
       title = "Annualized rolling 22-day volatility of the S&P 500")

## Plot log chart of S&P
sp %>% 
  ggplot(aes(date, price)) +
  geom_line(color = "blue") +
  scale_y_log10() +
  labs(x = "",
       y = "Index",
       title ="S&P 500 long-run chart with y-axis on log-scale")

## Table of volatility spike frequency
sp %>% 
  mutate(vol = rollapply(log_ret, width = 22, sd, align = "right", fill = NA)*sqrt(252)) %>% 
  mutate(Regime = ifelse(date < "2000-01-01", "Pre-2000", "Post-2000")) %>% 
  group_by(Regime) %>% 
  summarise("Frequency (%)" = round(sum(vol > 0.25, na.rm = TRUE)/n()*100,1)) %>% 
  knitr::kable(caption = "Frequency of volatility spikes above 25%")

## Find largest 20 moves
largest <- sp %>% 
  dplyr::filter(date >= "1980-01-01") %>% 
  top_n(abs(log_ret), n = 20)

## Plot volatility clustering
sp %>% 
  dplyr::filter(date >= "1980-01-01") %>% 
  ggplot(aes(date, log_ret)) + 
  geom_line(color = "blue") + 
  labs(title = "S&P 500 daily log returns",
       subtitle = "Red dots are the 20 largest absolute daily price changes",
       x = '',
       y = 'Returns',
       caption = 'Based on code written by MS Perlin') + 
  geom_point(data = largest, 
             aes(date, log_ret), 
             size = 3, color = 'red'  ) +
  geom_text(aes(x = as.Date("2015-01-01"), 
                y = -0.22, 
                label = "Large moves cluster"),
            color = "red") +
  geom_segment(aes(x = as.Date("2015-01-01"), xend = as.Date("2009-01-01"),
                   y = -0.2, yend = -0.12),
               arrow = arrow(length = unit(2, "mm")),
               color = "red") + 
  geom_segment(aes(x = as.Date("2015-06-01"), xend = as.Date("2020-02-01"),
                   y = -0.2, yend = -0.12),
               arrow = arrow(length = unit(2, "mm")),
               color = "red") + 
  scale_y_continuous(labels = scales::percent) +
  theme(plot.caption = element_text(hjust = 0))

## Run simulation
set.seed(20200320) 
n_sim <- 1000 
n_periods <- 72

# source functions
source('garch_fcts.R') # for sim_func. We altered this to include different simulation
# perods; e.g., days, months, quarters, etc.

df_sim <- sim_func(n_sim = n_sim, 
                 n_t = n_days_ahead, my_garch, "months", df_prices = sp_mon)

sp_temp <- sp_mon %>% 
  dplyr::filter(date >= "2010-01-03") #dplyr::filter masked by some other package!

ggplot() + 
  geom_line(data = sp_temp, 
            aes(x = date, y = price), color = 'blue', size = 0.75)  + 
  geom_line(data = df_sim, aes(x = ref_date, y = sim_price, group = i_sim),
            color = 'lightgrey', 
            size = 0.35) + 
  geom_hline(yintercept = max(sp_temp$price), color = "red", size = 0.35) + 
  labs(title = 'GARCH model projections for the S&P 500',
       subtitle = 'Based on 1000 simulations',
       x = '',
       y = 'Index',
       captioin = "Based on code written by MS Perlin") + 
  ylim(c(min(sp_temp$price), max(df_sim$sim_price))) +
  scale_y_continuous(labels = scales::comma) +
  theme(plot.caption = element_text(hjust = 0))


# ARCH test function from MS Perlin
source('garch_fcts.r')

# Run test up to six-month lag
max_lag <- 6

arch_test <- do_arch_test(x = sp_mon$log_ret, max_lag = max_lag)

arch_test %>% 
  mutate(statistic = round(statistic,1), 
         pvalue = format(round(pvalue,2),nsmall=1 )) %>% 
  rename("Lag" = lag,
         "Statistic" = statistic,
         "P-value" = pvalue) %>% 
  knitr::kable(caption = "S&P 500 Arch test", 
               align = c("l", "r", "r"))

tab_out <- out$tab_out

sp_long <- pivot_longer(data = tab_out %>%
                                 select(model_name, AIC, BIC),  cols = c('AIC', 'BIC'))

sp_best_models <- sp_long %>%
  group_by(name) %>%
  summarise(model_name = model_name[which.min(value)],
            value = value[which.min(value)])

models_names <- unique(sp_long$model_name)
best_model <- c(tab_out$model_name[which.min(tab_out$AIC)],
                 tab_out$model_name[which.min(tab_out$BIC)])

best_model %>% 
  mutate(value = round(value,2)) %>% 
  rename("Criterion" = name,
         "Model" = model_name,
         "Statistic" = value) %>% 
  knitr::kable(caption = "GARCH models with lowest information criterion")


df_sim %>% 
  group_by(i_sim) %>% 
  summarise(tot_ret = round((last(sim_price)/last(sp_mon$price)-1)*100,1)) %>% 
  arrange(tot_ret) %>% 
  slice(c(1:5, (n()-5):n())) %>% 
  rename("Simulation (#)" = i_sim,
         "Return (%)" = tot_ret) %>% 
  knitr::kable(caption = "Top and bottom five cumulative returns")

tab_prob <- df_sim %>%
  group_by(ref_date) %>%
  summarise(prob = mean(sim_price > max(sp$price)))


df_date <- tibble(idx = c(first(which(tab_prob$prob > 0.25)),
                          first(which(tab_prob$prob > 0.5)),
                          first(which(tab_prob$prob > 0.75))),
                  ref_date = tab_prob$ref_date[idx],
                  prob = tab_prob$prob[idx],
                  my_text = format(ref_date, '%b %Y'))

tab_prob %>% 
  ggplot(aes(x = ref_date, y = prob) ) + 
  geom_line(size = 1.25, color = "blue") + 
  labs(title = paste0('Probabilities of S&P 500 reaching former peak of ', 
                      round(max(sp$price),0)),
       subtitle = paste0('Calculations based on GARCH model simulations'),
       x = '',
       y = 'Probability') + 
  scale_y_continuous(labels = scales::percent) + 
  geom_point(data = df_date,
             aes(x = ref_date, y = prob), size = 3.5, color = 'red') + 
  geom_text(data = df_date, aes(x = ref_date, y = prob, 
                                label = my_text), 
            nudge_x = -200,
            # nudge_y = -0.05,
            color ='red', check_overlap = TRUE)