```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Time Series Data

### Lake Huron Time Series

Annual measurements of the level of Lake Huron in feet

```{r}
par(mar = c(3.2, 3.2, 0.5, 0.5), mgp = c(2, 0.5, 0), bty = "L")
data(LakeHuron)
plot(LakeHuron, ylab = "Depth in feet", xlab = "Year", las = 1)
grid()
```

### $\mathrm{CO}_{2}$ Concentration

Atmospheric concentrations of $\mathrm{CO}2$ are expressed in parts per million (ppm) and reported in the preliminary 1997 SIO manometric mole fraction scale.

```{r}
data(co2)

par(mar = c(3.8, 4, 0.8, 0.6))
plot(co2, las = 1, xlab = "", ylab = "")
mtext("Time (year)", side = 1, line = 2)
mtext(expression(paste("CO"[2], " Concentration (ppm)")), side = 2, line = 2.5)

# Seasonal and Trend decomposition using Loess (STL)
par(mar = c(4, 3.6, 0.8, 0.6))
stl <- stl(co2, s.window = "periodic")
plot(stl, las = 1)
```

### U.S. monthly unemployment rates

```{r}
library(quantmod)
getSymbols("UNRATE", src = "FRED")
head(UNRATE); tail(UNRATE)
chartSeries(UNRATE) 
```


## ARMA and ACF

```{r}
set.seed(123)
n = 100
WN <- rnorm(n)

par(mfrow = c(4, 2), mar = c(3.6, 3.6, 0.8, 0.6))
plot(1:n, WN, type = "l", las = 1, xlab = "", ylab = "")
mtext("WN", side = 2, line = 2)
acf(WN, xlab = "", ylab = "", main = "", las = 1)
mtext("ACF", side = 2, line = 2, cex = 0.8)
# AR(1) phi = 0.8
AR <- arima.sim(n = n, model = list(ar = 0.8))
plot(1:n, AR, type = "l", las = 1, xlab = "", ylab = "")
mtext("AR(1)", side = 2, line = 2)
acf(AR, xlab = "", ylab = "", main = "", las = 1)
mtext("ACF", side = 2, line = 2, cex = 0.8)
# MA(1) theta = 0.5
MA <- arima.sim(n = n, model = list(ma = 0.5))
plot(1:n, MA, type = "l", las = 1, xlab = "", ylab = "")
mtext("MA(1)", side = 2, line = 2)
acf(MA, xlab = "", ylab = "", main = "", las = 1)
mtext("ACF", side = 2, line = 2, cex = 0.8)
# ARMA(1, 1) phi = 0.8, theta = 0.5
ARMA <- arima.sim(n = n, model = list(ar = 0.8, ma = 0.5))
plot(1:n, ARMA, type = "l", las = 1, xlab = "", ylab = "")
mtext("ARMA(1,1)", side = 2, line = 2)
mtext("Time", side = 1, line = 2)
acf(ARMA, xlab = "", ylab = "", main = "", las = 1)
mtext("ACF", side = 2, line = 2, cex = 0.8)
mtext("Lag", side = 1, line = 2)
```

## Lake Huron Case Study

### Detrend 

```{r}
## Let us create a 'years' variable.
years <- time(LakeHuron)

plot(LakeHuron, ylab = "Depth (ft)", xlab = "Year", las = 1)
grid()
par(mfrow = c(2, 1), mar = c(4, 4, 1, 1))
acf(LakeHuron, xlab="Lag in years", ylab = "sample ACF", main = "")
pacf(LakeHuron, xlab="Lag in years", ylab = "sample PACF", main = "")

# Estimate the linear trend
lm <- lm(LakeHuron ~ years)
plot(LakeHuron, ylab = "Depth (ft)", xlab = "Year", las = 1); grid()
abline(lm, col = "blue")

deTrend <- resid(lm)

par(mfrow = c(2, 1), mar = c(4, 4, 1, 1))
acf(deTrend, xlab="Lag in years", ylab = "sample ACF", main = "")
pacf(deTrend, xlab="Lag in years", ylab = "sample PACF", main = "")

hist(deTrend, main = "", xlab = "", las = 1)
qqnorm(deTrend, main = "", pch = 16, cex = 0.8); qqline(deTrend, col = "blue")
```

### Model Selection/Fitting

```{r}
## AR(1)
ar1.model <- arima(deTrend, order = c(1, 0, 0))
ar1.model
ar1.resids <- resid(ar1.model)
plot(1875:1972, ar1.resids, type = "l", xlab = "year", ylab = "AR1 residuals")

## Sample ACF and PACF of the residuals
par(mfrow = c(1, 2))
acf(ar1.resids, ylab = "sample ACF", xlab = "lag (year)")
pacf(ar1.resids, ylab = "sample PACF", xlab = "lag (year)")
## Normal Q-Q plot for the residuals
qqnorm(ar1.resids, main = ""); qqline(ar1.resids, col = "blue")
## Test for time dependence for the residuals 
Box.test(ar1.resids, type = "Ljung-Box")

## AR(2)
ar2.model <- arima(deTrend, order = c(2, 0, 0))
## summarize the model
ar2.model
## calculate the residuals
ar2.resids <- resid(ar2.model)
## time series plot of the residuals
plot(1875:1972, ar2.resids, type = "l", xlab = "year", ylab = "AR2 residuals")
## Sample ACF and PACF of the residuals
par(mfrow=c(1,2))
acf(ar2.resids,  ylab = "sample ACF", xlab = "lag (year)")
pacf(ar2.resids,  ylab = "sample PACF", xlab = "lag (year)")

## Normal Q-Q plot for the residuals
qqnorm(ar2.resids, main = ""); qqline(ar2.resids, col = "blue")

## Test for time dependence for the residuals 
Box.test(ar2.resids, type = "Ljung-Box")

## Fit the ARMA(2, 1) model
arma21.model <- arima(deTrend, order = c(2, 0, 1))
## summarize the model
arma21.model
## calculate the residuals
arma21.resids <- resid(arma21.model)
## time series plot of the residuals
plot(1875:1972, arma21.resids, type = "l", xlab = "year", ylab = "ARMA(2, 1) residuals")
## Sample ACF and PACF of the residuals
par(mfrow=c(1,2))
acf(arma21.resids,  ylab = "sample ACF", xlab = "lag (year)")
pacf(arma21.resids,  ylab = "sample PACF", xlab = "lag (year)")
## Normal Q-Q plot for the residuals
qqnorm(arma21.resids, main = ""); qqline(arma21.resids, col = "blue")
## Test
Box.test(arma21.resids, type = "Ljung-Box")

arma21.model <- arima(deTrend, order=c(2, 0, 1))

AIC(ar1.model); AIC(ar2.model); AIC(arma21.model)
```

### AR(2) Fitting

```{r}
## Fit the full regression model plus AR(2) errors
full.model <- arima(LakeHuron, order = c(2, 0, 0),
                    xreg = cbind(rep(1, length(LakeHuron)), years), include.mean = FALSE)
full.model

## Examine the residuals of the model
par(mfrow = c(2, 2))
res <- full.model$residuals

plot(res, type = "l", xlab = "year", ylab = "AR(2) residuals", las = 1)
abline(h = 0, col = "blue")
qqnorm(res, main = ""); qqline(res, col = "blue")
acf(res,  ylab = "sample ACF", xlab = "lag (year)")
pacf(res,  ylab = "sample PACF", xlab = "lag (year)")
```

### Forecasting

```{r}
library(forecast)
(fit <- Arima(LakeHuron, order = c(2, 1, 0)))
checkresiduals(fit)
autoplot(forecast(fit, level = c(50, 95)))
```



