############# Autocorrelation Tests & Robust SE(Code #13) #############

####### Autocorrelation Tests & Robust SE

FX_da <- read.csv("http://www.bauer.uh.edu/rsusmel/4397/Stocks_FX_1973.csv",head=TRUE,sep=",")

x_date <- SFX_da$Date
x_cat <- SFX_da$CAT
x_dis <- SFX_da$DIS
x_ibm <- SFX_da$IBM
x_sp500 <- SFX_da$SP500
x_Mkt_RF <- SFX_da$Mkt_RF
x_SMB <- SFX_da$SMB
x_HML <- SFX_da$HML
x_RF <- SFX_da$RF

T <- length(x_cat)
T
 lr_cat <- log(x_cat[-1]/x_cat[-T])
 lr_dis <- log(x_dis[-1]/x_dis[-T])
 lr_ibm <- log(x_ibm[-1]/x_ibm[-T])
 lr_sp <- log(x_sp500[-1]/x_sp500[-T])
 date_1 <- x_date[-1]
 Mkt_RF <- x_Mkt_RF[-1]/100
 SMB <- x_SMB[-1]/100
 HML <- x_HML[-1]/100
 RF <- x_RF[-1]/100

ibm_x <- lr_ibm - RF
cat_x <- lr_cat - RF
dis_x <- lr_dis - RF

T <- length(ibm_x)

fit_ibm_ff3 <- lm (ibm_x ~ Mkt_RF + SMB + HML)
summary(fit_ibm_ff3)

fit_cat_ff3 <- lm (cat_x ~ Mkt_RF + SMB + HML)
summary(fit_cat_ff3)

fit_dis_ff3 <- lm (dis_x ~ Mkt_RF + SMB + HML)
summary(fit_dis_ff3)

## BG Test for IBM using bgtest from lmtest package (compared returns and excess returns)
library(lmtest)

bgtest(fit_ibm_ff3, order=12)			# IBM residuals with 12 lags
bgtest(fit_ibm_ff3, order=4)			# IBM residuals with 4 lags

bgtest(lr_ibm ~ Mkt_RF + SMB + HML, order=12)	# Using IBM returns (No excess returns) with 12 lagas
bgtest(lr_ibm ~ Mkt_RF + SMB + HML, order=4)	# Using IBM returns with 4 lags


## BG Test for CAT
bgtest(fit_cat_ff3, order=12)			# CAT residuals with 12 lags
bgtest(fit_cat_ff3, order=4)			# CAT residuals with 4 lags


## BG Test for DIS
bgtest(fit_dis_ff3, order=12)			# CAT residuals with 12 lags
bgtest(fit_dis_ff3, order=4)			# CAT residuals with 4 lags


### DW for IBM
RSS <- sum(e_ibm^2)
T <- length(e_ibm)
DW <- sum((e_ibm[1:(T-1)]-e_ibm[2:T])^2)/RSS		# DW stat
DW
2*(1 - cor(e_ibm[1:(T-1)],e_ibm[2:T]))			# approximate DW stat

# Same DW test done with dwtest from lmtest package
dwtest(fit_ibm_ff3)


### DW Test for CAT
dwtest(fit_cat_ff3)


### DW Test for DIS
dwtest(fit_dis_ff3)


### Exchange Rates Application
SF_da <- read.csv("http://www.bauer.uh.edu/rsusmel/4397/SpFor_prices.csv", head=TRUE, sep=",")
summary(SF_da)
x_date <- SF_da$Date
x_S <- SF_da$GBPSP 
cpi_uk <- SF_da$UK_CPI
cpi_us <- SF_da$US_CPI
T <- length(x_S)
i_us3 <- SF_da$Dep_USD3M
i_uk3 <- SF_da$Dep_UKP3M

lr_usdgbp <- log(x_S[-1]/x_S[-T])
I_us <- log(cpi_us[-1]/cpi_us[-T])
I_uk <- log(cpi_uk[-1]/cpi_uk[-T])
inf_dif <- (I_us - I_uk)
int_dif <- (i_us3[-1] - i_uk3[-1])/100
	
y <- lr_usdgbp
fit_gbp <- lm (y ~ inf_dif + int_dif)
summary(fit_gbp)


### BG test done with bgtest from lmtest package 
bgtest(fit_gbp, order=12)
bgtest(fit_gbp, order=4)


### DW Test for USD/GBP
dwtest(fit_gbp)


#### Q & LB Tests
y <- e_ibm
T <- length(y)
r_sum <- 0
lb_sum <- 0
p_lag <- 12
a <- 1

