chapter_5_lesson_4.qmd

---
title: "Transformations and Non-Linear Models"
subtitle: "Chapter 5: 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>
```


## Learning Outcomes

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




## Preparation

-   Read Sections 5.7-5.8



## 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







<!-- {{< include _packages.qmd >}} -->







## Class Activity: Simulate an Exponential Trend with a Seasonal Component (15 min)


<a id="Simulation">We</a> will simulate code that has a seasonal component and impose an exponential trend.

::: {.callout-note icon=false title="Figures" collapse="false"}
<!-- Begin hiding code for Figures -->

@fig-simTSplot shows the simulated time series and the time series after the natural logarithm is applied.

```{r}
#| label: fig-simTSplot
#| 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

set.seed(12345)

n_years <- 9 # Number of years to simulate
n_months <- n_years * 12 # Number of months
sigma <- .05 # Standard deviation of random term

z_t <- rnorm(n = n_months, mean = 0, sd = sigma)

dates_seq <- seq(floor_date(now(), unit = "year"), length.out=n_months + 1, by="-1 month") |>
  floor_date(unit = "month") |> sort() |> head(n_months)

sim_ts <- tibble(
  t = 1:n_months,
  dates = dates_seq,
  random = arima.sim(model=list(ar=c(.5,0.2)), n = n_months, sd = 0.02),
  x_t = exp(2 + 0.015 * t +
      0.03 * sin(2 * pi * 1 * t / 12) + 0.04 * cos(2 * pi * 1 * t / 12) + 
      0.05 * sin(2 * pi * 2 * t / 12) + 0.03 * cos(2 * pi * 2 * t / 12) +
      0.01 * sin(2 * pi * 3 * t / 12) + 0.005 * cos(2 * pi * 3 * t / 12) +
      random
    )
  ) |>
  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),
    sin6 = sin(2 * pi * 6 * t / 12)) |>
  mutate(std_t = (t - mean(t)) / sd(t)) |>
  as_tsibble(index = dates)

