_chapter_6_lesson_2.qmd

---
title: "Autoregressive Moving Average (ARMA) Models"
subtitle: "Chapter 6: Lesson 2"
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>
```


## Learning Outcomes

{{< include outcomes/_chapter_6_lesson_2_outcomes.qmd >}}




## Preparation

-   Read Sections 6.5.1 and 6.6-6.7



## Learning Journal Exchange (10 min)

-   Review another student's journal

-   What would you add to your learning journal after reading another student's?

-   What would you recommend the other student add to their learning journal?

-   Sign the Learning Journal review sheet for your peer




## Class Activity: Introduction to Autoregressive Moving Average (ARMA) Models (xx min)

### Autoregressive (AR) Models

In [Chapter 4, Lesson 3](https://byuistats.github.io/timeseries/chapter_4_lesson_3.html#ARdefinition), we learned the definition of an AR model:

::: {.callout-note icon=false title="Definition of an Autoregressive (AR) Model"}

The time series $\{x_t\}$ is an **autoregressive process of order $p$**, denoted as $AR(p)$, if
$$
  x_t = \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \alpha_3 x_{t-3} + \cdots + \alpha_{p-1} x_{t-(p-1)} + \alpha_p x_{t-p} + w_t ~~~~~~~~~~~~~~~~~~~~~~~ (4.15)
$$

where $\{w_t\}$ is white noise and the $\alpha_i$ are the model parameters with $\alpha_p \ne 0$.

:::

The $AR(p)$ model can be expressed as:
$$
  \underbrace{\left( 1 - \alpha_1 \mathbf{B} - \alpha_2 \mathbf{B}^2 - \cdots - \alpha_p \mathbf{B}^p \right)}_{\theta_p \left( \mathbf{B} \right)} x_t = w_t
$$

### Moving Average (MA) Models

The definition of an $MA(q)$ model is:

::: {.callout-note icon=false title="Definition of a Moving Average (MA) Model"}

We say that a time series $\{x_t\}$ is a **moving average process of order $q$**, denoted as $MA(q)$, if each term in the time series is a linear combination of the current white noise term and the $q$ most recent past white noise terms. 

It is given as:
$$
  x_t = w_t + \beta_1 w_{t-1} + \beta_2 w_{t-2} + \beta_3 w_{t-3} + \cdots + \beta_{q-1} w_{t-(q-1)} + \beta_q w_{t-q}
$$

where $\{w_t\}$ is white noise with zero mean and variance $\sigma_w^2$, and the $\beta_i$ are the model parameters with $\beta_q \ne 0$.

:::

Written in terms of the backward shift operator, we have 
$$
  x_t = \underbrace{\left( 1 + \beta_1 \mathbf{B} + \beta_2 \mathbf{B}^2  + \beta_3 \mathbf{B}^3  + \cdots + \beta_{q-1} \mathbf{B}^{q-1}  + \beta_q \mathbf{B}^{q} \right)}_{\phi_q(\mathbf{B})} w_t
$$

Putting the $AR$ and $MA$ models together, we get the $ARMA$ model. 

::: {.callout-note icon=false title="Definition of an Autogregressive Moving Average (ARMA) Model"}

A time series $\{ x_t \}$ follows an **autoregressive moving average (ARMA) model** of order $(p, q)$, which we write as $ARMA(p,q)$, if it can be written as:

$$
  x_t = 
    \underbrace{
      \alpha_1 x_{t-1} + \alpha_2 x_{t-2} 
      + \alpha_3 x_{t-3} 
      + \cdots 
      % + \alpha_{p-1} x_{t-(p-1)} + 
      \alpha_p x_{t-p}
    }_{
      AR(p) ~ \text{model}
    } 
    + \underbrace{
      w_t + \beta_1 w_{t-1} + \beta_2 w_{t-2} 
      + \beta_3 w_{t-3} 
      + \cdots 
      % + \beta_{q-1} w_{t-(q-1)} 
      + \beta_q w_{t-q}
    }_{
      MA(q) ~ \text{model}
    } 
$$
where $\{w_t\}$ is a white noise process.

:::

We can write this as:
$$
  \theta_p \left( \mathbf{B} \right) x_t
  =
  \phi_q \left( \mathbf{B} \right) w_t
$$

::: {.callout-caution icon=false title="Facts about ARMA Processes"}

The following facts are true for $ARMA(p,q)$ processes:

-   The ARMA process is stationary if all the roots of $\theta_p \left( \mathbf{B} \right)$ are greater than 1 in absolute value.
-   The ARMA process is invertible if all the roots of $\phi_q \left( \mathbf{B} \right)$ are greater than 1 in absolute value.
-   The special case $ARMA(p,0)$ is the $AR(p)$ model.
-   The special case $ARMA(0,p)$ is the $MA(q)$ model.
-   An $ARMA$ model will usually require fewer parameters than a single $MA$ or $AR$ model. This is called *parameter parsimony*.
-   If $\theta$ and $\phi$ have a common factor, a stationary model can be simplified. This is called *parameter redundancy*.
As an example, the model 
$$
  \left( 1 - \frac{1}{2} \mathbf{B} \right)\left( 1 - \frac{1}{3} \mathbf{B} \right) x_t 
  = 
  \left( 1-\frac{1}{2} \mathbf{B} \right)\left( 1 - \frac{1}{4} \mathbf{B} \right) w_t
$$ 
is the same as the model
$$
  \left( 1 - \frac{1}{3} \mathbf{B} \right) x_t 
  = 
  \left( 1 - \frac{1}{4} \mathbf{B} \right) w_t
$$

:::


## Class Activity: Model for the Residuals from the Rexburg Weather Model

### Review

We now review the model we built in [Chapter 5 Lesson 3](https://byuistats.github.io/timeseries/chapter_5_lesson_3.html#weather) for the monthly average of the daily high temperature in Rexburg, Idaho.

```{r}
#| label: weather1
#| code-fold: true
#| code-summary: "Show the code"
#| warning: false

