library(quantmod)
library(plotly)
library(PerformanceAnalytics)

# a. Use quantmod for downloading stock data.

saham <- c('BBCA.JK', 'ASII.JK', 'TLKM.JK', 'UNVR.JK', 'ICBP.JK')
sahamlabel = c("Bank Central Asia Tbk.", "Astra International Tbk.", "Telkom Indonesia (Persero)", "Unilever Indonesia Tbk.", "Indofood CBP Sukses Makmur Tbk.")

HargaPenutupan <- NULL                                                          
for (data in saham)
HargaPenutupan <- cbind(HargaPenutupan,getSymbols.yahoo(data, from="2022-01-01", auto.assign=FALSE)[,4])    # kombinasikan data harga saham (data frame)
HargaPenutupan <- HargaPenutupan[apply(HargaPenutupan,1,function(x) all(!is.na(x))),]                              # hapus semua data tanpa harga (jika ada) 
names(HargaPenutupan)[names(HargaPenutupan) ==
                 c("BBCA.JK.Close", "ASII.JK.Close", "TLKM.JK.Close", "UNVR.JK.Close", "ICBP.JK.Close")] <-
  
                 c("BBCA", "ASII", "TLKM", "UNVR", "ICBP")

# b. Define a function calculate_returns() to calculate daily returns from the closing prices of the stocks and apply this function to each stock's price series

calculate_returns <- function(hargasaham) {
  returns <- diff(log(hargasaham))              # Perhitungan return menggunakan metode log
  return(returns)
}
returnsaham <- lapply(HargaPenutupan, calculate_returns)

# c. Combine the returns of all stocks into a single data frame.

dfreturn <- do.call(cbind, returnsaham)

head(dfreturn)

# d. Calculate the correlation matrix of returns using cor().

korelasisaham = cor(dfreturn, use = "pairwise.complete.obs")    # Korelasi dengan metode Pairwise
data.frame(korelasisaham)

# e. Visualize the correlation matrix using corrplot().

colnames(korelasisaham) = rownames(korelasisaham) = sahamlabel # Penamaan pada data frame korelasi

mask = matrix(F, nrow=nrow(korelasisaham), ncol=ncol(korelasisaham))
mask[lower.tri(mask)] = T

heatmap = plot_ly(z=korelasisaham, x=sahamlabel, y=sahamlabel, type="heatmap")
heatmap = heatmap %>%
  layout(title = "Correlation Heatmap of Indonesian Stocks Return Price", 
         xaxis = list(title=""), 
         yaxis=list(title=""),
         colorscale= "Inferno",
         margin = list(l=60, r=80, t=100, b=60))

heatmap

# f. Calculate the portfolio's returns and standard deviation based on equal weights for each stock.

n <- length(HargaPenutupan)
weights <- rep(1/n, n)
  
portreturn <- colSums(na.omit(dfreturn * weights))   # na.omit untuk mengabaikan nilai NA
portsd <- sqrt(var(portreturn))
    
    
print(paste("Portfolio Returns:", mean(portreturn)))

print(paste("Portfolio Standard Deviation:", portsd)) 

return = na.omit(dfreturn)
colmean = apply(return, 2, mean, na.rm=T)
colsd <- apply(return, 2, sd, na.rm = TRUE)

print(colsd)

# g. Calculate the diversification ratio, which measures the extent to which the portfolio's risk is reduced due to diversification.
calculate_diversification_ratio <- function(portsd, colsd) {
weight_persaham <- sum(colsd) / length(colsd)
  
# Calculate diversification ratio
diversification_ratio <- portsd / weight_persaham
  
  return(diversification_ratio)
}

diversification_ratio <- calculate_diversification_ratio(portsd, colsd)
print(paste("Diversification Ratio:", diversification_ratio))

# h. Define the objective function and constraints for portfolio optimization. Here, we aim to maximize the portfolio return subject to the constraint that the sum of portfolio weights equals 1.
library(quadprog)

ekpektasireturn <- c(0.1, 0.1, 0.1, 0.1, 0.1)  # Return yang diinginkan dari masing-masing saham
cov_matrix <- matrix(korelasisaham, nrow=5)  # Covariance matrix

library(Matrix)

pd_D_mat <- nearPD(cov_matrix)

D_mat <- as.matrix(pd_D_mat$mat)
d_vec <- rep(0, length(ekpektasireturn))
A_mat <- cbind(rep(1, length(ekpektasireturn)), diag(length(ekpektasireturn)))
b_vec <- c(1, d_vec)

