load("dataset_2.RData")
load("FF3.RData")

NASDAQ100_EX_ret <- NASDAQ100_EX[,-c(1,2)]

#COMPUTE STOCKS EXPECTED RETURNS ACCORDING TO CAPM (USING MARKET EXCESS RETURN FROM FRENCH REPOSITORY)

#compute the betas
betas_stocks <- vector()
for (i in 1:ncol(NASDAQ100_EX_ret)) {
  model <- lm(NASDAQ100_EX_ret[,i] ~ FF3$MKT)
  betas_stocks[i] <- model$coefficients[2]
}
names(betas_stocks) <- names(NASDAQ100_EX_ret)

#we create a vector with the mean risk-free rate and one with the mean market excess return
rfree_mean <- rep(mean(NASDAQ100$Rfree),length(betas_stocks))
MKT_mean <- rep(mean(FF3$MKT),length(betas_stocks))

#then we apply the CAPM formula for the expected return
E_ret_capm <- rfree_mean + betas_stocks*MKT_mean

#COMPUTE STOCKS EXPECTED RETURNS ACCORDING TO FAMA-FRENCH THREE-FACTOR MODEL

#we fit the first stage regression, obtaining the loadings
LOADINGS_LIST <- list()
for (i in 1:ncol(NASDAQ100_EX_ret)) {
  model <- lm(NASDAQ100_EX_ret[,i] ~ FF3$MKT + FF3$SMB + FF3$HML)
  LOADINGS_LIST[[i]] <- as.vector(model$coefficients[-1])
}
names(LOADINGS_LIST) <- names(NASDAQ100_EX_ret)
LOADINGS <- data.frame(do.call(rbind,LOADINGS_LIST))
names(LOADINGS) <- c("b1","b2","b3")

#we fit the second stage regression, obtaining the factor premia
PREMIA_LIST <- list()
for (t in 1:nrow(NASDAQ100_EX_ret)) {
  model <- lm(unlist(NASDAQ100_EX_ret[t,]) ~ unlist(LOADINGS$b1) + unlist(LOADINGS$b2) + unlist(LOADINGS$b3))
  PREMIA_LIST[[t]] <- as.vector(model$coefficients[-1])
}
names(PREMIA_LIST) <- as.character(1:nrow(NASDAQ100_EX_ret))
PREMIA_T <- data.frame(do.call(rbind,PREMIA_LIST))
PREMIA <- colMeans(PREMIA_T) #compute the mean across the T periods
names(PREMIA) <- c("p1","p2","p3")

#now we can compute the expected returns according to the Fama-French three-factor model
E_ret_FF3 <- rfree_mean + LOADINGS$b1*PREMIA[1] + LOADINGS$b2*PREMIA[2] + LOADINGS$b3*PREMIA[3]
names(E_ret_FF3) <- names(E_ret_capm)

#COMPARE THE TWO SETS OF ESTIMATED EXPECTED RETURNS
cor(E_ret_capm,E_ret_FF3)
plot(E_ret_capm,E_ret_FF3, pch=19, xlab="CAPM estimates", ylab="FF3 estimates", main="Expected returns")

#COMPUTE EXPECTED RETURNS WITH CAPM AND FF3 USING A 10-YEAR ROLLING WINDOW AND USE THEM TO CREATE A LONG-SHORT PORTFOLIO

#estimation window
WIND <- 120

#asset returns to be used to compute out-of-sample portfolio returns
NASDAQ100_EX_ret_oos <- NASDAQ100_EX_ret[-(1:WIND),]
#risk free rate to be used to compute out-of-sample portfolio returns
RFree_oos <- NASDAQ100_EX[-(1:WIND),2]

#function to assign positive equal weights to the 20 stocks with highest estimated expected return,
#and negative equal weights to the 20 stocks with lowest estimated expected returns
ls_weights <- function(returns) {
  N <- length(returns)
  
  #find the index positions of the 20 highest and lowest values
  lowest_indices <- order(returns)[1:20]
  highest_indices <- order(returns)[(N-19):N]
  
  #replace the values at the identified indices
  returns[lowest_indices] <- -1/20
  returns[highest_indices] <- 1/20
  
  #set the other values to 0
  returns[-c(lowest_indices, highest_indices)] <- 0
  
  weights <- returns
  return(weights)
}

E_RET_CAPM_LIST <- list()
WEIGHTS_CAPM_LIST <- list()
#loop to compute the expected return according to CAPM for each of the out-of-sample periods
for (m in 1:(nrow(NASDAQ100_EX_ret)-WIND)){
  
  #select the data in the current estimation window
  NASDAQ100_m <- NASDAQ100[m:(WIND+m-1),]
  NASDAQ100_EX_ret_m <- NASDAQ100_EX_ret[m:(WIND+m-1),]
  FF3_m <- FF3[m:(WIND+m-1),]
  rfree_mean_m <- rep(mean(NASDAQ100_m$Rfree),ncol(NASDAQ100_EX_ret_m))
  MKT_mean_m <- rep(mean(FF3_m$MKT),ncol(NASDAQ100_EX_ret_m))
  
  #compute the betas
  betas_stocks_m <- vector()
  for (i in 1:ncol(NASDAQ100_EX_ret_m)) {
    model <- lm(NASDAQ100_EX_ret_m[,i] ~ FF3_m$MKT)
    betas_stocks_m[i] <- model$coefficients[2]
  }
  names(betas_stocks_m) <- names(NASDAQ100_EX_ret_m)
  
  #apply the CAPM formula for the expected return
  E_ret_capm_m <- rfree_mean_m + betas_stocks_m*MKT_mean_m
  E_RET_CAPM_LIST[[m]] <- E_ret_capm_m
  
  #assign weights for long-short portfolio
  weights_capm_m <- ls_weights(E_ret_capm_m)
  WEIGHTS_CAPM_LIST[[m]] <- weights_capm_m
}
#join the estimated expected returns in a matrix. each row is an out-of-sample period
E_RET_CAPM <- do.call(rbind,E_RET_CAPM_LIST)
#do the same for the weights
WEIGHTS_CAPM <- do.call(rbind,WEIGHTS_CAPM_LIST)

