---
title: "Lab 11"
author: "Name here"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE)
library(readr)
library(ggplot2)
library(dplyr)
library(lubridate)
library(ggfortify)
library(forecast)
library(tseries)
library(rugarch)
```
Please use D2L to turn in both the PDF/ Word output and your R Markdown file.


## Q1. Airline Passenger Data
For this question we will implement a GARCH model using with covariates for the airline passenger data. The intention in this lab is to be the final *teach you to teach yourself* moment from this class as we didn't directly cover this in lecture.

Assume you found the code example on stack overflow (re-created below) [https://stats.stackexchange.com/questions/93815/fit-a-garch-1-1-model-with-covariates-in-r](https://stats.stackexchange.com/questions/93815/fit-a-garch-1-1-model-with-covariates-in-r) and then looked at the vignette for the `rugarch` package [https://cran.r-project.org/web/packages/rugarch/vignettes/Introduction_to_the_rugarch_package.pdf](https://cran.r-project.org/web/packages/rugarch/vignettes/Introduction_to_the_rugarch_package.pdf).

```{r, eval = F}
## Model parameters
nb.period  <- 1000
omega      <- 0.00001
alpha      <- 0.12
beta       <- 0.87
lambda     <- c(0.001, 0.4, 0.2)

## Generate some covariates
set.seed(234)
ext.reg.1 <- 0.01 * (sin(2*pi*(1:nb.period)/nb.period))/2 + rnorm(nb.period, 0, 0.0001)
ext.reg.2 <- 0.05 * (sin(6*pi*(1:nb.period)/nb.period))/2 + rnorm(nb.period, 0, 0.001)
ext.reg   <- cbind(ext.reg.1, ext.reg.2)

## Generate some GARCH innovations
sim.spec    <- ugarchspec(variance.model     = list(model = "sGARCH", garchOrder = c(1,1)), 
                          mean.model         = list(armaOrder = c(0,0), include.mean = FALSE),
                          distribution.model = "norm", 
                          fixed.pars         = list(omega = omega, alpha1 = alpha, beta1 = beta))
path.sgarch <- ugarchpath(sim.spec, n.sim = nb.period, n.start = 1)
epsilon     <- as.vector(fitted(path.sgarch))

## Create the time series
y <- lambda[1] + lambda[2] * ext.reg[, 1] + lambda[3] * ext.reg[, 2] + epsilon

## Data visualization
par(mfrow = c(3,1))
plot(ext.reg[, 1], type = "l", xlab = "Time", ylab = "Covariate 1")
plot(ext.reg[, 2], type = "l", xlab = "Time", ylab = "Covariate 2")
plot(y, type = "h", xlab = "Time")
par(mfrow = c(1,1))

## Fit
fit.spec <- ugarchspec(variance.model     = list(model = "sGARCH",
                                                 garchOrder = c(1, 1)), 
                       mean.model         = list(armaOrder = c(0, 0),
                                                 include.mean = TRUE,
                                                 external.regressors = ext.reg), 
                       distribution.model = "norm")
fit      <- ugarchfit(data = y, spec = fit.spec)

## Results review
fit.val     <- coef(fit)
fit.sd      <- diag(vcov(fit))
true.val    <- c(lambda, omega, alpha, beta)
fit.conf.lb <- fit.val + qnorm(0.025) * fit.sd
fit.conf.ub <- fit.val + qnorm(0.975) * fit.sd
```


##### a.
What is the purpose of the `variance.model` argument in the `ugarchspec` model?

##### b.
What is the purpose of the `external.regressors` argument in the `mean.model` option with in the `ugarchspec` specification?

##### c.
Construct an `external.regressors` matrix for the airline data that contains a linear trend term and dummy variable for the month.

##### d.
In the code example from stack overflow interpret the output from the following command `print(fit.val)`.

##### e.
Update the code to fit the model specified for the simulated data using the airline passengers and summarize the results.

## Q2. Housing Data Variance
While this is not a time series question, recall the housing dataset from Lab 5.

```{r}
# create plot
housing <- read_csv('http://math.montana.edu/ahoegh/teaching/stat408/datasets/HousingSales.csv')
ggplot(aes(y = Closing_Price, x = Living_Sq_Ft), data = housing) + geom_point() + geom_smooth()
```

where the variance increased as a function of price. A fan shaped residual plot is common in scenarios like this.

Essentially, time series analysis boils down to creative regression modeling with complex variance structure. In a time series setting, variance can be modeled as a function of proximity in time, or proximity in a seasonal sense as in SARIMA. 

One solution might be to allow different variance terms within each zipcode. However, often the variance is related to a continuous covariate, such as the square footage of the house. In this situation, we need another creative solution to account for the variance structure. Write out a model that is structured so that the variance is a function covariates. How would you fit this model to estimate these parameters?

While not relevant to this particular question, note that the `ugarchspec()` function also permits covariates for the variance model.