---
title: "STAT 224 Chapter 8"
output:
  html_document:
    theme: cerulean
    css: styles.css
    code_folding: show
    df_print: kable
    number_sections: yes
    self_contained: yes
    toc: FALSE
    toc_depth: 4
    toc_float: FALSE
params:
  dotabs:  TRUE
---


```{r test-main, child = 'preamble.Rmd'}
```

# Chapter 8: Correlated Errors `r if(params$dotabs) paste('{.tabset}')`

### Diagnostics


### Another Problem to Solve: Autocorrelation

  + Regression assumption:   
    error terms $\epsilon_i$ and $\epsilon_j$ for observations $i$ and $j$ are uncorrelated. 

What does this mean?

  Now, we know that observed residuals are correlated (due to minimizing least squares)
    + We still expect apparent randomness (no systematic patterns) if regression model OK

  + If observations follow a natural sequential order (in terms of how measured, collected, etc),    
    then residuals may appear to have distinct serial patterns across time (order).
      + This suggest lack of independence. 
      + This is called **autocorrelation**.

Sketch some examples.

### Causes of autocorrelation

#### Time and space

Autocorrelation most often occurs often when 

  + The same observations are measured in adjacent time periods
    + stock market volume, weather parameters, etc, today and tomorrow.
    
  + Observations are spatially close, 
    + temperature measured  in Illinois, Indiana, and Michigan
      
  + Experiments are run on "batches" of animals, cells, at specific times, etc
  
  + Tall trees compete for sunlight, thus are spaced apart with small plants in between
    + spatially, tall plants spaced apart, short plants closer together

Autocorrelation occurs since information/phenomena change slowly through time/space,    
so factors influencing a variable are likely to be similar in adjacent periods/places.

#### Omission of a time/space-related predictor

  + This is usually a rectifiable problem   
    (like we solved non-linearity and non-constant spread via transformations)
  
  + For example, if experiments are run in batches, include the batch number as a predictor.

  + Autocorrelation can also be negative:    
    Tall trees are next to short trees, which are next to tall trees, which are next to short trees,...
 
#### So, what's the problem with autocorrelation?

<!--
  + The presence of autocorrelation implies heteroscedasticity    
    + variance of $Y$ at $x$ not same for all $x$) 
-->
    
Autocorrelation can affect OLS analysis in the several ways:

  +  OLS estimates are still unbiased but are not efficient (no longer have minimum variance)
  +  $\sigma^2$ and $\se(\beta)$s may be seriously under- (or over-) estimated. 
  + positive correlation, often under-estimation of $\sigma^2$ and $\se(\beta)$s
  + negative correlation (less common), often over-estimation
      
Underestimate of variability which leads to increased false positive decisions:   
rejecting the null hypothesis $\beta = 0$ when one should not

Whether positive or negative correlation,    
confidence intervals for $\beta$s may no longer be valid.

#### Detecting autocorrelation

Plots:

  + Standardized residuals vs. time (or "index" = order of data collection)
  + Examine the sequence plot: signs of residuals by time/index:    
    +  $+++++---++++-----+++$
  
Statistical tests:

  + The **runs test**: use a statistical test based on the sequence plot
  + Durbin-Watson statistic: Measure the correlation in the residuals
  
### Example: Consumer expenditure vs. money stock

Textbook Table 8.1 shows the data

Quarterly data from 1952 to 1956 on   
$y$ = U.S. consumer expenditure    
$x$ = U.S. money stock

Units: billions of current dollars

We regress the expenditure on money stock and plot vs. time   
(these data are naturally collected in time order)

```{r}
# expendData <- read.delim("http://statistics.uchicago.edu/~collins/data/RABE5/P211.txt")
expendData <- read.delim("./data/P211.txt")
# order data by time (year and quarter)
expendData <- arrange(expendData, Year, Quarter) 
head(expendData, 6)
```

```{r}
## fit linear regression: expenditures on money stock
lmfit <- lm(Expenditure ~ Stock, data = expendData)

tidy(lmfit)
glance(lmfit)[-c(7:10)]
```

### Plots and statistical tests for autocorrelation

#### Index plot

For  time series data (data collected where time order is relevant),    
a useful plot for analysis is the  index plot 

  + plot of the standardized residual versus 
    + time or 
    + some index of the collection order, 
    + whichever makes more sense. 
    
If errors are correlated,    
residuals of the same sign ($+$ or $-$) will occur in clusters or bunches called **runs**.

```{r}
## add residuals, standardized residuals, and an index variable to dataset
n <- nrow(expendData) ## sample size
expendData <- mutate(expendData, e = residuals(lmfit), r = rstandard(lmfit), index = 1:n)
```

```{r fig.height=4}
## index plot
gf_point(r ~ index, data = expendData, size = 2) %>% 
  gf_line() %>%
  gf_labs(x = "Quarter", y = "standardized residuals") %>%
  gf_lims(y = c(min(min(~r, data = expendData), -2.1), max(max(~r, data = expendData), 2.1))) %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")
```  

We connected the points from left to right to make it easier to see autocorrelation.

Here, the sequential residuals seem to follow an ordered pattern    
(except for a few points) with respect to magnitude and sign

#### One statistic = number of runs

**Number of runs:**

