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

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache=TRUE)
library(tidyverse)
library(fable)
library(tsibble)
library(feasts)
library(lubridate)
library(tsibbledata)
options(width=50)
```

# Time series in R


## Class packages
\fontsize{13}{15}\sf

```r
# Data manipulation and plotting functions
library(tidyverse)
# Forecasting functions
library(fable)
# Time series manipulation
library(tsibble)
# Time series graphics and statistics
library(feasts)
# Functions to work with date-times
library(lubridate)
# Tidy time series data
library(tsibbledata)
```

## `tsibbledata` datasets
\small

1. `anasett`: Passenger numbers on Anasett airline flights
2. `aus_livestock`: Meat production in Australia for human consumption from Q3 1965 to Q4 2018.
3. `aus_production`: Quarterly estimates of manufacturing production of selected commodities in Australia.
4. `aus_retail`: Australian retail trade turnover (total value of retail traded).
5. `gafa_stock`: GAFA stock prices.
6. `global_economy`: Global economic indicators.

## `tsibbledata` datasets
\small

7. `hh_budget`: Household budget characteristics.
8. `nyc_bikes `: A sample from NYC Citi Bike usage of 10 bikes throughout 2018.
9. `olympic_running`: Fastest running times for Olympic races.
10. `PBS`: Monthly Medicare Australia prescription data.
11. `pelt`: Pelt trading records.
12. `vic_elec`: Half-hourly electricity demand for Victoria, Australia


## `tsibble` objects

\fontsize{12}{13}\sf

A `tsibble` allows storage and manipulation of time series in R. 
A tsibble is a data- and model-oriented object.

It contains:

- Measured variable(s): numbers of interest
- Key variable(s): identifiers for each series
- An index: time information about the observation
- A tsibble is sorted by its key first and then index
- Key variable(s) together with the index uniquely identifies each record


## `tsibble` objects example
\fontsize{13}{15}\sf

```{r tstable, cache=TRUE}
#Creating a tsibble object
library(tsibble)
e1 <- tsibble(year = 2012:2016, x = c(1,2,1,2,2),
  y = c(123,39,78,52,110), index = year, key = x)
e1
```

## `tsibble` objects example
\fontsize{13}{15}\sf
```{r tstable2, cache=TRUE}
# yearquarter returns a numeric value that can be added
e2<- tsibble(
  qtr = rep(yearquarter("2010 Q1") + 0:9, 3),
  group = rep(c("x", "y", "z"), each = 10),
  value = rnorm(30),key = group, index = qtr)
e2
```



## The `tsibble` index

Common time index variables can be created with these functions:

```{r tstable4, echo=FALSE}
tribble(
  ~ `Frequency`, ~ Function,
  "Annual", "start year:end year",
  "Quarterly", "yearquarter()",
  "Monthly", "yearmonth()",
  "Weekly", "yearweek()",
  "Daily", "as_Date(), ymd()",
  "Sub-daily", "as_datetime()"
) %>%
  knitr::kable(booktabs=TRUE)
```

## Example for time index variables
\fontsize{13}{15}\sf
```{r tstable5, cache=TRUE}
2015:2020
yearquarter("2010 Q1")+0:3
yearmonth("2010 1")+0:3

```

## Example for time index variables
\fontsize{13}{15}\sf
```{r tstable6, cache=TRUE}
yearweek("2010 1")+0:3
as.Date("2020-01-22") + 0:3
```

## Example for time index variables
\fontsize{13}{15}\sf
```{r tstable7, cache=TRUE}
ymd("2020-01-22")+0:3
as_datetime("2020-01-22 00:50:50")+0:3
```


## Coerce a dataset to be an tsibbledata object

```{r key-tsibble}
olympic_running %>% as_tsibble(
  key = c(Length, Sex), index = Year)
```

## The `key` to many time series

```{r tstablekey1, eval=FALSE}
tsibbledata::olympic_running %>%
  group_by_key() %>%
  slice(1) %>%
  head(6) %>%
  knitr::kable(booktabs=TRUE)
```

## The `key` to many time series

```{r tstablekey2, echo=FALSE}
tsibbledata::olympic_running %>%
  group_by_key() %>%
  slice(1) %>%
  head(6) %>%
  knitr::kable(booktabs=TRUE)
```

## Australian GDP

\fontsize{12}{14}\sf

```{r, echo = TRUE}
#filter is to select a subset in rows
aus_economy <- global_economy %>%
  filter(Code == "AUS")