weather_df <- rio::import("https://byuistats.github.io/timeseries/data/rexburg_weather_monthly.csv") |>
  mutate(dates = my(date_text)) |>
  filter(dates >= my("1/2008") & dates <= my("12/2023")) |>
  rename(x = avg_daily_high_temp) |>
  mutate(TIME = 1:n()) |>
  mutate(
    cos1 = cos(2 * pi * 1 * TIME/12),
    cos2 = cos(2 * pi * 2 * TIME/12),
    cos3 = cos(2 * pi * 3 * TIME/12),
    cos4 = cos(2 * pi * 4 * TIME/12),
    cos5 = cos(2 * pi * 5 * TIME/12),
    cos6 = cos(2 * pi * 6 * TIME/12),
    sin1 = sin(2 * pi * 1 * TIME/12),
    sin2 = sin(2 * pi * 2 * TIME/12),
    sin3 = sin(2 * pi * 3 * TIME/12),
    sin4 = sin(2 * pi * 4 * TIME/12),
    sin5 = sin(2 * pi * 5 * TIME/12),
    sin6 = sin(2 * pi * 6 * TIME/12)) |>
  mutate(zTIME = (TIME - mean(TIME)) / sd(TIME)) |>
  as_tsibble(index = TIME)

weather_df |>
  as_tsibble(index = dates) |>
  autoplot(.vars = x) +
  geom_smooth(method = "lm", se = FALSE, color = "#F0E442") +
    labs(
      x = "Month",
      y = "Mean Daily High Temperature (Fahrenheit)",
      title = "Time Plot of Mean Daily Rexburg High Temperature by Month",
      subtitle = paste0("(", format(weather_df$dates %>% head(1), "%b %Y"), endash, format(weather_df$dates %>% tail(1), "%b %Y"), ")")
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0.5)
    )
```





<!-- Now, we compare all the models side-by-side. -->

<!-- ```{r} -->
<!-- #| label: weather17 -->
<!-- #| code-fold: true -->
<!-- #| code-summary: "Show the code" -->
<!-- #| output: false -->

<!-- model_combined <- weather_df |> -->
<!--   model( -->
<!--     full_cubic = TSLM(x ~ TIME + I(TIME^2) + I(TIME^3) + -->
<!--       sin1 + cos1 + sin2 + cos2 + sin3 + cos3  -->
<!--       + sin4 + cos4 + sin5 + cos5 + cos6), -->
<!--     full_quadratic = TSLM(x ~ TIME + I(TIME^2) + -->
<!--       sin1 + cos1 + sin2 + cos2 + sin3 + cos3  -->
<!--       + sin4 + cos4 + sin5 + cos5 + cos6), -->
<!--     full_linear = TSLM(x ~ TIME  + -->
<!--       sin1 + cos1 + sin2 + cos2 + sin3 + cos3  -->
<!--       + sin4 + cos4 + sin5 + cos5 + cos6 ), -->
<!--     reduced_quadratic_1  = TSLM(x ~ TIME + I(TIME^2) + sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + cos6), -->
<!--     reduced_quadratic_2  = TSLM(x ~ TIME + I(TIME^2) + sin1 + cos1 + sin2 + cos2 + sin3 + cos6), -->
<!--     reduced_quadratic_3  = TSLM(x ~ TIME + I(TIME^2) + sin1 + cos1 + sin2 + cos2 + sin3), -->
<!--     reduced_quadratic_4  = TSLM(x ~ TIME + I(TIME^2) + sin1 + cos1 + sin2 + cos2 + sin3 + cos3), -->
<!--     reduced_quadratic_5  = TSLM(x ~ TIME + I(TIME^2) + sin1 + cos1 + sin2 + cos2), -->
<!--     reduced_quadratic_6  = TSLM(x ~ TIME + I(TIME^2) + sin1 + cos1), -->
<!--     reduced_linear_1  = TSLM(x ~ TIME + sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + cos6), -->
<!--     reduced_linear_2  = TSLM(x ~ TIME + sin1 + cos1 + sin2 + cos2 + sin3 + cos6), -->
<!--     reduced_linear_3  = TSLM(x ~ TIME + sin1 + cos1 + sin2 + cos2 + sin3), -->
<!--     reduced_linear_4  = TSLM(x ~ TIME + sin1 + cos1 + sin2 + cos2 + sin3 + cos3), -->
<!--     reduced_linear_5  = TSLM(x ~ TIME + sin1 + cos1 + sin2 + cos2), -->
<!--     reduced_linear_6  = TSLM(x ~ TIME + sin1 + cos1) -->
<!--   ) -->

<!-- glance(model_combined) |> -->
<!--   select(.model, AIC, AICc, BIC) -->
<!-- ``` -->


<!-- ```{r} -->
<!-- #| label: weather18 -->
<!-- #| echo: false -->

<!-- combined_models <- glance(model_combined) |>  -->
<!--   select(.model, AIC, AICc, BIC) -->

<!-- minimum <- combined_models |> -->
<!--   summarize( -->
<!--     AIC = which(min(AIC)==AIC), -->
<!--     AICc = which(min(AICc)==AICc), -->
<!--     BIC = which(min(BIC)==BIC) -->
<!--   ) -->
<!-- combined_models |> -->
<!--   rename(Model = ".model") |> -->
<!--   round_df(1) |> -->
<!--   format_cells(rows = minimum$AIC, cols = 2, "bold") |> -->
<!--   format_cells(rows = minimum$AICc, cols = 3, "bold") |> -->
<!--   format_cells(rows = minimum$BIC, cols = 4, "bold") |> -->
<!--   display_table() -->
<!-- ``` -->


We chose the "Reduced Linear 5" model. For convenience, we reprint the coefficients here.

```{r}
#| label: weather22
#| code-fold: true
#| code-summary: "Show the code"

reduced5_linear_lm <- weather_df |>
  model(reduced_linear_5  = TSLM(x ~ zTIME + sin1 + cos1 + sin2 + cos2))

