# The effect of omitting long-run trends from factor models 

rm(list=ls()) #Removes all items in Environment!
#setwd(path)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# set output options
options(width = 70, digits=4)
listofpackages <- c("dygraphs", "dplyr","ellipse","reshape2","ggplot2","PerformanceAnalytics","zoo")

for (j in listofpackages){
  if(sum(installed.packages()[, 1] == j) == 0) {
    install.packages(j)
  }
  library(j, character.only = T)
}
install.packages("fEcofin", repos="http://R-Forge.R-project.org")
 
# load required packages
library(fEcofin)                # various data sets

################################################################################
# Data Loadings and Transform Descriptive Analysis 
################################################################################

# create data frame with dates as rownames
berndt.df = berndtInvest[, -1]
berndt.df$date <- as.Date(berndtInvest[, 1])
rownames(berndt.df) = as.character(berndtInvest[, 1])
colnames(berndt.df)
dimnames(berndt.df)[[2]] #command alternative to the previous one 

# transform the data and compute cumulative returns 

t0 <- which(berndt.df$date == "1978-01-01")
t1 <- which(berndt.df$date == "1987-12-01")

series_names <- c("CITCRP","CONED","CONTIL","DATGEN","DEC","DELTA","GENMIL","GERBER","IBM","MARKET","MOBIL","PANAM","PSNH","TANDY","TEXACO","WEYER","RKFREE")

for (name in series_names) {
  P_col_name <- paste0(name,"_P")
  LP_col_name<- paste0("L",P_col_name)
  berndt.df[t0, P_col_name] <- 1
  for (i in (t0+1):(t1)) {
    berndt.df[i, P_col_name][[1]] <- berndt.df[i-1, P_col_name][[1]] * (1+berndt.df[i, name][[1]] )
  } 
  berndt.df[, LP_col_name] <- log(berndt.df[, P_col_name]) 
}
# add a trend to the database 
berndt.df$TREND<- array(data = NA, dim = nrow(berndt.df))
berndt.df[t0, c("TREND")] <- 1 # don't need to repeat the value to make the array being assigned be of the same length. be careful though as it is one of the few cases of exception
for (i in (t0+1):(t1)) {
  berndt.df[i, "TREND"][[1]] <- berndt.df[i-1, "TREND"][[1]] +1
} 

################################################################################
# Descriptive Analysis 
################################################################################
ggplot(berndt.df, aes(x = date)) +
  geom_line(aes(y = LTEXACO_P, color = "TEXACO"), size = 1, linetype = "solid") +
  geom_line(aes(y = LMARKET_P, color = "MARKET"), size = 1, linetype = "solid") +
  labs(x = "Date", y = "Portfolios TEXACO and MKT") +
  ylim(-0.5, 2) +
  theme_minimal() +
  theme(
    legend.position = c(0.15, 0.95),  # Set the legend position (top-left)
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green"),
    labels = c("TEXACO", "MARKET")
  )


# plot returns
ggplot(berndt.df, aes(x = date)) +
  geom_line(aes(y = TEXACO, color = "TEXACO"), size = 1, linetype = "solid") +
  geom_line(aes(y = MARKET, color = "MARKET"), size = 1, linetype = "solid") +
  labs(x = "Date", y = "Returns TEXACO and MKT") +
  ylim(-0.45, 0.45) +
  theme_minimal() +
  theme(
    legend.position = c(0.15, 0.95),  # Set the legend position (top-left)
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8),
    plot.title = element_text(size = 12, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green"),
    labels = c("TEXACO", "MARKET")
  )

################################################################################
# Standard CAPM Factor Models
################################################################################
 
## use lm function to compute single index model regressions for each asset
##
returns.mat = as.matrix(berndt.df[, c(1:9,11:16)])
asset.names = colnames(returns.mat)
asset.names

# initialize list object to hold regression objects

reg.list = list()
# loop over all assets and estimate time series regression
for (i in asset.names) {
 reg.df = berndt.df[, c(i, "MARKET")]
 si.formula = as.formula(paste(i,"~", "MARKET", sep=" "))
 reg.list[[i]] = lm(si.formula, data=reg.df)
}

# examine the elements of reg.list  - they are lm objects!
names(reg.list)
class(reg.list$TEXACO)
reg.list$TEXACO
summary(reg.list$TEXACO)

# plot actual vs. fitted over time
# use chart.TimeSeries() function from PerformanceAnalytics package
dataToPlot = cbind(fitted(reg.list$TEXACO), berndt.df$TEXACO)
colnames(dataToPlot) = c("Fitted","Actual")
chart.TimeSeries(dataToPlot, main="Single Index Model for TEXACO",
                 colorset=c("black","blue"), legend.loc="bottomleft")

# scatterplot of the single index model regression
plot(berndt.df$MARKET, berndt.df$TEXACO, main="SI model for CITCRP",
     type="p", pch=16, col="blue",
     xlab="MARKET", ylab="TEXACO")
abline(h=0, v=0)
abline(reg.list$TEXACO, lwd=2, col="red")

## extract beta values, residual sd's and R2's from list of regression objects

reg.vals = matrix(0, length(asset.names), 3)
rownames(reg.vals) = asset.names
colnames(reg.vals) = c("beta", "residual.sd", "r.square")
for (i in names(reg.list)) {
    tmp.fit = reg.list[[i]]
    tmp.summary = summary(tmp.fit)
    reg.vals[i, "beta"] = coef(tmp.fit)[2]
    reg.vals[i, "residual.sd"] = tmp.summary$sigma
    reg.vals[i, "r.square"] = tmp.summary$r.squared
}
reg.vals

# print regression results