while (a <= p_lag) {
  za <- cor(y[(p_lag+1):T],y[a:(T-p_lag+a-1)])^2	 
  r_sum <- r_sum + za		#sum cor(e[(lag+1):T],e[a:(T-p_lag+a-1)])^2
  lb_sum <- lb_sum + za/(T-a) 
  a <- a + 1
}

qchisq(.95,p_lag)
Q <- T*r_sum
Q
1-pchisq(Q,p_lag)

LB <- T*(T-2)*lb_sum
LB
1- pchisq(LB,p_lag)


### Same Q & LB Tests done with Box.test from lmtest package
## Q for IBM
Box.test(e_ibm, lag = 12, type="Box-Pierce")

## LB for IBM
Box.test(e_ibm, lag = 12, type="Ljung-Box")


## Q & LB for CAT
e_cat <- fit_cat_ff3$residuals		# Extract OLS residuals from fit_cat_ff3
Box.test(e_cat, lag = 12, type="Box-Pierce")

Box.test(e_cat, lag = 12, type="Ljung-Box")


## Q & LB for DIS
e_dis <- fit_dis_ff3$residuals
Box.test(e_dis, lag = 12, type="Box-Pierce")

Box.test(e_dis, lag = 12, type="Ljung-Box")


## Q & LB Tests for USD/GBP
e_gbp <- fit_gbp$residuals
Box.test(e_gbp, lag = 12, type="Box-Pierce")

Box.test(e_gbp, lag = 12, type="Ljung-Box")


### Robust Q Tests, using vrtest
library(vrtest)
Auto.Q(e_ibm, 12)		#Maximum potential lag = 12

Auto.Q(e_dis, 12)		#Maximum potential lag = 12

Auto.Q(e_gbp, 12)		#Maximum potential lag = 12


### Time-varying Volatility: Q & LB Tests for squared returns done with Box.test
## Q & LB for IBM
e_ibm2 <- e_ibm^2
Box.test(e_ibm2, lag = 12, type="Box-Pierce")

Box.test(e_ibm2, lag = 12, type="Ljung-Box")


## Q & LB for CAT
e_cat2 <- e_cat^2
Box.test(e_cat2, lag = 12, type="Box-Pierce")

Box.test(e_cat2, lag = 12, type="Ljung-Box")

## Q & LB for DIS
e_dis2 <- e_dis^2
Box.test(e_dis2, lag = 12, type="Box-Pierce")

Box.test(e_dis2, lag = 12, type="Ljung-Box")

## Q & LB for USD/GBP
e_gbp2 <- e_gbp^2
Box.test(e_gbp2, lag = 12, type="Box-Pierce")

Box.test(e_gbp2, lag = 12, type="Ljung-Box")


#######  OLS with Robust SE (White & NW SEs)
b_ibm <- fit_ibm_ff3$coefficients		# Extract OLS coeff’s from fit_ibm_ff3
SE_OLS <- sqrt(diag(vcov(fit_ibm_ff3)))		# Extract OLS SE from fit_ibm_ff3
t_OLS <- b_ibm/SE_OLS

SE_OLS
t_OLS

### White SE
library(sandwich)
White <- vcovHC(fit_ibm_ff3, type = "HC0")		# HC0 = Standard White SE
SE_White <- sqrt(diag(White))				# White SE HC0
t_White <- b_ibm/SE_White

SE_White
t_White 

White3 <- vcovHC(fit_ibm_ff3, type = "HC3")		# HC3 refinement of White SE
SE_White3 <- sqrt(diag(White3))				# White SE HC0
t_White3 <- b_ibm/SE_White3

SE_White3
t_White3 

### NW SE
NW <- NeweyWest(fit_ibm_ff3, lag = 4, prewhite = FALSE)	# with 4 lags
SE_NW <- diag(sqrt(abs(NW)))
t_NW <- b_ibm/SE_NW

SE_NW
t_NW

NW <- NeweyWest(fit_ibm_ff3, lag = 8, prewhite = FALSE)
SE_NW <- diag(sqrt(abs(NW)))
t_NW <- b_ibm/SE_NW

SE_NW
t_NW


######## GLS & FGLS

fit_dis_ff3 <- lm (dis_x ~ Mkt_RF + SMB + HML)
summary(fit_dis_ff3)

#### GLS for DIS with Sigma2 = Mkt_RF^2
T <- length(dis_x)
Mkt_RF2 <- Mkt_RF^2					# (A3')
y_w <- dis_x/sqrt(Mkt_RF2)				# transformed y = y*
x0 <- matrix(1,T,1)
xx_w <- cbind(x0, Mkt_RF, SMB, HML)/sqrt(Mkt_RF2)	# transformed X = X*

fit_dis_wls <- lm(y_w ~ xx_w - 1)			# GLS (take automatic constant of lm with -1)
summary(fit_dis_wls)