The visual impression can be formally confirmed by counting      
the number of sequences, batches, clusters, (**runs**) of either positive or negative residuals

Statistic = number of runs

**Sequence Plot:**

One-dimensional plot of the signs of the residuals, taken in the order of the observations. 

For the expenditure vs. money stock example, the sequence plot is

```{r}
expendData <- mutate(expendData, positive = 1 * (r > 0), sign = factor(positive, labels = c("-", "+")))
```

$`r expendData$sign`$

<br>
If you are curious, here is what I have stored in the dataset up to this point:
```{r}
head(expendData, 6)
```

Now, back to looking at the sequence plot:

$`r expendData$sign`$

```{r}
## rle = run length encoding (in base R)
rle(as.vector(expendData$positive)) ## runs and their lengths
nruns <- length(rle(as.vector(expendData$sign))$lengths) ## number of runs

n <- nrow(expendData) ## sample size
n1 <- sum(expendData$positive) ## number of positive residuals
n2 <- n - n1 ## number of negative residuals
```

The sequence plot indicates that there are $`r nruns`$ runs. 

  + There are $n_1 = `r n1`$ positive residuals and $n_2 = `r n2`$ negative residuals.
  + Under the null hypothesis of no autocorrelation,    
    the expected number of runs $\mu$ and its variance $\sigma^2$ can be shown to be:
    
$\displaystyle \mu=\frac{2n_1n_2}{n_1+n_2}+1$

and

$\displaystyle \sigma^2=\frac{2n_1n_2(2n_1n_2-n_1-n_2)}{(n_1+n_2)^2(n_1+n_2-1)}$

```{r}
mu <- 1 + 2*n1*n2 / (n1 + n2)

sigmasqr <- 2*n1*n2 * (2*n1*n2 - n1 - n2) / ( (n1 + n2)^2 * (n1 + n2 -1) )
```

In our example, $n_1=`r n1`$ and $n_2=`r n2`$, so   
the expected number of runs is $\mu=`r mu`$ and    
the standard deviation is $\sigma = `r sqrt(sigmasqr)`$ 

We observed only $`r nruns`$ runs; which is `r round( (mu - nruns) / sqrt(sigmasqr), 1)` SDs below $\mu$.   

Is this unusual under the null hypothesis of no autocorrelation? 

#### The runs test

A well-known nonparametric test is called a **"runs test"**

In R, load the `tseries` (time series) package,    
then use the `runs.test` function  

The `runs.test` function expects you to give it a factor with two and only two levels.  
For example, 

`library(tseries)`  
`runs.test( as.factor( rstandard(lmfit) > 0) )`

```{r, message=FALSE, warning=FALSE}
library(tseries, warn.conflicts = FALSE)
data.frame(sample.size=n, number.of.runs=nruns, positive=n1, negative=n2)
runs.fit <- runs.test(expendData$sign)
runs.test.statistic <- runs.fit$statistic
runs.test.pvalue <- runs.fit$p.value
```

When I ran this code, I got a p-value of $`r round(runs.test.pvalue, 4)`$.   
This indicates we should reject the null hypothesis of no autocorrelation.

We kind of expected this given the appearance of the index plot.

#### The Durbin-Watson (DW) approach

A common approach to detecting serial correlation is known as the Durbin-Watson test. 
 
The test is based on assuming that successive errors are correlated.

The model is written this way:

\[ Y_t = \beta_0 + \beta_1 x_{1t} + \cdots + \beta_p x_{pt} + \epsilon_t\]

where $t$ is the time index and $\epsilon_t$ is the $t^{th}$ error in time.

If the errors are correlated, **perhaps** each residual ($\epsilon_t$):

  + depends **linearly** on the size and sign of the residual before it in time ($\epsilon_{t-1}$)m,
  + has correlation $\rho = \cor(\epsilon_t, \epsilon_{t-1})$, and
  + is **not** correlated with any earlier residuals ($\epsilon_{t-2}, \epsilon_{t-3}, \ldots$).
  
If this model is reasonable, then
\[
\epsilon_t=\gamma_0 + \gamma_1\epsilon_{t-1}+\omega_t = \rho\epsilon_{t-1}+\omega_t, \qquad |\rho|<1, 
\]
where $\rho$ is the correlation coefficient between successive errors   
and $\omega_t$ are independent, normally distributed, with constant variance across the errors ($\epsilon$s).

Why is the intercept ($\gamma_0$) equal to zero?
  
  + By definition, intercept = (mean of response) $-$ slope \times (mean of predictor).
  + mean of response = $E(\epsilon_t) = 0$
  + mean of predictor = $E(\epsilon_{t-1}) = 0$
  + $\gamma_0 = 0 - \gamma_1(0) = 0$
  
Why is the slope ($\gamma_1$) equal to $\rho$?

  + By definition, slope = correlation(response, predictor) $\times \displaystyle{\frac{\text{SD(response)}}{\text{SD(predictor)}}}$
    + correlation(response, predictor) = $\cor(\epsilon_t, \epsilon_{t-1}) = \rho$
    + SD(response) = SD($\epsilon_t$) $= \sigma$
    + SD(predictor) = SD($\epsilon_{t-1}$) $= \sigma$
  + $\gamma_1 = \rho \times (\sigma\,/\,\sigma) = \rho$
  
