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

closing = NULL
for (stock in tick) { closing = cbind(closing, Cl(get(stock)))}

return = Return.calculate(closing)
return = data.frame(return)
head (return)

Calculate the correlation matrix of returns using cor().

corstock = cor(na.omit(return))
corstock

##               AALI.JK.Close ABBA.JK.Close ABDA.JK.Close ABMM.JK.Close
## AALI.JK.Close    1.00000000    0.15500305    0.01175096    0.19766090
## ABBA.JK.Close    0.15500305    1.00000000   -0.01351559    0.07805656
## ABDA.JK.Close    0.01175096   -0.01351559    1.00000000    0.01949983
## ABMM.JK.Close    0.19766090    0.07805656    0.01949983    1.00000000
## ACES.JK.Close    0.11384187    0.17354554   -0.04318530   -0.01698310
##               ACES.JK.Close
## AALI.JK.Close     0.1138419
## ABBA.JK.Close     0.1735455
## ABDA.JK.Close    -0.0431853
## ABMM.JK.Close    -0.0169831
## ACES.JK.Close     1.0000000

Visualize the correlation matrix using corrplot().

library(corrplot)

## corrplot 0.92 loaded

corrplot(corstock, method = 'number')

colnames(corstock) = rownames(corstock) = ticklabel #Renaming 

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

heatmap = plot_ly(z=corstock, x=ticklabel, y=ticklabel, 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

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

 num_stocks <- length(closing)
  weights <- rep(1/num_stocks, num_stocks)
  
  portfolio_returns <- rowSums(closing * weights)
    portfolio_sd <- sqrt(var(portfolio_returns))
    
cat("Portfolio Returns:", mean(portfolio_returns), "\n")

## Portfolio Returns: 7.740259

cat("Portfolio Standard Deviation:", portfolio_sd, "\n")

## Portfolio Standard Deviation: 0.6390579

Calculate the diversification ratio, which measures the extent to which the portfolio's risk is reduced due to diversification.

return = na.omit(return)
column_mean = apply(return, 2, mean, na.rm=T)
column_sd <- apply(return, 2, sd, na.rm = TRUE)

# Print the standard deviation of each column
print(column_sd)

## AALI.JK.Close ABBA.JK.Close ABDA.JK.Close ABMM.JK.Close ACES.JK.Close 
##    0.01553983    0.03384390    0.01417629    0.03037228    0.02847660

calculate_diversification_ratio <- function(portfolio_sd, column_sd) {
  weighted_avg_individual_sd <- sum(column_sd) / length(column_sd)
  
  # Calculate diversification ratio
  diversification_ratio <- weighted_avg_individual_sd / portfolio_sd
  
  return(diversification_ratio)
}

diversification_ratio <- calculate_diversification_ratio(portfolio_sd, column_sd)
cat("Diversification Ratio:", diversification_ratio, "\n")

## Diversification Ratio: 0.03830917

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)
# Example data (replace with your actual data)
expected_returns <- c(0.08, 0.08, 0.08, 0.08, 0.08)  # Expected returns for assets
cov_matrix <- matrix(corstock, nrow=5)  # Covariance matrix

library(Matrix)
pd_D_mat <- nearPD(cov_matrix)

# Use the positive definite matrix for optimization
D_mat <- as.matrix(pd_D_mat$mat)
d_vec <- rep(0, length(expected_returns))
A_mat <- cbind(rep(1, length(expected_returns)), diag(length(expected_returns)))
b_vec <- c(1, d_vec)

# Solve quadratic programming problem
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)

## [1] 0.1535611 0.1746983 0.2540573 0.2028852 0.2147981

Visualize the optimal weights assigned to each stock in the portfolio.

pie <- plot_ly(labels = ticklabel, values = optimal_weights, type = "pie")

pie <- pie %>% layout(title = "Optimized Portfolio Weight", xaxis = list(title = ""), yaxis = list(title = ""))

pie

Calculate the expected return and volatility of the optimized portfolio using the optimal weights and

covariance matrix.

cat("Portfolio Volatility (Standard Deviation):", portfolio_volatility*portfolio_sd, "\n")

## Portfolio Volatility (Standard Deviation): 0.3183642

Soal 2

2. (CPL: KK3:PP1: PP2, Poin: 50 poin, Case Study) Assume you are working as an Actuary at an Insurance Company. Your job is to define insurance claims data and analyze the risk associated with different scenarios. Consider explaining your answer in details as the following instructions:

Generate insurance claims data for a hypothetical insurance company with 10000 policies over 1# years.

# Set the seed for reproducibility
set.seed(123)

# Number of policies
policies <- 10000

# Number of years
years <- 16

# 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(
  PolicyID = integer(),
  Year = integer(),
  Claims = numeric(),
  ClaimAmount = numeric()
)

# Generate claims data for each policy and year
for (policy_id in 1:policies) {
  for (year in 1:years) {
    # Generate number of claims for the policy and year
    num_claims <- rpois(1, lambda)
    
    if (num_claims > 0) {
    # Generate claim amounts for each claim
    claim_amounts <- rnorm(num_claims, 5006, 2006)
    
    # Append data to claims_data data frame
    claims_data <- rbind(claims_data, data.frame(
      PolicyID = rep(policy_id, num_claims),
      Year = rep(year, num_claims),
      Claims = num_claims,
      ClaimAmount = claim_amounts
    ))
  }
}
}
# View the first few rows of the claims data
head(claims_data)

Generate random policy premiums from a normal distribution with a mean of 10* and a standard deviation of 20*.

set.seed(123)
premiums <- abs(rnorm((nrow(claims_data)), mean = 1006, sd = 2006))

claims_data$Premiums = premiums

head(claims_data)

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

head(claims_data)

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$ClaimAmount, claims_data$PolicyID, sum)

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

loss_ratios <- total_claim_amounts / total_premiums

# Display the first few loss ratios
ratio = data.frame(loss_ratios)
head(ratio)

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 = "red"))

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

# Show the histogram
loss_ratio_histogram

boxplot(ratio)

Calculate summary statistics including the mean, median, and 95th percentile of the loss ratios.

Note: replace the * sign with the last two digits of your Student ID number. replace the # sign with the last one digit of your Student ID number.

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

## Mean Loss Ratio: 9.421615

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

## Median Loss Ratio: 2.954068

## Median Loss Ratio: 2.954068
cat("95th Percentile of Loss Ratios:", p95_loss_ratio, "\n")

## 95th Percentile of Loss Ratios: 16.92062