---
title: "STAT 436 / 536 - Lecture 14"
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)
knitr::opts_chunk$set(fig.align= 'center')
knitr::opts_chunk$set(fig.height= 6)
knitr::opts_chunk$set(fig.width = 6)
library(tidyverse)
library(gridExtra)
library(readr)
library(ggfortify)
scrape_BZNPM_Jul <- function(days){
  # Scrapes hourly PM2.5 readings in Bozeman for July 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=7/",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)
}
```


### State Space Models

- Most statistical analysis focuses on static parameter estimation; in other words, all of the values of the parameters are considered fixed across time. This is true for regression models, but what about ARIMA models and exponentially weighted moving average models?

\vfill

- The inherent assumption in these models is that the relationships are constant. Sometimes there are fundamental processes that might change in time such that the relationship between the variables and the response are not constant. 
\vfill
- State-space models and more specifically, Dynamic Linear Models (DLMs) can model parameter values with dynamics.
\vfill
A state space model is typically expressed using two levels, the observation equation and state equation:
\vfill
1. *Observation Equation:* 
\vfill
2. *State or Evolution Equation:* 
\vfill
Consider the example:
\begin{eqnarray*}
  y_t &=& \theta_t + \epsilon_t \; \; \; \; \; \text{(Observation Equation)} \\
  \theta_t &=& \theta_{t-1} + w_t \; \; \; \; \; \text{(Evolution Equation)},
\end{eqnarray*}
where $\epsilon_t \sim N(0, \sigma^2)$ and $w_t \sim N(0, \sigma^2_w)$. Note that $\nu$, $\sigma^2$, and $\sigma^2_w$ are all assumed known, but they can easily be estimated in an MCMC procedure.
\vfill

\newpage

#### Dynamic Linear Models

- Now consider the following representation:
\begin{eqnarray*}
  x_t &=& \tilde{F}^T \tilde{\theta_t} + v_t \\
  \tilde{\theta_t} &=& G \tilde{\theta_{t-1}} + w_t,
\end{eqnarray*}
where $\tilde{\theta_t}$ is known as the state vector for time $t,$ $\tilde{F}$ is vector of known constants or regressors, $v_t$ is observation noise, $G$ is the state evolution matrix which is typically assumed known, and $w_t$ is the innovation noise.
\vfill
- Dynamic linear models are a broad class of models that incoporate many existing models we have seen in this class.
\vfill
- What do we have if: $G$ is an identity matrix, and $w_t = 0$ $\forall t$?


\vfill
- Recall an autoregressive model can be written as
$$x_t = \alpha_1 x_{t-1} + \dots + \alpha_p x_{t-p} + w_t$$
\vfill
- How can we write an AR model in the DLM framework?
Starting from the last model, let $\tilde{F} = (x_{t-1}, \dots x_{t-p})$ and $\tilde{\theta} = (\alpha_1, \dots, \alpha_p)$
Thus $\tilde{\theta}_t = \tilde{\theta}_{t-1} = \tilde{\theta},$ so this can be re-written as
$$x_t = \tilde{F}^T \tilde{\theta} + v_t.$$
\vfill
- Now what would it mean if $w_t$ is non-zero in the two previous models?
\vfill
\newpage

#### Exponential Smoothing

- The exponential smoothing techniques (Holt-Winters) can be re-imagined in a state-space framework. In particular there is an entire book on this subject. *Forecasting with Exponential Smoothing: The State Space Approach* by Hyndman, Koehler, Ord and Snynder is a Spring book with a free e-link through the library. 
\vfill
- In this book the exponential smoothing techniques are referred to by the following three components:
    - T 
    - S 
    - E 
    
\vfill
- The methods are even generalized beyond what we have seen. For example consider the trend component, which is a combination of the level term (mean or intercept) and the growth term (slope). An additional term ($0 < \phi < 1$) is introduced to fit a dampening pattern. Let $T_h$ be the forecast trent over the next $h$ time periods.
\vfill
1. $T_h = l$
\vfill
2. $T_h = l + bh$ 
\vfill
3. $T_h = l + (\phi + \phi^2 + \dots + \phi^h)b$ 
\vfill
4. $T_h = lb^h$ 
\vfill
5. $T_h = l b^{(\phi + \phi^2 + \dots + \phi^h)}$ 
\vfill
- Furthermore, these can be generalized and combined with three levels of the seasonal component (None, Additive, Multiplicative) to create a set of 15 models. These 15 models would incorporate all of the exponential smoothing models we have seen thus far: No trend and no seasonal component would be the simple exponential smoothing method, additive and additive is the Holt-Winters 3-parameter model. 
\vfill
\newpage
- Recall that the exponential smoothing techniques do not have a probabilistic interpretation (For creating confidence intervals), but rather only produce point estimates:
$$\hat{y}_{t+1} = \alpha y_t + (1- \alpha) \hat{y}_t$$
or
\begin{eqnarray*}
  \hat{y}_{t+h} &=& l_t + b_t h\\
  l_t &=&  \alpha y_t + (1- \alpha) (l_{t-1} + b_{t-1})\\
  b_t &=& \beta (l_t - l_{t-1}) + (1-\beta) b_{t-1}
\end{eqnarray*}
\vfill
- We now consider the state-space representations for a one-step ahead prediction with Holt's Linear Method:
\begin{eqnarray*}
  \hat{y}_{t+1} &=& l_t + b_t + \epsilon_t \\
  l_t &=& l_{t-1} + b_{t-1} \alpha \epsilon_t\\
  b_t &=& b_{t-1} + \beta(l_t - l_{t-1} - b_{t-1}) = b_{t-1} + \alpha \beta \epsilon_t
\end{eqnarray*}
where $\epsilon_t = y_t - \mu_t$.
\vfill
- Using standard notation, let
\begin{eqnarray*}
 y_t &=& [1\;\; 1] x_{t-1} + \epsilon_t,\\
 x_t &=& \begin{bmatrix} 1 & 1\\ 0 & 1 \end{bmatrix} x_{t-1} + \begin{bmatrix} \alpha \\ \alpha \beta \end{bmatrix} \epsilon_t,
\end{eqnarray*}
where the state vector $x_t = (l_t, b_t)$.
\vfill
- Once a distribution is specified for $\epsilon_t$, the model is fully specified and can be estimated. Typically, $\epsilon_t \sim N (0, \sigma^2)$
\vfill
\newpage

##### State Space Model Estimation 
- The major benefit of the state-space model framework in this situation is a likelihood-based approach that permits uncertainty calculation and likelihood-based model comparison.
\vfill
- The state-space framework requires initial values for the state components, say $x_0$. From a Bayesian perspective, this is a prior distribution.
\vfill
- With only a single source of error, $\sigma_t \sim N(0, \sigma^2)$, Kalman filtering equations are not required in this setting.
\vfill
- The `forecast` package contains the `ets` function for fitting smoothing models with a state-space framework.
\vfill
- The `ets` function permits specifying a model directly using: A (additive), M (multiplicative), and N (none). Additionally, the model can be automatically selected.
\vfill

```{r, eval = F}
library(forecast)
ggtsdisplay(USAccDeaths)
fit <- ets(USAccDeaths); fit
plot(forecast(fit))
```

- Next, look at the Bozeman Air Quality data for July 3 - 31, 2018 and fit an ets model.
```{r, cache = T}
library(readr)
library(dplyr)
library(rvest)
air.quality <- scrape_BZNPM_Jul(3:31)
```

#### Exponential Smoothing Models with Regression Parameters

- One drawback of the exponential smoothing techniques is that regressor variables typically cannot be included in the model.
\vfill
- Regressor variables can take two forms: 
    - Explanatory variables such as the number of cups of coffee sold the previous day, or the weather the day before.
    - Intervention variables, which typically take the form of an indicator variables, for things like special holidays or extreme weather events. In the taxi data set some of the lowest use days were major snow storms in NYC.
    
\vfill

- The general state-space model for the exponential smoothing framework can be written as:
\begin{eqnarray*}
  y_t &=& w x_{t} + \epsilon_t\\
  x_t &=& F x_{t-1} + g \epsilon_t
\end{eqnarray*}
 now consider the extension
\begin{eqnarray*}
  y_t &=& w x_{t} + \;\;\;\; \;\;\;\;\;\; + \epsilon_t\\
  x_t &=& F x_{t-1} + g \epsilon_t,
\end{eqnarray*}
where $z$ are regression variables and $\theta$ are the associated covariates. The `ets` function does not permit variables, but instead users are encouraged to consider the DLM implementation of the ARIMA model framework.

    