```

<!--   * Object: "tsibble" -->
<!--   * Print and plotting methods available. -->

```{r, echo=FALSE}
aus_economy
```


# Time plots

## Australian GDP

```{r, echo=TRUE, fig.height=3}
aus_economy %>% autoplot(GDP)
```

## Time plots

\small

```{r, echo=TRUE, fig.height=3}
ansett %>%
  filter(Airports=="MEL-SYD", Class=="Economy") %>%
  autoplot(Passengers)
```

## Time plots

\small
```{r a11, echo = TRUE}
# Taking a subset of the time series according to time
a10.subset <- PBS %>% filter(ATC2 == "A10")%>%
  filter_index("2008 Jan")
a10.subset
```

## Time plots
```{r a10, echo = TRUE}
a10 <- PBS %>%
  filter(ATC2 == "A10") %>%
  summarise(Cost = sum(Cost)/1e6)
```


## Time plots

\small

```{r a10-plot, echo=TRUE, dependson="a10", fig.height=3}
a10 %>% autoplot(Cost) +
  ylab("$ million") + xlab("Year") +
  ggtitle("Antidiabetic drug sales")
```

## Your turn

- Create plots of the following time series: Bricks from `aus_production`, Lynx from `pelt`, Google from `gafa_stock`
- Use `help()` to find out about the data in each series.
- For the last plot, modify the axis labels and title.

## Are time plots the best?

\fontsize{12}{14}\sf

```{r maxtemp, echo=FALSE,fig.height=3}
maxtemp <- vic_elec %>%
  index_by(Day = lubridate::date(Time)) %>%
  summarise(
    Temperature = max(Temperature),
    Demand = sum(Demand),
    Holiday = any(Holiday)
  )
```

```{r maxtemp1, warning=FALSE, message=FALSE, dependson="maxtemp", fig.height=3}
maxtemp %>%
  autoplot(Temperature) +
  xlab("Week") + ylab("Max temperature")
```

## Are time plots the best?

\fontsize{12}{14}\sf

```{r maxtemp2, warning=FALSE, message=FALSE, dependson="maxtemp", fig.height=3}
maxtemp %>%
  ggplot(aes(x = Day, y = Temperature)) + geom_point() +
  xlab("Week") + ylab("Max temperature")
```

## Are time plots the best?

```{r maxtemp3, warning=FALSE, message=FALSE, echo=FALSE, dependson="maxtemp", fig.height=3}
maxtemp %>%
  ggplot(aes(x = Day, y = 1)) +
  geom_tile(aes(fill = Temperature)) +
  scale_fill_gradient2(
    low = "navy",
    mid = "yellow",
   high = "red", midpoint=28) +
  ylab("") + scale_y_discrete(expand=c(0,0))
```

## Are time plots the best?

\fullheight{TemperatureBlanket}

# Seasonal plots

## Seasonal plots

\footnotesize

```{r, echo=TRUE, dependson="a10"}
a10 %>% gg_season(Cost, labels = "both") +
  ylab("$ million") + ggtitle("Seasonal plot: antidiabetic drug sales")
```

## Seasonal plots

  * Data plotted against the individual "seasons" in which the data were observed.  (In this case a "season" is a month.)
  * Something like a time plot except that the data from each season are overlapped.
  * Enables the underlying seasonal pattern to be seen more clearly, and also allows any substantial departures from the seasonal pattern to be easily identified.
  * In R: `gg_season()`

<!-- ## Seasonal polar plots -->

<!-- \small -->

<!-- ```{r, fig.height=6, out.width="7.5cm"} -->
<!-- a10 %>% gg_season(polar = TRUE) + -->
<!--   ylab("$ million") -->
<!-- ``` -->

## Seasonal subseries plots

\small

```{r, echo=TRUE, dependson="a10",fig.height=3}
a10 %>%
  gg_subseries(Cost) + ylab("$ million") +
  ggtitle("Subseries plot: antidiabetic drug sales")
```

## Seasonal subseries plots

  * Data for each season collected together in time plot as separate time series.
  * Enables the underlying seasonal pattern to be seen clearly, and changes in seasonality over time to be visualized.
  * In R: `gg_subseries()`

## Quarterly Australian Beer Production
\fontsize{13}{15}\sf

```{r, fig.height=3}
beer <- aus_production %>%
  select(Quarter, Beer) %>%
  filter(year(Quarter) >= 1992)