reduced5_linear_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)

r5lin_coef_unrounded <- reduced5_linear_lm |>
  tidy() |>
  select(term, estimate, std.error)

r5lin_coef <- r5lin_coef_unrounded |>
  round_df(3)

stats_unrounded <- weather_df |>
  as_tibble() |>
  dplyr::select(TIME) |>
  summarize(mean = mean(TIME), sd = sd(TIME))

stats <- stats_unrounded |>
  round_df(3)
```

The fitted model is:

\begin{align*}
  x_t 
      &= \hat \beta_0 + \hat \beta_1 \left( \frac{t - \bar t}{s_t} \right)  \\
      & ~~~~~~~~~~ + \hat \beta_2 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \hat \beta_3 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~ + \hat \beta_4 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \hat \beta_5 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) 
      \\
      &= `r r5lin_coef$estimate[1]` 
            + `r r5lin_coef$estimate[2]` \left( \frac{t - `r stats$mean`}{`r stats$sd`} \right) \\ 
      & ~~~~~~~~~~~~~~~~~ + (`r r5lin_coef$estimate[3]`) \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + (`r r5lin_coef$estimate[4]`) \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~ + `r r5lin_coef$estimate[5]` \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + (`r r5lin_coef$estimate[6]`) \cos \left( \frac{2\pi \cdot 2 t}{12} \right) 
      \\
\end{align*}

### Fiiting an ARMA(p,q) Model

First, we create an acf and pacf plot of the residuals from the model above.

```{r}
#| code-fold: true
#| code-summary: "Show the code"
#| label: fig-acfResidWeather
#| fig-cap: "ACF Plot of the Residuals from the Rexburg Weather Model"

reduced5_linear_lm |>
  residuals() |>
  ACF() |>
  autoplot(var = .resid)
```


```{r}
#| code-fold: true
#| code-summary: "Show the code"
#| label: fig-pacfResidWeather
#| fig-cap: "PACF Plot of the Residuals from the Rexburg Weather Model"

reduced5_linear_lm |> 
  residuals() |>
  PACF() |>
  autoplot(var = .resid)
```


<!-- Check Your Understanding -->

::: {.callout-tip icon=false title="Check Your Understanding"}

-   What ARMA model does the combined information in the acf and pacf plots suggest would be appropriate for the residuals of the Rexburg weather model?

:::

Here are summaries of some ARMA models that could be constructed.

```{r}
#| output: false

model_resid <- reduced5_linear_lm |>
  residuals() |>
  select(-.model) |>
  model(
    auto = ARIMA(.resid ~ 1 + pdq(0:2,0,0:2) + PDQ(0, 0, 0)),
    a000 = ARIMA(.resid ~ 1 + pdq(0,0,0) + PDQ(0, 0, 0)),
    a001 = ARIMA(.resid ~ 1 + pdq(0,0,1) + PDQ(0, 0, 0)),
    a002 = ARIMA(.resid ~ 1 + pdq(0,0,2) + PDQ(0, 0, 0)),
    a100 = ARIMA(.resid ~ 1 + pdq(1,0,0) + PDQ(0, 0, 0)),
    a101 = ARIMA(.resid ~ 1 + pdq(1,0,1) + PDQ(0, 0, 0)),
    a102 = ARIMA(.resid ~ 1 + pdq(1,0,2) + PDQ(0, 0, 0)),
    a200 = ARIMA(.resid ~ 1 + pdq(2,0,0) + PDQ(0, 0, 0)),
    a201 = ARIMA(.resid ~ 1 + pdq(2,0,1) + PDQ(0, 0, 0)),
    a202 = ARIMA(.resid ~ 1 + pdq(2,0,2) + PDQ(0, 0, 0))) 

model_resid |>
  glance() 
```


```{r}
#| echo: false

# model_best |>
#   residuals() |>
#   ACF() |>
#   autoplot()


# new_data <- tibble(
#   date = seq(
#     max(cbe_ts$date) + months(1),
#       max(cbe_ts$date) + months(36),
#       by = "1 months"),
#   month = tsibble::yearmonth(date),  
#   time = seq(nrow(cbe_ts), length = 36),
#   imth = rep(1:12, 3)) |>
#   as_tsibble(index = month)
# 
# plot_dat <- new_data |>
#   mutate(
#     residuals = forecast(
#       model_best,
#       h = "3 years") |>
#       pull(.mean),
#     expected_log = forecast(
#       elec_lm,
#       new_data = new_data) |>
#       pull(.mean),
#     expected = exp(expected_log + residuals)  
#   )
# 
# ggplot() +
#   geom_line(data = cbe_ts, aes(x = date, y = elec)) +
#   geom_line(data = plot_dat, aes(x = date, y = expected), linetype = 2)
```




```{r}
#| label: tbl-ModelComparison
#| tbl-cap: "Comparison of the AIC, AICc, and BIC values for the models fitted to the logarithm of the simulated time series."
#| echo: false

combined_models <- glance(model_resid) |> 
  select(.model, AIC, AICc, BIC)

minimum <- combined_models |>
  reframe(
    AIC = which(min(AIC)==AIC),
    AICc = which(min(AICc)==AICc),
    BIC = which(min(BIC)==BIC)
  )

combined_models |>
  rename(Model = ".model") |>
  round_df(1) |>
  format_cells(rows = minimum$AIC, cols = 2, "bold") |>
  format_cells(rows = minimum$AICc, cols = 3, "bold") |>
  format_cells(rows = minimum$BIC, cols = 4, "bold") |>
  display_table()
```

```{r}
model_resid |>
  tidy() |>
  filter(.model == "a101") 

model_resid |>
  select(a101) |>
  residuals() |>
  ACF() |>
  autoplot()

model_resid |>
  select(a101) |>
  residuals() |>
  PACF() |>
  autoplot()