Thus, if $\epsilon_t$ depends linearly on $\epsilon_{t-1}$ and    
is not correlated with the size of earlier errors, then
this model is reasonable, then

\[
\epsilon_t=\rho\epsilon_{t-1}+\omega_t, \qquad |\rho|<1, 
\]

The relationship is likely more complex than this, but it's a nice, simple starting point.  
We call this a **"first-order autoregressive model"**

#### The Durbin-Watson statistic

The DW statistic (often labeled $d$) is defined as:
\[
d\;=\; \frac{\sum_{t=2}^n (e_t-e_{t-1})^2}{\sum_{t=1}^n e_t^2} \;=\; \frac{\sum_{t=2}^n (e_t-e_{t-1})^2}{\SSE},
\]
where $e_t$ is the $t^{th}$ observed residual with the data arranged in time order. 

**NOTE:** Yes, the denominator of $d$ sums from $t=1$, the but numerator cannot.    
The numerator has $(n-1)$ terms:
$\;e_t\; (e_2, e_3, \ldots, e_n)\;$ and $\;e_{t-1}\; (e_1, e_2, \ldots, e_{n-1})$

The DW statistic  is used to test $H_0:\rho=0$ vs $H_1:\rho >0$.    
When $\rho=0$ cannot be rejected, we conclude that the errors are uncorrelated.

**NOTE:** $H_1: \rho < 0$ is not common to test or observe,   
but we can by using the statistic ($4-d$) instead of $d$ and proceeding as below.

#### Properties of the DW statistic

We can show that the $d$ statistic can be approximated by $\widehat{\rho}$:

\[d \;\approx\; 2(1-\widehat{\rho}),  
\qquad\textrm{ where }\qquad \widehat{\rho} 
\;=\; \frac{\sum_{t=2}^n e_te_{t-1}}{\sum_{t=1}^n e_t^2}
\;=\; \frac{\sum_{t=2}^n e_te_{t-1}}{\SSE}.
 \]

**NOTE:** Yes, the denominator of $\widehat{\rho}$ sums from $t=1$, the but numerator cannot.    
The numerator has $(n-1)$ terms:
$\;e_t\; (e_2, e_3, \ldots, e_n)\;$ and $\;e_{t-1}\; (e_1, e_2, \ldots, e_{n-1})$

When $\rho = 0$ (no autocorrelation), $d \approx 2$.   
When $\rho = 1$ (perfect positive correlation), $d \approx 0$.    
When $\rho = -1$ (perfect negative correlation), $d \approx 4$.

The range of $d$ is approximately $0 \le d \le 4$.

So, when $d$ is close to 2, there is not much autocorrelation.   
When $d$ is small, it suggests evidence of positive autocorrelation ($\rho > 0$).    
When $(4-d)$ is small, it suggests evidence of negative autocorrelation ($\rho < 0$).    

If your alternative hypothesis is $H_1: \rho < 0$  ($4-d$) as your statistic.   
Usually, the one-sided hypothesis is proposed prior to looking at the data (based on context).

#### Using the DW statistic

How close is considered "close to 2"?   
For positive autocorrelations, ($d$ is below 2 leads us to reject null), use the following process:

 If $d<d_L$, we reject the null hypothesis that there is no positive autocorrelation.     
 
 If $d>d_U$, we fail to reject the null hypothesis. There is not much evidence of positive autocorrelation.   
 
 If $d_L < d < d_U$, the test is inconclusive.
 
 
Like many nonparametric tests, values are tabled for reference.   
See the textbook Appendix: Tables A.6 and A.7.

In R, there are a couple choices.

`lmfit <- lm(y ~ x, data = mydata)`   

`library(car)`   
`durbinWatsonTest(lmfit)`

`library(lmtest)`   
`dwtest(lmfit)`

**SHOW R CODE TO SEE RESULTS**

```{r}
library(car, warn.conflicts = FALSE)
durbinWatsonTest(lmfit, alternative = "positive")
```



```{r,echo=TRUE}
library(lmtest, warn.conflicts = FALSE)
dwtest(lmfit, alternative = "greater")
```

Not surprisingly (given our earlier look at the residuals by time and the runs test),    
the DW statistic also indicates strong autocorrelation.

### Transformation

### Removal of autocorrelation: transformation

When there is autocorrelation present in the errors,     
one way to remove it is to transform the variables. 

#### Cochrane-Orcutt transformation

The Cochrane-Orcutt transformation (textbook Section 8.4) is one method.

To see how it works, can write mathematically the errors for two adjacent observations:

$$\begin{align}
\epsilon_t & \;=\; y_t-\beta_0-\beta_1 x_t \\
\epsilon_{t-1} & \;=\; y_{t-1}-\beta_0-\beta_1 x_{t-1}
\end{align}$$

When errors are first-order autocorrelated with correlation coefficient $\rho$, we have 

$\epsilon_t=\rho\epsilon_{t-1}+\omega_t$

where the errors $\omega_t$ are independent, normally distributed, with constant variance across $\epsilon_t$s


Plugging in the above expressions for two adjacent residuals gives 