# i. Find the optimal portfolio weights that maximize the portfolio return given the covariance matrix and constraints.
# Metode optimasi maksimum
library(quadprog)
output <- solve.QP(Dmat = D_mat, dvec = d_vec, Amat = A_mat, bvec = b_vec, meq = 1)

# Optimal portfolio weights
optimal_weights <- output$solution
print(optimal_weights)

# j. Visualize the optimal weights assigned to each stock in the portfolio.
library(plotly)

colors <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd")
plot_data <- data.frame(Symbol = saham, Weight = optimal_weights, Color = colors)

bar_plot <- plot_ly(plot_data, x = ~Symbol, y = ~Weight, type = 'bar', 
                    marker = list(color = ~Color)) %>%
  layout(title = "Optimal Weights of Stocks in Portfolio",
         xaxis = list(title = "Stocks"),
         yaxis = list(title = "Weight"))

bar_plot

# k. Calculate the expected return and volatility of the optimized portfolio using the optimal weights and covariance matrix.
portekpektasireturn <- sum(optimal_weights * ekpektasireturn)

# Calculate portfolio volatility (standard deviation)
portvolatility <- sqrt(t(optimal_weights) %*% cov_matrix %*% optimal_weights)

# Print the results

print(paste("Expected Return of the Portfolio:", portekpektasireturn))

print(paste("Portfolio Volatility (Standard Deviation):", portvolatility*portsd))

# a. Generate insurance claims data for a hypothetical insurance company with 10000 policies over 17 years.
set.seed(863)
policies <- 10000
years <- 17

# Mean number of claims per policy per year (assuming Poisson distribution)
lambda <- 0.1


# Initialize empty data frame to store claims data
claims_data <- data.frame(
  Policy_ID = integer(),
  Year = integer(),
  Claims = numeric(),
  Claim_Amount = numeric()
)

# Generate claims data for each policy and year
for (policy_id in 1:policies) {
  for (year in 1:years) {
# c. Simulate claim frequencies for each policy from a Poisson distribution with a mean of 0.1 (i.e., on average, each policy has 0.1 claims per year).
    num_claims <- rpois(1, lambda)
    
    if (num_claims > 0) {
# d. Simulate claim amounts for each policy and each year from a normal distribution with a mean of 5007 and a standard deviation of 2007.
    claim_amounts <- round(rnorm(num_claims, 5007, 2007),2)
    
    claims_data <- rbind(claims_data, data.frame(
      Policy_ID = rep(policy_id, num_claims),
      Year = rep(year, num_claims),
      Claims = num_claims,
      Claim_Amount = claim_amounts
    ))
  }
}
}

claims_data

# b. Generate random policy premiums from a normal distribution with a mean of 1007 and a standard deviation of 2007.
set.seed(863)
policy_premiums <- abs(round(rnorm((nrow(claims_data)), mean = 1007, sd = 2007), 2))
claims_data$Premiums = policy_premiums

# e. Calculate total claim amounts per policy and then calculate loss ratios (total claims divided by premiums) to assess the risk associated with each policy.
total_claim_amounts <- tapply(claims_data$Claim_Amount, claims_data$Policy_ID, sum)

total_premiums <- tapply(claims_data$Premiums, claims_data$Policy_ID, sum)

loss_ratios <- total_claim_amounts / total_premiums

loss = data.frame(loss_ratios)
loss

# f. Plot a histogram of the loss ratios to visualize the distribution of risk.
loss_ratio_histogram <- plot_ly(x = loss_ratios, type = "histogram", 
                                name = "Loss Ratios", 
                                xbins = list(size = 0.1), 
                                marker = list(color = "yellow"))

loss_ratio_histogram <- loss_ratio_histogram %>% 
  layout(title = "Histogram of Loss Ratios",
         xaxis = list(title = "Loss Ratio"),
         yaxis = list(title = "Frequency"))

loss_ratio_histogram

# g. Calculate summary statistics including the mean, median, and 95th percentile of the loss ratios.
mean_loss_ratio <- mean(loss_ratios)
median_loss_ratio <- median(loss_ratios)
p95_loss_ratio <- quantile(loss_ratios, probs = 0.95)

# Display summary statistics
cat("Mean Loss Ratio:", mean_loss_ratio, "\n")

cat("Median Loss Ratio:", median_loss_ratio, "\n")

cat("95th Percentile of Loss Ratios:", p95_loss_ratio, "\n")