plot_raw <- sim_ts |>
    autoplot(.vars = x_t) +
    labs(
      x = "Month",
      y = "x_t",
      title = "Simulated Time Series"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

plot_log <- sim_ts |>
    autoplot(.vars = log(x_t)) +
    labs(
      x = "Month",
      y = "log(x_t)",
      title = "Logarithm of Simulated Time Series"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

plot_raw | plot_log
```

We will compute the (natural) logarithm of the time series values before fitting any linear models. So, our response variable will be $\log(x_t)$, rather than $x_t$.

Even though there is no visual evidence of curvature in the trend for the logarithm of the time series, we will start by fitting a model that allows for a cubic trend. (In practice, we would probably not fit this model. However, there are some things that will occur that are helpful to understand the underlying process.)

:::
<!-- End figures -->


::: {.callout-note icon=false title="Cubic Model" collapse="false"}
<!-- Begin hiding code for cubic model -->

#### Cubic Model

After taking the (natural) logarithm of $x_t$, we fit a cubic model to the log-transformed time series.

##### Full Cubic Model

\begin{align*}
  \log(x_t) &= \beta_0 
            + \beta_1 \left( \frac{t - \mu_t}{\sigma_t} \right) 
            + \beta_2 \left( \frac{t - \mu_t}{\sigma_t} \right)^2 
            + \beta_3 \left( \frac{t - \mu_t}{\sigma_t} \right)^3 \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_4 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \beta_5 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_6 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \beta_7 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_8 \sin \left( \frac{2\pi \cdot 3 t}{12} \right)
            + \beta_9 \cos \left( \frac{2\pi \cdot 3 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_{10} \sin \left( \frac{2\pi \cdot 4 t}{12} \right)
            + \beta_{11} \cos \left( \frac{2\pi \cdot 4 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_{12} \sin \left( \frac{2\pi \cdot 5 t}{12} \right)
            + \beta_{13} \cos \left( \frac{2\pi \cdot 5 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            \phantom{+ \beta_{15} \sin \left( \frac{2\pi \cdot 6 t}{12} \right)}
            + \beta_{14} \cos \left( \frac{2\pi \cdot 6 t}{12} \right) 
      + z_t
\end{align*}


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

# Cubic model with standardized time variable

full_cubic_lm <- sim_ts |>
  model(full_cubic = TSLM(log(x_t) ~ 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
full_cubic_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)
```


Note that neither the quadratic nor the cubic terms are statistically significant in this model. We will eliminate the cubic term and fit a model with a quadratic trend.

:::
<!-- End of cubic model -->



::: {.callout-note icon=false title="Quadratic Model" collapse="false"}
<!-- Begin hiding code for quadratic model -->

#### Quadratic Model

We now fit a quadratic model to the log-transformed time series.

##### Full Quadratic Model

The full model with a quadratic trend is written as:

\begin{align*}
  \log(x_t) &= \beta_0 
            + \beta_1 \left( \frac{t - \mu_t}{\sigma_t} \right) 
            + \beta_2 \left( \frac{t - \mu_t}{\sigma_t} \right)^2 \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_3 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \beta_4 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_5 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \beta_6 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_7 \sin \left( \frac{2\pi \cdot 3 t}{12} \right)
            + \beta_8 \cos \left( \frac{2\pi \cdot 3 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_9 \sin \left( \frac{2\pi \cdot 4 t}{12} \right)
            + \beta_{10} \cos \left( \frac{2\pi \cdot 4 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_{11} \sin \left( \frac{2\pi \cdot 5 t}{12} \right)
            + \beta_{12} \cos \left( \frac{2\pi \cdot 5 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           \phantom{+ \beta_{14} \sin \left( \frac{2\pi \cdot 6 t}{12} \right)}
            + \beta_{13} \cos \left( \frac{2\pi \cdot 6 t}{12} \right) 
      + z_t
\end{align*}

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

full_quad_lm <- sim_ts |>
  model(full_quad = TSLM(log(x_t) ~ std_t + I(std_t^2) + 
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
    + sin4 + cos4 + sin5 + cos5 + cos6 )) # Note sin6 is omitted
full_quad_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05) 
```

Now, note that the quadratic term is statistically significant. It was not significant when the cubic term was included in the model, but it is now.

##### Reduced Quadratic Trend: Model 1

This model omits all the Fourier (sine and cosine) terms that are not significant in the previous model.

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

reduced_quadratic_lm1 <- sim_ts |>
  model(reduced_quadratic1 = TSLM(log(x_t) ~ std_t + I(std_t^2) + 
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3))
reduced_quadratic_lm1 |>
  tidy() |>
  mutate(sig = p.value < 0.05)
```

All the terms are statistically significant in this model.

##### Reduced Quadratic Trend: Model 2

This model only includes the Fourier series terms where $i=1$ or $i=2$.

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

reduced_quadratic_lm2 <- sim_ts |>
  model(reduced_quadratic2 = TSLM(log(x_t) ~ std_t + I(std_t^2) + 
    sin1 + cos1 + sin2 + cos2))
reduced_quadratic_lm2 |>
  tidy() |>
  mutate(sig = p.value < 0.05) 
```

As we would expect, all the terms are statistically significant. (They were all significant in the previous model, so it is not surprising that they are still significant.)

##### Reduced Quadratic Trend: Model 3

This model is reduced to include only the Fourier series terms for $i=1$.

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

reduced_quadratic_lm3 <- sim_ts |>
  model(reduced_quadratic3 = TSLM(log(x_t) ~ std_t + I(std_t^2) + 
    sin1 + cos1)) 
reduced_quadratic_lm3 |>
  tidy() |>
  mutate(sig = p.value < 0.05) 
```

All the terms in this parsimonious model are statistically significant.

:::
<!-- End of quadratic model -->



::: {.callout-note icon=false title="Linear Model" collapse="false"}
<!-- Begin hiding code for linear model -->

#### Linear Model

Even though the quadratic terms were statistically significant, there is no visual indication that there is a quadratic trend in the time series after taking the logarithm.
Hence, we will now fit a linear model to the log-transformed time series. We want to be able to compare the fit of models with a linear trend to the models with quadratic trends.

##### Full Linear Model

First, we fit a full model with a linear trend. We can express this model as:

\begin{align*}
  \log(x_t) &= \beta_0 
            + \beta_1 \left( \frac{t - \mu_t}{\sigma_t} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_2 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \beta_3 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_4 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \beta_5 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_6 \sin \left( \frac{2\pi \cdot 3 t}{12} \right)
            + \beta_7 \cos \left( \frac{2\pi \cdot 3 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_8 \sin \left( \frac{2\pi \cdot 4 t}{12} \right)
            + \beta_9 \cos \left( \frac{2\pi \cdot 4 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_{10} \sin \left( \frac{2\pi \cdot 5 t}{12} \right)
            + \beta_{11} \cos \left( \frac{2\pi \cdot 5 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            \phantom{+ \beta_{13} \sin \left( \frac{2\pi \cdot 6 t}{12} \right)}
            + \beta_{12} \cos \left( \frac{2\pi \cdot 6 t}{12} \right) 
      + z_t
\end{align*}


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

full_linear_lm <- sim_ts |>
  model(full_linear = TSLM(log(x_t) ~ std_t + 
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
    + sin4 + cos4 + sin5 + cos5 + cos6 )) # Note sin6 is omitted
full_linear_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)
```

##### Reduced Linear Trend: Model 1

This model excludes the terms that were not significant in the full model with a linear trend. 

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

reduced_linear_lm1 <- sim_ts |>
  model(reduced_linear1 = TSLM(log(x_t) ~ std_t + 
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3))
reduced_linear_lm1 |>
  tidy() |>
  mutate(sig = p.value < 0.05)
```

All of the terms are significant in this model.

##### Reduced Linear Trend: Model 2

We reduce the model to see if a more parsimonious model will suffice. This model contains a linear trend and the Fourier series terms associated with $i=1$ and $i=2$.

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

reduced_linear_lm2 <- sim_ts |>
  model(reduced_linear2 = TSLM(log(x_t) ~ std_t + 
    sin1 + cos1 + sin2 + cos2))
reduced_linear_lm2 |>
  tidy() |>
  mutate(sig = p.value < 0.05) 
```

All the terms in this model are statistically significant.

##### Reduced Linear Trend: Model 3

Finally, we fit a more reduced model that includes a linear trend and only the Fourier terms associated with $i=1$.

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

reduced_linear_lm3 <- sim_ts |>
  model(reduced_linear3 = TSLM(log(x_t) ~ std_t + 
    sin1 + cos1)) 
reduced_linear_lm3 |>
  tidy() |>
  mutate(sig = p.value < 0.05) 

```

Each of these terms is significant.

:::
<!-- End of linear model -->




::: {.callout-note icon=false title="Model Comparison" collapse="false"}
<!-- Begin hiding code for model comparison -->

#### Comparison of Fitted Models

##### AIC, AICc, and BIC

We will now compare the models we fitted above. #tbl-ModelComparison gives the AIC, AICc, and BIC of the models fitted above.

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

model_combined <- sim_ts |>
  model(
    full_cubic = TSLM(log(x_t) ~ std_t + I(std_t^2) + I(std_t^3) +
                    sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                  + sin4 + cos4 + sin5 + cos5 + cos6 ),
    full_quad = TSLM(log(x_t) ~ std_t + I(std_t^2) + 
                       sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                     + sin4 + cos4 + sin5 + cos5 + cos6 ),
    reduced_quadratic1 = TSLM(log(x_t) ~ std_t + I(std_t^2) + sin1 + cos1 + sin2 + cos2 + sin3 + cos3),
    reduced_quadratic2 = TSLM(log(x_t) ~ std_t + I(std_t^2) + sin1 + cos1 + sin2 + cos2),
    reduced_quadratic3 = TSLM(log(x_t) ~ std_t + I(std_t^2) + sin1 + cos1),
    full_linear = TSLM(log(x_t) ~ std_t + 
                         sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                       + sin4 + cos4 + sin5 + cos5 + cos6 ),
    reduced_linear1 = TSLM(log(x_t) ~ std_t + sin1 + cos1 + sin2 + cos2 + sin3 + cos3),
    reduced_linear2 = TSLM(log(x_t) ~ std_t + sin1 + cos1 + sin2 + cos2),
    reduced_linear3 = TSLM(log(x_t) ~ std_t + sin1 + cos1) 
  )

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

```{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_combined) |> 
  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()
```

The model with the lowest AIC and AICc values is the reduced quadratic trend model 1.
This can be written as:

\begin{align*}
  \log(x_t) &= \beta_0 
            + \beta_1 \left( \frac{t - \mu_t}{\sigma_t} \right) 
            + \beta_2 \left( \frac{t - \mu_t}{\sigma_t} \right)^2 \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_3 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \beta_4 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_5 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \beta_6 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_7 \sin \left( \frac{2\pi \cdot 3 t}{12} \right)
            + \beta_8 \cos \left( \frac{2\pi \cdot 3 t}{12} \right)
      + z_t
\end{align*}

The model with the lowest BIC is the reduced linear trend model 1. We express this model as:

\begin{align*}
  \log(x_t) &= \beta_0 
            + \beta_1 \left( \frac{t - \mu_t}{\sigma_t} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_2 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \beta_3 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_4 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \beta_5 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_6 \sin \left( \frac{2\pi \cdot 3 t}{12} \right)
            + \beta_7 \cos \left( \frac{2\pi \cdot 3 t}{12} \right) 
      + z_t
\end{align*}

Both of these models have very low AIC, AICc, and BIC values. 

We will take a deeper look at the residuals of these two models to assess if there is evidence of autocorrelation in the random terms. We compute the autocorrelation function and the partial autocorrelation function of the residuals for both. If there is evidence of autocorrelation, we will use the GLS algorithm to fit the models, since it will take into account the autocorrelation in the terms.

::: 
<!-- End of model comparison -->



::: {.callout-note icon=false title="Autocorrelation of the Random Component" collapse="false"}
<!-- Begin hiding code for checking autocorrelation of the random component -->

#### Investigating Autocorrelation of the Random Component

Recall that if there is autocorrelation in the random component, the standard error of the parameter estimates tends to be underestimated. We can account for this autocorrelation using Generalized Least Squares, GLS, if needed.

##### Reduced Quadratic Trend: Model 1

@fig-QuadACF illustrates the ACF of the reduced model 1 with quadratic trend.

```{r}
#| label: fig-QuadACF
#| fig-cap: "ACF of reduced model 1 with a quadratic trend"
#| code-fold: true
#| code-summary: "Show the code"

resid_q1_ts <- reduced_quadratic_lm1 |>
  residuals() 

acf(resid_q1_ts$.resid, plot=TRUE, lag.max = 25)
```

Notice that the residual correlogram indicates a positive autocorrelation in the values. This suggests that the standard errors of the regression coefficients will be underestimated, which means that some predictors can appear to be statistically significant when they are not. 


@fig-QuadPACF illustrates the PACF of the reduced model 1 with quadratic trend.

```{r}
#| label: fig-QuadPACF
#| fig-cap: "PACF of reduced model 1 with a quadratic trend"
#| code-fold: true
#| code-summary: "Show the code"

pacf(resid_q1_ts$.resid, plot=TRUE, lag.max = 25)
```

Only the first partial autocorrelation is statistically significant. The partial autocorrelation plot indicates that an $AR(1)$ model could adequately model the random component of the logarithm of the time series.
Recall that in [Chapter 5, Lesson 1](https://byuistats.github.io/timeseries/chapter_5_lesson_1.html#FittingRegModelWithPACF), we fitted a linear regression model using the value of the partial autocorrelation function for $k=1$. This helps account for the autocorrelation in the residuals.

The first few partial autocorrelation values are:

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

pacf(resid_q1_ts$.resid, plot=FALSE, lag.max = 10)
```



```{r}
#| echo: false

alphas_quad <- pacf(resid_q1_ts$.resid, plot=FALSE, lag.max = 25) 
```

The partial autocorrelation when $k=1$ is approximately `r alphas_quad$acf[1] |> round(3)`. We will use this value as we recompute the regression coefficients.

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

# Load additional packages
pacman::p_load(tidymodels, multilevelmod,
  nlme, broom.mixed)

temp_spec <- linear_reg() |>
  set_engine("gls", correlation = nlme::corAR1(0.479))

temp_gls <- temp_spec |>
  fit(log(x_t) ~ std_t + I(std_t^2) + sin1 + cos1 + sin2 + cos2 + sin3 + cos3, data = sim_ts)

tidy(temp_gls) |>
  mutate(
    lower = estimate + qnorm(0.025) * std.error,
    upper = estimate + qnorm(0.975) * std.error
  ) 
```

The quadratic term is not statistically significant in this model! After accounting for the autocorrelation in the random component, the quadratic component of the trend is not statistically significant. This is a great example of an instance where ordinary linear regression leads to errant results.

We now consider the reduced model 1 where the trend is linear.


##### Reduced Linear Trend: Model 1

@fig-LinearACF illustrates the ACF of the reduced model 1 with linear trend.

```{r}
#| label: fig-LinearACF
#| fig-cap: "ACF of reduced model 1 with a linear trend"

resid_lin1_ts <- reduced_linear_lm1 |>
  residuals() 

acf(resid_lin1_ts$.resid, plot=TRUE, lag.max = 25)
```

We observe evidence of autocorrelation in the random terms.

@fig-LinearPACF illustrates the PACF of the reduced model 1 with linear trend.

```{r}
#| label: fig-LinearPACF
#| fig-cap: "ACF of reduced model 1 with a linear trend"
#| code-fold: true
#| code-summary: "Show the code"

alphas_lin <- pacf(resid_lin1_ts$.resid, plot=FALSE, lag.max = 25) 
pacf(resid_lin1_ts$.resid, plot=TRUE, lag.max = 25)
```

We will use the PACF when $k=1$ to apply the GLS algorithm.
The first few partial autocorrelation values are:

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

pacf(resid_lin1_ts$.resid, plot=FALSE, lag.max = 10)
```

```{r}
#| echo: false

alphas_lin <- pacf(resid_lin1_ts$.resid, plot=FALSE, lag.max = 25) 
```

The partial autocorrelation when $k=1$ is approximately `r alphas_lin$acf[1] |> round(3)`. We will use this value as we recompute the regression coefficients.

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

# Load additional packages
pacman::p_load(tidymodels, multilevelmod,
  nlme, broom.mixed)

temp_spec <- linear_reg() |>
  set_engine("gls", correlation = nlme::corAR1(0.497))

temp_gls <- temp_spec |>
  fit(log(x_t) ~ std_t + sin1 + cos1 + sin2 + cos2 + sin3 + cos3, data = sim_ts)

tidy(temp_gls) |>
  mutate(
    lower = estimate + qnorm(0.025) * std.error,
    upper = estimate + qnorm(0.975) * std.error
  ) 
```


@fig-finalFittedPlot illustrates the original time series (in black) and the fitted model (in blue). 
For reference, a dotted line illustrating the simple least squares line is plotted on this figure for reference. It helps highlight the exponential shape of the trend.

```{r}
#| label: fig-finalFittedPlot
#| fig-cap: "Time plot of the time series and the fitted linear regression model"
#| code-fold: true
#| code-summary: "Show the code"

forecast_df <- reduced_linear_lm1 |> 
  forecast(sim_ts) |> # computes the anti-log of the predicted values and returns them as .mean
  as_tibble() |> 
  dplyr::select(std_t, .mean) |> 
  rename(pred = .mean)

sim_ts |>
  left_join(forecast_df, by = "std_t") |>
  as_tsibble(index = dates) |>
  autoplot(.vars = x_t) +
  geom_smooth(method = "lm", formula = 'y ~ x', se = FALSE, color = "#E69F00", linewidth = 0.5, linetype = "dotted") +
  geom_line(aes(y = pred), color = "#56B4E9", alpha = 0.75) +
    labs(
      x = "Month",
      y = "Simulated Time Series",
      title = "Time Plot of Simulated Time Series with an Exponential Trend",
      subtitle = "Predicted values based on the full cubic model are given in blue"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0.5)
    )
```


:::
<!-- End of Checking Autocorrelation of Residuals -->




::: {.callout-note icon=false title="Comparison of the Fitted Coefficients to the Simulation Parameters" collapse="false"}
<!-- Begin hiding code for the comparison -->

#### Comparison of Model Coefficients

```{r}
#| echo: false

coeffs <- tidy(temp_gls) |> 
  select(term, estimate) 
```

Note that the mean of the $x_t$ values is $\bar x_t = `r mean(sim_ts$x_t) |> round(3)`$, and the standard deviation is $s_t = `r sd(sim_ts$x_t) |> round(3)`$.

The fitted model is:
\begin{align*}
  \hat x_t 
      &= e^{
          `r coeffs |> filter(term == "(Intercept)") |> select(estimate) |> pull() |> round(3)` 
          + `r coeffs |> filter(term == "std_t") |> select(estimate) |> pull() |> round(3)` 
                ~ \left( \frac{t - `r mean(sim_ts$x_t) |> round(3)`}{`r sd(sim_ts$x_t) |> round(3)`} \right)
          + `r coeffs |> filter(term == "sin1") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos1") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin2") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos2") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin3") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 3 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos3") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 3 t }{ 12 } \right)
     }
     \\
     &= e^{
          `r coeffs |> filter(term == "(Intercept)") |> select(estimate) |> pull() |> round(3)` 
          + `r coeffs |> filter(term == "std_t") |> select(estimate) |> pull() |> round(3)` 
                ~ \left( \frac{t}{`r sd(sim_ts$x_t) |> round(3)`} \right)
          - `r coeffs |> filter(term == "std_t") |> select(estimate) |> pull() |> round(3)` 
                ~ \left( \frac{`r mean(sim_ts$x_t) |> round(3)`}{`r sd(sim_ts$x_t) |> round(3)`} \right)
          + `r coeffs |> filter(term == "sin1") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos1") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin2") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos2") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin3") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 3 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos3") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 3 t }{ 12 } \right)
     }
     \\
     &= e^{
          `r coeffs |> filter(term == "(Intercept)") |> select(estimate) |> pull() |> round(3)` 
          + `r ((coeffs |> filter(term == "std_t") |> select(estimate) |> pull() ) / sd(sim_ts$x_t)) |> round(3)` ~ t
          - `r ((coeffs |> filter(term == "std_t") |> select(estimate) |> pull()) * (mean(sim_ts$x_t) / sd(sim_ts$x_t))) |> round(3)`
          + `r coeffs |> filter(term == "sin1") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos1") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin2") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos2") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin3") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 3 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos3") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 3 t }{ 12 } \right)
     }
     \\
     &= e^{
          `r ((coeffs |> filter(term == "(Intercept)") |> select(estimate) |> pull()) - ((coeffs |> filter(term == "std_t") |> select(estimate) |> pull()) * (mean(sim_ts$x_t) / sd(sim_ts$x_t)))) |> round(3)`
          + `r ((coeffs |> filter(term == "std_t") |> select(estimate) |> pull() ) / sd(sim_ts$x_t)) |> round(3)` ~ t
          + `r coeffs |> filter(term == "sin1") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos1") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin2") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos2") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "sin3") |> select(estimate) |> pull() |> round(3)` 
                ~ \sin \left( \frac{2 \pi \cdot 3 t }{ 12 } \right) 
          + `r coeffs |> filter(term == "cos3") |> select(estimate) |> pull() |> round(3)` 
                ~ \cos \left( \frac{2 \pi \cdot 3 t }{ 12 } \right)
     }
\end{align*}

Notice how well this matches the original model used to simulate the data.

\begin{align*}
  x_t &= e^{
          2 + 0.015 t +
       0.03 ~ \sin \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) + 0.04 ~ \cos \left( \frac{2 \pi \cdot 1 t }{ 12 } \right) + 
       0.05 ~ \sin \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) + 0.03 ~ \cos \left( \frac{2 \pi \cdot 2 t }{ 12 } \right) +
       0.01 ~ \sin \left( \frac{2 \pi \cdot 3 t }{ 12 } \right) + 0.005 ~ \cos \left( \frac{2 \pi \cdot 3 t }{ 12 } \right) +
        z_t
     }
\end{align*}
where $z_t = 0.5 z_{t-1} + 0.2 z_{t-2} + w_t$ and $w_t$ is a white noise process with standard deviation of 0.02.

:::
<!-- End of the big reveal -->



## Small-Group Activity: Retail Sales (20 min)

@fig-RetailSalesGeneralMerch gives the total sales (in millions of U.S. dollars) for the category "all other general merchandise stores (45299)."

```{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 == 45299) |>
  filter(as_date(month) >= my("Jan 1998")) |>
  as_tsibble(index = month)

retail_ts |>
  autoplot(.vars = sales_millions) +
    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))

```


<!-- Check Your Understanding -->

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

Use @fig-RetailSalesGeneralMerch to explain the following questions to a partner.

-   What is the shape of the trend of this time series?
-   Which decomposition would be more appropriate: additive or multiplicative? Justify your answer.
-   Apply the appropriate transformation to the time series. 
-   Fit appropriate  models utilizing the Fourier terms for seasonal components.
-   Determine the "best" model for these data. Justify your decision.
<!-- -   Check for autocorrelation in the random terms. -->
<!-- -   Use GLS to fit the model, if there is evidence of autocorrelation in the random component. -->
<!-- -   Refine your choice of the "best" model, if necessary. Justify your actions. -->
-   Plot the fitted values and the time series on the same figure.

Do the following to explore a cubic trend model with all the seasonal Fourier terms.

-   Fit a full model with a cubic trend and the logarithm of the time series for the response, including all the Fourier seasonal terms. (Be sure to include the linear and quadratic components of the trend.)
-   Fit a full model with a cubic trend and the logarithm of the time series for the response, but use indicator seasonal variables.  (Be sure to include the linear and quadratic components of the trend.)
-   Compare the coefficients on the linear, quadratic, and cubic terms between the two models above. Why would we observe this result?
-   Why does the coefficient on the intercept term differ in the two models?

:::








## Class Activity: Simulated Non-Linear Series (10 min)

In Section 5.8 of the textbook, we are introduced to the possibility that a time series could have an exponential trend but also exhibit negative values. In this case, the logarithm of the negative values is not defined, so we need a different approach to fit the model.

We will fit the non-linear time series

$$
x_t = -c_0 + e^{\alpha_0 + \alpha_1 t} + z_t
$$

where the time series can assume negative values. 

We present the simulation from the textbook here:

```{r}
#| label: fig-TextbokSimulationTimePlot
#| fig-cap: "Time plot of a simulated non-linear series"
#| code-fold: true
#| code-summary: "Show the code"

set.seed(1)
dat <- tibble(w = rnorm(100, sd = 10)) |>
    mutate(
        Time = 1:n(),
         z = purrr::accumulate2(
            lag(w), w, 
            \(acc, nxt, w) 0.7 * acc + w,
            .init = 0)[-1],
          x = exp(1 + 0.05 * Time) + z) |>
    tsibble::as_tsibble(index = Time)

autoplot(dat, .var = x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "Time",
    y = "Value",
    title = "Time Plot of Simulated Non-Linear Series"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

```

We will fit the non-linear model given above.

```{r}
equation <- x ~ exp(alp0 + alp1 * Time)

dat_nls <- nls(equation, data = dat,
  start = list(alp0 = 0.1, alp1 = 0.5))

summary(dat_nls)$parameters
```

```{r}
#| echo: false

params <- summary(dat_nls)$parameters |> as_tibble() |> select(Estimate) |> pull() |> round(3)
```

The model we fitted was 

$$
  x_t = e^{1 + 0.05 t} + z_t
$$

where 

$$
  z_t = 0.7 z_{t-1} + w_t
$$
and $w_t$ is a white noise process. 

Our fitted model is:
$$
  \hat x_t = e^{`r params[1]` + `r params[2]` t}
$$

Notice that our estimates are quite close to the parameters used in the simulation.



<!-- Check your Understanding -->

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

-   Modify the time series and repeat the simulation.
-   Fit your model using the `nls` function in R.
-   Compare the parameter estimates with the actual values from your simulation.

:::





## 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_5_4.qmd" download="homework_5_4.qmd"> homework_5_4.qmd </a>

:::





















<!-- Solutions at the bottom of the page -->

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

::: {#Solutions style="display:none;"}

#### Small-Group Activity: Retail Sales Code

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


::: {.callout-caution icon=false title="Figures" collapse="false"}

@fig-RetailSalesGeneralMerch2 gives the total sales (in millions of U.S. dollars) for the category "all other general merchandise stores (45299)," beginning with January 1998.

```{r}
#| label: fig-RetailSalesGeneralMerch2
#| 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 == 45299) |>
  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)
  ) |>
  filter(as_date(month) >= my("Jan 1998")) |>
  as_tsibble(index = month)

retail_ts |>
  autoplot(.vars = sales_millions) +
    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}
#| 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

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)
    )

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)
    )

plot_raw | plot_log
```

:::
<!-- Figures -->


::: {.callout-caution icon=false title="Cubic Trend" collapse="false"}

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

# Cubic model with standardized time variable

full_cubic_lm <- retail_ts |>
  model(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
full_cubic_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)
```

Note that each of the Fourier terms is statistically significant. 

:::
<!-- Cubic -->


::: {.callout-caution icon=false title="Quadratic Trend" collapse="false"}

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

full_quad_lm <- retail_ts |>
  model(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
full_quad_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05) 
```


:::
<!-- Quadratic -->


::: {.callout-caution icon=false title="Linear Trend" collapse="false"}

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

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

:::
<!-- Linear -->


::: {.callout-caution icon=false title="Model Comparison" collapse="false"}

We will now compare the models we fitted above. #tbl-ModelComparison2 gives the AIC, AICc, and BIC of the models fitted above. In addition, other models with a reduced number of Fourier terms are included. For example, the model labeled `reduced_quadratic_fourier_i5` includes linear and quadratic trend terms but also the Fourier terms where $i \le 5$.

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

model_combined <- retail_ts |>
  model(
    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 ),
    full_quadratic = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
                       sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                     + sin4 + cos4 + sin5 + cos5 + cos6 ),
    reduced_quadratic_fourier_i5 = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
                       sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                     + sin4 + cos4 + sin5 + cos5),
    reduced_quadratic_fourier_i4 = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
                       sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                     + sin4 + cos4),
    # reduced_quadratic_fourier_i3 = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
    #                    sin1 + cos1 + sin2 + cos2 + sin3 + cos3),
    # reduced_quadratic_fourier_i2 = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
    #                    sin1 + cos1 + sin2 + cos2),
    # reduced_quadratic_fourier_i1 = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
    #                    sin1 + cos1),
    full_linear = TSLM(log(sales_millions) ~ std_t + 
                         sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                       + sin4 + cos4 + sin5 + cos5 + cos6 ),
    reduced_linear_fourier_i5 = TSLM(log(sales_millions) ~ std_t + 
                       sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                     + sin4 + cos4 + sin5 + cos5),
    reduced_linear_fourier_i4 = TSLM(log(sales_millions) ~ std_t + 
                       sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
                     + sin4 + cos4),
    # reduced_linear_fourier_i3 = TSLM(log(sales_millions) ~ std_t + 
    #                    sin1 + cos1 + sin2 + cos2 + sin3 + cos3),
    # reduced_linear_fourier_i2 = TSLM(log(sales_millions) ~ std_t + 
    #                    sin1 + cos1 + sin2 + cos2),
    # reduced_linear_fourier_i1 = TSLM(log(sales_millions) ~ std_t + 
    #                    sin1 + cos1)
  )

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

```{r}
#| label: tbl-ModelComparison2
#| tbl-cap: "Comparison of the AIC, AICc, and BIC values for the models fitted to the logarithm of the retail sales time series."
#| 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()
```

The model with the lowest AIC, AICc, and BIC values is the full quadratic model.
This can be written as:

\begin{align*}
  \log(x_t) &= \beta_0 
            + \beta_1 \left( \frac{t - \mu_t}{\sigma_t} \right) 
            + \beta_2 \left( \frac{t - \mu_t}{\sigma_t} \right)^2 \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_3 \sin \left( \frac{2\pi \cdot 1 t}{12} \right) 
            + \beta_4 \cos \left( \frac{2\pi \cdot 1 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_5 \sin \left( \frac{2\pi \cdot 2 t}{12} \right) 
            + \beta_6 \cos \left( \frac{2\pi \cdot 2 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_7 \sin \left( \frac{2\pi \cdot 3 t}{12} \right)
            + \beta_8 \cos \left( \frac{2\pi \cdot 3 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_9 \sin \left( \frac{2\pi \cdot 4 t}{12} \right)
            + \beta_{10} \cos \left( \frac{2\pi \cdot 4 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            + \beta_{11} \sin \left( \frac{2\pi \cdot 5 t}{12} \right)
            + \beta_{12} \cos \left( \frac{2\pi \cdot 5 t}{12} \right) \\
      & ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            \phantom{+ \beta_{13} \sin \left( \frac{2\pi \cdot 6 t}{12} \right)}
            + \beta_{13} \cos \left( \frac{2\pi \cdot 6 t}{12} \right)
      + z_t
\end{align*}

:::
<!-- model comparison -->


::: {.callout-caution icon=false title="Autocorrelation of the Random Component" collapse="false"}

We check for autocorrelation in the random component to determine if using GLS is warranted.

@fig-FullQuadACF illustrates the ACF of the full model with a quadratic trend.

```{r}
#| label: fig-FullQuadACF
#| fig-cap: "ACF of the full model with a quadratic trend"

full_quad_ts <- full_quad_lm |>
  residuals() 

acf(full_quad_ts$.resid, plot=TRUE, lag.max = 25)
```

We observe evidence of autocorrelation in the random terms. In fact, there is something happening on an annual 

@fig-FullQuadPACF illustrates the PACF of the reduced model 1 with linear trend.

```{r}
#| label: fig-FullQuadPACF
#| fig-cap: "PACF of the full model with a quadratic trend"
#| code-fold: true
#| code-summary: "Show the code"

alphas_quad <- pacf(full_quad_ts$.resid, plot=FALSE, lag.max = 25) 
pacf(full_quad_ts$.resid, plot=TRUE, lag.max = 25)
```

:::
<!-- Autocorrelation of the random component -->


::: {.callout-caution icon=false title="Applying Genearlized Least Squares, GLS" collapse="false"}

Recall that in [Chapter 5, Lesson 1](https://byuistats.github.io/timeseries/chapter_5_lesson_1.html#FittingRegModelWithPACF), we fitted a linear regression model using the value of the partial autocorrelation function for $k=1$. This helps account for the autocorrelation in the residuals.

We will use the PACF when $k=1$ to apply the GLS algorithm.
The first few partial autocorrelation values are:

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

pacf(full_quad_ts$.resid, plot=FALSE, lag.max = 10)
```

The partial autocorrelation when $k=1$ is approximately `r alphas_quad$acf[1] |> round(3)`. We will use this value as we recompute the regression coefficients.

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

# Load additional packages
pacman::p_load(tidymodels, multilevelmod,
  nlme, broom.mixed)

temp_spec <- linear_reg() |>
  set_engine("gls", correlation = nlme::corAR1(0.497))

temp_gls <- temp_spec |>
  fit(log(sales_millions) ~ std_t + sin1 + cos1 + sin2 + cos2 + sin3 + cos3, data = retail_ts)

tidy(temp_gls) |>
  mutate(
    lower = estimate + qnorm(0.025) * std.error,
    upper = estimate + qnorm(0.975) * std.error
  ) 
```


@fig-finalFittedPlot2 illustrates the original time series (in black) and the fitted model (in blue). 
For reference, a dotted line illustrating the simple least squares line is plotted on this figure for reference. It helps highlight the exponential shape of the trend.

```{r}
#| label: fig-finalFittedPlot2
#| 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"

forecast_df <- full_quad_lm |> 
  forecast(retail_ts) |>  # computes the anti-log of the predicted values and returns them as .mean
  as_tibble() |> 
  dplyr::select(std_t, .mean) |> 
  rename(pred = .mean)

retail_ts |>
  left_join(forecast_df, by = "std_t") |>
  as_tsibble(index = month) |>
  autoplot(.vars = sales_millions) +
  geom_smooth(method = "lm", formula = 'y ~ x', se = FALSE, color = "#E69F00", linewidth = 0.5, linetype = "dotted") +
  geom_line(aes(y = pred), color = "#56B4E9", alpha = 0.75) +
    labs(
      x = "Month",
      y = "Simulated Time Series",
      title = "Time Plot of Simulated Time Series with an Exponential Trend",
      subtitle = "Predicted values based on the full cubic model are given in blue"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0.5)
    )
```


:::
<!-- GLS -->

::: {.callout-caution icon=false title="Cubic Trend Revisited" collapse="false"}

Recall the full cubic model:

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

# Cubic model with standardized time variable

full_cubic_lm <- retail_ts |>
  model(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
full_cubic_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)
```


We now fit a model where the Fourier terms are replaced with indicator variables.

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

full_cubic_lm2 <- retail_ts |>
  model(full_quad = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + I(std_t^3) + 
    factor(month(month)))) # factor variable for the month number
full_cubic_lm2 |>
  tidy() |>
  mutate(sig = p.value < 0.05) 
```

Note that the coefficients on the linear, quadratic, and cubic terms are identical between the two models. The other terms have different coefficients. 

The constant (intercept) term differs in the two models. This is because it is affected by the choice of how the seasonal variables are parameterized.

:::
<!-- Cubic2 -->

:::
<!-- End of check your understanding colored code -->

:::
<!-- End of folded solution code -->