chapter_3_lesson_4_handout.qmd

---
title: "Holt-Winters Method (Additive Models) - Part 2"
subtitle: "Chapter 3: Lesson 4"
format: html
editor: source
sidebar: false
---

```{r}
#| include: false
source("common_functions.R")
```

```{=html}
<script type="text/javascript">
 function showhide(id) {
    var e = document.getElementById(id);
    e.style.display = (e.style.display == 'block') ? 'none' : 'block';
 }
 
 function openTab(evt, tabName) {
    var i, tabcontent, tablinks;
    tabcontent = document.getElementsByClassName("tabcontent");
    for (i = 0; i < tabcontent.length; i++) {
        tabcontent[i].style.display = "none";
    }
    tablinks = document.getElementsByClassName("tablinks");
    for (i = 0; i < tablinks.length; i++) {
        tablinks[i].className = tablinks[i].className.replace(" active", "");
    }
    document.getElementById(tabName).style.display = "block";
    evt.currentTarget.className += " active";
 }    
</script>
```


```{r}
#| echo: false
#| warning: false

nat_gas <- rio::import("https://byuistats.github.io/timeseries/data/natural_gas_res.csv") |>
  mutate(date = my(month)) |>
  filter(date >= my("Jan 2017"))|>
  mutate(quarter = yearquarter(date)) |>
  group_by(quarter) |> 
  summarize(
    gas_use_mmcf = sum(residential_nat_gas_consumption),
    n = n()
  ) |>
  filter(n == 3) |>  # Eliminate partial quarter(s)
  dplyr::select(-n) |>
  mutate(gas_billion_cf = round(gas_use_mmcf / 10^3))
nat_gas_ts <- nat_gas %>%
  as_tsibble(index = quarter)
```

#### Figure 1: Quarterly U.S. natural gas consumption (Bcf)

```{r}
#| echo: false
#| warning: false

# This is a variation on the function hw_additive_slope_additive_seasonal(), but it rounds the data to the nearest one unit at each step
hw_additive_slope_additive_seasonal_gas <- function(df, date_var, value_var, p = 12, predict_periods = 18, alpha = 0.2, beta = 0.2, gamma = 0.2, s_initial = rep(0,p)) {
  
  # Get expanded data frame
  df <- df |> expand_holt_winters_df(date_var, value_var, p, predict_periods)
  
  # Fill in prior belief about s_t
  for (t in 1:p) {
    df$s_t[t] <- s_initial[t]
  }
  
  # Fill in first row of values
  offset <- p # number of header rows to skip
  df$a_t[1 + offset] <- round( df$x_t[1 + offset] )
  df$b_t[1 + offset] <- round( (1 / p) * mean(df$x_t[(p + 1 + offset):(2 * p + offset)] - df$x_t[(1 + offset):(p + offset)]) )
  df$s_t[1 + offset] <- round( (1 - gamma) * df$s_t[1] )
  df$xhat_t[1 + offset] <- round( df$x_t[1 + offset] )

  # Fill in remaining rows of body of df with values
  for (t in (2 + offset):(nrow(df) - predict_periods) ) {
    df$a_t[t] = round( alpha * (df$x_t[t] - df$s_t[t-p]) + (1 - alpha) * (df$a_t[t-1] + df$b_t[t-1]) )
    df$b_t[t] = round( beta * (df$a_t[t] - df$a_t[t-1]) + (1 - beta) * df$b_t[t-1] )
    df$s_t[t] = round( gamma * (df$x_t[t] - df$a_t[t]) + (1 - gamma) * df$s_t[t-p] )
    df$xhat_t[t] = round( (df$a_t[t] + df$s_t[t]) )
  }
  
  df <- df |>
    mutate(k = ifelse(row_number() >= nrow(df) - predict_periods, row_number() - (nrow(df) - predict_periods), NA))
  
  # Fill in forecasted values
  offset <- nrow(df) - predict_periods
  for (t in offset:nrow(df)) {
    df$s_t[t] = round( df$s_t[t - p] )
    df$xhat_t[t] = round( df$a_t[offset] + df$k[t] * df$b_t[offset] + df$s_t[t - p] )
  }
  
  # Delete temporary variable k
  df <- df |> select(-k)
  
  return(df)
}

nat_gas_ts <- nat_gas |>
  hw_additive_slope_additive_seasonal_gas("quarter", "gas_billion_cf", p = 4, predict_periods = 9, s_initial = c(1000, -1000, -1000, 1000)) |>
  as_tsibble(index = date) %>%  
  mutate(t = (1 - 4):(nrow(.) - 4)) |> # Create column t
  dplyr::select(date, t, x_t, a_t, b_t, s_t, xhat_t)

nat_gas_ts |>
  filter(t > 0) |>
  filter(!is.na(x_t)) |> ###### Comment out to get the solution
  ggplot(aes(x = date)) +
    geom_line(aes(y = x_t), color = "black", size = 1) +
    ################### Uncomment these lines to show the solution ############################
    # geom_line(aes(y = a_t + s_t, color = "Combined", alpha=0.5), size = 1) +
    # geom_line(aes(y = xhat_t, color = "Combined", alpha=0.5), linetype = "dashed", size = 1) +
    coord_cartesian(ylim = c(0,2500)) +
    labs(
      x = "Date",
      y = "Natural Gas Use (Billions of Cubic Feet)",
      title = "U.S. Residential Natural Gas Consumption, by Quarter",
      color = "Components"
    ) +
    theme_minimal() +
    theme(legend.position = "none") +
    theme(
      plot.title = element_text(hjust = 0.5)
    )
```

$\ $
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>

#### Table 1: Holt-Winters filter for the U.S. quarterly natural gas consumption (in Bcf)
```{r}
#| echo: false
#| warning: false

nat_gas_ts |>
  # Hide values from students
  mutate(
    xhat_t = 
      case_when(
        t %in% c(1:4, 26:36) ~ NA,
        t %in% c(7:25) ~ a_t + s_t,
        TRUE ~ xhat_t
      ),
    a_t = ifelse(t %in% c(1:4), NA, a_t),
    b_t = ifelse(t %in% c(1:4), NA, b_t),
    s_t = ifelse(t %in% c(1:4), NA, s_t)
    # xhat_t = ifelse(t %in% c(24:27), NA, xhat_t)
  ) |>
  replace_na_with_char() |>
  # Add emdashes to make it clear where students should not put values
  mutate(
    x_t = ifelse(t %in% c(-3:0, 28:36), emdash, x_t),
    a_t = ifelse(t %in% c(-3:0, 28:36), emdash, a_t),
    b_t = ifelse(t %in% c(-3:0, 28:36), emdash, b_t),
    xhat_t = ifelse(t %in% c(-3:0), emdash, xhat_t)
  ) |>
  rename(
    "$$Quarter$$" = date,
    "$$t$$" = t,
    "$$x_t$$" = x_t,
    "$$a_t$$" = a_t,
    "$$b_t$$" = b_t,
    "$$s_t$$" = s_t,
    "$$\\hat x_t$$" = xhat_t,
  ) |>
  display_table("0.75in")
```