Load Necessary Libraries

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## 
## The following object is masked from 'package:graphics':
## 
##     legend
library(xts)
library(quadprog)
library(readr)

Load Data

etf_data <- read.csv("myetf4.csv")

# Convert date column to Date type
colnames(etf_data)[1] <- "Date"  # Rename 'Index' to 'Date' if necessary
etf_data$Date <- as.Date(etf_data$Date, format="%Y/%m/%d")

Select and Process ETFs

etfs <- etf_data %>%
  filter(Date >= "2015-12-14" & Date <= "2018-12-28") %>%
  select(Date, `tw0050`, `tw0056`, `tw006205`, `tw00646`)

# Convert to xts for time-series operations
etfs_xts <- xts(etfs[,-1], order.by=etfs$Date)

Compute Daily Log Returns

returns <- na.omit(Return.calculate(etfs_xts, method = "log"))

Compute GMVP (Daily Returns)

cov_matrix <- cov(returns)
mean_returns <- colMeans(returns)

n <- ncol(returns)
Dmat <- 2 * cov_matrix
dvec <- rep(0, n)
Amat <- cbind(rep(1, n))
bvec <- 1
solution <- solve.QP(Dmat, dvec, Amat, bvec, meq=1)
weights_gmvp <- solution$solution

gmvp_return <- sum(weights_gmvp * mean_returns)
gmvp_sd <- sqrt(t(weights_gmvp) %*% cov_matrix %*% weights_gmvp)

list(Weights = weights_gmvp, GMVP_Return = gmvp_return, GMVP_SD = gmvp_sd)
## $Weights
## [1] -0.2241314  0.7298726  0.1080760  0.3861827
## 
## $GMVP_Return
## [1] 0.0002263853
## 
## $GMVP_SD
##             [,1]
## [1,] 0.005921611

Compute GMVP (Monthly Returns)

monthly_returns <- apply.monthly(returns, colSums)
cov_matrix_monthly <- cov(monthly_returns)
mean_returns_monthly <- colMeans(monthly_returns)

Dmat <- 2 * cov_matrix_monthly
dvec <- rep(0, n)
solution <- solve.QP(Dmat, dvec, Amat, bvec, meq=1)
weights_gmvp_monthly <- solution$solution

gmvp_return_monthly <- sum(weights_gmvp_monthly * mean_returns_monthly)
gmvp_sd_monthly <- sqrt(t(weights_gmvp_monthly) %*% cov_matrix_monthly %*% weights_gmvp_monthly)

list(Weights_Monthly = weights_gmvp_monthly, GMVP_Return_Monthly = gmvp_return_monthly, GMVP_SD_Monthly = gmvp_sd_monthly)
## $Weights_Monthly
## [1] -0.0009382175  0.4936002513 -0.0105293618  0.5178673280
## 
## $GMVP_Return_Monthly
## [1] 0.006042594
## 
## $GMVP_SD_Monthly
##            [,1]
## [1,] 0.02499232

Compute Tangency Portfolio

w_tangency <- solve.QP(Dmat, mean_returns_monthly, Amat, bvec, meq=1)$solution

tangency_return <- sum(w_tangency * mean_returns_monthly)
tangency_sd <- sqrt(t(w_tangency) %*% cov_matrix_monthly %*% w_tangency)

list(Weights_Tangency = w_tangency, Tangency_Return = tangency_return, Tangency_SD = tangency_sd)
## $Weights_Tangency
## [1]  4.925598 -1.410903 -3.657954  1.143259
## 
## $Tangency_Return
## [1] 0.05775268
## 
## $Tangency_SD
##           [,1]
## [1,] 0.1627257