############# Heteroscedasticity and Autocorrelation (Code #12) #############

####### Heteroscedasticity and Autocorrelation

SFX_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_RMW <- SFX_da$RMW
x_CMA <- SFX_da$CMA
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
 RMW <- x_RMW[-1]/100
 CMA <- x_CMA[-1]/100
 RF <- x_RF[-1]/100

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

y <- ibm_x; 
x1 <-  Mkt_RF
x2 <-  SMB
x3 <- HML
T <- length(x1)
x0 <- matrix(1,T,1)
x <- cbind(x0,x1,x2,x3)
k <- ncol(x)

T <- length(y)

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


### Model Selection - Best Subset
library(olsrr)
ff_step_data <- data.frame(Mkt_RF, SMB, HML) 
fit_ibm_ff3 <- lm(ibm_x ~ ., data = ff_step_data)	# default p-value (penter) is 0.3 ols_step_best_subset(fit_ibm_ff3 , details = TRUE)	#  long final output

#### GQ Test

library(lmtest)
gqtest(fit_ibm_ff3 , fraction = .20)

gqtest(cat_x ~ Mkt_RF + SMB + HML, fraction = .20)


#### LM-BP- test

## IBM -- Checking for Mkt_RF^2 & SMB^2
e_ibm <- fit_ibm_ff3$residuals
e2 <- e_ibm^2
Mkt_RF_2 <- Mkt_RF^2
SMB_2 <- SMB^2
fit_ibm_BP <- lm (e2 ~ Mkt_RF_2 + SMB_2)
Re_e2 <- summary(fit_ibm_BP)$r.squared
Re_e2 
LM_BP_test <- Re_e2 * T
LM_BP_test

 p_val <- 1 - pchisq(LM_BP_test, df = 2) 	 # p-value of LM_test 
 p_val

## IBM -- Checking for Mkt_RF_2 + SMB_2 + HML^2
HML_2 <- HML^2
fit_BP <- lm (e2 ~ Mkt_RF_2 + SMB_2 + HML_2)
Re_e2 <- summary(fit_BP)$r.squared
Re_e2 
LM_BP_test2 <- Re_e2 * T
LM_BP_test2

p_val <- 1 - pchisq(LM_BP_test2, df = 2) 	 # p-value of LM_test 
p_val


## IBM -- Using bptest only checks for variables in model (Mkt_RF, SMB, HML)
bptest(ibm_x ~ Mkt_RF + SMB + HML)


## CAT-- Checking for Mkt_RF^2 & SMB^2 & HML^2
fit_cat_ff3 <- lm (cat_x ~ Mkt_RF + SMB + HML)
summary(fit_cat_ff3)

e_cat <- fit_cat_ff3$residuals
e_cat2 <- e_cat^2
fit_cat_BP <- lm(e_cat2 ~ Mkt_RF_2 + SMB_2 + HML_2)
Re_e2 <- summary(fit_cat_BP)$r.squared
LM_BP_test <- Re_e2 * T
LM_BP_test

 p_val <- 1 - pchisq(LM_BP_test, df = 1) 	 # p-value of LM_test 
 p_val

## CAT -- Using bptest only checks for variables in model (Mkt_RF, SMB, HML)
bptest(fit_cat_ff3)



#### LM-White- test

## IBM 
Mkt_HML <- Mkt_RF*HML
Mkt_SMB <- Mkt_RF*SMB
SMB_HML <- SMB*HML
xx2 <- cbind(Mkt_RF_2, SMB_2, HML_2, Mkt_HML,Mkt_SMB, SMB_HML)
fit_ibm_W <- lm (e2 ~ Mkt_RF + SMB + HML + xx2)
R2_e2 <- summary(fit_ibm_W)$r.squared
LM_W_test <- R2_e2 * T
LM_W_test
df_lm <- ncol(xx2)
df_lm
qchisq(.95, df = df_lm)
p_val <- 1 - pchisq(LM_W_test, df = df_lm) 	 # p-value of LM_test 
p_val


## CAT
fit_cat_W <- lm (e_cat2 ~ Mkt_RF + SMB + HML + xx2)
R2_e2 <- summary(fit_cat_W)$r.squared
LM_W_test <- R2_e2 * T
LM_W_test

