---
title: "Hypothesis testing in the CER Model"
author: "Eric Zivot"
date: "July 30, 2015"
output: slidy_presentation
---

## Set options and load packages

```{r}
options(digits=3, width=70)
# install IntroCompFinR package from R-forge
# use install.packages("IntroCompFinR", repos="http://R-Forge.R-project.org")
library(IntroCompFinR)
library(mvtnorm)
library(PerformanceAnalytics)
library(tseries)
library(zoo)
Sys.setenv(TZ="UTC")
```

## Get data from package IntroCompfinR

```{r}
# get data from IntroCompFin package
data(msftDailyPrices, sp500DailyPrices, sbuxDailyPrices)
msftPrices = to.monthly(msftDailyPrices, OHLC=FALSE)
sp500Prices = to.monthly(sp500DailyPrices, OHLC=FALSE)
sbuxPrices = to.monthly(sbuxDailyPrices, OHLC=FALSE)
# set sample to match book chapter
smpl = "1998-01::2012-05"
msftPrices = msftPrices[smpl]
sp500Prices = sp500Prices[smpl]
sbuxPrices = sbuxPrices[smpl]
# calculate returns
msftRetS = Return.calculate(msftPrices, method="simple")
sp500RetS = Return.calculate(sp500Prices, method="simple")
sbuxRetS = Return.calculate(sbuxPrices, method="simple")
msftRetS = msftRetS[-1]
sp500RetS = sp500RetS[-1]
sbuxRetS = sbuxRetS[-1]
msftRetC = log(1 + msftRetS)
sp500RetC = log(1 + sp500RetS)
sbuxRetC = log(1 + sbuxRetS)
# merged data set
cerRetS = merge(msftRetS, sbuxRetS, sp500RetS)
cerRetC = merge(msftRetC, sbuxRetC, sp500RetC)
colnames(cerRetS) = colnames(cerRetC) = c("MSFT", "SBUX", "SP500")
head(cerRetC, n=3)

```

## Testing $H_0: \mu = 0$ vs. $H_1: \mu \ne 0$

Brute force calculation of t-statistic

1. Get estimates of $\mu$

```{r}
n.obs = nrow(cerRetC)
muhat.vals = apply(cerRetC, 2, mean)
muhat.vals
```

2. Calculate SEs

```{r}
# calculate standard errors
sigmahat.vals = apply(cerRetC, 2, sd)
se.muhat = sigmahat.vals/sqrt(n.obs)
se.muhat
```

3. Calculate t-statistics

```{r}
t.stats = muhat.vals/se.muhat
abs(t.stats)
```

4. Apply decision rule:  all t-statistics have absolute value less than 2, so do not reject null hypothesis at 5% level of significance.

## P-value calculation

P-value = significance level at which null hypothesis is just rejected

```{r}
# calculate 2-sided p-value
2*(1-pnorm(abs(t.stats)))
```

All p-values are less than 5%, so do not reject null hypothesis at 5% level.

## Use 95% confidence interval to perform hypothesis test

$H_0: \mu = 0$ vs. $H_1: \mu \ne 0$ is rejected at the 5% level if $\mu = 0$ is not in the 95% confidence interval.

```{r}
# calculate 95% confidence interval
lower = muhat.vals - 2*se.muhat
upper = muhat.vals + 2*se.muhat
cbind(lower, upper)
```

Here, $\mu = 0$ is in all of the 95% confidence intervals so we do not reject $H_0: \mu = 0$ vs. $H_1: \mu \ne 0$ at the 5% level for any asset.

## Test for normal distribution

Graphical descriptive statistics suggest that the distribution of MSFT returns have fatter tails than the normal distribution.

```{r}
fourPanelPlot(cerRetC[, "MSFT"])
```

However, graphical diagnostics are not a formal statistical test.

## Jarque-Bera test for normality

Excess kurtosis value indicates non-normality.

```{r}
msft.skew = skewness(cerRetC[, "MSFT"])
msft.ekurt = kurtosis(cerRetC[, "MSFT"])
msft.skew
msft.ekurt
```

JB statistic confirms non-mormality

```{r}
JB = n.obs*(msft.skew^2 + 0.25*msft.ekurt^2)/6
JB
```

JB > 6 so reject $H_0: r_t \sim N(\mu, \sigma^2)$ at 5% level.

## R package **tseries** function <code>jarque.bera.test()</code>

```{r}
library(tseries)
jarque.bera.test(cerRetC[, "MSFT"])
```

Here, "X-squared" denotes the JB statistic. The small p-value indicates rejection of null hypothesis at any reasonable significance level.

## Test for no serial correlation

The R function <code>acf()</code> computes the sample autocorrelations $\hat{\rho_j}$, plots them and shows the critical values for rejecting the null hypothesis $H_0: \rho_j = 0$

