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

### You can compute White SE using function coeftest, from lmtest package (uses sandwich package too)
library(lmtest)
coeftest(fit_ibm_ff3, vcov = vcovHC(fit_ibm_ff3, type = 'HC3'))

# Same result using vcov = White3
coeftest(fit_ibm_ff3, vcov = 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

### Again, you can compute NW SE using function coeftest, with vcov = MV 
library(lmtest)
coeftest(fit_ibm_ff3, vcov = 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)

#### FGLS for DIS with Sigma2 = gamma_0 + gamma_1 * Mkt_RF^2 + gamma_2 * SMB^2
fit_dis_ff3 <- lm(dis_x ~ Mkt_RF + SMB + HML)
e_dis <- fit_dis_ff3$residuals
e_dis2 <- e_dis^2
SMB2 <- SMB^2
fit_dis2 <- lm(e_dis2 ~ Mkt_RF2 + SMB2)
summary(fit_dis2)

var_dis2 <- fit_dis2$fitted 			# Estimated variance vector, with elements Sigma2_hat
w_fgls <- sqrt(var_dis2)			# to be used as 1/w_fgls = weights	

y_fw <- dis_x/w_fgls				# transformed y
xx_fw <- cbind(x0, Mkt_RF, SMB, HML)/w_fgls	# transformed x
fit_dis_fgls <- lm(y_fw ~ xx_fw - 1)		# OLS with transformed y & x
summary(fit_dis_fgls)
	

#### FGLS for AR(1) - Using Cochrane-Orcutt

# C.O. function requires Y, X (with constant), OLS b.
c.o.proc <- function(Y,X,b_0,tol){
  T <- length(Y)
  e <- Y - X%*%b_0 	# OLS residuals
rss <- sum(e^2)	# Initial RSS of model, RSS9 		
rss_1 <- rss	# RSS_1 will be used to reset RSS after each iteration
d_rss = rss  	# initialize d_rss: difference between RSSi & RSSi-1
  e2 <- e[-1]	# adjust sample size for et	   
	e3 <- e[-T]	# adjust sample size for et-1	 
  ols_e0 <- lm(e2 ~ e3 - 1) 	# OLS to estimate rho	
  rho <- ols_e0$coeff[1]	# initial value for rho, 0
 i<-1
while (d_rss > tol) {	# tolerance of do loop. Stop when diff in RSS < tol
rss <- rss_1	# RSS at iter (i-1)
   YY <- Y[2:T] - rho * Y[1:(T-1)]	# pseudo-diff Y
   XX <- X[2:T, ] - rho * X[1:(T-1), ] 	# pseudo-diff X
ols_yx <- lm(YY ~ XX - 1)	# adjust if constant included in X
b <- ols_yx$coef	# updated OLS b at iteration i
#  b[1] <- b[1]/(1-rho)	# If constant not pseudo-differenced remove tag #
e1 <- Y - X%*%b	# updated residuals at iteration i
	 e2 <- e1[-1]	# adjust sample size for updated et		
    e3 <- e1[-T]	# adjust sample size for updated e_t-1 (lagged et)
    ols_e1 <- lm(e2~e3-1)	# updated regression to value for rho at iteration i
    rho <- ols_e1$coeff[1]	# updated value of rho at iteration i, i
 rss_1 <- sum(e1^2)	# updated value of RSS at iteration i, RSSi
 d_rss <- abs(rss_1 - rss)	# diff in RSS (RSSi - RSSi-1)
i <- i+1
    }

 result <-list()
 result$Cochrane_Orc.Proc <- summary(ols_yx)
 result$rho.regression <- summary(ols_e1)
 #  result$Corrected.b_1 <- b[1]
 result$Iterations <- i-1
 return(result)
}

### FGLS for Mexican interest rates

FX_da <- read.csv("http://www.bauer.uh.edu/rsusmel/4397/FX_USA_MX.csv", head=TRUE, sep=",")
names(FX_da)
us_CPI <- FX_da$US_CPI					# Read US CPI (US_CPI)
us_i <- FX_da$US_int					# Read US 3-mo Interest rate (i_us) data from file
mx_CPI <- FX_da$MX_CPI					# Read MEX CPI (MX_CPI)
mx_GDP <- FX_da$MX_GDP					# Read MEX CPI (MX_GDP)
mx_i <- FX_da$MX_int					# Read MEX 3-mo Interest rate (i_mx) data from file
S <- FX_da$MXN_USD					# Read changes in MXN/USD (e)


T <- length(S)
S_T <- S[T]
e_mx <- log(S[-1]/S[-T])				# Tranform S into log changes (e_mx) 
us_I <- log(us_CPI[-1]/us_CPI[-T])			# Tranform US CPI into Inflation (us_I) 
us_i_1 <- us_i[-1]/100	 				# Adjust sample size (we lose one observation from log changes)
mx_I <- log(mx_CPI[-1]/mx_CPI[-T])			# Tranform MX CPI into Inflation (mx_I) 
mx_y <- log(mx_GDP[-1]/mx_GDP[-T])			# Tranform MX GDP into growth rate (mx_y) 

mx_i_1 <- mx_i[-1]/100	
y <- mx_i_1
T_mx <- length(mx_i_1)
xx_i <- cbind(us_i_1, e_mx, mx_I, mx_y)		
x0 <- matrix(1,T_mx,1)
X <- cbind(x0,xx_i)					# X matrix
fit_i <- lm(mx_i_1 ~ us_i_1 + e_mx + mx_I + mx_y)	# OLS estimation
b_i <-fit_i$coefficients				# extract coefficients from lm	
summary(fit_i)

c.o.proc(y, X, b_i, .0001)				# Cochrane-Orcutt


### Cochrane-Orcutt with R package Orcutt for Mexican interest rates
install.packages("orcutt")
library(orcutt)
coch_i <- cochrane.orcutt(fit_i, convergence = 8, max.iter=100)
coch_i
t_coch_i <- coch_i$t.value	# Extract t-values from coch_i
t_coch_i