beer %>% autoplot(Beer)
```

## Quarterly Australian Beer Production

```{r}
beer %>% gg_season(Beer, labels="right")
```

## Quarterly Australian Beer Production

```{r}
beer %>% gg_subseries(Beer)
```

## Your turn

Look at the quarterly tourism data for the Snowy Mountains

```r
snowy <- filter(tourism,
  Region == "Snowy Mountains",
  Purpose == "Holiday")
```

  - Use `autoplot()`, `gg_season()` and `gg_subseries()` to explore the data.
  - What do you learn?

# Seasonal or cyclic?

## Time series patterns

Trend
  : pattern exists when there is a long-term increase or decrease in the data.

Seasonal
  : pattern exists when a series is influenced by seasonal factors (e.g., the quarter of the year, the month, or day of the week).

Cyclic
  : pattern exists when data exhibit rises and falls that are \emph{not of fixed period}. 

## Time series components

### Differences between seasonal and cyclic patterns:

* seasonal pattern constant length; cyclic pattern variable length
* average length of cycle longer than length of seasonal pattern
* magnitude of cycle more variable than magnitude of seasonal pattern

## Time series patterns
\fontsize{13}{15}\sf

```{r, fig.height=3}
as_tsibble(fma::elec) %>%
  filter(index >= 1980) %>% # or filter_index("1980 Jan" ~.)
  autoplot(value) + xlab("Year") + ylab("GWh") +
  ggtitle("Australian electricity production")
```

## Time series patterns
```{r, warning=FALSE, fig.height=3}
as_tsibble(fma::elec) %>%
  filter(index >= 1980) %>% gg_subseries(value)+ 
  xlab("Year") + ylab("GWh") +ggtitle("Australian electricity production")
```

## Time series patterns
```{r, warning=FALSE, fig.height=3}
as_tsibble(fma::elec) %>%
  filter(index >= 1980) %>% gg_season(value)+ 
  xlab("Year") + ylab("GWh") +ggtitle("Australian electricity production")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, warning=FALSE, fig.height=3}
aus_production %>%
  autoplot(Bricks) +
  ggtitle("Australian clay brick production") +
  xlab("Year") + ylab("million units")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, warning=FALSE, fig.height=3}
aus_production %>%
  gg_subseries(Bricks) +
  ggtitle("Australian clay brick production") +
  xlab("Year") + ylab("million units")
```


## Time series patterns
\fontsize{13}{15}\sf

```{r, warning=FALSE, fig.height=3}
aus_production %>%
  gg_season(Bricks) +
  ggtitle("Australian clay brick production") +
  xlab("Year") + ylab("million units")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, fig.height=3}
as_tsibble(fma::hsales) %>%
  autoplot(value) +
  ggtitle("Sales of new one-family houses, USA") +
  xlab("Year") + ylab("Total sales")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, fig.height=3}
as_tsibble(fma::hsales) %>%
  gg_subseries(value) +
  ggtitle("Sales of new one-family houses, USA") +
  xlab("Year") + ylab("Total sales")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, fig.height=3}
as_tsibble(fma::hsales) %>%
  gg_season(value) +
  ggtitle("Sales of new one-family houses, USA") +
  xlab("Year") + ylab("Total sales")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, fig.height=3}
as_tsibble(fma::ustreas) %>%
  autoplot(value) +
  ggtitle("US Treasury Bill Contracts") +
  xlab("Day") + ylab("price")
```

## Time series patterns
\fontsize{13}{15}\sf

```{r, fig.height=3}
pelt %>%
  autoplot(Lynx) +
  ggtitle("Annual Canadian Lynx Trappings") +
  xlab("Year") + ylab("Number trapped")
```


## Seasonal or cyclic?

\alert{Differences between seasonal and cyclic patterns:}

  * seasonal pattern constant length; cyclic pattern variable length
  * average length of cycle longer than length of seasonal pattern
  * magnitude of cycle more variable than magnitude of seasonal pattern

\pause

\begin{alertblock}{}
The timing of peaks and troughs is predictable with seasonal data, but unpredictable in the long term with cyclic data.
\end{alertblock}

# Lag plots and autocorrelation

## Example: Beer production
\fontsize{13}{15}\sf

```{r}
new_production <- aus_production %>%
  filter(year(Quarter) >= 1992) # or filter_index("1992 Q1"~.)