```{r}
acf(cerRetC[, "MSFT"], lwd=2)
```

Here, $\hat{rho_1}$ and $\hat{rho_3}$ look to be statistically differnet from zero at the 5% level.

## Compute rolling means

Use the **zoo** function <code>rollapply()</code> to compute rolling mean estimates.

```{r}
# compute rolling means over 24 month windows
roll.muhat = rollapply(cerRetC[, "MSFT"], width=24,
                       FUN=mean, align="right")
class(roll.muhat)
# first 24 values are NA
roll.muhat[23:25]
```

## Plot rolling means with returns

```{r}
plot.zoo(merge(roll.muhat,cerRetC[, "MSFT"]), plot.type="single",
     main="24 month rolling means for MSFT",ylab="returns",
     lwd=c(2,2), col=c("blue","orange"))
abline(h=0)
grid()
legend(x="bottomleft",legend=c("Rolling mean","Monthly returns"),
       lwd=c(2,2), col=c("blue","orange"))
```

* Rolling means start out positive during dot-com boom
* Drop negative during dot-com bust
* close to zero during run-up to financial crisis
* Dip negative during financial crisis
* become positive after financial crisis

## 24-month rolling means for S&P500 index 

```{r}
roll.muhat = rollapply(cerRetC[, "SP500"], width=24,
                       FUN=mean, align="right")
plot.zoo(merge(roll.muhat,cerRetC[, "SP500"]), plot.type="single",
     main="24-month rolling means for SP500",ylab="returns",
     lwd=c(2,2), col=c("blue","orange"))
abline(h=0)
grid()
legend(x="bottomleft",legend=c("Rolling mean","Monthly returns"),
       lwd=c(2,2), col=c("blue","orange"))
```

## Compute rolling volatility ($\sigma$) estimates

Use the **zoo** function <code>rollapply()</code> to compute rolling $\sigma$ estimates.

```{r}
# compute rolling sd over 24 month windows
roll.sigmahat = rollapply(cerRetC[, "MSFT"], width=24,
                       FUN=sd, align="right")
class(roll.sigmahat)
# first 24 values are NA
roll.sigmahat[23:25]
```

## Plot rolling volatility with returns

```{r}
plot.zoo(merge(roll.sigmahat,cerRetC[, "MSFT"]), plot.type="single",
     main="24-month rolling volatility for MSFT",ylab="returns",
     lwd=c(2,2), col=c("blue","orange"))
abline(h=0)
grid()
legend(x="bottomleft",legend=c("Rolling sd","Monthly returns"),
       lwd=c(2,2), col=c("blue","orange"))
```

* Rolling means start out large during dot-com boom/bust
* Decrease during run-up to financial crisis
* increase slightly during financial crisis


## 24-month rolling volatility for S&P500 index 

```{r}
roll.sigmahat = rollapply(cerRetC[, "SP500"], width=24,
                       FUN=sd, align="right")
plot.zoo(merge(roll.sigmahat,cerRetC[, "SP500"]), plot.type="single",
     main="24-month rolling volatility for SP500",ylab="returns",
     lwd=c(2,2), col=c("blue","orange"))
abline(h=0)
grid()
legend(x="bottomleft",legend=c("Rolling sd","Monthly returns"),
       lwd=c(2,2), col=c("blue","orange"))
```

## Compute rolling correlations

To compute rolling correlations, you need to write a function to return the pairwise correlation etween two returns.

```{r}
rhohat = function(x) {
	cor(x)[1,2]
}
```

Then use <code>rollapply()</code> with this function.

```{r}
roll.rhohat = rollapply(cerRetC[,c("MSFT","SP500")],
                       width=24,FUN=rhohat, by.column=FALSE,
                       align="right")
class(roll.rhohat)
roll.rhohat[23:25]
```

## Plot rolling correlations ($\rho_{ij}$)

```{r}
plot.zoo(roll.rhohat, main="24-month rolling correlations b/w MSFT and SP500",
     ylim=c(0,1), lwd=2, col="blue", ylab="rho.hat")
abline(h=0)   
grid()
```

* Correlations start out around 0.6
* Drop close to zero in mid 2005
* Increase during run up to financial crisis
* Reach 0.8 after financial crisis

## Summary of Hypothesis Testing in CER Model

* Hypothesis tests about $\mu$ are not very powerful because $SE(\hat{\mu})$ is typically very large
* Can ofter reject the null hypothesis that monthly returns are normally distributed
* Typically cannot reject the null hypothesis that monthly returns are serially uncorrelated
* Rolling window estimates indicate that $\mu$, $\sigma$ and $\rho_{ij}$ are non constant over time (covariance stationary assumption is suspect over long time spans)
