```{r ts-setup, include=FALSE}
knitr::opts_knit$set(unnamed.chunk.label = "ts-")
knitr::opts_chunk$set(echo = TRUE, comment=NA, cache=TRUE, tidy.opts=list(width.cutoff=60), tidy=TRUE, fig.align='center', out.width='80%')
```

# Introduction to time series {#chap-ts}
\chaptermark{Introduction to time series}

At a very basic level, a time series is a set of observations taken sequentially in time. It is different than non-temporal data because each data point has an order and is, typically, related to the data points before and after by some process.

A script with all the R code in the chapter can be downloaded  [here](./Rcode/intro-to-ts.R).  The Rmd for this chapter can be downloaded [here](./Rmds/intro-to-ts.Rmd).


<center><img src="https://www.esrl.noaa.gov/gmd/webdata/ccgg/trends/co2_data_mlo.png"></center>


## Examples of time series {#sec-ts-examples}

```{r ts-plot-www, fig.cap = "Number of users connected to the internet"}
data(WWWusage, package="datasets")
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(WWWusage, ylab = "", las = 1, col = "blue", lwd = 2)
```

```{r ts-plot-lynx, fig.cap = "Number of lynx trapped in Canada from 1821-1934"}
data(lynx, package="datasets")
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(lynx, ylab = "", las = 1, col = "blue", lwd = 2)
```


## Classification of time series {#sec-ts-classification}

A ts can be represented as a set

$$
\{ x_1,x_2,x_3,\dots,x_n \} 
$$

For example,
$$
\{ 10,31,27,42,53,15 \}
$$
It can be further classified.

### By some _index set_ 

Interval across real time; $x(t)$    

* begin/end: $t \in [1.1,2.5]$    

Discrete time; $x_t$    

* Equally spaced: $t = \{1,2,3,4,5\}$  
* Equally spaced w/ missing value: $t = \{1,2,4,5,6\}$  
* Unequally spaced: $t = \{2,3,4,6,9\}$  


### By the _underlying process_    