new_production
```


## Lagged scatterplots

 - Each graph shows $y_t$ plotted against $y_{t-k}$ for
different values of $k$.
 - Vertical axis: lagged observation
 - Horizontal axis: current observation
 - Colors: points are colored by the current quarter


## Example: Beer production
\fontsize{13}{15}\sf


```{r, fig.height=6.5, fig.width=6.5, out.width="8cm"}
new_production %>% gg_lag(Beer, geom='point')
```



## Lagged scatterplots

  * The autocorrelations are the correlations associated
with these scatterplots.

## Autocorrelation

**Covariance** and **correlation**: measure extent of **linear relationship** between two variables ($y$ and $X$).\pause

**Autocovariance** and **autocorrelation**: measure linear relationship between **lagged values** of a time series $y$.\pause

We measure the relationship between:

  * $y_{t}$ and $y_{t-1}$
  * $y_{t}$ and $y_{t-2}$
  * $y_{t}$ and $y_{t-3}$
  * etc.

## Autocorrelation

We denote the sample autocovariance at lag $k$ by $c_k$ and the sample autocorrelation at lag $k$ by $r_k$.  Then define

\begin{block}{}
\begin{align*}
c_k &= \frac{1}{T}\sum_{t=k+1}^T (y_t-\bar{y})(y_{t-k}-\bar{y}) \\[0.cm]
\text{and}\qquad
r_{k} &= c_k/c_0
\end{align*}
\end{block}\pause\small

  * $r_1$ indicates how successive values of  $y$  relate to each other
  * $r_2$ indicates how  $y$ values two periods apart relate to each other
  * $r_k$ is \textit{almost} the same as the sample correlation between $y_t$ and $y_{t-k}$.

## Autocorrelation


Results for first 9 lags for beer data:

\fontsize{11}{13}\sf

```{r, echo=TRUE}
new_production %>% ACF(Beer, lag_max = 9)
```

## Autocorrelation

Results for first 9 lags for beer data:

\fontsize{11}{13}\sf

```{r beeracf, fig.height=2.5}
new_production %>% ACF(Beer, lag_max = 9) %>% autoplot()
```

\vspace*{10cm}

## Autocorrelation

  * $r_{4}$  higher than for the other lags. This is due to **the seasonal pattern in the data**: the peaks tend to be **4 quarters** apart and the troughs tend to be **2 quarters** apart.
  * $r_2$ is more negative than for the other lags because troughs tend to be 2 quarters behind peaks.
  * Together, the autocorrelations at lags 1, 2, \dots, make up the \emph{autocorrelation} or ACF.
  * The plot is known as a **correlogram**

## ACF

```{r, fig.height=3, echo=TRUE}
new_production %>% ACF(Beer) %>% autoplot()
```

## Trend and seasonality in ACF plots

- When data have a trend, the autocorrelations for small lags tend to be large and positive.
- When data are seasonal, the autocorrelations will be larger at the seasonal lags (i.e., at multiples of the seasonal frequency)
- When data are trended and seasonal, you see a combination of these effects.

## Aus monthly electricity production

```{r,  fig.height=3}
elec2 <- as_tsibble(fma::elec) %>%
  filter(year(index) >= 1980)
elec2 %>% autoplot(value)
```

## Aus monthly electricity production

```{r, fig.height=3, warning=FALSE}
elec2 %>% ACF(value, lag_max=48) %>%
  autoplot()
```

## Aus monthly electricity production

Time plot shows clear trend and seasonality.

The same features are reflected in the ACF.

  * The slowly decaying ACF indicates trend.
  * The ACF peaks at lags 12, 24, 36, \dots, indicate seasonality of length 12.

## Google stock price
\fontsize{12}{14}\sf

```{r google-2015, warning=FALSE}
google_2015 <- gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) == 2015) %>%
  select(Date, Close)
google_2015
```

## Google stock price

```{r}
google_2015 %>% autoplot(Close)
```

## Google stock price
\fontsize{12}{16}\sf

```r
google_2015 %>%
  ACF(Close, lag_max=100)
