############# Review of Testing AND Maximum Likelihood Estimation & Data Problems (Code #5) #############


#### Review: Testing CAPM (constant = alpha = 0) 
SFX_da <- read.csv("http://www.bauer.uh.edu/rsusmel/4397/Stocks_FX_1973.csv", head=TRUE, sep=",")
names(SFX_da)
summary(SFX_da)

## Extract variables from imported data
x_ibm <- SFX_da$IBM				# extract IBM price data
x_Mkt_RF <- SFX_da$Mkt_RF			# extract Market excess returns (in %)
x_RF <- SFX_da$RF				# extract Risk-free rate (in %)
x_SMB <- SFX_da$SMB
x_HML <- SFX_da$HML

# Define log returns & adjust size of variables accordingly
T <- length(x_ibm)				# sample size
lr_ibm <- log(x_ibm[-1]/x_ibm[-T])		# create IBM log returns (in decimal returns)
Mkt_RF <- x_Mkt_RF[-1]/100			# Adjust sample size to ( T-1) by removing 1st obs 
RF <- x_RF[-1]/100				# Adjust sample size and use decimal returns.
SMB <- x_SMB[-1]/100
HML <- x_HML[-1]/100

y <- lr_ibm - RF

### Test of H0: alpha = 0, using  R's lm function
fit_ibm_capm <- lm(y ~ Mkt_RF)
summary(fit_ibm_capm)
b_ibm <- fit_ibm_capm$coefficients		# Extract from lm function OLS coefficients
SE_ibm <- sqrt(diag(vcov(fit_ibm_capm)))	# SE from fit_ibm (also a kx1 vector)
t_alpha <- b_ibm[1]/SE_ibm[1]			# t-stat for H0: Alpha = 0 (same as in 3rd column of summary)
t_alpha						# Reject H0 if |t_alpha| > 1.96 (at 5% level)
p_val_alpha <- (1- pnorm(abs(t_alpha	))) * 2	# p-value for t_alpha
p_val_alpha					# (if p-value < .05 => reject)


##### Testing CAPM with F-M two pass regression with 25 ME/BM Fama-French portfolios (1926:July - 2023:July)

FF_p_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/FF_25_portfolios.csv", head=TRUE, sep=",")
FF_f_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/FF_3_factors.csv", head=TRUE, sep=",")
names(FF_p_da)
names(FF_f_da)
summary(FF_p_da)
summary(FF_f_da)

# Extract variables from imported data
Mkt_RF_fm <- FF_f_da$Mkt_RF		# extract Market excess returns (in %)
HML_fm <- FF_f_da$HML			# extract HML returns (in %)
SMB_fm <- FF_f_da$SMB			# extract HML returns (in %)
RF_fm <- FF_f_da$RF			# extract Risk-free rate (in %)

Y_p <- FF_p_da[,2:26]	- RF_fm		# Extract returns of 25 FF portfolios 


T <- length(HML_fm)			# Number of observations (1926:July on) 
x0 <- matrix(1,T,1)
x <- cbind(x0, Mkt_RF_fm)
summary(Mkt_RF_fm)
summary(Y_p)

k <- ncol(Y_p)
## First Pass
Allbs = NULL 				# Initialize empty (a space to put betas)		
for (i in seq(1,k,1)){
  y <- Y_p[,i]     # select Y (portfolio)
  b <- solve(t(x)%*% x)%*% t(x)%*%y 	# OLS regression = (X'X)^(-1) X'y
  Allbs  =cbind(Allbs,b)  		# accumulate b as rows
} 

mb <- rowMeans(Allbs)
mb
varb <- colVars(t(Allbs))
sd_b <- sqrt(varb)
sd_b

beta_ret <- cbind(colMeans(Y_p),t(Allbs))
cor(beta_ret[,1], beta_ret[,3])

plot(beta_ret[,1], beta_ret[,3], main="Scatterplot: Portfolio Returns & CAPM Beta",
     xlab="Mean Portfolio Returns ", ylab="Market Beta", pch=19)

## Second Pass (CAPM)
fit_fm_capm_25 <- lm(beta_ret[,1] ~ beta_ret[,3])
summary(fit_fm_capm_25)



#### Using F-F 3-factor Model
x <- cbind(x0, Mkt_RF_fm, SMB_fm, HML_fm)


## First Pass
Allbs=NULL 				# Initialize empty (a space to put betas)		
for (i in seq(1,k,1)){
  y <- Y_p[,i]     # select Y (portfolio)
  b <-solve(t(x)%*% x)%*% t(x)%*%y 	# OLS regression
  Allbs=cbind(Allbs,b)  		# accumulate b as rows
} 

beta_ret <- cbind(colMeans(Y_p),t(Allbs))
cor(beta_ret[,1], beta_ret[,3])

plot(beta_ret[,1], beta_ret[,3], main="Scatterplot: Portfolio Returns & Market Beta",
     xlab="Mean Portfolio Returns ", ylab="Market Beta", pch=19)

## Second Pass (3 FF) - OLS regression mean portfolio returns against betas
fit_fm_ff3_25 <- lm(beta_ret[,1] ~ beta_ret[,3:5])	# OLS mean portfolio returns vs betas (Mkt, SMB, HML)
summary(fit_fm_ff3_25)


#### Using F-F 3-factor Model starting in 1973:Jan (T1 = 559)
T1 <- 559				# T1 = 871 if data starts in 1990:Jan
x_T1 <- x[T1:T,]
Y_p_T1 <- Y_p[T1:T,]


## First Pass
Allbs=NULL 				# Initialize empty (a space to put betas)		
for (i in seq(1,k,1)){
  y <- Y_p_T1[,i]     # select Y (portfolio)
  b <-solve(t(x_T1)%*% x_T1)%*% t(x_T1)%*%y 	# OLS regression
  Allbs=cbind(Allbs,b)  		# accumulate b as rows
} 

beta_ret <- cbind(colMeans(Y_p),t(Allbs))
cor(beta_ret[,1], beta_ret[,3])
plot(beta_ret[,1], beta_ret[,3], main="Scatterplot: Portfolio Returns & Market Beta (1973-2023)",
     xlab="Mean Portfolio Returns ", ylab="Market Beta", pch=19)

## Second Pass (3 FF) - OLS regression mean portfolio returns against betas
fit_fm_ff3_25 <- lm(beta_ret[,1] ~ beta_ret[,3:5])	# OLS mean portfolio returns vs betas (Mkt, SMB, HML)
summary(fit_fm_ff3_25)



### Test of H0: Market Beta = 1 using  R's lm function
y <- lr_ibm - RF

fit_ibm_ff3 <- lm(y ~ Mkt_RF + SMB + HML)			# CAPM regression
summary(fit_ibm_ff3)
b_ibm <- fit_ibm_ff3$coefficients		# Extract from lm function OLS coefficients
SE_ibm <- sqrt(diag(vcov(fit_ibm_ff3)))		# SE from fit_ibm (also a kx1 vector)
t_beta_mkt <- (b_ibm[2] - 1)/SE_ibm[2]		# t-stat for H0: Market Beta = 1 (same as in 3rd column of summary)
t_beta_mkt					# Reject H0 if |t_beta_mkt| > 1.96 (at 5% level)
p_val_beta_mkt <- (1- pnorm(abs(t_beta_mkt))) * 2	# pvalue for t_beta
p_val_beta_mkt					# (if p-value < .05 => reject)