Data Preparation
library(tidyverse)
library(lubridate)
library(quadprog)
library(knitr)
library(kableExtra)
# Load the data
myetf4 <- read.csv("C:\\Users\\Admin\\Downloads\\myetf4.csv", stringsAsFactors = FALSE)
colnames(myetf4) <- c("Date", "tw0050", "tw0056", "tw006205", "tw00646")
myetf4$Date <- as.Date(myetf4$Date)
# Filter in-sample: 2015/12/14 - 2018/12/28
insample <- myetf4 %>%
filter(Date >= as.Date("2015-12-14") & Date <= as.Date("2018-12-28")) %>%
arrange(Date)
cat("In-sample rows:", nrow(insample), "\n")
## In-sample rows: 751
cat("Date range:", format(min(insample$Date)), "to", format(max(insample$Date)), "\n")
## Date range: 2015-12-14 to 2018-12-28
head(insample) %>%
kable(caption = "First 6 rows of in-sample data") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
First 6 rows of in-sample data
|
Date
|
tw0050
|
tw0056
|
tw006205
|
tw00646
|
|
2015-12-14
|
53.29
|
18.25
|
31.06
|
19.61
|
|
2015-12-15
|
53.33
|
18.38
|
31.59
|
19.63
|
|
2015-12-16
|
54.14
|
18.56
|
31.60
|
19.89
|
|
2015-12-17
|
54.77
|
18.81
|
32.23
|
20.05
|
|
2015-12-18
|
54.50
|
18.95
|
32.18
|
19.85
|
|
2015-12-21
|
54.41
|
19.02
|
33.00
|
19.64
|
Compute Daily Returns
prices_daily <- insample[, c("tw0050", "tw0056", "tw006205", "tw00646")]
ret_daily <- as.data.frame(apply(prices_daily, 2, function(x) diff(x) / head(x, -1)))
ret_daily <- ret_daily * 100
cat("Number of daily return observations:", nrow(ret_daily), "\n")
## Number of daily return observations: 750
summary(ret_daily) %>%
kable(caption = "Summary of Daily Returns (%)") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Summary of Daily Returns (%)
|
|
tw0050
|
tw0056
|
tw006205
|
tw00646
|
|
|
Min. :-7.03406
|
Min. :-5.43041
|
Min. :-6.83477
|
Min. :-4.10000
|
|
|
1st Qu.:-0.40630
|
1st Qu.:-0.28541
|
1st Qu.:-0.53456
|
1st Qu.:-0.30140
|
|
|
Median : 0.07183
|
Median : 0.08187
|
Median : 0.00000
|
Median : 0.04148
|
|
|
Mean : 0.04632
|
Mean : 0.03846
|
Mean :-0.02118
|
Mean : 0.02554
|
|
|
3rd Qu.: 0.50863
|
3rd Qu.: 0.41763
|
3rd Qu.: 0.48307
|
3rd Qu.: 0.41841
|
|
|
Max. : 3.73230
|
Max. : 2.48353
|
Max. : 6.00775
|
Max. : 5.63087
|
Compute GMVP Weights (Daily)
The Global Minimum Variance Portfolio (GMVP)
minimizes portfolio variance:
\[\min_{\mathbf{w}} \; \mathbf{w}^\top
\Sigma \mathbf{w} \quad \text{s.t.} \quad \mathbf{1}^\top \mathbf{w} =
1, \quad \mathbf{w} \geq 0\]
Sigma_daily <- cov(ret_daily)
cat("Covariance Matrix (Daily Returns):\n")
## Covariance Matrix (Daily Returns):
round(Sigma_daily, 6) %>%
kable(caption = "Covariance Matrix – Daily Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Covariance Matrix – Daily Returns
|
|
tw0050
|
tw0056
|
tw006205
|
tw00646
|
|
tw0050
|
0.783706
|
0.455916
|
0.446726
|
0.366339
|
|
tw0056
|
0.455916
|
0.452641
|
0.267367
|
0.235354
|
|
tw006205
|
0.446726
|
0.267367
|
1.304184
|
0.291037
|
|
tw00646
|
0.366339
|
0.235354
|
0.291037
|
0.590289
|
n <- ncol(ret_daily)
Dmat <- 2 * Sigma_daily
dvec <- rep(0, n)
Amat <- cbind(rep(1, n), diag(n))
bvec <- c(1, rep(0, n))
meq <- 1
sol_daily <- solve.QP(Dmat, dvec, Amat, bvec, meq = meq)
w_daily <- sol_daily$solution
names(w_daily) <- colnames(ret_daily)
cat("\nGMVP Weights (Daily):\n")
##
## GMVP Weights (Daily):
data.frame(ETF = names(w_daily),
Weight = paste0(round(w_daily * 100, 4), "%")) %>%
kable(caption = "GMVP Optimal Weights – Daily Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
GMVP Optimal Weights – Daily Returns
|
ETF
|
Weight
|
|
tw0050
|
0%
|
|
tw0056
|
57.183%
|
|
tw006205
|
8.3682%
|
|
tw00646
|
34.4487%
|
GMVP Return and Standard Deviation (Daily)
mu_daily <- colMeans(ret_daily)
gmvp_ret_d <- sum(w_daily * mu_daily)
gmvp_var_d <- as.numeric(t(w_daily) %*% Sigma_daily %*% w_daily)
gmvp_sd_d <- sqrt(gmvp_var_d)
cat(sprintf("GMVP Daily Return : %.6f%%\n", gmvp_ret_d))
## GMVP Daily Return : 0.029021%
cat(sprintf("GMVP Daily Std Dev : %.6f%%\n", gmvp_sd_d))
## GMVP Daily Std Dev : 0.601901%
data.frame(
Metric = c("Expected Daily Return (%)", "Daily Std Dev (%)"),
Value = round(c(gmvp_ret_d, gmvp_sd_d), 6)
) %>%
kable(caption = "GMVP Performance – Daily Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
GMVP Performance – Daily Returns
|
Metric
|
Value
|
|
Expected Daily Return (%)
|
0.029021
|
|
Daily Std Dev (%)
|
0.601901
|
Q2 – GMVP Using Monthly Returns
Aggregate to Monthly Returns
insample_monthly <- insample %>%
mutate(YM = format(Date, "%Y-%m")) %>%
group_by(YM) %>%
slice_tail(n = 1) %>%
ungroup() %>%
arrange(Date)
prices_monthly <- insample_monthly[, c("tw0050", "tw0056", "tw006205", "tw00646")]
ret_monthly <- as.data.frame(apply(prices_monthly, 2, function(x) diff(x) / head(x, -1)))
ret_monthly <- ret_monthly * 100
cat("Number of monthly return observations:", nrow(ret_monthly), "\n")
## Number of monthly return observations: 36
summary(ret_monthly) %>%
kable(caption = "Summary of Monthly Returns (%)") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Summary of Monthly Returns (%)
|
|
tw0050
|
tw0056
|
tw006205
|
tw00646
|
|
|
Min. :-10.760
|
Min. :-7.6505
|
Min. :-17.30709
|
Min. :-8.7947
|
|
|
1st Qu.: -1.085
|
1st Qu.:-1.4821
|
1st Qu.: -2.97001
|
1st Qu.:-0.5192
|
|
|
Median : 1.195
|
Median : 0.5376
|
Median : -0.07805
|
Median : 1.1126
|
|
|
Mean : 0.882
|
Mean : 0.7087
|
Mean : -0.53555
|
Mean : 0.4511
|
|
|
3rd Qu.: 2.539
|
3rd Qu.: 2.3660
|
3rd Qu.: 2.48567
|
3rd Qu.: 2.3882
|
|
|
Max. : 6.047
|
Max. : 7.8421
|
Max. : 8.55715
|
Max. : 4.0964
|
Compute GMVP Weights (Monthly)
Sigma_monthly <- cov(ret_monthly)
cat("Covariance Matrix (Monthly Returns):\n")
## Covariance Matrix (Monthly Returns):
round(Sigma_monthly, 6) %>%
kable(caption = "Covariance Matrix – Monthly Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Covariance Matrix – Monthly Returns
|
|
tw0050
|
tw0056
|
tw006205
|
tw00646
|
|
tw0050
|
11.751458
|
8.661004
|
8.472189
|
3.928466
|
|
tw0056
|
8.661004
|
9.080806
|
5.553289
|
3.572509
|
|
tw006205
|
8.472189
|
5.553289
|
24.412877
|
6.736296
|
|
tw00646
|
3.928466
|
3.572509
|
6.736296
|
8.605161
|
Dmat_m <- 2 * Sigma_monthly
dvec_m <- rep(0, n)
Amat_m <- cbind(rep(1, n), diag(n))
bvec_m <- c(1, rep(0, n))
sol_monthly <- solve.QP(Dmat_m, dvec_m, Amat_m, bvec_m, meq = 1)
w_monthly <- sol_monthly$solution
names(w_monthly) <- colnames(ret_monthly)
cat("\nGMVP Weights (Monthly):\n")
##
## GMVP Weights (Monthly):
data.frame(ETF = names(w_monthly),
Weight = paste0(round(w_monthly * 100, 4), "%")) %>%
kable(caption = "GMVP Optimal Weights – Monthly Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
GMVP Optimal Weights – Monthly Returns
|
ETF
|
Weight
|
|
tw0050
|
0.3184%
|
|
tw0056
|
47.4049%
|
|
tw006205
|
0.1204%
|
|
tw00646
|
52.1563%
|
GMVP Return and Standard Deviation (Monthly)
mu_monthly <- colMeans(ret_monthly)
gmvp_ret_m <- sum(w_monthly * mu_monthly)
gmvp_var_m <- as.numeric(t(w_monthly) %*% Sigma_monthly %*% w_monthly)
gmvp_sd_m <- sqrt(gmvp_var_m)
cat(sprintf("GMVP Monthly Return : %.6f%%\n", gmvp_ret_m))
## GMVP Monthly Return : 0.573367%
cat(sprintf("GMVP Monthly Std Dev : %.6f%%\n", gmvp_sd_m))
## GMVP Monthly Std Dev : 2.490441%
data.frame(
Metric = c("Expected Monthly Return (%)", "Monthly Std Dev (%)"),
Value = round(c(gmvp_ret_m, gmvp_sd_m), 6)
) %>%
kable(caption = "GMVP Performance – Monthly Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
GMVP Performance – Monthly Returns
|
Metric
|
Value
|
|
Expected Monthly Return (%)
|
0.573367
|
|
Monthly Std Dev (%)
|
2.490441
|
Q3 – Tangency Portfolio (Monthly Returns, Rf = 0)
The Tangency Portfolio maximizes the Sharpe
Ratio:
\[\max_{\mathbf{w}} \;
\frac{\mathbf{w}^\top \boldsymbol{\mu} - r_f}{\sqrt{\mathbf{w}^\top
\Sigma \mathbf{w}}} \quad \text{s.t.} \quad \mathbf{1}^\top \mathbf{w} =
1, \quad \mathbf{w} \geq 0\]
With \(r_f = 0\), the closed-form
solution is:
\[\mathbf{z} = \Sigma^{-1}
\boldsymbol{\mu}, \quad \mathbf{w}^* = \frac{\mathbf{z}}{\mathbf{1}^\top
\mathbf{z}}\]
rf <- 0
excess_mu <- mu_monthly - rf
Sigma_inv <- solve(Sigma_monthly)
z_vec <- Sigma_inv %*% excess_mu
w_tan <- z_vec / sum(z_vec)
names(w_tan) <- colnames(ret_monthly)
cat("Tangency Portfolio Weights:\n")
## Tangency Portfolio Weights:
data.frame(ETF = names(w_tan),
Weight = paste0(round(w_tan * 100, 4), "%")) %>%
kable(caption = "Tangency Portfolio Weights – Monthly Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Tangency Portfolio Weights – Monthly Returns
|
ETF
|
Weight
|
|
tw0050
|
130.5054%
|
|
tw0056
|
-15.7681%
|
|
tw006205
|
-84.7532%
|
|
tw00646
|
70.0159%
|
if (any(w_tan < -1e-8)) {
cat("Negative weights detected -> applying constrained Tangency via quadprog.\n\n")
Dmat_t <- 2 * Sigma_monthly
dvec_t <- rep(0, n)
Amat_t <- cbind(mu_monthly, diag(n))
bvec_t <- c(1, rep(0, n))
sol_tan <- solve.QP(Dmat_t, dvec_t, Amat_t, bvec_t, meq = 1)
w_tan_raw <- sol_tan$solution
w_tan_c <- w_tan_raw / sum(w_tan_raw)
names(w_tan_c) <- colnames(ret_monthly)
cat("Constrained Tangency Portfolio Weights:\n")
data.frame(ETF = names(w_tan_c),
Weight = paste0(round(w_tan_c * 100, 4), "%")) %>%
kable(caption = "Constrained Tangency Weights") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
w_tan_final <- w_tan_c
} else {
w_tan_final <- as.vector(w_tan)
names(w_tan_final) <- colnames(ret_monthly)
}
## Negative weights detected -> applying constrained Tangency via quadprog.
##
## Constrained Tangency Portfolio Weights:
Tangency Portfolio Return, Std Dev, and Sharpe Ratio
tan_ret <- sum(w_tan_final * mu_monthly)
tan_var <- as.numeric(t(w_tan_final) %*% Sigma_monthly %*% w_tan_final)
tan_sd <- sqrt(tan_var)
tan_sharpe <- (tan_ret - rf) / tan_sd
cat(sprintf("Tangency Monthly Return : %.6f%%\n", tan_ret))
## Tangency Monthly Return : 0.758289%
cat(sprintf("Tangency Monthly Std Dev : %.6f%%\n", tan_sd))
## Tangency Monthly Std Dev : 2.860317%
cat(sprintf("Tangency Sharpe Ratio : %.6f\n", tan_sharpe))
## Tangency Sharpe Ratio : 0.265107
data.frame(
Metric = c("Expected Monthly Return (%)", "Monthly Std Dev (%)", "Sharpe Ratio (Rf=0)"),
Value = round(c(tan_ret, tan_sd, tan_sharpe), 6)
) %>%
kable(caption = "Tangency Portfolio Performance – Monthly Returns") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Tangency Portfolio Performance – Monthly Returns
|
Metric
|
Value
|
|
Expected Monthly Return (%)
|
0.758289
|
|
Monthly Std Dev (%)
|
2.860317
|
|
Sharpe Ratio (Rf=0)
|
0.265107
|
Summary Comparison
data.frame(
Portfolio = c("GMVP (Daily)", "GMVP (Monthly)", "Tangency (Monthly)"),
tw0050 = paste0(round(c(w_daily[1], w_monthly[1], w_tan_final[1]) * 100, 2), "%"),
tw0056 = paste0(round(c(w_daily[2], w_monthly[2], w_tan_final[2]) * 100, 2), "%"),
tw006205 = paste0(round(c(w_daily[3], w_monthly[3], w_tan_final[3]) * 100, 2), "%"),
tw00646 = paste0(round(c(w_daily[4], w_monthly[4], w_tan_final[4]) * 100, 2), "%"),
Return = paste0(round(c(gmvp_ret_d, gmvp_ret_m, tan_ret), 6), "%"),
Std_Dev = paste0(round(c(gmvp_sd_d, gmvp_sd_m, tan_sd), 6), "%")
) %>%
kable(caption = "Portfolio Comparison Summary") %>%
kable_styling(bootstrap_options = c("striped", "hover", "bordered"), full_width = TRUE)
Portfolio Comparison Summary
|
Portfolio
|
tw0050
|
tw0056
|
tw006205
|
tw00646
|
Return
|
Std_Dev
|
|
GMVP (Daily)
|
0%
|
57.18%
|
8.37%
|
34.45%
|
0.029021%
|
0.601901%
|
|
GMVP (Monthly)
|
0.32%
|
47.4%
|
0.12%
|
52.16%
|
0.573367%
|
2.490441%
|
|
Tangency (Monthly)
|
60.49%
|
18.07%
|
0%
|
21.44%
|
0.758289%
|
2.860317%
|
Efficient Frontier Plot (Monthly)
set.seed(123)
N_sim <- 5000
sim_ret <- numeric(N_sim)
sim_sd <- numeric(N_sim)
for (i in seq_len(N_sim)) {
w_rand <- runif(n)
w_rand <- w_rand / sum(w_rand)
sim_ret[i] <- sum(w_rand * mu_monthly)
sim_sd[i] <- sqrt(as.numeric(t(w_rand) %*% Sigma_monthly %*% w_rand))
}
plot(sim_sd, sim_ret,
pch = 16, cex = 0.3, col = "lightblue",
xlab = "Monthly Std Dev (%)", ylab = "Monthly Return (%)",
main = "Efficient Frontier – Taiwan ETFs (Monthly Returns)")
points(gmvp_sd_m, gmvp_ret_m, pch = 17, col = "red", cex = 1.8)
text(gmvp_sd_m, gmvp_ret_m, labels = "GMVP", pos = 4, col = "red", cex = 0.9)
points(tan_sd, tan_ret, pch = 15, col = "darkgreen", cex = 1.8)
text(tan_sd, tan_ret, labels = "Tangency", pos = 4, col = "darkgreen", cex = 0.9)
for (j in seq_len(n)) {
points(sqrt(Sigma_monthly[j, j]), mu_monthly[j], pch = 19, col = "orange", cex = 1.5)
text(sqrt(Sigma_monthly[j, j]), mu_monthly[j],
labels = colnames(ret_monthly)[j], pos = 3, cex = 0.8, col = "orange")
}
legend("topleft",
legend = c("Random Portfolios", "GMVP", "Tangency", "Individual ETFs"),
pch = c(16, 17, 15, 19),
col = c("lightblue", "red", "darkgreen", "orange"),
cex = 0.85, bty = "n")
