---
title: "STAT 436 / 536 - Lecture 5"
date: September 12, 2018
output: pdf_document
---

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

### Forecasting Strategies

Consider an hourly time series of particulate matter less than 2.5 micrometers from Bozeman, MT in August of 2018.  

\vfill

```{r echo=F}
library(rvest)
library(dplyr)
scrape_BZNPM_AUG <- function(days){
  # Scrapes hourly PM2.5 readings in Bozeman for August 2018
  # inputs: sequence of days
  # outputs: data frame that contains day, hour, and hourly average PM2.5 concentration.
  smoke.df <- data.frame(day = NULL, hour = NULL, conc = NULL)
  for (d in days){
    bzn_aq <- read_html(paste("http://svc.mt.gov/deq/todaysair/AirDataDisplay.aspx?siteAcronym=BH&targetDate=8/",d,"/2018", sep=''))
    daily.smoke <- cbind(rep(d,24),html_table(bzn_aq)[[1]][,1:2])
    colnames(daily.smoke) <- c('day','hour','conc')
    daily.smoke$conc[daily.smoke$conc == 'DU'] <- 'NA'
    daily.smoke$conc <- as.numeric(daily.smoke$conc)
    smoke.df <- bind_rows(smoke.df,daily.smoke)
  }
  return(smoke.df)
}
air.quality <- scrape_BZNPM_AUG(1:31)
aq <- ts(air.quality$conc, freq = 24, start = 1)
library(ggfortify)
autoplot(aq) + ggtitle(label = 'PM2.5 Measurements in Bozeman, MT for August 2018', subtitle = 'Source: MT DEQ') + ylab('µg/m3') + xlab('Day')
```

- What information would be useful to predict the PM2.5 measurments on September 1?  How about the PM2.5 measurements for tomorrow?
\vfill
    -  \vfill
    -  \vfill
    -  \vfill
    -  \vfill
    
\newpage

### Leading Variables and Associated Variables

- In some situations another variable may be useful informing our forecasts. A common example of this is in the housing industry where building approvals would be a *leading indicator* of building activity. \vfill
- Consider natural gas consumption and heating degree days (HDD).  \vfill
- What do you assume is the relationship between the heating degree days during a given month and the reported natural gas consumption? \vfill
```{r, echo = F, fig.height=6}
library(readr)
library(dplyr)
library(gridExtra)
ng.consumption <- read_csv("/Volumes/GoogleDrive/My Drive/Teaching/STAT436_536/data/MT_NG.csv")
ng.consumption.small <- ng.consumption[as.numeric(substr(ng.consumption$Month,5,6)) < 18,]
ts.consumption <- ts(rev(ng.consumption.small$MT_NG), frequency = 12, start = c(2000,1))
weather <- read_csv("/Volumes/GoogleDrive/My Drive/Teaching/STAT436_536/data/BZN_weather.csv")
sac.weather <- weather %>% filter(NAME == 'SACAJAWEA, MT US')
ts.weather <- ts(sac.weather$HTDD[-c(1:3)], frequency = 12, start = c(2000,1))
library(ggfortify)
f1 <- autoplot(ts.consumption) + ggtitle(label = "Monthly Natural Gas Consumption in Montana", subtitle = 'Source: US DOE') + ylim(0,4500) + ylab('million cubic feet')
f2 <- autoplot(ts.weather) + ggtitle(label = "Monthly Heating Degree Days at Sacajawea School, Bozeman, MT", subtitle = 'Source: NOAA') + ylim(0,1500) + ylab('Heating Degree Days')
grid.arrange(f1,f2,nrow=2)
```
\vfill

#### Cross-Correlation

- Suppose two times series models, say x and y, are stationary with a constant mean and variance. The variables, x and y, may be serially correlated, perhaps with different time lags.
\vfill
- The cross covariance function can be written as:  \vfill

\newpage
- Q: Is the cross-correlation function symmetric? In other words, does $\gamma_k(x,y) = \gamma_k(y,x)$?
\vfill
- The sample cross-correlation function (ccf) can be written as 
\vfill
- Similarly to single variable autocorrelation calculations, sample based computations can be made:
    - Sample cross covariance:  \vfill
    - Sample cross correlation:  \vfill
    
- Note the cross correlation and cross covariance functions are defined on the de-seasonalized and de-trended time series, so it will not capture similar patterns in those structures. See HW2 Q1 for a demonstration.
\vfill

```{r, fig.cap = 'No Decomposition of the Time Series'}
acf(ts.union(ts.weather,ts.consumption))
```

```{r, fig.cap = 'Using the decomposed random component'}
hdd.random <- decompose(ts.weather)$random %>% window(start=c(2001,1), end = c(2016,12))
gas.random <- decompose(ts.consumption)$random %>% window(start=c(2001,1), end = c(2016,12))
acf(ts.union(hdd.random,gas.random))
```

- 
\vfill