$$\begin{align}
\epsilon_t & \;=\; \rho\epsilon_{t-1}+\omega_t\\
(y_t-\beta_0-\beta_1 x_t) & \;=\; \rho \left(y_{t-1}-\beta_0-\beta_1 x_{t-1} \right)+\omega_t.
\end{align}$$


Rearranging the above equation we have
\[
(y_t-\rho y_{t-1}) \;=\; \beta_0(1-\rho)  + \beta_1( x_t-\rho x_{t-1}) + \omega_t.
\]

This is equivalent to the linear model 
\[ y_t^\star=\beta_0^\star+\beta_1^\star x_t^\star +\omega_t \]

with    
$y_t^\star \;=\; (y_t-\rho y_{t-1})$,    
$x_t^\star \;=\; (x_t-\rho x_{t-1})$,    
$\beta_0^\star \;=\; \beta_0(1-\rho)$, and    
$\beta_1^\star \;=\; \beta_1$

This linear model (with error term $\omega_t$, which are uncorrelated by definition) is then fit to the data

#### Implementing Cochrane-Orcutt transformation

  1. Compute $\b0$ and $\b1$ by OLS as usual: $Y_t = \beta_0 + \beta_1 x_t + \epsilon_t$
  
  2. Calculate the residuals from Step 1 and estimate $\rho$: $\displaystyle \widehat{\rho} \;=\; \frac{\sum_{t=2}^n e_te_{t-1}}{\sum_{t=1}^n e_t^2}$
  
  3. Compute OLS with response $(y_t - \widehat{\rho} y_{t-1})$ and predictor $(x_t-\widehat{\rho} x_{t-1})$
  
  Let $\b0 = \displaystyle{\frac{\b0^*}{1-\widehat{\rho}}}$ and $\b1 = \b1^*$.
  
  $\se(\b1)$ can be obtained from the OLS estimate: $\se(\b1^*)$.    
  $\se(\b0)$ can be esimated from $\widehat{\rho}$ and the OLS estimate: $\se(\b0^*)$.
  
  That is, $\displaystyle \se(\b0) \approx \frac{\se(\b0^*)}{1-\widehat{\rho}}$
  
  4. If the new residuals from  Step 3 still show autocorrelation, go back to Step 2 using $\b0$ and $\b1$ from Step 3.
  
Repeat until residuals after Step 3 show no evidence of autocorrelation.

#### Example: Expenditures vs. money stock

**Step 1:**

```{r}
lmfit <- lm(Expenditure ~ Stock, data = expendData)

head(expendData, 6)
```

**Step 2:**

```{r}
# estimate rho (correlation coefficient)
# "by hand"
e <- residuals(lmfit)
et <- e[2:n]
et1 <- e[1:(n-1)]
rhohat <- sum(et * et1) / sum(e^2)
rhohat
## or grab rhohat from the Durbin-Watson test
durbinWatsonTest(lmfit, alternative = "positive")
objects(durbinWatsonTest(lmfit, alternative = "positive"))
rhohat <- durbinWatsonTest(lmfit, alternative = "positive")$r
rhohat
```

$\widehat{\rho} = `r rhohat`$

**Step 3:**

```{r}
y <- with(expendData, Expenditure - rhohat * lag(Expenditure))
x <- with(expendData, Stock - rhohat * lag(Stock))

lmfit.CO1 <- lm(y ~ x)
tidy(lmfit.CO1)

b0star <- tidy(lmfit.CO1)$estimate[1]
b1star <- tidy(lmfit.CO1)$estimate[2]
```

$\b0^* = `r b0star`$    
$\b1^* = `r b1star`$    

```{r}
b0hat <- b0star / (1 - rhohat)
b1hat <- b1star
data.frame(b0hat, b0star, b1hat, b1star)
```

**Step 4:**

**SHOW R CODE AND OUTPUT**:

```{r}
dwtest(lmfit.CO1, alternative = "greater")

## grab the p-value from the Durbin-Watson test
objects(dwtest(lmfit.CO1, alternative = "greater"))
dw.CO1.stat <- dwtest(lmfit.CO1, alternative = "greater")$statistic
dw.CO1.pvalue <- dwtest(lmfit.CO1, alternative = "greater")$p.value
```
Durbin-Watson statistic is now closer to 2: $\; d = `r dw.CO1.stat`$

But, the pvalue is still pretty small: p-value = $`r dw.CO1.pvalue`$.   
I would like to make one more iteration.   

So, we go back to Step 2 with $\b0 = `r b0hat`$ and $\b1 = `r b1hat`$ calculated in Step 3 above.

**Step 2:** Estimate $\rho$

The residuals are now 
$\quad e_t = y_t - (\b0 + \b1 x_t) = y_t - (`r b0hat` + `r b1hat` x_t)$

```{r}
## Calculate rhohat "manually"
e <- with(expendData, Expenditure - (b0hat + b1hat * Stock))
et <- e[2:n]
et1 <- e[1:(n-1)]
rhohat1 <- sum(et * et1) / sum(e^2)
rhohat1
```

I found $\widehat{\rho} = `r rhohat1`$.

**Step 3:** Cochrane-Orcutt-transform $y$ and $x$ and re-estimate $\b0$ and $\b1$

