---
title: "STAT481/581: Introduction to Time Series Analysis"
author: "Dynamic regression models"
date: "OTexts.org/fpp3/"
toc: true
colortheme: monashwhite
output:
  binb::monash:
    fig_width: 7
    fig_height: 3.5
    includes:
      in_header: header.tex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, cache=TRUE, warning=FALSE, message=FALSE)
options(digits=4, width=60)
library(tidyverse)
library(fable)
library(tsibble)
library(tsibbledata)
library(feasts)
library(lubridate)
us_change <- read_csv("data/uschange.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time) %>%
  select(Income, Consumption)

vic_elec_daily <- vic_elec %>%
  filter(year(Time) == 2014) %>%
  index_by(Date = date(Time)) %>%
  summarise(
    Demand = sum(Demand)/1e3,
    Temperature = max(Temperature),
    Holiday = any(Holiday)
  ) %>%
  mutate(Day_Type = case_when(
    Holiday ~ "Holiday",
    wday(Date) %in% 2:6 ~ "Weekday",
    TRUE ~ "Weekend"
  ))
```

# Regression with ARIMA errors

## Regression with ARIMA errors
\fontsize{13}{15}\sf

\begin{block}{Regression models}\vspace*{-0.2cm}
\[
  y_t = \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \varepsilon_t,
\]
\end{block}

  * $y_t$ modeled as function of $k$ explanatory variables
$x_{1,t},\dots,x_{k,t}$.
  * In regression, we assume that $\varepsilon_t$ was  WN.
  * Now we want to allow $\varepsilon_t$ to be autocorrelated.
\vspace*{0.3cm}
 
\begin{alertblock}{Example: ARIMA(1,1,1) errors}\vspace*{-0.2cm}
\begin{align*}
  y_t &= \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t,\\
      & (1-\phi_1B)(1-B)\eta_t = (1+\theta_1B)\varepsilon_t,
\end{align*}
\end{alertblock}
\rightline{where $\varepsilon_t$ is white noise.}

## Residuals and errors

\begin{alertblock}{Example: $\eta_t$ = ARIMA(1,1,1)}\vspace*{-0.2cm}
\begin{align*}
  y_t &= \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t,\\
      & (1-\phi_1B)(1-B)\eta_t = (1+\theta_1B)\varepsilon_t,
\end{align*}\end{alertblock} 

  * Be careful in distinguishing $\eta_t$ from $\varepsilon_t$.
  * Only the errors $\varepsilon_t$ are assumed to be white noise.
  * In ordinary regression, $\eta_t$ is assumed to be white noise and so $\eta_t = \varepsilon_t$.

## Estimation
\fontsize{14}{14}\sf

If we minimize $\sum \eta_t^2$ (by using ordinary regression):

  1. Estimated coefficients $\hat{\beta}_0,\dots,\hat{\beta}_k$ are no longer optimal as some information ignored;
  2. Statistical tests associated with the model (e.g., t-tests on the coefficients) are incorrect.
  3. $p$-values for coefficients usually too small (``spurious regression'').
  4. AIC of fitted models misleading.

 

 * Minimizing $\sum \varepsilon_t^2$ avoids these problems.
 * Maximizing likelihood similar to minimizing $\sum \varepsilon_t^2$.

## Stationarity

\begin{block}{Regression with ARMA errors}
\[
  y_t = \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t,
\]
where $\eta_t$ is an ARMA process.
\end{block}

  * All variables in the model must be stationary.
  * If we estimate the model while any of these are non-stationary, the estimated coefficients can be incorrect.
  * Difference variables until all stationary.
  * If necessary, apply same differencing to all variables.

## Stationarity

\begin{block}{Model with ARIMA(1,1,1) errors}\vspace*{-0.2cm}
\begin{align*}
  y_t &= \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t,\\
      & (1-\phi_1B)(1-B)\eta_t = (1+\theta_1B)\varepsilon_t,
\end{align*}
\end{block} 

\begin{block}{Equivalent to model with ARIMA(1,0,1) errors}\vspace*{-0.2cm}
\begin{align*}
  y'_t &= \beta_1 x'_{1,t} + \dots + \beta_k x'_{k,t} + \eta'_t,\\
       & (1-\phi_1B)\eta'_t = (1+\theta_1B)\varepsilon_t,
\end{align*}
\end{block}
where $y'_t=y_t-y_{t-1}$, $x'_{t,i}=x_{t,i}-x_{t-1,i}$ and  $\eta'_t=\eta_t-\eta_{t-1}$.

## Regression with ARIMA errors
\fontsize{13}{15}\sf

Any regression with an ARIMA
error can be rewritten as a regression with an ARMA error by differencing all
variables with the same differencing operator as in the ARIMA model. 

\begin{block}{Original data}\vspace*{-0.2cm}
\begin{align*}
  y_t & = \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t\\
  \mbox{where}\quad
      & \phi(B)(1-B)^d\eta_t = \theta(B)\varepsilon_t
\end{align*}
\end{block} \vspace*{-0.1cm}
\begin{block}{After differencing all variables}\vspace*{-0.2cm}
\begin{align*}
  y'_t & = \beta_1 x'_{1,t} + \dots + \beta_k x'_{k,t} + \eta'_t.\\
  \mbox{where}\quad
       & \phi(B)\eta_t = \theta(B)\varepsilon_t \\
\text{and}\quad & y_t' = (1-B)^dy_t
\end{align*}
\end{block}

## Model selection
\fontsize{13}{15}\sf

  * Fit regression model with automatically selected ARIMA errors. (R will take care of differencing before estimation.)
  * Check that $\varepsilon_t$ series looks like white noise.

 \vspace*{0.5cm}

### Selecting predictors

  * AICc can be calculated for final model.
  * Repeat procedure for all subsets of predictors to be considered, and select model with lowest AICc value.

## US personal consumption and income

```{r usconsump, fig.height=6, fig.width=8.5}
us_change %>%
  gather(key='variable', value='value') %>%
  ggplot(aes(y=value, x=Time, group=variable)) +
  geom_line() + facet_grid(variable ~ ., scales='free_y') +
  xlab("Year") + ylab("") +
  ggtitle("Quarterly changes in US consumption and personal income")
```

## US personal consumption and income

```{r, fig.height=6, fig.width=8.5}
us_change %>%
  ggplot(aes(x=Income,y=Consumption)) +
  geom_point() +
  ggtitle("Quarterly changes in US consumption and personal income")
```

## US personal consumption and income

  * No need for transformations or further differencing.
  *  Increase in income does not necessarily translate into instant increase in consumption (e.g., after the loss of a job, it may take a few months for expenses to be reduced to allow for the new circumstances). We will ignore this for now.

## US personal consumption and income
\fontsize{11}{12}\sf

```{r usconsump2, echo=TRUE, fig.height=3}
fit <- us_change %>% model(ARIMA(Consumption ~ Income))
report(fit)
```

 \begin{alertblock}{}
Write down the equations for the fitted model.
\end{alertblock}

## US personal consumption and income
\fontsize{11}{14}\sf

```{r , echo=TRUE, fig.height=3.7}
residuals(fit, type='regression') %>%
  gg_tsdisplay(.resid) + ggtitle("Regression errors")
```

## US personal consumption and income
\fontsize{11}{14}\sf

```{r , echo=TRUE, fig.height=3.7}
residuals(fit, type='response') %>%
  gg_tsdisplay(.resid) + ggtitle("ARIMA errors")
```

## US personal consumption and income

\fontsize{9}{12}\sf

```{r , echo=TRUE, fig.height=3.7}
augment(fit) %>%
  features(.resid, ljung_box, dof = 5, lag = 12)
```

## US personal consumption and income
\fontsize{9}{12}\sf

```{r usconsump3, echo=TRUE, fig.height=3}
us_change_future <- new_data(us_change, 8) %>%
  mutate(Income = mean(us_change$Income))
forecast(fit, new_data = us_change_future) %>%
  autoplot(us_change) +
  labs(x = "Year", y = "Percentage change",
       title = "Forecasts from regression with ARIMA(1,0,2) errors")
```

## Forecasting

  * To forecast a regression model with ARIMA errors, we need to forecast the
regression part of the model and the ARIMA part of the model and combine the
results.
  * Some predictors are known into the future (e.g., time, dummies).
  * Separate forecasting models may be needed for other predictors.
  * Forecast intervals ignore the uncertainty in forecasting the predictors.

## Daily electricity demand
\fontsize{12}{13}\sf

Model daily electricity demand as a function of temperature using quadratic regression with ARMA errors.

\fontsize{9}{10}\sf

```{r, echo=TRUE, fig.height=3}
vic_elec_daily %>%
  ggplot(aes(x=Temperature, y=Demand, colour=Day_Type)) +
  geom_point() +
  labs(x = "Maximum temperature", y = "Electricity demand (GW)")
```

## Daily electricity demand
\fontsize{11}{12}\sf

```{r, echo=TRUE, fig.height=3.6}
vic_elec_daily %>%
  gather("var", "value", Demand, Temperature) %>%
  ggplot(aes(x = Date, y = value)) + geom_line() +
  facet_grid(vars(var), scales = "free_y")
```

## Daily electricity demand
\fontsize{9}{10}\sf

```{r, echo=TRUE}
fit <- vic_elec_daily %>%
  model(ARIMA(Demand ~ Temperature + I(Temperature^2) +
                (Day_Type=="Weekday")))
report(fit)
```

## Daily electricity demand
\fontsize{9}{10}\sf

```{r, echo=TRUE}
augment(fit) %>%
  gg_tsdisplay(.resid, plot_type = "histogram")
```

## Daily electricity demand
\fontsize{9}{10}\sf

```{r, echo=TRUE}
augment(fit) %>%
  features(.resid, ljung_box, dof = 8, lag = 14)
```

## Daily electricity demand
\fontsize{10}{13}\sf

```{r, echo=TRUE}
# Forecast one day ahead
vic_next_day <- new_data(vic_elec_daily, 1) %>%
  mutate(Temperature = 26, Day_Type = "Holiday")
forecast(fit, vic_next_day)
```

## Daily electricity demand
\fontsize{10}{13}\sf
```{r, echo=TRUE}
vic_elec_future <- new_data(vic_elec_daily, 14) %>%
  mutate(
    Temperature = 26,
    Holiday = c(TRUE, rep(FALSE, 13)),
    Day_Type = case_when(
      Holiday ~ "Holiday",
      wday(Date) %in% 2:6 ~ "Weekday",
      TRUE ~ "Weekend"
    )
  )
```

## Daily electricity demand
\fontsize{10}{13}\sf
```{r, echo = TRUE}
forecast(fit, vic_elec_future) %>%
  autoplot(vic_elec_daily) + ylab("Electricity demand (GW)")
```

# Stochastic and deterministic trends

## Stochastic \& deterministic trends

\begin{block}{Deterministic trend}
\[ y_t = \beta_0 + \beta_1 t + \eta_t \]
where $\eta_t$ is ARMA process. 
\end{block}

\begin{block}{Stochastic trend}
\[ y_t = \beta_0 + \beta_1 t + \eta_t \]
where $\eta_t$ is ARIMA process with $d\ge1$. 

Difference both sides until $\eta_t$ is stationary:
\[ y'_t = \beta_1 + \eta'_t \]
where $\eta'_t$ is ARMA process.
\end{block}

## International visitors

```{r}
aus_visitors <- as_tsibble(fpp2::austa)
aus_visitors %>%
  autoplot(value) +
  labs(x = "Year", y = "millions of people",
       title = "Total annual international visitors to Australia")
```

## International visitors
\fontsize{10}{10}\sf

**Deterministic trend**

```{r, echo=TRUE}
fit_deterministic <- aus_visitors %>%
  model(Deterministic = ARIMA(value ~ trend() + pdq(d = 0)))
report(fit_deterministic)
```
 \vspace*{-0.3cm}

```{r austaparams, echo=FALSE}
coef <- set_names(tidy(fit_deterministic)$estimate, tidy(fit_deterministic)$term)
phi1 <- coef['ar1']
phi2 <- coef['ar2']
intercept <- coef['intercept']
slope <- coef['trend()']
sigma2 <- glance(fit_deterministic)$sigma2
```

## International visitors
\vspace*{-.2cm}
\begin{align*}
  y_t &= `r format(intercept,digits=2)` + `r format(slope, digits=2)` t + \eta_t \\
  \eta_t &= `r format(phi1, digits=3)` \eta_{t-1} `r format(phi2, digits=2)` \eta_{t-2} + \varepsilon_t\\
  \varepsilon_t &\sim \text{NID}(0,`r format(sigma2, digits=3)`).
\end{align*}

## International visitors
\fontsize{10}{10}\sf

**Stochastic trend**

```{r, echo=TRUE}
fit_stochastic <- aus_visitors %>%
  model(Stochastic = ARIMA(value ~ pdq(d=1)))
report(fit_stochastic)
```
 \vspace*{-0.3cm}

```{r austaparams2, cache=TRUE, echo=FALSE}
coef <- set_names(tidy(fit_stochastic)$estimate, tidy(fit_stochastic)$term)
drift <- coef['constant']
theta1 <- coef['ma1']
sigma2 <- glance(fit_stochastic)$sigma2
```

###
\vspace*{-.2cm}
\begin{align*}
y_t-y_{t-1} &= `r format(drift, digits=2)` + \varepsilon_t \\
  y_t &= y_0 + `r format(drift, digits=2)` t + \eta_t \\
  \eta_t &= \eta_{t-1} + `r format(theta1,digits=2, nsmall=2)` \varepsilon_{t-1} + \varepsilon_{t}\\
  \varepsilon_t &\sim \text{NID}(0,`r format(sigma2, digits=3)`).
\end{align*}

## International visitors

```{r}
fc_deterministic <- forecast(fit_deterministic, h = 10)
fc_stochastic <- forecast(fit_stochastic, h = 10)
rbind(fc_deterministic, fc_stochastic) %>%
  autoplot(aus_visitors) +
  facet_grid(vars(.model)) +
  labs(x = "Year", y = "Visitors to Australia (millions)",
       title = "Forecasts from trend models") +
  guides(colour = FALSE)
```

## International visitors

```{r echo=FALSE, fig.height=5.5, fig.width=8.5}
rbind(fc_deterministic, fc_stochastic) %>%
  autoplot(aus_visitors, alpha = 0.5) +
  guides(colour=guide_legend(title="Forecast")) +
  labs(x = "Year", y = "Visitors to Australia (millions)",
       title = "Forecasts from trend models")
```

## Forecasting with trend

  * Point forecasts are almost identical, but prediction intervals differ.
  * Stochastic trends have much wider prediction intervals because the errors are non-stationary.
  * Be careful of forecasting with deterministic trends too far ahead.

# Dynamic harmonic regression

## Dynamic harmonic regression

 **Combine Fourier terms with ARIMA errors**

\fontsize{13}{14}\sf

### Advantages
   * it allows any length seasonality;
   * for data with more than one seasonal period, you can include Fourier terms of different frequencies;
   * the seasonal pattern is smooth for small values of $K$ (but more wiggly seasonality can be handled by increasing $K$);
   * the short-term dynamics are easily handled with a simple ARMA error.

### Disadvantages
 * seasonality is assumed to be fixed

## Eating-out expenditure

\fontsize{10}{11}\sf

```{r cafe, echo=TRUE, fig.height=3.8, fig.width=8}
aus_cafe <- aus_retail %>% filter(
    Industry == "Cafes, restaurants and takeaway food services",
    year(Month) %in% 2004:2018
  ) %>% summarise(Turnover = sum(Turnover))
aus_cafe %>% autoplot(Turnover)
```

## Eating-out expenditure

\fontsize{9}{11}\sf

```{r cafefit, dependson='cafe', fig.height=5, echo=TRUE, results='hide'}
fit <- aus_cafe %>% model(
    `K = 1` = ARIMA(log(Turnover) ~ fourier(K = 1) + PDQ(0,0,0)),
    `K = 2` = ARIMA(log(Turnover) ~ fourier(K = 2) + PDQ(0,0,0)),
    `K = 3` = ARIMA(log(Turnover) ~ fourier(K = 3) + PDQ(0,0,0)),
    `K = 4` = ARIMA(log(Turnover) ~ fourier(K = 4) + PDQ(0,0,0)),
    `K = 5` = ARIMA(log(Turnover) ~ fourier(K = 5) + PDQ(0,0,0)),
    `K = 6` = ARIMA(log(Turnover) ~ fourier(K = 6) + PDQ(0,0,0)))
glance(fit)
```
```{r, echo = FALSE}
glance(fit) %>%select(.model, sigma2, AIC, AICc, BIC) %>% knitr::kable()
```

## Eating-out expenditure

```{r}
cafe_plot <- function(...){
  fit %>%
    select(...) %>%
    forecast() %>% autoplot(aus_cafe) +
    ggtitle(sprintf("Log transformed %s, fourier(K = %s)", model_sum(select(fit,...)[[1]][[1]]), deparse(..1))) +
    geom_label(
      aes(x = yearmonth("2007 Jan"), y = 4250, label = paste0("AICc = ", format(AICc))),
      data = glance(select(fit,...))
    ) +
    geom_line(aes(y = .fitted), colour = "red", augment(select(fit, ...))) +
    ylim(c(1500, 5100))
}
```

```{r cafe1, dependson='cafe', fig.height=5, echo=FALSE}
cafe_plot(K = 1)
```

## Eating-out expenditure

```{r cafe2, dependson='cafe', fig.height=5, echo=FALSE}
cafe_plot(K = 2)
```

## Eating-out expenditure

```{r cafe3, dependson='cafe', fig.height=5, echo=FALSE}
cafe_plot(K = 3)
```

## Eating-out expenditure

```{r cafe4, dependson='cafe', fig.height=5, echo=FALSE}
cafe_plot(K = 4)
```

## Eating-out expenditure

```{r cafe5, dependson='cafe', fig.height=5, echo=FALSE}
cafe_plot(K = 5)
```

## Eating-out expenditure

```{r cafe6, dependson='cafe', fig.height=5, echo=FALSE}
cafe_plot(K = 6)
```

## Example: weekly gasoline products
\fontsize{8}{8}\sf
```{r, echo = FALSE}
options(width = 70)
```

```{r gasmodel, echo=TRUE}
gasoline <- as_tsibble(fpp2::gasoline)
fit <- gasoline %>% model(ARIMA(value ~ fourier(K = 13) + PDQ(0,0,0)))
report(fit)
```

## Example: weekly gasoline products
\fontsize{11}{12}\sf

```{r gasf, echo=TRUE, fig.height=3.5}
forecast(fit, h = "3 years") %>%
  autoplot(gasoline)
```

## 5-minute call centre volume

\fontsize{9}{10}\sf

```{r calls, echo=TRUE, fig.height=4}
(calls <- read_tsv("http://robjhyndman.com/data/callcenter.txt") %>%
  gather("date", "volume", -X1) %>% transmute(
    time = X1, date = as.Date(date, format = "%d/%m/%Y"),
    datetime = as_datetime(date) + time, volume) %>%
  as_tsibble(index = datetime))
```

## 5-minute call centre volume

```{r calls-plot, echo=TRUE, fig.height=4}
calls %>% fill_gaps() %>% autoplot(volume)
```

## 5-minute call centre volume

\fontsize{12}{13}\sf

```{r calls-season, echo=TRUE, fig.height=3}
calls %>% fill_gaps() %>%
  gg_season(volume, period = "day", alpha = 0.1) +
  guides(colour = FALSE)
```

## 5-minute call centre volume

\fontsize{10}{11}\sf
```{r calls-cal, echo=TRUE, fig.height=3}
library(sugrrants)
calls %>% filter(month(date, label = TRUE) == "Apr") %>%
  ggplot(aes(x = time, y = volume)) +
  geom_line() + facet_calendar(date)
```

## 5-minute call centre volume
\fontsize{8}{8}\sf

```{r callsmodel, echo=TRUE}
calls_mdl <- calls %>%
  mutate(idx = row_number()) %>%
  update_tsibble(index = idx)
fit <- calls_mdl %>%
  model(ARIMA(volume ~ fourier(169, K = 10) + pdq(d=0) + PDQ(0,0,0)))
report(fit)
```

## 5-minute call centre volume
\fontsize{10}{13}\sf

```{r callsres, echo=TRUE}
augment(fit) %>%
  gg_tsdisplay(.resid, plot_type = "histogram", lag_max = 338)
```

## 5-minute call centre volume
\fontsize{11}{12}\sf

```{r callsf, echo=TRUE, fig.height=4}
fit %>% forecast(h = 1690) %>%
  autoplot(calls_mdl)
```

# Lagged predictors

## Lagged predictors

\structure{Sometimes a change in $x_t$ does not affect $y_t$ instantaneously} 
\begin{block}{}
\begin{itemize}
  \item $y_t=$ sales, $x_t=$ advertising.
  \item $y_t=$ stream flow, $x_t=$ rainfall.
  \item $y_t=$ size of herd, $x_t=$ breeding stock.
\end{itemize}
\end{block}
 

  * These are dynamic systems with input ($x_t$) and output $(y_t)$.
  * $x_t$ is often a leading indicator.
  * There can be multiple predictors.

## Lagged predictors

The model include present and past values of predictor: $x_t,x_{t-1},x_{t-2},\dots.$
\begin{block}{}
\centerline{$
y_t = a + \nu_0x_t + \nu_1x_{t-1} + \dots + \nu_kx_{t-k} + \eta_t$}
\end{block}
where $\eta_t$ is an ARIMA process. 

## Lagged predictors

**Rewrite model as **
\begin{align*}
y_{t} & = a+ (\nu_{0} + \nu_{1} B + \nu_{2} B^{2} + \dots + \nu_{k} B^{k}) x_{t} +\eta_t \\
      & = a+ \nu(B) x_{t} +\eta_t.
\end{align*} \vspace*{-0.6cm}

  * $\nu(B)$ is called a \textit{transfer function} since it describes how
change in $x_t$ is transferred to $y_t$.
  * $x$ can influence $y$, but $y$ is not allowed to influence $x$.

## Example: Insurance quotes and TV adverts

```{r tvadvert}
insurance <- as_tsibble(fpp2::insurance, pivot_longer = FALSE)

insurance %>% gather("key", "value", Quotes, TV.advert) %>%
  ggplot(aes(x = index, y = value)) + geom_line() +
  facet_grid(vars(key), scales = "free_y") +
  labs(x = "Year", y = NULL, title = "Insurance advertising and quotations")
```

## Example: Insurance quotes and TV adverts
\fontsize{10}{10}\sf

```{r, echo=TRUE}
fit <- insurance %>%
  # Restrict data so models use same fitting period
  mutate(Quotes = c(NA,NA,NA,Quotes[4:40])) %>%
  # Estimate models
  model(
    ARIMA(Quotes ~ pdq(d = 0) + TV.advert),
    ARIMA(Quotes ~ pdq(d = 0) + TV.advert + lag(TV.advert)),
    ARIMA(Quotes ~ pdq(d = 0) + TV.advert + lag(TV.advert) +
            lag(TV.advert, 2)),
    ARIMA(Quotes ~ pdq(d = 0) + TV.advert + lag(TV.advert) +
            lag(TV.advert, 2) + lag(TV.advert, 3))
  )
```

## Example: Insurance quotes and TV adverts
\fontsize{10}{10}\sf

```{r, echo=TRUE, results = 'hide'}
glance(fit)
```
```{r, echo = FALSE}
glance(fit) %>%
  transmute(`Lag order` = 0:3, sigma2, AIC, AICc, BIC) %>%
  knitr::kable()
```

## Example: Insurance quotes and TV adverts

\fontsize{8}{9}\sf

```{r tvadvertagain, echo=TRUE}
fit <- insurance %>%
  model(ARIMA(Quotes ~ pdq(3, 0, 0) + TV.advert + lag(TV.advert)))
report(fit)
```

 
\vspace{-1em}

```{r tvadvertparam, echo=FALSE}
# Store coefficients
coef <- set_names(tidy(fit)$estimate, tidy(fit)$term)
phi1 <- coef['ar1']
phi2 <- coef['ar2']
phi3 <- coef['ar3']
intercept <- coef['intercept']
gamma0 <- coef['TV.advert']
gamma1 <- coef['lag(TV.advert)']
```

###
\begin{align*}
  y_t &= `r format(intercept, digits=3)` +
         `r format(gamma0, digits=3)` x_t +
         `r format(gamma1, digits=2)` x_{t-1} + \eta_t,\\
  \eta_t &= `r format(phi1, digits=3)` \eta_{t-1}
        `r format(phi2, digits=2)` \eta_{t-2} +
        `r format(phi3, digits=2)` \eta_{t-3} + \varepsilon_t,
\end{align*}

## Example: Insurance quotes and TV adverts
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3}
advert_a <- new_data(insurance, 20) %>%
  mutate(TV.advert = 10)
forecast(fit, advert_a) %>% autoplot(insurance)
```

## Example: Insurance quotes and TV adverts
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3}
advert_b <- new_data(insurance, 20) %>%
  mutate(TV.advert = 8)
forecast(fit, advert_b) %>% autoplot(insurance)
```

## Example: Insurance quotes and TV adverts
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3}
advert_c <- new_data(insurance, 20) %>%
  mutate(TV.advert = 6)
forecast(fit, advert_c) %>% autoplot(insurance)
```