par(mfrow=c(1,2))
barplot(reg.vals[,1], horiz=T, main="Beta values", col="blue", cex.names = 0.75, las=1)
barplot(reg.vals[,3], horiz=T, main="R-square values", col="blue", cex.names = 0.75, las=1)
par(mfrow=c(1,1))


################################################################################
# CAPM in levels 
################################################################################

## use lm function to compute single index model regressions for each asset
##
selected_columns <- c("LCITCRP_P","LCONED_P","LCONTIL_P","LDATGEN_P","LDEC_P","LDELTA_P","LGENMIL_P","LGERBER_P","LIBM_P","LMOBIL_P","LPANAM_P","LPSNH_P","LTANDY_P","LTEXACO_P","LWEYER_P")

# Extract the specified columns and store them in a matrix
lprices.mat <- as.matrix(berndt.df[, selected_columns])

asset.names = colnames(lprices.mat)
asset.names

# initialize list object to hold regression objects

reg1.list = list()
# loop over all assets and estimate time series regression
for (i in asset.names) {
  #reg.df = berndt.df[, c(i, "LMARKET_P")]
  si.formula = as.formula(paste(i,"~", "LMARKET_P+TREND", sep=" "))
  reg1.list[[i]] = lm(si.formula, data=berndt.df)
}

# examine the elements of reg.list  - they are lm objects!
names(reg1.list)
class(reg1.list$LTEXACO_P)
reg1.list$LTEXACO_P
summary(reg1.list$LTEXACO_P)

# plot actual vs. fitted over time
# use chart.TimeSeries() function from PerformanceAnalytics package
dataToPlot = cbind(fitted(reg1.list$LTEXACO_P), berndt.df$LTEXACO_P)
colnames(dataToPlot) = c("Fitted","Actual")
chart.TimeSeries(dataToPlot, main="Single Index Model for price TEXACO",
                 colorset=c("black","blue"), legend.loc="bottomleft")

# scatterplot of the single index model regression
plot(berndt.df$LMARKET_P, berndt.df$LTEXACO_P, main="SI model for LTEXACO_P",
     type="p", pch=16, col="blue",
     xlab="MARKET", ylab="TEXACO")
abline(h=0, v=0)
abline(reg.list$TEXACO, lwd=2, col="red")

## extract beta values, residual sd's and R2's from list of regression objects

reg.vals1 = matrix(0, length(asset.names), 3)
rownames(reg.vals1) = asset.names
colnames(reg.vals1) = c("beta", "residual.sd", "r.square")
for (i in names(reg1.list)) {
  tmp.fit = reg1.list[[i]]
  tmp.summary = summary(tmp.fit)
  reg.vals1[i, "beta"] = coef(tmp.fit)[2]
  reg.vals1[i, "residual.sd"] = tmp.summary$sigma
  reg.vals1[i, "r.square"] = tmp.summary$r.squared
}
reg.vals1

# print regression results

par(mfrow=c(1,2))
barplot(reg.vals1[,1], horiz=T, main="Beta values", col="blue", cex.names = 0.75, las=1)
barplot(reg.vals1[,3], horiz=T, main="R-square values", col="blue", cex.names = 0.75, las=1)
par(mfrow=c(1,1))

########################################
# a Single Factor Model with Predictability : an illustration with TEXACO 
#######################################

#Log Level linear model between LCITCRP_P TREND an MARKET 
model_TEXACO_P=lm(berndt.df$LTEXACO_P ~ berndt.df$LMARKET_P+berndt.df$TREND)
summary(model_TEXACO_P)

ggplot(berndt.df, aes(x = date)) +
  geom_line(aes(y = LTEXACO_P, color = "TEXACO"), size = 1, linetype = "solid") +
  geom_line(aes(y = fitted(model_TEXACO_P), color = "Fitted"), size = 1, linetype = "solid") +
  labs(x = "Date", y = "Actual and Fitted") +
  ylim(-0.5, 2) +
  theme_minimal() +
  theme(
    legend.position = c(0.15, 0.95),  # Set the legend position (top-left)
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green"),
    labels = c("TEXACO", "Fitted")
  )


#store log level residuals as u and test for their stationarity 
u_TEXACO=as.matrix(model_TEXACO_P$residuals)
DuTEXACO=diff(u_TEXACO,lag=1)
model_DuTEXACO=lm(DuTEXACO ~  u_TEXACO[1:(nrow(u_TEXACO)-1)]-1)
summary(model_DuTEXACO)
D12uTEXACO=diff(u_TEXACO,lag=12)
model_D12uTEXACO=lm(D12uTEXACO ~  u_TEXACO[1:(nrow(u_TEXACO)-12)]-1)
summary(model_D12uTEXACO)
########################################################################

#Compute Log Returns for 1M ahead timespan (1 months)
logret1M_TEXACO=diff(berndt.df$LTEXACO_P,lag=1)
logret1M_MARKET=diff(berndt.df$LMARKET_P,lag=1)

#Compute Log Returns for 1Y ahead timespan (12 months)
logret1Y_TEXACO=diff(berndt.df$LTEXACO_P,lag=12)
logret1Y_MARKET=diff(berndt.df$LMARKET_P,lag=12)


########################################################################

#model the regression on log returns appending the u residuals as another variable
model_d_TEXACO_1M=lm(logret1M_TEXACO ~ logret1M_MARKET+u_TEXACO[1:(nrow(u_TEXACO)-1)])
summary(model_d_TEXACO_1M)

model_d_TEXACO_1Y=lm(logret1Y_TEXACO ~ logret1Y_MARKET+u_TEXACO[1:(nrow(u_TEXACO)-12)])
summary(model_d_TEXACO_1Y)

model_d_TEXACO_1Y_CAPM=lm(logret1Y_TEXACO ~ logret1Y_MARKET)
summary(model_d_TEXACO_1Y_CAPM)