```{r}
y <- with(expendData, Expenditure - rhohat1 * lag(Expenditure))
x <- with(expendData, Stock - rhohat1 * lag(Stock))

lmfit.CO2 <- lm(y ~ x)
tidy(lmfit.CO2)

b0star <- tidy(lmfit.CO2)$estimate[1]
b1star <- tidy(lmfit.CO2)$estimate[2]
```

$\b0^* = `r b0star`$    
$\b1^* = `r b1star`$    

```{r}
b0hat <- b0star / (1 - rhohat1)
b1hat <- b1star
data.frame(b0hat, b0star, b1hat, b1star)
```

**Step 4:**

**SHOW R CODE AND OUTPUT**:

```{r}
dwtest(lmfit.CO2, alternative = "greater")

## grab the p-value from the Durbin-Watson test
objects(dwtest(lmfit.CO2, alternative = "greater"))
dw.CO2.stat <- dwtest(lmfit.CO2, alternative = "greater")$statistic
dw.CO2.pvalue <- dwtest(lmfit.CO2, alternative = "greater")$p.value
```

Durbin-Watson statistic is even closer to 2: $\; d = `r dw.CO2.stat`$

But, the pvalue is still pretty small: p-value = $`r dw.CO2.pvalue`$.   
I would like to make one more iteration.   

So, we go back to Step 2 with $\b0 = `r b0hat`$ and $\b1 = `r b1hat`$ calculated in Step 3 above.

**Step 2:** Estimate $\rho$

The residuals are now 
$\quad e_t = y_t - (\b0 + \b1 x_t) = y_t - (`r b0hat` + `r b1hat` x_t)$

```{r}
## Calculate rhohat "manually"
e <- with(expendData, Expenditure - (b0hat + b1hat * Stock))
et <- e[2:n]
et1 <- e[1:(n-1)]
rhohat2 <- sum(et * et1) / sum(e^2)
rhohat2
```

I found $\widehat{\rho} = `r rhohat2`$.

**Step 3:** Cochrane-Orcutt-transform $y$ and $x$ and re-estimate $\b0$ and $\b1$

```{r}
y <- with(expendData, Expenditure - rhohat2 * lag(Expenditure))
x <- with(expendData, Stock - rhohat2 * lag(Stock))

lmfit.CO3 <- lm(y ~ x)
tidy(lmfit.CO3)

b0star <- tidy(lmfit.CO3)$estimate[1]
b1star <- tidy(lmfit.CO3)$estimate[2]
```

$\b0^* = `r b0star`$    
$\b1^* = `r b1star`$    

```{r}
b0hat <- b0star / (1 - rhohat2)
b1hat <- b1star
data.frame(b0hat, b0star, b1hat, b1star)
```

**Step 4:**

**SHOW R CODE AND OUTPUT**:

```{r}
dwtest(lmfit.CO3, alternative = "greater")

## grab the p-value from the Durbin-Watson test
objects(dwtest(lmfit.CO3, alternative = "greater"))
dw.CO3.stat <- dwtest(lmfit.CO3, alternative = "greater")$statistic
dw.CO3.pvalue <- dwtest(lmfit.CO3, alternative = "greater")$p.value
```

Durbin-Watson statistic is even closer to 2: $\; d = `r dw.CO3.stat`$

p-value = $`r dw.CO2.pvalue`$.   

We fail to reject the null of no autocorrelation, so we stop here after three iterations.

**NOTE:** I'll show you below how we can let the computer do iterations automatically.

$\b0 = \displaystyle{\frac{\b0^*}{1-\widehat{\rho}}} = `r b0hat`$

$\b1 = \b1^* = `r b1hat`$

Thus, our estimated regression equation is

$y_t = `r b0hat` + `r b1hat` x_t$


#### How does C-O compare to OLS regression?

Regression coefficient estimates
```{r}
tidy(lmfit)
tidy(lmfit.CO3)
```

The slopes are similar and $\se(\b1)$ is larger after transformation.

The intercept and it's SE do not get appropriately transformed in the standard regression output above.   
We can transform the output ourselves:

```{r}
se.b0hat <- tidy(lmfit.CO3)$std.error[1] / (1 - rhohat2)
b0hat <- b0star / (1 - rhohat2)
data.frame(b0hat, se.b0hat)
```

$\b0 = \displaystyle{\frac{\b0^*}{1-\widehat{\rho}}} = `r b0hat`$

$\se(\b0) \approx \displaystyle{\frac{\se(\b0^*)}{1-\widehat{\rho}}} = `r se.b0hat`$

t-statistic $\; = \displaystyle{\frac{\b0 - 0}{\se(\b0)}} = `r b0hat / se.b0hat`$ as above.


Regression goodness-of-fit statistics

```{r}
glance(lmfit)[-c(7:10)]
glance(lmfit.CO3)[-c(7:10)]
```

Durbin-Watson statistic  

**SHOW R CODE AND OUTPUT**

```{r}
durbinWatsonTest(lmfit, alternative = "positive")
```

**SHOW R CODE AND OUTPUT**

```{r}
durbinWatsonTest(lmfit.CO3, alternative = "positive")
```

After the model is corrected for autocorrelation by Cochrane-Orcutt transformaion:

  + DW statistic much closer to 2 (weaker evidence of autocorrelation)
  
  + Autocorrelation is much smaller.
  
  + Slope coefficients are similar, but $\se(\b1)$ is more than 3 times larger!
  
  + $R^2$ is smaller.
  