qchisq(.95, df = df_lm) 

 p_val <- 1 - pchisq(LM_W_test, df = 6) 	 # p-value of LM_test 
 p_val


### 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)


#### GQ Test
gqtest(fit_gbp, fraction = .20)


## BP Test
e_gbp <- fit_gbp$residuals
e_gbp2 <- e_gbp^2
int_dif_2 <- int_dif^2
inf_dif_2 <- inf_dif^2
fit_gbp_BP <- lm (e_gbp2 ~ int_dif_2 + inf_dif_2)
Re_e2 <- summary(fit_gbp_BP)$r.squared
Re_e2

LM_BP_test <- Re_e2 * T
LM_BP_test

p_val <- 1 - pchisq(LM_BP_test, df = 1)  		# p-value of LM_test 
p_val


## White Test
int_inf <- inf_dif * int_dif 
fit_gbpt_W <- lm(e_gbp2 ~ int_dif + inf_dif + int_dif_2 + inf_dif_2 + int_inf)
R2_e2 <- summary(fit_cat_W)$r.squared
LM_W_test <- R2_e2 * T
LM_W_test

 p_val <- 1 - pchisq(LM_W_test, df = 6) 	 # p-value of LM_test 
 p_val


#### Autocorrelation tests

# BG Test for IBM
e_ibm <- fit_ibm_ff3$residuals
p_lag <- 12
e_lag <- matrix(0,T-p_lag,p_lag)
a <- 1
while (a <= p_lag) {
  za <- e_ibm[a:(T-p_lag+a-1)]
  e_lag[,a] <- za 
a <- a+1
}
e_lag

Mkt_RF_p <- Mkt_RF[(p_lag+1):T]
SMB_p <- SMB[(p_lag+1):T]
HML_p <- HML[(p_lag+1):T]

fit_ibm_p <- lm(e_ibm[(p_lag+1):T] ~ e_lag + Mkt_RF_p + SMB_p + HML_p)
r2_e1 <- summary(fit_ibm_p)$r.squared
lm_t <- (T-p_lag)*r2_e1
lm_t
df <- ncol(e_lag)
1-pchisq(lm_t,df)


fit_ibm_1 <- lm(e_ibm[(p_lag+1):T] ~ e_lag)
r2_e1 <- summary(fit_ibm_1)$r.squared
lm_t_1 <- (T-p_lag)*r2_e1
lm_t_1
df <- ncol(e_lag)
1-pchisq(lm_t_1,df)

# Same BG test done with bgtest from lmtest package (compared returns and excess returns)
bgtest(ibm_x ~ Mkt_RF + SMB + HML, order=12)

bgtest(lr_ibm ~ Mkt_RF + SMB + HML, order=12)

bgtest(ibm_x ~ Mkt_RF + SMB + HML, order=4)
bgtest(ibm_x ~ Mkt_RF + SMB + HML, order=12)


# 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)



#### BG tests for CAT & DIS in terms of returns and excess returns

bgtest(lr_cat ~ Mkt_RF + SMB + HML, order=4)
bgtest(lr_cat ~ Mkt_RF + SMB + HML, order=12)

bgtest(cat_x ~ Mkt_RF + SMB + HML, order=4)
bgtest(cat_x ~ Mkt_RF + SMB + HML, order=12)

bgtest(dis_x ~ Mkt_RF + SMB + HML, order=4)
bgtest(dis_x ~ Mkt_RF + SMB + HML, order=12)

bgtest(lr_dis ~ Mkt_RF + SMB + HML, order=4)
bgtest(lr_dis ~ Mkt_RF + SMB + HML, order=12)



#### Q & LB Tests
r_sum <- 0
lb_sum <- 0
p_lag <- 12
a <- 1
while (a <= p_lag) {
  za <- cor(e_ibm[(p_lag+1):T], e_ibm[a:(T-p_lag+a-1)])	
  r_sum <- r_sum + (za)^2 			#sum cor(e[(lag+1):T],e[a:(T-p_lag+a-1)])^2
  lb_sum <- lb_sum + (za)^2/(T-a) 
  a <- a + 1
}
Q <- T*r_sum
LB <- T*(T-2)*lb_sum
qchisq(.95,p_lag)
Q
LB
1-pchisq(Q,p_lag)
1- pchisq(LB,p_lag)


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

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


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

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