############# Regression with lm function & Testing single H0 (Code #4) #############

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


### Regression with Linear Algebra
y <- lr_ibm - RF
T <- length(y)
x0 <- matrix(1,T,1)
x <- cbind(x0, Mkt_RF, SMB, HML)
k <- ncol(x)					# number of regressors (including constant)
k
b <- solve(t(x)%*% x)%*% t(x)%*%y	 	# b = (X'X)^(-1) X' y  	(OLS regression)
b

e <- y - x%*%b					# OLS "errors" = residuals
e2 <- e^2					# squared residual vector
RSS <- sum(e2)					# sum of squared residuals (RSS)
sum_e2_1 <- t(e)%*%e				# sum of squared residuals using dot product
Sigma2 <- as.numeric(RSS/(T-k))			# Estimated σ2 = s2 
SE_reg <- sqrt(Sigma2) 				# Estimated σ – Regression stand error
Var_b <- Sigma2*solve(t(x)%*% x)		# Estimated Var[b|X] = s2 * (X′X)^(-1) 
SE_b <- sqrt(diag(Var_b))			# SE[b|X] 
SE_b
t_b <- b/SE_b					# t-values
t_b
y_hat <- x%*%b					# fitted values


# Scatter plot of fitted values 
plot(y_hat, type="l", main="Plot: Fitted IBM Returns",
   xlab="IBM returns ", ylab="Fitted IBM returns ", pch=19)
lines(y, type="l", col="red")
legend("topleft",                                       		# Add legend to plot
       legend = c("Actual", "Fitted"),        col = c("blue", "red"),        lty = 1)


### Test of H0: beta = 1
t_b2 <- (b[2]-1)/SE_b[2]
t_b2 
p_val <- (1 - pnorm(abs(t_b2))) * 2		# pvalue for t_b2 (adjust b/c two sided test)
p_val

## Interpretation of p_value and test result using cat in R
if (p_val < 0.05) {
  cat("p_values is:", p_val, "=> beta is significantly different from 1.\n")
} else {
  cat("p_values is:", p_val, "=> beta is not ignificantly different from 1. \n")
}


### Test of H0 using  R's lm function
fit_ibm_ff3 <- lm(y ~ Mkt_RF + SMB + HML)
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_beta1 <- (b_ibm[2] - 1)/SE_ibm[2]		# t-stat for H0: Beta1 - 1
p_val_beta1 <- (1- pnorm(abs(t_beta1))) * 2	# pvalue for t_beta
p_val_beta1


#### F-M two pass regression with 25 Fama-French portfolios

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 <- FF_f_da$Mkt_RF		# extract Market excess returns (in %)
HML <- FF_f_da$HML			# extract HML returns (in %)
SMB <- FF_f_da$SMB			# extract HML returns (in %)
RF <- FF_f_da$RF			# extract Risk-free rate (in %)

Y_p <- FF_p_da[,2:26]			# Extract 25 FF portfolios an


T <- length(HML)
x0 <- matrix(1,T,1)
x <- cbind(x0, Mkt_RF, SMB, HML)
summary(Mkt_RF)
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
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])
cor(beta_ret)

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

plot(beta_ret[,1], beta_ret[,3], main="Scatterplot: Portfolio Returns & Market 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)