```





<!-- ## SMALL GROUP ACTIVITY -->

<!-- ```{r} -->
<!-- #| label: fig-RetailSalesGeneralMerch -->
<!-- #| fig-cap: "Time plot of the total monthly retail sales for all other general merchandise stores (45299)" -->
<!-- #| code-fold: true -->
<!-- #| code-summary: "Show the code" -->

<!-- # Read in retail sales data for "all other general merchandise stores" -->
<!-- retail_ts <- rio::import("https://byuistats.github.io/timeseries/data/retail_by_business_type.parquet") |> -->
<!--   filter(naics == 442) |> -->
<!--   # filter(as_date(month) >= my("Jan 1998")) |> -->
<!--   mutate(month = as_date(month)) |> -->
<!--   mutate(t = 1:n()) |>  -->
<!--   mutate(std_t = (t - mean(t)) / sd(t)) |> -->
<!--   mutate( -->
<!--     cos1 = cos(2 * pi * 1 * t / 12), -->
<!--     cos2 = cos(2 * pi * 2 * t / 12), -->
<!--     cos3 = cos(2 * pi * 3 * t / 12), -->
<!--     cos4 = cos(2 * pi * 4 * t / 12), -->
<!--     cos5 = cos(2 * pi * 5 * t / 12), -->
<!--     cos6 = cos(2 * pi * 6 * t / 12), -->
<!--     sin1 = sin(2 * pi * 1 * t / 12), -->
<!--     sin2 = sin(2 * pi * 2 * t / 12), -->
<!--     sin3 = sin(2 * pi * 3 * t / 12), -->
<!--     sin4 = sin(2 * pi * 4 * t / 12), -->
<!--     sin5 = sin(2 * pi * 5 * t / 12) -->
<!--   ) |> -->
<!--   as_tsibble(index = month) -->

<!-- retail_ts |> -->
<!--   autoplot(.vars = sales_millions) + -->
<!--   stat_smooth(method = "lm",  -->
<!--               formula = y ~ x,  -->
<!--               geom = "smooth", -->
<!--               se = FALSE, -->
<!--               color = "#E69F00", -->
<!--               linetype = "dotted") + -->
<!--     labs( -->
<!--       x = "Month", -->
<!--       y = "Sales (Millions of U.S. Dollars)", -->
<!--       title = paste0(retail_ts$business[1], " (", retail_ts$naics[1], ")") -->
<!--     ) + -->
<!--     theme_minimal() + -->
<!--     theme(plot.title = element_text(hjust = 0.5)) -->

<!-- ``` -->

<!-- <!-- @fig-retailSideBySidePlot shows the "All other general merchandise" retail sales data. --> -->
<!-- ```{r} -->
<!-- #| include: false -->
<!-- #| label: fig-retailSideBySidePlot -->
<!-- #| fig-cap: "Time plot of the time series (left) and the natural logarithm of the time series (right)" -->
<!-- #| code-fold: true -->
<!-- #| code-summary: "Show the code" -->
<!-- #| results: asis -->
<!-- #| fig-height: 3.5 -->

<!-- retail_plot_raw <- retail_ts |> -->
<!--     autoplot(.vars = sales_millions) + -->
<!--     labs( -->
<!--       x = "Month", -->
<!--       y = "sales_millions", -->
<!--       title = "Other General Merchandise Sales" -->
<!--     ) + -->
<!--     theme_minimal() + -->
<!--     theme( -->
<!--       plot.title = element_text(hjust = 0.5) -->
<!--     ) -->

<!-- retail_plot_log <- retail_ts |> -->
<!--     autoplot(.vars = log(sales_millions)) + -->
<!--     labs( -->
<!--       x = "Month", -->
<!--       y = "log(sales_millions)", -->
<!--       title = "Logarithm of Other Gen. Merch. Sales" -->
<!--     ) + -->
<!--     theme_minimal() + -->
<!--     theme( -->
<!--       plot.title = element_text(hjust = 0.5) -->
<!--     ) -->

<!-- retail_plot_raw | retail_plot_log -->
<!-- ``` -->



<!-- <!-- Check Your Understanding --> -->

<!-- ::: {.callout-tip icon=false title="Check Your Understanding"} -->

<!-- Use the retail sales data to do the following. -->

<!-- -   Select an appropriate fitted model using the AIC, AICc, or BIC critera. -->
<!-- -   Use the residuals to determine the appropriate correction for the data. -->
<!-- -   Forecast the data for the next 5 years. -->
<!-- -   Apply the appropriate correction to the forecasted values. -->
<!-- -   Plot the fitted (forecasted) values along with the time series. -->

<!-- ::: -->












### Small-Group Activity: Industrial Electricity Consumption in Texas


These data represent the amount of electricity used each month for industrial applications in Texas.

```{r}
#| code-fold: true
#| code-summary: "Show the code"

elec_ts <- rio::import("https://byuistats.github.io/timeseries/data/electricity_tx.csv") |>
  dplyr::select(-comments) |>
  mutate(month = my(month)) |>
  mutate(
    t = 1:n(),
    std_t = (t - mean(t)) / sd(t)
  ) |>
  mutate(
    cos1 = cos(2 * pi * 1 * t / 12),
    cos2 = cos(2 * pi * 2 * t / 12),
    cos3 = cos(2 * pi * 3 * t / 12),
    cos4 = cos(2 * pi * 4 * t / 12),
    cos5 = cos(2 * pi * 5 * t / 12),
    cos6 = cos(2 * pi * 6 * t / 12),
    sin1 = sin(2 * pi * 1 * t / 12),
    sin2 = sin(2 * pi * 2 * t / 12),
    sin3 = sin(2 * pi * 3 * t / 12),
    sin4 = sin(2 * pi * 4 * t / 12),
    sin5 = sin(2 * pi * 5 * t / 12)
  ) |>
  as_tsibble(index = month)

