```{r setup, cache=FALSE, echo=FALSE}
library("RColorBrewer")    # brewer.pal
library("knitr")           # opts_chunk

# terminal output
options(width = 80)

# color palette
palette(brewer.pal(6, "Set1"))

# set random number generator seed to ensure consistent runs
set.seed(0)

# code chunk options
opts_chunk$set(cache=TRUE, fig.align="center", comment=NA, echo=TRUE,
               tidy=FALSE)
```

---

```{r}
library("forecast")
```

Annual U.S. Housing Starts
--------------------------

```{r}
data <- read.csv("http://ptrckprry.com/course/forecasting/data/housingstarts.csv")
housing.starts <- data$housing.starts
date <- as.Date(data$date)
n <- length(date)
time <- 1:n
```


Time Series Plot, ACF, and PACF
-------------------------------

```{r}
plot(date, housing.starts, type="l", col=2)
Acf(housing.starts)
Pacf(housing.starts)
```


AICC for Candidate Models
-------------------------

```{r}
# Without constant:
fit.00 <- Arima(housing.starts, c(0, 0, 0), include.constant=FALSE)
fit.01 <- Arima(housing.starts, c(0, 0, 1), include.constant=FALSE)
fit.02 <- Arima(housing.starts, c(0, 0, 2), include.constant=FALSE)
fit.10 <- Arima(housing.starts, c(1, 0, 0), include.constant=FALSE)
fit.11 <- Arima(housing.starts, c(1, 0, 1), include.constant=FALSE)
fit.12 <- Arima(housing.starts, c(1, 0, 2), include.constant=FALSE)
fit.20 <- Arima(housing.starts, c(2, 0, 0), include.constant=FALSE)
fit.21 <- Arima(housing.starts, c(2, 0, 1), include.constant=FALSE)
fit.22 <- Arima(housing.starts, c(2, 0, 2), include.constant=FALSE)

# With constant:
fit.00c <- Arima(housing.starts, c(0, 0, 0), include.constant=TRUE)
fit.01c <- Arima(housing.starts, c(0, 0, 1), include.constant=TRUE)
fit.02c <- Arima(housing.starts, c(0, 0, 2), include.constant=TRUE)
fit.10c <- Arima(housing.starts, c(1, 0, 0), include.constant=TRUE)
fit.11c <- Arima(housing.starts, c(1, 0, 1), include.constant=TRUE)
fit.12c <- Arima(housing.starts, c(1, 0, 2), include.constant=TRUE)
fit.20c <- Arima(housing.starts, c(2, 0, 0), include.constant=TRUE)
fit.21c <- Arima(housing.starts, c(2, 0, 1), include.constant=TRUE)
fit.22c <- Arima(housing.starts, c(2, 0, 2), include.constant=TRUE)

# Summarize Results
models <- data.frame(p = rep(c(0, 0, 0, 1, 1, 1, 2, 2, 2), 2),
                     d = rep(0, 18),
                     q = rep(c(0, 1, 2), 6),
                     include.constant = c(rep(FALSE, 9), rep(TRUE, 9)),
                     loglik = c(fit.00$loglik, fit.01$loglik, fit.02$loglik,
                                fit.10$loglik, fit.11$loglik, fit.12$loglik,
                                fit.20$loglik, fit.21$loglik, fit.22$loglik,
                                fit.00c$loglik, fit.01c$loglik, fit.02c$loglik,
                                fit.10c$loglik, fit.11c$loglik, fit.12c$loglik,
                                fit.20c$loglik, fit.21c$loglik, fit.22c$loglik),
                     aicc = c(fit.00$aicc, fit.01$aicc, fit.02$aicc,
                                fit.10$aicc, fit.11$aicc, fit.12$aicc,
                                fit.20$aicc, fit.21$aicc, fit.22$aicc,
                                fit.00c$aicc, fit.01c$aicc, fit.02c$aicc,
                                fit.10c$aicc, fit.11c$aicc, fit.12c$aicc,
                                fit.20c$aicc, fit.21c$aicc, fit.22c$aicc)
                     )
print(models, digits=6)
```

Selected Model
--------------

```{r}
fit <- Arima(housing.starts, c(1, 0, 2), include.constant=TRUE)
print(fit)
```

Fitted model is

    [x(t) - 1442.5] = 0.425 * [x(t-1) - 1442.5]
                    + eps(t) + 0.943 * eps(t-1) + 0.464 * eps(t-2)


Residuals
---------

```{r}
resid <- residuals(fit)
plot(date, resid, type="l", col=2)
Acf(resid)
Pacf(resid)
```


Ljung-Box Goodness-of-Fit Test
------------------------------

```{r}
# fitdf = 4 (p = 1, q = 2, model includes constant)
Box.test(resid, lag=12, type = "Ljung-Box", fitdf=4)
Box.test(resid, lag=24, type = "Ljung-Box", fitdf=4)
Box.test(resid, lag=36, type = "Ljung-Box", fitdf=4)
Box.test(resid, lag=48, type = "Ljung-Box", fitdf=4)
```



Forecasts From Selected Model
-----------------------------
```{r}
plot(forecast(fit, h=30, level=95), col=2)
```


Model Without Constant
----------------------

```{r}
fit.nc <- Arima(housing.starts, c(1, 0, 2), include.constant=FALSE)
print(fit.nc)
```


Residuals (Model Without Constant)
----------------------------------

```{r}
resid.nc <- residuals(fit.nc)
plot(date, resid.nc, type="l", col=2)
Acf(resid.nc)
Pacf(resid.nc)
```


Ljung-Box Goodness-of-Fit Test (Model Without Constant)
-------------------------------------------------------

```{r}
# fitdf = 3 (p = 1, q = 2, no constant)
Box.test(resid.nc, lag=12, type = "Ljung-Box", fitdf=3)
Box.test(resid.nc, lag=24, type = "Ljung-Box", fitdf=3)
Box.test(resid.nc, lag=36, type = "Ljung-Box", fitdf=3)
Box.test(resid.nc, lag=48, type = "Ljung-Box", fitdf=3)
```



Forecasts From Model Without Constant
-------------------------------------
```{r}
plot(forecast(fit.nc, h=100, level=95), col=2)
```