**WARNING:**  
Ignoring positive autocorrelation often leads to   

  + underestimated errors,   
  + too many false positives,   
  + too narrow and invalid confidence intervals.


Index plots of residuals:

```{r, fig.width=10}
## index plot of regular OLS residuals
p1 <- gf_point(r ~ index, data = expendData, size = 2) %>% 
  gf_line() %>%
  gf_lims(y = c(-2.1, 2.1)) %>%
  gf_labs(x = "Quarter", y = "OLS residuals") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")

## index plot of Cochrane-Orcutt residuals after three iterations
p2 <- gf_point(rstandard(lmfit.CO3) ~ 2:n, size = 2) %>%
  gf_line() %>%
  gf_lims(y = c(-2.1, 2.1)) %>%
  gf_labs(x = "Quarter", y = "Cochrane-Orcutt residuals (3 iterations)") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")  

grid.arrange(p1, p2, ncol=2)
```

#### Using R to iterate for Cochrane-Orcutt transformation

`library(orcutt)`   
`cochrane.orcutt(lmfit)`

```{r}
library(orcutt, warn.conflicts = FALSE)
cochrane.orcutt.fit <- cochrane.orcutt(lmfit)
```

OLS
```{r}
tidy(lmfit)
```

Several iterations of Cochrane-Orcutt

```{r}
tidy(cochrane.orcutt.fit)
```

The slope coefficients are similar, but $\se(\b1)$ is almost 4 times larger!
 
Comparison of goodness-of-fit statistics: 

```{r}
glance(lmfit)[-c(7:10)]
glance(cochrane.orcutt.fit)
```

Index plots of residuals:

```{r, fig.height=4, fig.width=10}
## index plot of regular OLS residuals
p1 <- gf_point(r ~ index, data = expendData, size = 2) %>% 
  gf_line() %>%
  gf_lims(y = c(-2.1, 2.1)) %>%
  gf_labs(x = "Quarter", y = "OLS residuals") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")

## index plot of Cochrane-Orcutt residuals after three iterations

p2 <- gf_point(rstandard(lmfit.CO3) ~ 2:n, size = 2) %>%
  gf_line() %>%
  gf_lims(y = c(-2.1, 2.1)) %>%
  gf_labs(x = "Quarter", y = "Cochrane-Orcutt residuals (3 iterations)") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")  

## index plot of Cochrane-Orcutt "standardized" residuals after several iterations
rhohat.orcutt <- cochrane.orcutt.fit$rho
y <- with(expendData, Expenditure - rhohat.orcutt * lag(Expenditure))[-1]
x <- with(expendData, Stock - rhohat.orcutt * lag(Stock))[-1]
b0hat <- cochrane.orcutt.fit$coeff[1] * (1 - rhohat.orcutt)
b1hat <- cochrane.orcutt.fit$coeff[2]
e <- y - b0hat - b1hat * x 
## simple standardization of residuals (z-score)
z <- (e - mean(e)) / sd(e)

p3 <- gf_point(z ~ 2:n, size = 2) %>%
  gf_line() %>%
  gf_lims(y = c(-2.1, 2.1)) %>%
  gf_labs(x = "Quarter", y = paste("Cochrane-Orcutt residuals (", cochrane.orcutt.fit$number.interaction, " iterations)", sep="")) %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")  

grid.arrange(p1, p2, p3, ncol=3)
```

### Omitted Variables

### Correlated Errors due to Omitted Variables

Correlated errors can occur when a variable **not** in the model is related   
to the response and related to the time or order in which the data were collected

#### Example: Relationship between housing starts and population of potential buyers

From textbook Table 8.4:    
A midwestern construction industry association wants to know the relationship between housing starts and population growth. They obtain annual data on regional housing starts and try to relate it to number of potential home buyers in the region.

```{r}
# houseData <- read.delim("http://statistics.uchicago.edu/~collins/data/RABE5/P219.txt")
houseData <- read.delim("./data/P219.txt")
names(houseData) <- c("houseStarts", "popnSize", "mortgageMoney")

head(houseData)
```

```{r fig.height=4}
gf_point(houseStarts ~ popnSize, data = houseData, size = 2) %>% 
  gf_labs(y = "Housing Starts") %>%
  gf_labs(x="Population")
```  


```{r}
lmfit <- lm(houseStarts ~ popnSize, data=houseData)
#plot(fit$residuals)
tidy(lmfit)
glance(lmfit)[-c(7:10)]
```


```{r fig.height=4}
## add residuals, standardized residuals, and an index variable to dataset
n <- nrow(houseData) ## sample size
houseData <- mutate(houseData, e = residuals(lmfit), r = rstandard(lmfit), index = 1:n)
## index plot
p1 <- gf_point(r ~ index, data = houseData, size = 2) %>% 
  gf_line() %>%
    gf_labs(title='Fit with population only') %>%
  gf_labs(y = "standardized residuals") %>%
  gf_lims(y = c(min(min(~r, data = houseData), -2.1), max(max(~r, data = houseData), 2.1))) %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")
p1
```  

This residual plot shows clear evidence of autocorrelation! 

We can confirm this using a runs test:


```{r}
houseData <- mutate(houseData, positive = 1 * (r > 0), sign = factor(positive, labels = c("-", "+")))
```

Signs of the residuals in time order:

`r houseData$sign`

```{r}
runs.test(factor(sign(houseData$r)))
runs.fit <- runs.test(factor(sign(houseData$r)))
runs.test.statistic <- runs.fit$statistic
runs.test.pvalue <- runs.fit$p.value
```
p-value = $`r round(runs.test.pvalue, 4)`$.       
This indicates we should reject the null hypothesis of no autocorrelation.

We can also calculate the Durbin-Watson statistic:

```{r}
durbinWatsonTest(lmfit)
```

The residual plot, runs test, and Durbin-Watson statistic   
all indicate the presence of autocorrelation.


The researchers further conjectured that the availability of mortgage money for the region might also be a predictor for housing starts, and decided to include it in the model.

Available mortgage money also varies by time.

```{r}
fit.mortgage <- lm(houseStarts ~ popnSize + mortgageMoney, data = houseData)
tidy(fit.mortgage)
glance(fit.mortgage)[-c(7:10)]
```

```{r fig.height=4}
# add residuals and standardized residuals to dataset
n <- nrow(houseData) # sample size
houseData <- mutate(houseData, e.mortgage = residuals(fit.mortgage), r.mortgage = rstandard(fit.mortgage))
## index plot
p2 <- gf_point(r.mortgage ~ index, data = houseData, size = 2) %>% 
  gf_line() %>%
  gf_labs(title='Fit with population and mortgage money') %>%
  gf_labs(y = "standardized residuals") %>%
  gf_lims(y = c(min(min(~ r.mortgage, data = houseData), -2.1), max(max(~ r.mortgage, data = houseData), 2.1))) %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")
p2
```  

This looks as though the autocorrelation has been fixed.   

Let's check first with a runs test:

```{r}
houseData <- mutate(houseData, positive.mortgage = 1 * (r.mortgage > 0), sign.mortgage = factor(positive.mortgage, labels = c("-", "+")))
```

Signs of the residuals in time order:

`r houseData$sign.mortgage`

```{r}
runs.fit <- runs.test(factor(sign(houseData$r.mortgage)))
runs.test.statistic <- runs.fit$statistic
runs.test.pvalue <- runs.fit$p.value
```

p-value = $`r round(runs.test.pvalue, 4)`$.   
This indicates we should not reject the null hypothesis of no autocorrelation.

We can also calculate the Durbin-Watson statistic:

```{r}
durbinWatsonTest(fit.mortgage)
```

```{r fig.width=10}
grid.arrange(p1, p2, ncol=2)
```

Can we understand what happened? 

First, make an added variable plot     
to see how the mortgage money affects housing starts,    
after controlling both mortgage money and housing for population:

```{r}
avPlots(fit.mortgage, terms = ~ mortgageMoney, pch=16, xlab="Residuals of mortgage money after regressing on population", ylab="Resid. of housing after regressing on population")
```

We see that there is a positive relationship:   
more mortgage money means more housing, all else being equal.

We can look at how residuals of mortgage money after regressing on population change over time:

```{r, fig.width=10}
fit.m <- lm(mortgageMoney ~ popnSize, data=houseData) # Regress mortgage money on residuals
houseData <- mutate(houseData, m.on.h.resids = fit.m$residuals)
p3 <- gf_point(m.on.h.resids ~ index, data = houseData, size = 2) %>% 
  gf_line() %>%
  gf_labs(y = "Residuals of mortgage money after regressing on population")%>%
 # gf_lims(y = c(min(min(~r.mortgage, data = house), -2.1), max(max(~r.mortgage, data = house), 2.1))) %>%
  gf_hline(yintercept = c(0), linetype = "dashed")
grid.arrange(p3 , p1, ncol=2)
```

It seems that in the original fit, with mortgage money omitted,     
the model underestimated housing starts at those times when     
mortgage money was higher than would have been predicted from population alone.

### Seasonality

### Autocorrelation: One More Example

The analytic methods examined thus far (DW statistic, runs test) deal with     
first order autocorrelation,    
meaning serial correlation of adjacent values. 

There are cases where autocorrelation may indeed be present,    
but with other, perhaps more complex, patterns of time-dependence.

From textbook Table 8.8,    
we consider the exports of a company that produces and markets ski equipment in the US.

We hope to summarize a relationship of quarterly sales to some leading economic indicators,   
like personal disposable income (PDI) in billions of current dollars. 

The initial model is

$$ S_t = \B0 + \B1 PDI_t + \epsilon_t $$

where $S_t$ is the ski Sales in period $t$.     

The sales period naturally follows time order, and cycles through seasons.

```{r}
# skiData <- read.delim("http://statistics.uchicago.edu/~collins/data/RABE5/P224.txt")
skiData <- read.delim("./data/P224.txt")
head(select(skiData, Sales, PDI), 8)
```

Let's look at the scatterplot of Sales vs PDI:

```{r}
gf_point(Sales ~ PDI, data = skiData) %>% 
  gf_labs(title = 'Ski Sales') %>%
  gf_labs(y = "Sales", x = "PDI")
```

This looks pretty good.  A clear linear relationship.    



