rm(list=ls())
#Check the working directory before importing else provide full path
#setwd(path)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# packages used
listofpackages <- c("dygraphs", "tidyverse","ellipse","reshape2","highcharter","xts","xlsx","readxl","quantmod","quadprog")

for (j in listofpackages){
  if(sum(installed.packages()[, 1] == j) == 0) {
    install.packages(j)
  }
  library(j, character.only = T)
}

raw_data           = read_xlsx("../data/2023_monthly_stocks.xlsx") 
names(raw_data)[1] = 'Date'
typeof(raw_data)
typeof(raw_data$Date)
typeof(raw_data$AXP)
typeof(raw_data$CSCO)

dates <-seq(as.Date("1985-02-01"),length=462, by="months")
params <- c("Date","AXP","AMGN","AAPL","BA","CAT","CSCO","CVX","GS",
            "HD","HON","IBM","INTC","JNJ","KO","JPM")
data <- raw_data[, c(params)]
data<- na.omit(data)
data <- data %>% 
  mutate(Date = as.Date(Date, format = "%Y-%m-%d"))

t1<-nrow(data)
series_names <- c("AXP","AMGN","AAPL","BA","CAT","CSCO","CVX","GS",
                  "HD","HON","IBM","INTC","JNJ","KO","JPM")

for (name in series_names) {
  return_col_name <- paste0(name, "_ret")
  data[, return_col_name] <- array(data = NA, dim = t1)
  for (i in 2:nrow(data)) {
    data[i, return_col_name][[1]] <- (data[i, name][[1]] - data[i - 1, name][[1]]) / data[i - 1, name][[1]]
  } 
}

params1 <- c("AXP_ret","AMGN_ret","AAPL_ret","BA_ret","CAT_ret","CSCO_ret","CVX_ret","GS_ret",
            "HD_ret","HON_ret","IBM_ret","INTC_ret","JNJ_ret","KO_ret","JPM_ret")
returns_data <- data[, c(params1)]
returns_data <- na.omit(returns_data)

returns_matrix<-as.matrix(returns_data)

# Calculate the mean returns and covariance matrix
mean_returns <- colMeans(returns_matrix)
cov_matrix <- cov(returns_matrix)

# Define the objective function to minimize risk (standard deviation)
portfolio.objective <- function(weights, cov_matrix) {
  portfolio_return <- sum(weights * mean_returns)
  portfolio_volatility <- sqrt(t(weights) %*% cov_matrix %*% weights)
  return(portfolio_volatility)
}

# Set up constraints (sum of weights = 1)
#A_eq <- matrix(1, nrow = 1, ncol = ncol(returns_matrix))
#b_eq <- matrix(1, nrow = 1)

# Generate a range of target returns
target_returns <- seq(min(mean_returns), max(mean_returns), length.out = 1000)

# Initialize vectors to store results
portfolio_returns <- numeric(length(target_returns))
portfolio_volatilities <- numeric(length(target_returns))

# Optimize for each target return
for (i in 1:length(target_returns)) {
  target_return <- target_returns[i]
  
  # Set up the quadratic programming problem
  Dmat <- 2*cov_matrix
  dvec <- matrix(rep(0, ncol(returns_matrix)),ncol=1)
  a1mat<- matrix(mean_returns, nrow =ncol(returns_matrix))
  a2mat<-matrix(rep(1, ncol(returns_matrix)), nrow =ncol(returns_matrix))
  Amat <- cbind(a1mat, a2mat)
  bvec <- matrix(c(target_return, 1),ncol=1)
  
  # Solve the optimization problem
  weights <- solve.QP(Dmat, dvec, Amat, bvec, meq = 2)$solution
  
  # Calculate portfolio risk (volatility)
  portfolio_volatility <- portfolio.objective(weights, cov_matrix)
  
  # Store results
  portfolio_returns[i] <- target_return
  portfolio_volatilities[i] <- portfolio_volatility
}

# Create a data frame for efficient frontier points
efficient_frontier <- data.frame(Return = portfolio_returns, Volatility = portfolio_volatilities)

# Plot the efficient frontier
ggplot(efficient_frontier, aes(x = Volatility, y = Return)) +
  geom_point() +
  labs(x = "Standard Deviation (Risk)", y = "Mean Return") +
  ggtitle("Efficient Frontier") +
  theme_minimal()


# Your efficient frontier plot...
gg <- ggplot(efficient_frontier, aes(x = Volatility, y = Return)) +
  geom_point(color = "blue", size = 0.5) +
  labs(x = "Standard Deviation (Risk)", y = "Mean Return") +
  ggtitle("Efficient Frontier") +
  theme_minimal()

# Add points for individual stocks with labels and adjust font size
stocks_data <- data.frame(Return = mean_returns, Volatility = sqrt(diag(cov_matrix)), Stock = colnames(returns_matrix))
gg <- gg + geom_point(data = stocks_data, aes(x = Volatility, y = Return), color = "red", size = 3) +
  geom_text(data = stocks_data, aes(x = Volatility, y = Return, label = Stock), hjust = 0, vjust = 0, nudge_x = 0.000, nudge_y = 0.000, size = 2)  # Adjust the size here

# Print the plot
print(gg)