Discrete (eg, total # of fish caught per trawl)

Continuous (eg, salinity, temperature)


### By the _number of values recorded_    

Univariate/scalar (eg, total # of fish caught)    

Multivariate/vector (eg, # of each spp of fish caught)


### By the _type of values recorded_    

Integer (eg, # of fish in 5 min trawl = 2413)

Rational (eg, fraction of unclipped fish = 47/951)

Real (eg, fish mass = 10.2 g)

Complex (eg, cos(2 $\pi$ 2.43) + _i_ sin(2 $\pi$ 2.43))


## Statistical analyses of time series {#sec-ts-stat-analysis}

Most statistical analyses are concerned with estimating properties of a population from a sample.  For example, we use fish caught in a seine to infer the mean size of fish in a lake.  Time series analysis, however, presents a different situation:

* Although we could vary the _length_ of an observed time series, it is often impossible to make multiple observations at a _given_ point in time

For example, one can’t observe today’s closing price of Microsoft stock more than once. Thus, conventional statistical procedures, based on large sample estimates, are inappropriate.

```{r ts-load-quantmod, message=FALSE, warning=FALSE, echo=FALSE, results='hide'}
if (!require("quantmod")) {
    install.packages("quantmod")
    library(quantmod)
}
start <- as.Date("2016-01-01")
end <- as.Date("2016-10-01")
getSymbols("MSFT", src = "yahoo", from = start, to = end)
plot(MSFT[, "MSFT.Close"], main = "MSFT")
```

## What is a time series model? {#sec-ts-definition}

We use a time series model to analyze time series data.  A _time series model_ for $\{x_t\}$ is a specification of the joint distributions of a sequence of random variables $\{X_t\}$, of which $\{x_t\}$ is thought to be a realization.  

Here is a plot of many realizations from a time series model.

```{r ts-plot-joint-dist, echo=FALSE, fig.cap="Distribution of realizations"}
set.seed(123)
nn <- 50
tt <- 40
ww <- matrix(rnorm(nn*tt), tt, nn)
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
matplot(ww, type="l", lty="solid",  las = 1,
        ylab = expression(italic(X[t])), xlab = "Time",
        col = gray(0.5, 0.4))
```

These lines represent the distribution of possible realizations.  However, we have only one realization. The time series model allows us to use the one realization we have to make inferences about the underlying joint distribution from whence our realization came.

```{r ts-plot-joint-dist-2, echo=FALSE, fig.cap="Blue line is our one realization."}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
matplot(ww, type="l", lty="solid",  las = 1,
        ylab = expression(italic(X[t])), xlab = "Time",
        col = gray(0.5, 0.4))
lines(ww[,1], col = "blue", lwd = 2)
```


## Two simple and classic time series models {#sec-ts-two-examples}

White noise: $x_t \sim N(0,1)$

```{r ts-ex-WN}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
matplot(ww, type="l", lty="solid",  las = 1,
        ylab = expression(italic(x[t])), xlab = "Time",
        col = gray(0.5, 0.4))
```

Random walk: $x_t = x_{t-1} + w_t,~\text{with}~w_t \sim N(0,1)$

```{r ts-ex-RW}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
matplot(apply(ww, 2, cumsum), type="l", lty="solid",  las = 1,
        ylab = expression(italic(x[t])), xlab = "Time",
        col = gray(0.5, 0.4))
```


## Classical decomposition {#sec-ts-classical-decomposition}

Model time series $\{x_t\}$ as a combination of

1. trend ($m_t$)  
2. seasonal component ($s_t$)  
3. remainder ($e_t$)

$x_t = m_t + s_t + e_t$


### 1. The trend ($m_t$)

We need a way to extract the so-called _signal_.  One common method is via "linear filters"

$$
m_t = \sum_{i=-\infty}^{\infty} \lambda_i x_{t+1}
$$

For example, a moving average

$$
m_t = \sum_{i=-a}^{a} \frac{1}{2a + 1} x_{t+i}
$$

If $a = 1$, then

$$
m_t = \frac{1}{3}(x_{t-1} + x_t + x_{t+1})
$$


### Example of linear filtering

Here is a time series.

```{r ts-plot-airpass, echo=FALSE, fig.cap = "Monthly airline passengers from 1949-1960"}
xx <- AirPassengers
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(xx, las = 1, ylab = "")
```

A linear filter with $a=3$ closely tracks the data.

```{r ts-plot-airpass-fltr1, echo=FALSE, fig.cap = "Monthly airline passengers from 1949-1960 with a low filter."}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(xx, las = 1, ylab = "")
## weights for moving avg
fltr <- c(1,1,1)/3
trend <- filter(xx, filter=fltr, method="convo", sides=2)
lines(trend, col = "blue", lwd = 2)
text(x = 1949, y = max(trend, na.rm = TRUE),
     labels = expression(paste(lambda, " = 1/3")),
     adj = c(0,0), col = "blue")
```

As we increase the length of data that is averaged from 1 on each side ($a=3$) to 4 on each side ($a=9$), the trend line is smoother.

```{r ts-plot-airpass-fltr2, echo=FALSE, fig.cap = "Monthly airline passengers from 1949-1960 with a medium filter."}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(xx, las = 1, ylab = "")
## weights for moving avg
fltr2 <- rep(1,9)/9
trend2 <- filter(xx, filter=fltr2, method="convo", sides=2)
lines(trend, col = "blue", lwd = 2)
lines(trend2, col = "darkorange", lwd = 2)
text(x = 1949, y = max(trend, na.rm = TRUE),
     labels = expression(paste(lambda, " = 1/3")),
     adj = c(0,0), col = "blue")
text(x = 1949, y = max(trend, na.rm = TRUE)*0.9,
     labels = expression(paste(lambda, " = 1/9")),
     adj = c(0,0), col = "darkorange")
```

When we increase up to 13 points on each side ($a=27$), the trend line is very smooth.

```{r ts-plot-airpass-fltr3, echo=FALSE, fig.cap = "Monthly airline passengers from 1949-1960 with a high filter."}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(xx, las = 1, ylab = "")
## weights for moving avg
fltr3 <- rep(1,27)/27
trend3 <- filter(xx, filter=fltr3, method="convo", sides=2)
lines(trend, col = "blue", lwd = 2)
lines(trend2, col = "darkorange", lwd = 2)
lines(trend3, col = "darkred", lwd = 2)
text(x = 1949, y = max(trend, na.rm = TRUE),
     labels = expression(paste(lambda, " = 1/3")),
     adj = c(0,0), col = "blue")
text(x = 1949, y = max(trend, na.rm = TRUE)*0.9,
     labels = expression(paste(lambda, " = 1/9")),
     adj = c(0,0), col = "darkorange")
text(x = 1949, y = max(trend, na.rm = TRUE)*0.8,
     labels = expression(paste(lambda, " = 1/27")),
     adj = c(0,0), col = "darkred")
```


### 2. Seasonal effect ($s_t$)

Once we have an estimate of the trend $m_t$, we can estimate $s_t$ simply by subtraction:

$$
s_t = x_t - m_t
$$

```{r ts-plot-airpass-decomp-seas, echo=FALSE, fig.cap = ""}
seas <- trend2 - xx
  
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(seas, las = 1, ylab = "")
# text(x = 1949, y = max(trend, na.rm = TRUE)*0.9,
#      labels = expression(paste(lambda, " = 1/9")),
#      adj = c(0,0), col = "darkorange")
```

This is the seasonal effect ($s_t$), assuming $\lambda = 1/9$, 
but, $s_t$ includes the remainder $e_t$ as well.  Instead we can estimate the mean seasonal effect ($s_t$).

```{r ts-mean-seasonal-effects, fig.cap="Mean seasonal effect."}
seas_2 <- decompose(xx)$seasonal
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(seas_2, las = 1, ylab = "")
```


### 3. Remainder ($e_t$)

Now we can estimate $e_t$ via subtraction:

$$
e_t = x_t - m_t - s_t
$$


```{r ts-errors, fig.cap="Errors."}
ee <- decompose(xx)$random
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(ee, las = 1, ylab = "")
```

## Decomposition on log-transformed data {#sec-ts-decomposition-log-data}

Let's repeat the decomposition with the log of the airline data.

```{r ts-plot-ln-airpass, fig.cap = "Log monthly airline passengers from 1949-1960"}
lx <- log(AirPassengers)
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(lx, las = 1, ylab = "")
```

### The trend ($m_t$)

```{r ts-plot-lin-trend, echo=FALSE}
tt <- as.vector(time(xx))
cc <- coef(lm(lx ~ tt))
pp <- cc[1] + cc[2] * tt
  
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot(tt, lx, type="l", las = 1,
     xlab = "Time", ylab = "")
lines(tt, pp, col = "blue", lwd = 2)
```

### Seasonal effect ($s_t$) with error ($e_t$)

```{r seas_ln_dat, echo=FALSE}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(lx-pp)
```

### Mean seasonal effect ($s_t$)

```{r mean_seas_effects, echo=FALSE}
## length of ts
ll <- length(lx)
## frequency (ie, 12)
ff <- frequency(lx)
## number of periods (years); %/% is integer division
periods <- ll %/% ff
## index of cumulative month
index <- seq(1,ll,by=ff) - 1
## get mean by month
mm <- numeric(ff)
for(i in 1:ff) {
  mm[i] <- mean(lx[index+i], na.rm=TRUE)
}
## subtract mean to make overall mean=0
mm <- mm - mean(mm)
seas_2 <- ts(rep(mm, periods+1)[seq(ll)],
               start=start(lx), 
               frequency=ff)
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(seas_2, las = 1, ylab = "")
```

### Remainder ($e_t$)

```{r ln_errors}
le <- lx - pp - seas_2
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(le, las = 1, ylab = "")
```

