---
title: 'Lab 5: ARIMA models and decomposition methods'
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. The `fma::plastics` data set consists of the monthly sales (in thousands) of product A for a plastics manufacturer for five years.
 * Plot the time series of sales of product A. Can you identify seasonal fluctuations and/or a trend-cycle?
 * Use a classical multiplicative decomposition to calculate the trend-cycle and seasonal indices.
 * Do the results support the graphical interpretation from part a?
 * Compute and plot the seasonally adjusted data.
 * Change one observation to be an outlier (e.g., add 500 to one observation), and recompute the seasonally adjusted data. What is the effect of the outlier?
 * Does it make any difference if the outlier is near the end rather than in the middle of the time series?
 
2. Monthly Australian retail data is provided in `aus_retail`. Select one of the time series as follows (but choose your own seed value). Decompose the series using X11. Does it reveal any outliers, or unusual features. 

```{r, echo=TRUE}
set.seed(12345678)
myseries <- aus_retail %>%
 filter(`Series ID` == sample(aus_retail$`Series ID`,1))

myseries %>% model(feasts:::X11(Turnover,type="additive"))%>% 
  components()%>%autoplot() +
  ggtitle("X11 decomposition of myseries")
```

3. For the canadian_gas data (monthly Canadian gas production in billions of cubic metres, January 1960 – February 2005),

 * Plot the data using `autoplot()`, `gg_subseries()` and `gg_season()` to look at the effect of the changing seasonality over time. What do you think is causing it to change so much?
 
 * Do an STL decomposition of the data. You will need to choose a seasonal window to allow for the changing shape of the seasonal component.
 
 * How does the seasonal shape change over time?
 
4. Figure 1 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.
```{r pressure, echo=FALSE, fig.cap="ACF plots", out.width = '80%'}
knitr::include_graphics("lab6-acf.png")
```
 * Explain the differences among these figures. Do they all indicate that the data are white noise?
 * Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?


5. For the `expsmooth::enplanements` dataset, find an appropriate differencing (after transformation if necessary) to obtain stationary data. Write down the differences you chose above using backshift operator notation. 

6. Use R to simulate and plot some data from simple ARIMA models. As an example, the following R code generates data from an AR(1) model with $\phi_1 = 0.6$ and $\sigma^2 = 1$. The process starts with $y_1=0$.
 
```{r, echo=TRUE}
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
```

 * Produce a time plot for the above series. How does the plot change as you change $\phi_1$?

 * Write your own code to generate data from an MA(1) model with $\theta_1 = 0.6$ and $\sigma^2 = 1$. You may start with $\epsilon_0 = 0$. Produce a time plot for the series. How does the plot change as you change $\theta_1$?
 * Graduate students only: Generate data from an ARMA(1,1) model with $\phi_1 = 0.6, \theta_1 = 0.6$ and $\sigma^2 = 1$.
 * Graduate students only: Generate data from an AR(2) model with $\phi_1 = -0.8, \phi_2 = 0.3$ and $\sigma^2 = 1$. Note that these parameters will give a non-stationary series. For the AR(2) model, the characteristic function is $\phi(B) = (1+0.8B-0.2B^2)$. You may find the mode of the complex roots for $\phi(B)=0$ using the following codes. Are the modulus of all the complex roots greater than 1?
 
```{r, echo=TRUE}
# polyroot for finding complex roots
# Mod for finding the modulus of the number
Mod(polyroot (c(1,0.8,-0.3))) 
```
 * Graduate students only: Graph the latter two series and compare them.

7. Undergraduate students only: for an MA(1), prove that the autocorrelation function is $$\rho_1 = \frac{\theta_1}{1+\theta_1^2}$$ and $\rho_h = 0$ for $h \geq 2$. 

8. Graduate students only: for an MA(q), prove that the autocorrelation function is $$\rho_h = \frac{\sum_{j=0}^{q-h}\theta_j\theta_{j+h}}{1+\theta_1^2+\dots+\theta_q^2} $$ for $1\leq h \leq q$ and $\rho_h = 0$ for $h>q$.
