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

####### REVIEW MODEL SELECTION

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


## GETS with package olsrr
library(olsrr)

## We need to put all the explanatory variables in a data frame
## January Dummy
Jan <- rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (round(T)/12+1))	# Create January dummy
T <- length(ibm_x)
T2 <- T+1
Jan_1 <- Jan[2:T2]

## Financial Crisis Dummy
t_sb <- 428							# Structural break date (End of 1st-regime)
T_s_1 <- T - t_sb
d_0 <- matrix(0, t_sb, 1)
d_1 <- matrix(1, T_s_1, 1)
Fin_cris <- rbind(d_0,d_1)
length(Fin_cris )

## Interactive terms
Mkt_RF2 <- Mkt_RF^2
SMB2 <- SMB^2
HML2 <- HML^2
Mkt_SMB <- Mkt_RF * SMB
Mkt_HML <- Mkt_RF * HML
SMB_HML <- SMB * HML
Mkt_Jan <- Mkt_RF * Jan_1
SMB_Jan <- SMB * Jan_1
HML_Jan <- HML * Jan_1
Mkt_FC <- Mkt_RF * Fin_cris 
HML_FC <- HML * Fin_cris 
SMB_FC <- SMB * Fin_cris 

ff_step_data <- data.frame(Mkt_RF, SMB, HML, Jan_1, Fin_cris, Mkt_RF2, HML2, SMB2, Mkt_HML, Mkt_SMB, SMB_HML, Mkt_Jan, SMB_Jan, HML_Jan, Mkt_FC, HML_FC, SMB_FC) 

## Fit GUM with lm
ibm_ff_gum <- lm(ibm_x ~ ., data = ff_step_data)
summary(ibm_ff_gum)	

## Use GETS (backward search, in this package). We can specify the p-value (prem) to eliminate variables.
ols_step_backward_p(ibm_ff_gum, prem = .05, details = TRUE)	# long final output & prem = p-value	


### Model Selection - Best Subset
ff_step_data_bs <- data.frame(Mkt_RF, SMB, HML, CMA, RMW, Jan_1)  # Define data set to be selected/reduced by best subset
fit_ibm_ff_bs <- lm(ibm_x ~ ., data = ff_step_data_bs)	# default p-value (penter) is 0.3 ols_step_best_subset(fit_ibm_ff3 , details = TRUE)	#  long final output

ibm_fit_best_subset <- ols_step_best_subset(fit_ibm_ff_bs, metric ="adjr")
ibm_fit_best_subset								# print result

# S3 plot for ols_step_best_subset
plot(ibm_fit_best_subset)


####### Heteroscedasticity and Autocorrelation

#### Fit IBM 3-factor FF Model
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)

#### 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			# Extract IBM residuals
e2 <- e_ibm^2					# Squared IBM residuals
Mkt_RF_2 <- Mkt_RF^2				# Potential Driver of IBM residual variance	
SMB_2 <- SMB^2					# Potential Driver of IBM residual variance
fit_ibm_BP <- lm (e2 ~ Mkt_RF_2 + SMB_2)	# Aux Regression: Regress squared IBM residuals on drivers
Re_e2 <- summary(fit_ibm_BP)$r.squared		# Extract R^2 from Aux Regression, Re_e2
Re_e2 
LM_BP_test <- Re_e2 * T				# Compute LM Test: Re_e2 * length of data
LM_BP_test

 p_val <- 1 - pchisq(LM_BP_test, df = 2) 	 # p-value of LM_test & df= 2 variables in Aux Reg
 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 = 3) 	 # p-value of LM_test & df= 3 variables in Aux Reg
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 = 3) 	 # p-value of LM_test & df= 3 variables in Aux Reg
 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) + 3				# df = 3 original variables + extra xx2 variables
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 = df_lm) 	 # 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=",")
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 = 2)  		# 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 = 5) 	 # 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 ( excess returns, with different lags & simple returns)
bgtest(ibm_x ~ Mkt_RF + SMB + HML, order=12)
bgtest(ibm_x ~ Mkt_RF + SMB + HML, order=4)

## In general, results are not that different from returns (lr_ibm) 
bgtest(ibm_x ~ Mkt_RF + SMB + HML, order=12)

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