---
title: 'Lab 3: Simple forecasting tools and linear regression'
author: "Introduction to Time Series Analysis"
output:
  pdf_document: default
  word_document: default
  html_document:
    df_print: paged
header-includes:
- \usepackage{amsmath}
---

\newcommand{\ba}{{\mathbf a}}

Name:  
This lab is to be done in class (completed outside of class if need be). You can collaborate with your classmates, but you must identify their names above, and you must submit **your own** lab as an knitted pdf file. To answer the questions, display the results and write your answers if asked.


```{r, include= FALSE}
knitr::opts_chunk$set(cache=TRUE, autodep=TRUE, cache.comments=TRUE)
library(tidyverse)
library(fable)
library(tsibble)
library(feasts)
library(lubridate)
library(tsibbledata)
```


1.  For each of the following time series, calculate the residuals from a naive or seasonal naïve forecast. Test if the residuals are white noise and normally distributed. What do you conclude?

  * Quarterly Australian beer production data from 1992
  * Australian Exports series from `global_economy `
  * Bricks series from `aus_production`

```{r}
# Quarterly Australian beer production data from 1992
recent_production <- aus_production %>% filter(year(Quarter) >= 1992)
fit <- recent_production %>% model(SNAIVE(Beer))
fit %>% gg_tsresiduals()
fit %>% forecast() %>% autoplot(recent_production)
```


2.  For the turnover time series in the following dataset, conduct the following tasks:

```{r}
   set.seed(12345678)
   myts <- aus_retail %>%filter(`Series ID` == sample(aus_retail$`Series ID`,1))
   
   myts_train <- myts %>%filter(Month <= yearmonth("2010 Dec"))
   myts%>%autoplot(Turnover) +autolayer(myts_train,Turnover, colour = "red")
   fit <- myts_train %>% model(SNAIVE(Turnover))
   fc <- fit %>% forecast()
   fit %>% accuracy()
   fc %>% accuracy(myts)
   fit %>% gg_tsresiduals()
```

  * Create a training dataset consisting of observations before 2011
  * Calculate seasonal naïve forecasts using `SNAIVE()` applied to your training data 
  * Compare the accuracy (e.g. using RMSE, MAE, MASE) of your forecasts against the actual values.
  * Check the residuals. Do the residuals appear to be uncorrelated and normally distributed?
  * How sensitive are the accuracy measures to the amount of training data used?
  
  
3. Consider the number of pigs slaughted in New South Wales (data set `aus_livestock`).

   * Produce some plots of the data in order to become familiar with it.
   * Create a training set of 486 observations, witholding a test set of 72 observations (6 years).
   * Try using various simple methods to forecast the training set and compare the results on the test set. Which method did best?
   * Check the residuals of your preferred method. Do they resemble white noise?

4. Half-hourly electricity demand for Victoria, Australia is contained in vic_elec. Extract the January 2014 electricity demand, and aggregate this data to daily with daily total demands and maximum temperatures.

```{r}
 jan14_vic_elec <- vic_elec %>%
  filter(yearmonth(Time) == yearmonth("2014 Jan")) %>%
  index_by(Date = as_date(Time)) %>%
  summarise(Demand = sum(Demand), Temperature = max(Temperature))
```

   * Plot the data and find the regression model for Demand with temperature as an explanatory variable. Why is there a positive relationship?
   * Produce a residual plot. Is the model adequate? Are there any outliers or influential observations?
   * Use the model to forecast the electricity demand that you would expect for the next day if the maximum temperature was $15^o$ and compare it with the forecast if the with maximum temperature was $35^o$.
   * Comment on these forecasts? Give prediction intervals for your forecasts. 
   
```{r}
jan14_vic_elec %>%
  gather("Variable", "Value", Demand, Temperature) %>%
  ggplot(aes(x = Date, y = Value)) +
  geom_line() +
  facet_grid(vars(Variable), scales = "free_y")
jan14_vic_elec %>%
  ggplot(aes(x=Temperature, y=Demand)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE)
#> `geom_smooth()` using formula 'y ~ x'
fit <- jan14_vic_elec %>% model(TSLM(Demand ~ Temperature))
fit %>% gg_tsresiduals()
fit %>% forecast(new_data(jan14_vic_elec, 1) %>% mutate(Temperature = 15)) %>%
  mutate(interval = hilo(.distribution, 80))
```
5. The `us_gasoline` series from package `fpp3` consists of weekly data for supplies of US finished motor gasoline product, from 2 February 1991 to 20 January 2017. The units are in ``million barrels per day''. Consider only the data to the end of 2004. 

   * Fit a harmonic regression with trend to the data. Experiment with changing the number Fourier terms. Plot the observed gasoline and fitted values and comment on what you see.
   * Select the appropriate number of Fourier terms to include by minimising the AICc or CV value.
   * Check the residuals of the final model using the gg_tsresiduals() function. Use a Ljung-Box test to check for correlation in the residuals. Even though the residuals fail the correlation tests, the results are probably not severe enough to make much difference to the forecasts and prediction intervals. (Note that the correlations are relatively small, even though they are significant.)
   * Forecast the next year of data.
   * Plot the forecasts along with the actual data for 2005. What do you find?
   
6. Graduate students only:   
```{r pressure, echo=FALSE, out.width = '80%'}
knitr::include_graphics("lmderiv.png")
```

   