####### Lecture GLS & FGLS

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

x_dis <- SFX_da$DIS
x_ibm <- SFX_da$IBM
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_dis)

 lr_dis <- log(x_dis[-1]/x_dis[-T])
 lr_ibm <- log(x_ibm[-1]/x_ibm[-T])
 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
dis_x <- lr_dis - RF

T <- length(dis_x)

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


######## White & NW SE

library(sandwich)
b_dis <- fit_dis_ff3$coefficients
White <- vcovHC(fit_dis_ff3, type = "HC0")		# HC0 = Standard White SE
SE_White <- sqrt(diag(White))				# White SE HC0
t_White <- b_dis/SE_White

t_White 

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

t_White3 

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

t_NW

NW <- NeweyWest(fit_dis_ff3, lag = 12, prewhite = FALSE)
SE_NW <- diag(sqrt(abs(NW)))
t_NW <- b_dis/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)


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