#compute out-of-sample returns of the long-short portfolio using weights suggested by CAPM estimates
EXR_CAPM <- vector()
for (m in 1:nrow(WEIGHTS_CAPM)){
  EXR_CAPM[m] <- t(unlist(WEIGHTS_CAPM[m,]))%*%unlist(NASDAQ100_EX_ret_oos[m,])
}
#add the risk-free rate to the excess returns to get the returns
R_CAPM <- EXR_CAPM + RFree_oos
#compute mean and standard deviation
PERF_CAPM <- c(mean(R_CAPM), sd(R_CAPM), mean(EXR_CAPM)/sd(EXR_CAPM))
names(PERF_CAPM) <- c("mean","sd","SR")

E_RET_FF3_LIST <- list()
WEIGHTS_FF3_LIST <- list()
#loop to compute the expected return according to FF3 for each of the out-of-sample periods
for (m in 1:(nrow(NASDAQ100_EX_ret)-WIND)){
  
  #select the data in the current estimation window
  NASDAQ100_m <- NASDAQ100[m:(WIND+m-1),]
  NASDAQ100_EX_ret_m <- NASDAQ100_EX_ret[m:(WIND+m-1),]
  FF3_m <- FF3[m:(WIND+m-1),]
  rfree_mean_m <- rep(mean(NASDAQ100_m$Rfree),ncol(NASDAQ100_EX_ret_m))
  MKT_mean_m <- rep(mean(FF3_m$MKT),ncol(NASDAQ100_EX_ret_m))
  
  #we fit the first stage regression, obtaining the loadings for the current period
  LOADINGS_LIST_M <- list()
  for (i in 1:ncol(NASDAQ100_EX_ret_m)) {
    model <- lm(NASDAQ100_EX_ret_m[,i] ~ FF3_m$MKT + FF3_m$SMB + FF3_m$HML)
    LOADINGS_LIST_M[[i]] <- as.vector(model$coefficients[-1])
  }
  names(LOADINGS_LIST_M) <- names(NASDAQ100_EX_ret_m)
  LOADINGS_M <- data.frame(do.call(rbind,LOADINGS_LIST_M))
  names(LOADINGS_M) <- c("b1","b2","b3")
  
  #we fit the second stage regression, obtaining the factor premia
  PREMIA_LIST_M <- list()
  for (t in 1:nrow(NASDAQ100_EX_ret_m)) {
    model <- lm(unlist(NASDAQ100_EX_ret_m[t,]) ~ unlist(LOADINGS_M$b1) + unlist(LOADINGS_M$b2) + unlist(LOADINGS_M$b3))
    PREMIA_LIST_M[[t]] <- as.vector(model$coefficients[-1])
  }
  names(PREMIA_LIST_M) <- as.character(1:nrow(NASDAQ100_EX_ret_m))
  PREMIA_T_M <- data.frame(do.call(rbind,PREMIA_LIST_M))
  PREMIA_M <- colMeans(PREMIA_T_M) #compute the mean across the T periods
  names(PREMIA_M) <- c("p1","p2","p3")
  
  #compute the expected returns according to the Fama-French three-factor model
  E_ret_FF3_m <- rfree_mean_m + LOADINGS_M$b1*PREMIA_M[1] + LOADINGS_M$b2*PREMIA_M[2] + LOADINGS_M$b3*PREMIA_M[3]
  names(E_ret_FF3_m) <- names(NASDAQ100_EX_ret_m)
  E_RET_FF3_LIST[[m]] <- E_ret_FF3_m
  
  #assign weights for long-short portfolio
  weights_FF3_m <- ls_weights(E_ret_FF3_m)
  WEIGHTS_FF3_LIST[[m]] <- weights_FF3_m
}
#join the estimated expected returns in a matrix. each row is an out-of-sample period
E_RET_FF3 <- do.call(rbind,E_RET_FF3_LIST)
#do the same for the weights
WEIGHTS_FF3 <- do.call(rbind,WEIGHTS_FF3_LIST)

#compute out-of-sample returns of the long-short portfolio using weights suggested by FF3 estimates
EXR_FF3 <- vector()
for (m in 1:nrow(WEIGHTS_FF3)){
  EXR_FF3[m] <- t(unlist(WEIGHTS_FF3[m,]))%*%unlist(NASDAQ100_EX_ret_oos[m,])
}
#add the risk-free rate to the excess returns to get the returns
R_FF3 <- EXR_FF3 + RFree_oos
#compute mean and standard deviation
PERF_FF3 <- c(mean(R_FF3), sd(R_FF3), mean(EXR_FF3)/sd(EXR_FF3))
names(PERF_FF3) <- c("mean","sd","SR")

#compare performance
print(round(PERF_CAPM,4))
print(round(PERF_FF3,4))

#compute and plot the wealth
DATE <- seq(as.Date("2009-01-01"),length=nrow(NASDAQ100_EX_ret_oos)+1,by="months")-1
CAPM_wealth <- c(1, cumprod(1 + R_CAPM))
plot(CAPM_wealth~DATE, type="l", lwd=2, col="red", ylab="Wealth")
FF3_wealth <- c(1, cumprod(1 + R_FF3))
lines(FF3_wealth~DATE, lwd=2, col="blue")