# Error: Can't handle tsibble of irregular interval.
```

\pause

```{r}
google_2015
```

## Google stock price
\fontsize{12}{14}\sf

```{r}
#mutate is to create a new variable
google_2015 <- google_2015 %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index=trading_day, regular=TRUE)
google_2015
```

## Google stock price

```{r, fig.height=3,warning=FALSE}
google_2015 %>%ACF(Close, lag_max=100) %>% autoplot()
```

## Your turn

We have introduced the following functions:

  - `gg_lag`
  - `ACF`

Explore the following time series using these functions. Can you spot any seasonality, cyclicity and trend? What do you learn about the series?

  - Bricks from `aus_production`
  - Lynx from `pelt`
  - Victorian Electricity Demand from `aus_elec`

## Which is which?

```{r, fig.height=6, fig.width=12, echo=FALSE, warning=FALSE, out.width="11.5cm"}
cowtemp <- as_tsibble(fma::cowtemp)
USAccDeaths <- as_tsibble(USAccDeaths)
AirPassengers <- as_tsibble(AirPassengers)
mink <- as_tsibble(fma::mink)
tp1 <- autoplot(cowtemp, value) + xlab("") + ylab("chirps per minute") +
  ggtitle("1. Daily temperature of cow")
tp2 <- autoplot(USAccDeaths, value) + xlab("") + ylab("thousands") +
  ggtitle("2. Monthly accidental deaths")
tp3 <- autoplot(AirPassengers, value) + xlab("") + ylab("thousands") +
  ggtitle("3. Monthly air passengers")
tp4 <- autoplot(mink, value) + xlab("") + ylab("thousands") +
  ggtitle("4. Annual mink trappings")
acfb <- ACF(cowtemp, value) %>% autoplot() + xlab("") + ggtitle("B") + ylim(-0.4,1)
acfa <- ACF(USAccDeaths, value) %>% autoplot() + xlab("") + ggtitle("A") + ylim(-0.4,1)
acfd <- ACF(AirPassengers, value) %>% autoplot() + xlab("") + ggtitle("D") + ylim(-0.4,1)
acfc <- ACF(mink, value) %>% autoplot() + xlab("") + ggtitle("C") + ylim(-0.4,1)
gridExtra::grid.arrange(tp1,tp2,tp3,tp4,
                        acfa,acfb,acfc,acfd,nrow=2)
```

# White noise

## Example: White noise
\fontsize{13}{15}\sf

```{r, fig.height=3}
wn <- tsibble(t = seq_len(36), y = rnorm(36),
              index = t)
wn %>% autoplot(y)
```

## Example: White noise
\fontsize{10}{15}\sf\tabcolsep=0.1cm

```{r, echo=TRUE, warning=FALSE}
wn %>% ACF(y, lag_max = 10) %>%
  as_tibble() %>%
  tidyr::spread(lag, acf) %>%
  rename_all(function(x){paste("$r_{",x,"}$",sep="")}) %>%
  knitr::kable(booktabs=TRUE,
               escape=FALSE, align="c", digits=3,
               format.args=list(nsmall=3))
```

## Example: White noise
\fontsize{10}{15}\sf\tabcolsep=0.1cm

```{r, echo=FALSE, fig.height=3}
wn %>% ACF(y) %>% autoplot()
```

## \large Sampling distribution of autocorrelations

Sampling distribution of $r_k$ for white noise data is asymptotically N(0,$1/T$).\pause

  *  95% of all $r_k$ for white noise must lie within $\pm 1.96/\sqrt{T}$.
  * If this is not the case, the series is probably not WN.
  * Common to plot lines at $\pm 1.96/\sqrt{T}$ when plotting ACF.
These are the \alert{critical values}.

## Example: Pigs slaughtered

\fontsize{12}{13}\sf

```{r, fig.height=2.5}
pigs <- aus_livestock %>%
  filter(State == "Victoria", Animal == "Pigs",
         year(Month) >= 2014)
pigs %>% autoplot(Count/1e3) +
  xlab("Year") + ylab("Thousands") +
  ggtitle("Number of pigs slaughtered in Victoria")
```

## Example: Pigs slaughtered

```{r,fig.height=2.5,warning=FALSE}
pigs %>% ACF(Count) %>% autoplot()
```

## Example: Pigs slaughtered

Monthly total number of pigs slaughtered
in the state of Victoria, Australia, from January 2014 through December 2018
(Source: Australian Bureau of Statistics.)\pause

  * Difficult to detect pattern in time plot.
  * ACF shows significant autocorrelation for lag 2 and 12.
  * Indicate some slight seasonality.

\pause

These show the series is **not a white noise series**.

## Your turn

You can compute the daily changes in the Google stock price in 2018 using

\fontsize{11.5}{15}\sf


```{r, eval = FALSE}
dgoog <- gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) >= 2018) %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index=trading_day, regular=TRUE) %>%
  mutate(diff = difference(Close))
```

\fontsize{14}{16}\sf

Does `diff` look like white noise?
