capm <- read.csv("D:/data/capm.csv")
#View(capm)
head(capm)
tail(capm)

ibm.r <- diff(log(capm[,5])) # IBM log returns
sp500.r <- diff(log(capm[,10])) # S&P 500 log returns

rf <- capm[,8] # Risk-free rate

ibm.ex <- ibm.r-rf[-1] # IBM excess return, drop first observation
sp500.ex <- sp500.r-rf[-1] # S&P excess return

#data <- as.data.frame(cbind(ibm.ex,sp500.ex))


m1 <- lm(ibm.ex~sp500.ex) # IBM capm
summary(m1)

names(m1)
param <- m1$coefficients
param

y.hat <- m1$fitted.values
x <- ts(cbind(ibm.ex,y.hat),start=c(1990,4),end=c(2004,7),frequency=12)
ts.plot(x,col=1:2,lwd=2,main="IBM excess returns and fitted values")

# Residuals
u <- ts(m1$residuals,start=c(1990,4),end=c(2004,7),frequency=12)
plot(u)
abline(h=0)

hist(u)
library(MASS)
truehist(u,main="Residuals")

### Tests

# H0: b=1; H1:b<>1
test.statistic <- (param[2]-1)/0.141260
test.statistic # H0 not rejected

# Using the car package
# See, the R_guide.pdf
# install.packages("car") # Install only once
library(car) # Load the package

# H0: b=1 using F test
linearHypothesis(m1,c("sp500.ex = 1"))
# The F-statistic of 2.10 with a corresponding p-value of 0.1486 implies that the null hypothesis of the
# CAPM beta of the IBM stock being 1 is convincingly not rejected and hence the estimated beta of 1.2050 is
# not significantly different from 1.

# H0: a=0 AND b=1
# F test
linearHypothesis(m1, c("(Intercept) = 0", "sp500.ex = 1"))
# The F-statistic of 1.4091 with a corresponding p-value of 0.2472 implies that the null hypothesis of
# the alpha = 0 AND beta = 1 is clearly not rejected.

### Gauss-Markov assumptions

# Assumption (2) var(u) = constant (homoskedasticity)

#install.packages("lmtest")
library(lmtest)

# Breusch-Pagan test
# H0: constant variance; homoskedasticity
bptest(formula(m1),studentize = FALSE)
# The null of constant variance of residuals is not 
# rejected with test statistic equal to 2.7801 and p-value 0.09544
bptest(formula(m1),studentize = TRUE)
# This result is robust to the alternative distributional
# assumption, since the test statistic and p-value also suggest
# that there is not a problem of heteroscedastic errors
# for the CAPM model

# Goldfeld-Quandt test

# library(lmtest)
gqtest(m1, order.by = ~sp500.ex, fraction = 34)
# H0: Homoskedasticity is present
# H1: Heteroskedasticity is present
# Since the p-value is not less than 0.05, we fail to reject 
# the null hypothesis. We do not have sufficient evidence to 
# say that heteroscedasticity is present in the regression model.

# Assumption (3): Test for autocorrelation in the residuals

# Durbin-Watson test
# library(car)
# H0 (null hypothesis): There is no correlation among the residuals.
# H1 (alternative hypothesis): The residuals are autocorrelated.
durbinWatsonTest(m1)
# From the output we can see that the test statistic is 2.128 and the 
# corresponding p-value is 0.43. Since this p-value is larger than 0.05, 
# we cannot reject the null hypothesis and conclude that the residuals 
# in this regression model are not autocorrelated, they are homoskedastic

# Breusch-Godfrey test
# library(lmtest)
# Test for high order autocorrelation
bgtest(m1, order = 6, order.by = NULL, type = c("Chisq"))
bgtest(m1, order = 6, order.by = NULL, type = c("F"))
# The p-values are 0.60 and 0.61 -> there is no autocorrelation
# in lags 1-6.

# Compute the Box-Pierce or Ljung-Box test statistic for examining 
# the null hypothesis of independence in a given time series. These 
# are sometimes known as 'portmanteau' tests. 
Box.test(u, lag = 6, type = c("Box-Pierce"), fitdf = 0)
Box.test(u, lag = 6, type = c("Ljung-Box"), fitdf = 0)
# The p-values are 0.65 and 0.63 -> there is no autocorrelation 
# in lags 1-6

# Test for normality
# Jarque-Bera test
library(tseries)
# The Jarque-Bera test is a goodness-of-fit test that determines 
# whether or not sample data have skewness and kurtosis that matches 
# a normal distribution.
jarque.bera.test(u)
# This tells us that the test statistic is 35.435 and the p-value of 
# the test is 0. In this case, we would fail not reject the null 
# hypothesis that the data is normally distributed.

library(DistributionUtils)

skewness(u, na.rm = FALSE)
kurtosis(u, na.rm = FALSE)

#y + b*x + u
#y + b*x + c*x^2 + u
#y + b*x + d*x^3 + u

# Ramsey's RESET test for functional form (null: linear model)
resettest(m1, power = 2:4, type = c("regressor"))
resettest(m1, power = 2:4, type = c("fitted"))
# Do not reject the null of a linear model

# Performs Quandt Likelihood Ratio Test for structural breaks
# with unknown breakdate. 
# install.packages("strucchange")
library(strucchange)

data <- list(ibm.ex=ibm.ex,sp500.ex=sp500.ex)
ibm1 <- efp(m1, type="OLS-CUSUM",data=data)
ibm2 <- efp(m1, type="ME",data=data,h=0.2)

#model <- lm(ibm.ex~sp500.ex,data=data)
qlr <- Fstats(ibm.ex~sp500.ex,data=data)
plot(qlr,alpha=0.05)
# If the F statistics cross their boundary, there is evidence for
# a structural change (at the level alpha=0.05).
# There is no evidence of a structural change in the CAPM

### Robust standard errors and t-statistics

vcov(m1) # covariance matrix of parameter estimates

vcovHC(m1) # heteroskedasticity consistent covariance matrix

vcovHAC(m1)# heteroskedasticity and autocorrelationconsistent covariance matrix (Newey-West)

summary(m1)

#library(lmtest)

coeftest(m1, vcov = vcovHC)

coeftest(m1, vcov = vcovHAC) # Newey-West

t(sapply(c("const", "HC0", "HC1", "HC2", "HC3", "HC4"),
         function(x) sqrt(diag(vcovHC(m1, type = x)))))