h

library(quadprog)

returns_data <- matrix(c(
  0.01, 0.02, -0.01, 0.005,  # Daily returns for ETF "0050" (complete the sequence)
  0.015, 0.018, -0.012, 0.008,  # Daily returns for ETF "0056" (complete the sequence)
  0.012, 0.025, -0.015, 0.007,  # Daily returns for ETF "006205" (complete the sequence)
  0.009, 0.021, -0.008, 0.006  # Daily returns for ETF "00646" (complete the sequence)
), nrow = 4, byrow = TRUE)

colnames(returns_data) <- c("0050", "0056", "006205", "00646")

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

# Add small positive value to the diagonal of covariance matrix for regularization
epsilon <- 1e-10
cov_matrix_reg <- cov_matrix + epsilon * diag(ncol(returns_data))

# Define objective function for optimization (minimize portfolio risk)
objective_function <- function(w, cov_matrix) {
  portfolio_risk <- sqrt(t(w) %*% cov_matrix %*% w)
  return(portfolio_risk)
}

# Constraints for optimization (sum of weights equals 1)
constraints <- matrix(1, nrow = 1, ncol = ncol(returns_data))
bounds <- list(lower = rep(0, ncol(returns_data)), upper = rep(1, ncol(returns_data)))

# Perform optimization to find optimal weights
optimal_weights <- solve.QP(Dmat = 2 * cov_matrix_reg, dvec = rep(0, ncol(cov_matrix)),
                            Amat = t(constraints), bvec = 1, meq = 1)$solution

# Calculate return and standard deviation of the GMVP
GMVP_return <- sum(optimal_weights * mean_returns)
GMVP_std_dev <- sqrt(t(optimal_weights) %*% cov_matrix %*% optimal_weights)

print(optimal_weights)

## [1]  0.7336424  0.4253957  0.4679659 -0.6270040