elec_plot_raw <- elec_ts |>
    autoplot(.vars = megawatthours) +
    labs(
      x = "Month",
      y = "Megawatt-hours",
      title = "Texas' Industrial Electricity Use"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

elec_plot_log <- elec_ts |>
    autoplot(.vars = log(megawatthours)) +
    labs(
      x = "Month",
      y = "log(Megwatt-hours)",
      title = "Log of Texas' Industrial Electricity Use"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

elec_plot_raw | elec_plot_log
```

<!-- Check Your Understanding -->

::: {.callout-tip icon=false title="Check Your Understanding"}

Use the Texas industrial electricity consumption data to do the following.

-   Select an appropriate fitted model using the AIC, AICc, or BIC critera.
-   Use the residuals to determine the appropriate correction for the data.
<!-- -   Forecast the data for the next 5 years. -->
<!-- -   Apply the appropriate correction to the forecasted values. -->
<!-- -   Plot the fitted (forecasted) values along with the time series. -->

:::
















## Homework Preview (5 min)

-   Review upcoming homework assignment
-   Clarify questions




::: {.callout-note icon=false}

## Download Homework

<a href="https://byuistats.github.io/timeseries/homework/homework_6_2.qmd" download="homework_6_2.qmd"> homework_6_2.qmd </a>

:::





<!-- <a href="javascript:showhide('Solutions1')" -->
<!-- style="font-size:.8em;">Small-Group Activity: Retail Sales</a> -->

<!-- ::: {#Solutions1 style="display:none;"} -->

<!-- <!-- Check Your Understanding --> -->

<!-- ::: {.callout-tip icon=false title="Check Your Understanding"} -->

<!-- Use the retail sales data to do the following. -->

<!-- -   Select an appropriate fitted model using the AIC, AICc, or BIC critera. -->
<!-- -   Use the residuals to determine the appropriate correction for the data. -->
<!-- -   Forecast the data for the next 5 years. -->
<!-- -   Apply the appropriate correction to the forecasted values. -->
<!-- -   Plot the fitted (forecasted) values along with the time series. -->

<!-- ::: -->

<!-- ```{r} -->
<!-- #| label: ExponentialQuadraticFull1a -->
<!-- #| code-fold: true -->
<!-- #| code-summary: "Show the code" -->

<!-- retail_full_quad_lm <- retail_ts |> -->
<!--   model(retail_full_quad = TSLM(log(sales_millions) ~ std_t + I(std_t^2) +  -->
<!--     sin1 + cos1 + sin2 + cos2 + sin3 + cos3  -->
<!--     + sin4 + cos4 + sin5 + cos5 + cos6 )) # Note sin6 is omitted -->
<!-- retail_full_quad_lm |> -->
<!--   tidy() |> -->
<!--   mutate(sig = p.value < 0.05)  -->

<!-- retail_resid_df <- retail_full_quad_lm |>  -->
<!--   residuals() |>  -->
<!--   as_tibble() |>  -->
<!--   dplyr::select(.resid) |> -->
<!--   rename(x = .resid)  -->

<!-- retail_resid_df |> -->
<!--   mutate(density = dnorm(x, mean(retail_resid_df$x), sd(retail_resid_df$x))) |> -->
<!--   ggplot(aes(x = x)) + -->
<!--     geom_histogram(aes(y = after_stat(density)), -->
<!--         color = "white", fill = "#56B4E9", binwidth = 0.02) + -->
<!--     geom_line(aes(x = x, y = density)) + -->
<!--     theme_bw() + -->
<!--     labs( -->
<!--       x = "Values", -->
<!--       y = "Frequency", -->
<!--       title = "Histogram of Residuals from the Full Quadratic Model" -->
<!--     ) + -->
<!--     theme( -->
<!--       plot.title = element_text(hjust = 0.5) -->
<!--     ) -->

<!-- skewness(retail_resid_df$x) -->
<!-- ``` -->

<!-- ## Small-Group Activity: Fitting ARMA Models xxxxxxxxxxxxXXXXXXXXXXxxxxxxxxx (xxx min) -->



<!-- <!-- Check your Understanding --> -->

<!-- ::: {.callout-tip icon=false title="Check Your Understanding"} -->

<!-- -   Write Equation (6.1) in terms of the backward shift operator. Your answer will be of the form: -->

<!-- $$ -->
<!--   x_t  -->
<!--     = (\text{some}~q^{th}~\text{degree polynomial in}~\mathbf{B}) w_t -->
<!--     = \phi_q(\mathbf{B}) w_t -->
<!-- $$ -->

<!-- ::: -->


<!-- ::: {.callout-caution icon=false title="Note"} -->

<!-- An $MA(q)$ process is comprised of a finite summation of stationary white noise terms. Hence, an $MA(q)$ process will be stationary with a time-invariante mean and autocovariance. -->

<!-- The mean and variance of $\{x_t\}$ are easily derived. The mean must be zero, because each term is a sum of scaled white noise terms with mean zero. -->

<!-- The variance of an $MA(q)$ process is ${ \sigma_w^2 \left( 1 + \beta_1^2 + \beta_2^2 + \beta_3^2 + \cdots + \beta_{q-1}^2 + \beta_q^2 \right) }$. This can be seen, because the white noise terms are independent with the same variance. -->

<!-- So, the autocorrelation function is -->

<!-- $$ -->
<!--   \rho(k) = -->
<!--   cor(x_t, x_{t+k}) = -->
<!--     \begin{cases} -->
<!--       1, & k=0 \\ -->
<!--       ~\\ -->
<!--       \dfrac{ \sum\limits_{i=0}^{q-k} \beta_i \beta_{i+k} }{ \sum\limits_{i=0}^q \beta_i^2 }, & k = 1, 2, \ldots, q \\ -->
<!--       ~\\ -->
<!--       0, & k > q -->
<!--     \end{cases} -->
<!-- $$ -->
<!-- where $\beta_0 = 1$. -->

<!-- Note that the autocorrelation function is zero if $k>q$, because $x_t$ and $x_{t+k}$ would be independent weighted summations of white noise processes and hence the covariance between them would be zero. -->

<!-- ::: -->

<!-- We now define an invertible $MA$ process. -->

<!-- ::: {.callout-note icon=false title="Definition of an Invertible $MA$ Process"} -->

<!-- An $MA$ process is said to be **invertible** if it can be expressed as a stationary autoregressive process (of possibly infinite order) with no error term. -->

<!-- ::: -->

<!-- #### Example of an Invertible MA Process -->

<!-- Recall that  -->
<!-- $$ -->
<!--   (1-x)(1 + x + x^2 + x^3 + \cdots) = 1 -->
<!-- $$ -->

<!-- or,  -->

<!-- $$ -->
<!--   (1-x)^{-1} = (1 + x + x^2 + x^3 + \cdots) -->
<!-- $$ -->
<!-- if $|x|<1$. -->

<!-- Now, note that the $MA$ process  -->

<!-- $$ -->
<!--   x_t = \left( 1 - \beta \mathbf{B} \right) w_t -->
<!-- $$ -->

<!-- can be written as: -->

<!-- \begin{align*} -->
<!--   \left( 1 - \beta \mathbf{B} \right)^{-1} x_t &= w_t \\ -->
<!--   \left( 1 + \beta \mathbf{B} + \beta^2 \mathbf{B}^2 + \beta^3 \mathbf{B}^3 + \cdots \right) x_t &= w_t \\ -->
<!--   x_t + \beta x_{t-1} + \beta^2 x_{t-2} + \beta^3 x_{t-3} + \cdots &= w_t \\ -->
<!--   x_t &= \left( -\beta x_{t-1} - \beta^2 x_{t-2} - \beta^3 x_{t-3} - \cdots \right) + w_t -->
<!-- \end{align*} -->

<!-- assuming that $|\beta|<1$. Note that this series will not converge unless $|\beta|<1$.  -->

<!-- We have just shown that the $MA$ process  -->
<!-- $$ -->
<!--   x_t = \left( 1 - \beta \mathbf{B} \right) w_t -->
<!-- $$ -->
<!-- is invertible. -->


<!-- ::: {.callout-note icon=false title="Theorem: Invertibility of an $MA(q)$ Process"} -->

<!-- The $MA(q)$ process  -->
<!-- $$ -->
<!--   x_t = \phi_q(\mathbf{B}) w_t -->
<!-- $$ -->
<!-- will be invertible if the solutions to the equation -->
<!-- $$ -->
<!--   \phi_q(\mathbf{B}) = 0 -->
<!-- $$ -->
<!-- are all greater than 1 in absolute value. -->

<!-- ::: -->

<!-- <a id="FittedModelWillBeInvertible">Does</a> this remind you of the test for the stationarity of an $AR(p)$ model? -->

<!-- Note that the autocovariance function (acvf) will identify a unique $MA(q)$ process only if the process is invertible. Fortunately, the algorithm R uses to estimate an $MA(q)$ process always leads to an invertible model. -->


<!-- ## Class Activity: Simulating an $MA(q)$ Model (5 min) -->

<!-- The textbook gives a simulation of an $MA(3)$ process: -->

<!-- $$ -->
<!--   x_t = w_t + \beta_1 w_{t-1} + \beta_2 w_{t-2} + \beta_3 w_{t-3} -->
<!-- $$ -->

<!-- where $\beta_1 = 0.7$, $\beta_1 = 0.5$, and $\beta_3 = 0.2$. This shiny app allows you to simulate from this process. -->







<!-- ![Shiny App goes here](images/simulation_goes_here.png) -->





<!-- <!-- The code below is replaced by the shiny app above. --> -->



<!-- <!-- ```{r} --> -->
<!-- <!-- #| code-fold: true --> -->
<!-- <!-- #| code-summary: "Show the code" --> -->

<!-- <!-- pacman::p_load("tsibble", "fable", "feasts", --> -->
<!-- <!--     "tsibbledata", "fable.prophet", "tidyverse", --> -->
<!-- <!--     "patchwork", "slider", "urca") --> -->

<!-- <!-- # define the parameters of the simulation --> -->
<!-- <!-- beta1 <- 0.7 --> -->
<!-- <!-- beta2 <- 0.5 --> -->
<!-- <!-- beta3 <- 0.2 --> -->

<!-- <!-- # function to compute the autocorrelation --> -->
<!-- <!-- rho <- function(k, beta) { --> -->
<!-- <!--   q <- length(beta) - 1 --> -->
<!-- <!--   if (k > q) ACF <- 0 else { --> -->
<!-- <!--     s1 <- 0; s2 <- 0 --> -->
<!-- <!--     for (i in 1:(q-k+1)) s1 <- s1 + beta[i] * beta[i+k] --> -->
<!-- <!--     for (i in 1:(q+1)) s2 <- s2 + beta[i]^2 --> -->
<!-- <!--     ACF <- s1 / s2} --> -->
<!-- <!--   ACF --> -->
<!-- <!-- } --> -->

<!-- <!-- # create the tibble --> -->
<!-- <!-- acf_dat <- tibble( --> -->
<!-- <!--   order = 0:10, --> -->
<!-- <!--   betas = list(c(1, beta1, beta2, beta3)), --> -->
<!-- <!--   rho.k = map2_dbl(order, betas, ~rho(.x, .y))) --> -->

<!-- <!-- # generate the autocorrelation plot --> -->
<!-- <!-- acf_dat |> --> -->
<!-- <!--   ggplot(aes(x = order, y = rho.k)) + --> -->
<!-- <!--   geom_hline(yintercept = 0, color = "darkgrey") + --> -->
<!-- <!--   geom_point() + --> -->
<!-- <!--   labs(y = expression(rho[k]), x = "lag k") + --> -->
<!-- <!--     labs( --> -->
<!-- <!--       x = "Time", --> -->
<!-- <!--       y = "ACF", --> -->
<!-- <!--       title = "Theoretical ACF for the Simulated MA(3) Process" --> -->
<!-- <!--     ) + --> -->
<!-- <!--     theme_bw() + --> -->
<!-- <!--     theme( --> -->
<!-- <!--       plot.title = element_text(hjust = 0.5) --> -->
<!-- <!--     ) --> -->
<!-- <!-- ``` --> -->

<!-- <!-- Now, we simulate data from this process. --> -->

<!-- <!-- ```{r} --> -->
<!-- <!-- #| code-fold: true --> -->
<!-- <!-- #| code-summary: "Show the code" --> -->

<!-- <!-- set.seed(1234) --> -->
<!-- <!-- dat <- tibble( --> -->
<!-- <!--   w = rnorm(1000), --> -->
<!-- <!--   betas = list(c(beta1, beta2, beta3))) |> --> -->
<!-- <!--   mutate( --> -->
<!-- <!--     w_lag = slide(w, ~.x, .before = 3, .after = -1), --> -->
<!-- <!--     w_lag = map(w_lag, ~rev(.x)), --> -->
<!-- <!--     t = 1:n()) |> --> -->
<!-- <!--   slice(-c(1:3)) |> --> -->
<!-- <!--   mutate( --> -->
<!-- <!--     lag_betas = map2_dbl( --> -->
<!-- <!--       w_lag, --> -->
<!-- <!--       betas, --> -->
<!-- <!--       \(.x, .y) sum(.x *.y)), --> -->
<!-- <!--     x = w + lag_betas) |> --> -->
<!-- <!--   tsibble::as_tsibble(index = t) --> -->
<!-- <!-- autoplot(dat, .var = x) --> -->
<!-- <!-- ``` --> -->

<!-- <!-- Here is the acf function computed from the simulated data. --> -->

<!-- <!-- ```{r} --> -->
<!-- <!-- #| code-fold: true --> -->
<!-- <!-- #| code-summary: "Show the code" --> -->

<!-- <!-- dat |> --> -->
<!-- <!--   ACF(y = x) |> --> -->
<!-- <!--   autoplot() --> -->
<!-- <!-- ``` --> -->



<!-- <!-- Check Your Understanding --> -->

<!-- ::: {.callout-tip icon=false title="Check Your Understanding"} -->

<!-- Use the simulation above to do the following: -->

<!-- -   Generate the theoretical acf plot for the $MA(3)$ model -->
<!-- $$ -->
<!--   x_t = w_t - 0.7 w_{t-1} + 0.5 w_{t-2} - 0.2 w_{t-3} -->
<!-- $$ -->
<!-- -   How does the value of the $\beta$'s affect the acf? -->
<!-- -   Simulate 1000 observations from this $MA(3)$ process. -->
<!--     -   Give the time plot of the simulated data -->
<!--     -   Plot the acf of the simulated data. -->
<!-- -   Compare the acf from the simulated data with the theoretical acf. -->

<!-- ::: -->


<!-- ## Class Activity: Identifying AR and MA Models from the ACF and PACF (5 min) -->

<!-- #### AR Process -->

<!-- Recall that on page 81, the textbook states that in general, the partial autocorrelation at lag $k$ is the $k^{th}$ coefficient of a fitted $AR(k)$ model. -->
<!-- This implies that if the underlying process is $AR(p)$, then all the coefficients $\alpha_k$ will equal 0 whenever $k>p$. So, an $AR(p)$ process will result in partial correlations that are zero after lag $p$. So, we can look at the correlogram of partial autocorrelations to determine the order of an appropriate $AR$ process to model a time series. -->

<!-- #### MA Process -->

<!-- Similarly, for an $MA(q)$ process, the coefficients $\beta_k$ will equal 0 whenever $k > q$. Hence, an $MA(q)$ process will demonstrate autocorrelations that are 0 after lag $q$. So, considering the correlogram of autocorrelations, we can assess if an $MA(q)$ model would be appropriate.  -->

<!-- Bless their hearts, the textbook authors give a bad example in Section 6.4.2. They even state that it is "not a realistic realisation." MA processes naturally arise in ratios of observed data. Multi-period asset returns (i.e. ratios of some previous term's value) tend to follow an MA process. -->

<!-- For example, if there are 252 trading days in a year, then the daily series of year-over-year returns (this year's value divided by last year's value) follows an $MA(252-1)$ process. If we are comparing values observed to those from one week ago, we would have an $MA(7-1)$ process. -->

<!-- #### Comparison  -->
<!-- ::: {.callout-note icon=false title="ACF and PACF of an $AR(p)$ Process"} -->

<!-- We can use the pacf and acf plots to assess if an $AR(p)$ or $MA(q)$ model is appropriate.  -->
<!-- For an $AR(p)$ or $MA(q)$ process, we observe the following: -->

<!-- <center> -->

<!-- |      | AR(p)                  | MA(q)                  | -->
<!-- |------|------------------------|------------------------| -->
<!-- | ACF  | Tails off              | Cuts off after lag $q$ | -->
<!-- | PACF | Cuts off after lag $p$ | Tails off              | -->

<!-- </center> -->

<!-- <!-- https://people.cs.pitt.edu/~milos/courses/cs3750/lectures/class16.pdf --> -->

<!-- <!-- |      | AR(p)                  | MA(q)                  | ARMA(p,q)                | --> -->
<!-- <!-- |------|------------------------|------------------------|--------------------------| --> -->
<!-- <!-- | ACF  | Tails off              | Cuts off after lag $q$ | Tails off                | --> -->
<!-- <!-- | PACF | Cuts off after lag $p$ | Tails off              | Tails off                | --> -->

<!-- ::: -->








<!-- # xxxyyyzzz -->

<!-- ```{r} -->
<!-- # Read in retail sales data for "all other general merchandise stores" -->
<!-- retail_ts <- rio::import("https://byuistats.github.io/timeseries/data/retail_by_business_type.parquet") |> -->
<!--   filter(naics == 445) |> # 45299 -->
<!--   filter(as_date(month) >= my("Jan 1996")) |> -->
<!--   filter(as_date(month) < my("Jan 2020")) |> -->
<!--   mutate(t = 1:n()) |> -->
<!--   mutate(std_t = (t - mean(t)) / sd(t)) |> -->
<!--   mutate( -->
<!--     cos1 = cos(2 * pi * 1 * t / 12), -->
<!--     cos2 = cos(2 * pi * 2 * t / 12), -->
<!--     cos3 = cos(2 * pi * 3 * t / 12), -->
<!--     cos4 = cos(2 * pi * 4 * t / 12), -->
<!--     cos5 = cos(2 * pi * 5 * t / 12), -->
<!--     cos6 = cos(2 * pi * 6 * t / 12), -->
<!--     sin1 = sin(2 * pi * 1 * t / 12), -->
<!--     sin2 = sin(2 * pi * 2 * t / 12), -->
<!--     sin3 = sin(2 * pi * 3 * t / 12), -->
<!--     sin4 = sin(2 * pi * 4 * t / 12), -->
<!--     sin5 = sin(2 * pi * 5 * t / 12) -->
<!--   ) |> -->
<!--   as_tsibble(index = month) -->

<!-- retail_ts |> -->
<!--   autoplot(.vars = log(sales_millions)) + -->
<!--   stat_smooth(method = "lm",  -->
<!--               formula = y ~ x,  -->
<!--               geom = "smooth", -->
<!--               se = FALSE, -->
<!--               color = "#E69F00", -->
<!--               linetype = "dotted") + -->
<!--     labs( -->
<!--       x = "Month", -->
<!--       y = "Sales (Millions of U.S. Dollars)", -->
<!--       title = paste0(retail_ts$business[1], " (", retail_ts$naics[1], ")") -->
<!--     ) + -->
<!--     theme_minimal() + -->
<!--     theme(plot.title = element_text(hjust = 0.5)) -->
<!-- ``` -->

<!-- ```{r} -->
<!-- model_combined <- retail_ts |> -->
<!--   model( -->
<!--     retail_full_cubic = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + I(std_t^3) +  -->
<!--         sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + -->
<!--         sin4 + cos4 + sin5 + cos5        + cos6 ), # Note sin6 is omitted -->
<!--     retail_full_quad = TSLM(log(sales_millions) ~ std_t + I(std_t^2) +   -->
<!--         sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + -->
<!--         sin4 + cos4 + sin5 + cos5        + cos6 ), # Note sin6 is omitted -->
<!--     retail_full_linear = TSLM(log(sales_millions) ~ std_t +  -->
<!--         sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + -->
<!--         sin4 + cos4 + sin5 + cos5        + cos6 ), # Note sin6 is omitted -->
<!--   ) -->

<!-- model_combined |> -->
<!--   select(retail_full_linear) |> -->
<!--   tidy() |> -->
<!--   mutate(sig = p.value < 0.05)  -->

<!-- glance(model_combined) |> -->
<!--   select(.model, AIC, AICc, BIC) -->

<!-- resid_ts <- model_combined |> -->
<!--   select(retail_full_linear) |> -->
<!--   residuals() |> -->
<!--   rename(.model0 = .model) -->

<!-- resid_ts |> -->
<!--   autoplot(.vars = .resid) + -->
<!--   stat_smooth(method = "lm",  -->
<!--               formula = y ~ x,  -->
<!--               geom = "smooth", -->
<!--               se = FALSE, -->
<!--               color = "#E69F00", -->
<!--               linetype = "dotted") + -->
<!--     labs( -->
<!--       x = "Month", -->
<!--       y = "Residual", -->
<!--       title = "Residual from the Linear Model", -->
<!--       subtitle = "Food and Beverage Sales (in Millions of US Dollars)" -->
<!--     ) + -->
<!--     theme_minimal() + -->
<!--     theme( -->
<!--       plot.title = element_text(hjust = 0.5), -->
<!--       plot.subtitle = element_text(hjust = 0.5) -->
<!--     ) -->

<!-- resid_ts |> -->
<!--   select(.resid) |> -->
<!--   acf() -->

<!-- resid_ts |> -->
<!--   select(.resid) |> -->
<!--   pacf() -->

<!-- gdp_ma <- resid_ts |> -->
<!--   model(arima = ARIMA(.resid ~ 1 + pdq(3,0,4) + PDQ(0, 0, 0))) -->

<!-- tidy(gdp_ma) -->

<!-- gdp_ma |> -->
<!--   residuals() |> -->
<!--   ACF() |> -->
<!--   autoplot() -->
<!-- ``` -->







<!-- # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -->

<!-- ## Small Group Activity: SectionTitle (xxx min) -->
<!-- ## Class Activity: SectionTitle (xxx min) -->


<!-- <!-- Check your Understanding --> -->

<!-- ::: {.callout-tip icon=false title="Check Your Understanding"} -->

<!-- -   Question1 -->
<!-- -   Question2 -->

<!-- ::: -->





<!-- ## Homework Preview (5 min) -->

<!-- -   Review upcoming homework assignment -->
<!-- -   Clarify questions -->




<!-- ::: {.callout-note icon=false} -->

<!-- ## Download Homework -->

<!-- <a href="https://byuistats.github.io/timeseries/homework/homework_6_2.qmd" download="homework_6_2.qmd"> homework_6_2.qmd </a> -->

<!-- ::: -->





<!-- <a href="javascript:showhide('Solutions1')" -->
<!-- style="font-size:.8em;">Class Activity</a> -->

<!-- ::: {#Solutions1 style="display:none;"} -->


<!-- ::: -->




<!-- <a href="javascript:showhide('Solutions')" -->
<!-- style="font-size:.8em;">Class Activity</a> -->

<!-- ::: {#Solutions2 style="display:none;"} -->


<!-- ::: -->




<!-- <a href="javascript:showhide('Solutions3')" -->
<!-- style="font-size:.8em;">Class Activity</a> -->

<!-- ::: {#Solutions3 style="display:none;"} -->


<!-- ::: -->