Since the data were collected over time,    
we should look for autocollinearity.


Start with an index plot of the `Sales` variable:

```{r}
n <- nrow(skiData)
skiData <- mutate(skiData, index = 1:n)
gf_point(Sales ~ index, data = skiData, size = 2) %>% 
  gf_line(color='black') %>%
  gf_labs(title='Ski Sales') %>%
  gf_labs(y = "Sales", x="Index")
``` 

The Sales values go up-and-down and up-and-down over time.

We go ahead and fit the simple linear regression model for now
to check further for autocorrelation.

```{r}
ski.fit <- lm(Sales ~ PDI, data = skiData)
skiData <- mutate(skiData, r = rstandard(ski.fit))

tidy(ski.fit)
glance(ski.fit)[-c(7:10)]
```

Look at an index plot of the residuals:

```{r}
gf_point(r ~ index , data = skiData) %>% 
  gf_line(color = "black") %>%
  gf_labs(title = 'Ski Sales') %>%
  gf_labs(y = "r", x = "index") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")
```

This makes clear the up-and-down and up-and-down pattern.

What do we see from the Durbin-Watson statistic?

```{r}
durbinWatsonTest(ski.fit)
```

What?   
The autocorrelation is basically zero.   
The Durbin-Watson statistic is near 2.   

There is no evidence of **first-order** autocorrelation.
There must be something else that explains the pattern.

Were any other variables collected over time?

```{r}
skiData <- mutate(skiData, theSeason = factor(Season, labels = c('Spring/Summer', 'Fall/Winter')))

head(skiData, 8)
```

What is the pattern in the data by season?
```{r}
gf_point(Sales ~ PDI, color = ~ theSeason, data = skiData) %>% 
  gf_labs(title = 'Ski Sales') %>%
  gf_labs(y = "Sales", x = "PDI") %>%
  gf_theme(legend.position = "bottom")
```

Perhaps two parallel lines are needed (same slope, different intercept by season).

```{r}
n <- nrow(skiData)
skiData <- mutate(skiData, index = 1:n)
gf_point(Sales ~ index, color = ~ theSeason, data = skiData, size = 2) %>% 
  gf_line(color='black') %>%
  gf_labs(title='Ski Sales') %>%
  gf_labs(y = "Sales", x="Index")
``` 

It's clear Sales are higher in Fall/Winter versus lower Sales in Spring/Summer    
and this information not accounted for in our model.


Here is the residual plot with seasons indicated:

```{r}
gf_point(r ~ index , color = ~ theSeason, data = skiData) %>% 
  gf_line(color = "black") %>%
  gf_labs(title = 'Ski Sales') %>%
  gf_labs(y = "r", x = "index") %>%
  gf_theme(legend.position = "bottom")%>%
#  gf_labs(legend=c('Fall and Winter','Spring and summer')) %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")
```

The higher Sales in Fall/Winter lead to consistently higher residuals in these seasons.    

This is a form of autocorrelation:    
high Sales in Fall are paired with high Sales in the next Winter,     
while low Sales in Spring are adjacent to low Sales in Summer.

Reminder: This autocorrelation was **NOT** detected by the Durbin-Watson test?

```{r}
durbinWatsonTest(ski.fit)
```

The positive correlation within adjacent Fall/Winter and Spring/Summer quarters     
is balanced by a **negative** correlation between adjacent Winter/Spring and Summer/Fall seasons,   
leading to a $\widehat\rho$ that is close to zero.

Also, the linear relationship between `Sales` and `PDI`   
is different by season (same slope, different intercept).

Let's try a different approach.     
We'll include `Season` as an indicator variable, with value 1 in the Fall/Winter.

```{r fig.width=10}
ski.fit2 <- lm(Sales ~ PDI + Season, data = skiData)
skiData <- mutate(skiData, r2 = rstandard(ski.fit2))

tidy(ski.fit2)
glance(ski.fit2)[-c(7:10)]

p1 <- gf_point(r2 ~ index, color = ~theSeason, data = skiData) %>% 
  gf_labs(title = 'Ski Sales') %>%
  gf_labs(y = "r", x="index") %>%
  gf_theme(legend.position="bottom") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")

p2 <- gf_point(r2 ~ index, color = ~theSeason, data = skiData) %>% 
  gf_line(color = "black") %>%
  gf_labs(title = 'Ski Sales') %>%
  gf_labs(y = "r", x="index") %>%
  gf_theme(legend.position="bottom") %>%
  gf_hline(yintercept = c(-2, 0, 2), linetype = "dashed")

grid.arrange(p1, p2, ncol=2)
```


This looks much better. 

I'm a little concerned that there may be some other autocorrelation going on:   
that residuals are consistently higher in observations 15-30 or so. 

But this autocorrelation, if it exists, would be first-order,     
and I should be able to capture it with the Durbin-Watson test:

```{r}
durbinWatsonTest(ski.fit2)
```

OK, this looks good!  No evidence of autocorrelation.

Last, I will consider including an interaction term,    
in case the slope differs for Fall/Winter versus Spring/Summer:

```{r}
ski.fit3 <- lm(Sales ~ PDI * Season, data=skiData)

tidy(ski.fit3)
```

The interaction term isn't significant, so I can stick with the simpler model.

