# Load Required Libraries
library(xts)
library(quantmod)
library(zoo)
library(quadprog)
library(DEoptim)
library(copula)
# Create the Data Files
start.date <- "2014-04-01"
end.date <- "2019-04-01"
tickers <- c("NVDA","MSFT","JPM","BAC","AMT","PLD")
for (tk in tickers) {
getSymbols(tk, from = start.date, to = end.date, src = "yahoo", auto.assign = TRUE)
}
data.NVDA <- NVDA; names(data.NVDA) <- paste0("NVDA.", c("Open","High","Low","Close","Volume","Adjusted"))
data.MSFT <- MSFT; names(data.MSFT) <- paste0("MSFT.", c("Open","High","Low","Close","Volume","Adjusted"))
data.JPM <- JPM; names(data.JPM) <- paste0("JPM.", c("Open","High","Low","Close","Volume","Adjusted"))
data.BAC <- BAC; names(data.BAC) <- paste0("BAC.", c("Open","High","Low","Close","Volume","Adjusted"))
data.AMT <- AMT; names(data.AMT) <- paste0("AMT.", c("Open","High","Low","Close","Volume","Adjusted"))
data.PLD <- PLD; names(data.PLD) <- paste0("PLD.", c("Open","High","Low","Close","Volume","Adjusted"))
# Subset to Fifth Year
hw4.start <- "2018-04-01"
hw4.end <- "2019-03-31"
NVDA.hw4 <- window(data.NVDA, start = hw4.start, end = hw4.end)
MSFT.hw4 <- window(data.MSFT, start = hw4.start, end = hw4.end)
JPM.hw4 <- window(data.JPM, start = hw4.start, end = hw4.end)
BAC.hw4 <- window(data.BAC, start = hw4.start, end = hw4.end)
AMT.hw4 <- window(data.AMT, start = hw4.start, end = hw4.end)
PLD.hw4 <- window(data.PLD, start = hw4.start, end = hw4.end)
# Merge adjusted prices and compute daily log returns
prices <- merge(NVDA.hw4$NVDA.Adjusted,
MSFT.hw4$MSFT.Adjusted,
JPM.hw4$JPM.Adjusted,
BAC.hw4$BAC.Adjusted,
AMT.hw4$AMT.Adjusted,
PLD.hw4$PLD.Adjusted,
join = "inner")
prices <- na.omit(prices)
colnames(prices) <- tickers
ret <- diff(log(prices))
ret <- na.omit(ret)
ret.mat <- as.matrix(ret)
mu.vec <- colMeans(ret.mat)
vcov.mat <- cov(ret.mat)
n <- ncol(ret.mat)
cat("Dimensions of return matrix:", dim(ret.mat), "\n")
## Dimensions of return matrix: 250 6
cat("Date range:", as.character(index(ret)[1]),
"to", as.character(index(ret)[nrow(ret)]), "\n")
## Date range: 2018-04-03 to 2019-03-29
# Part A - Section 1a: Variance-Covariance Matrix
cat("=== Variance-Covariance Matrix ===\n")
## === Variance-Covariance Matrix ===
print(round(vcov.mat, 8))
## NVDA MSFT JPM BAC AMT PLD
## NVDA 0.00107213 0.00030116 0.00014893 0.00017694 0.00002460 0.00006576
## MSFT 0.00030116 0.00026564 0.00010345 0.00012199 0.00004033 0.00007461
## JPM 0.00014893 0.00010345 0.00016103 0.00016965 0.00001410 0.00003221
## BAC 0.00017694 0.00012199 0.00016965 0.00022979 0.00000731 0.00003456
## AMT 0.00002460 0.00004033 0.00001410 0.00000731 0.00012975 0.00008136
## PLD 0.00006576 0.00007461 0.00003221 0.00003456 0.00008136 0.00015975
cat("\n=== Annualized Std Deviations ===\n")
##
## === Annualized Std Deviations ===
ann.sd <- sqrt(diag(vcov.mat) * 252) * 100
print(round(ann.sd, 2))
## NVDA MSFT JPM BAC AMT PLD
## 51.98 25.87 20.14 24.06 18.08 20.06
cat("\n=== Correlation Matrix ===\n")
##
## === Correlation Matrix ===
print(round(cor(ret.mat), 4))
## NVDA MSFT JPM BAC AMT PLD
## NVDA 1.0000 0.5643 0.3584 0.3565 0.0660 0.1589
## MSFT 0.5643 1.0000 0.5002 0.4938 0.2172 0.3622
## JPM 0.3584 0.5002 1.0000 0.8819 0.0975 0.2008
## BAC 0.3565 0.4938 0.8819 1.0000 0.0423 0.1804
## AMT 0.0660 0.2172 0.0975 0.0423 1.0000 0.5651
## PLD 0.1589 0.3622 0.2008 0.1804 0.5651 1.0000
# Part A - Section 1b: Minimum Variance Portfolio
Dmat <- 2 * vcov.mat
dvec <- rep(0, n)
Amat <- cbind(rep(1, n), diag(n))
bvec <- c(1, rep(0, n))
mvp.sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
mvp.w <- mvp.sol$solution
mvp.w <- mvp.w / sum(mvp.w)
names(mvp.w) <- tickers
mvp.ret <- sum(mvp.w * mu.vec) * 252
mvp.sd <- sqrt(t(mvp.w) %*% vcov.mat %*% mvp.w) * sqrt(252)
cat("=== Minimum Variance Portfolio ===\n")
## === Minimum Variance Portfolio ===
cat("Weights:\n"); print(round(mvp.w, 4))
## Weights:
## NVDA MSFT JPM BAC AMT PLD
## 0.0000 0.0232 0.3826 0.0000 0.4283 0.1658
cat("Annualized Return:", round(mvp.ret * 100, 2), "%\n")
## Annualized Return: 16.93 %
cat("Annualized Std Dev:", round(as.numeric(mvp.sd) * 100, 2), "%\n")
## Annualized Std Dev: 13.78 %
# Part A - Section 1c: Tangency Portfolio
rf.daily <- 0.02 / 252
excess <- mu.vec - rf.daily
vcov.inv <- solve(vcov.mat)
z <- vcov.inv %*% excess
tan.w <- pmax(z, 0)
tan.w <- tan.w / sum(tan.w)
names(tan.w) <- tickers
tan.ret <- sum(tan.w * mu.vec) * 252
tan.sd <- sqrt(t(tan.w) %*% vcov.mat %*% tan.w) * sqrt(252)
tan.sr <- (tan.ret - 0.02) / as.numeric(tan.sd)
cat("=== Tangency Portfolio ===\n")
## === Tangency Portfolio ===
cat("Weights:\n"); print(round(tan.w, 4))
## Weights:
## [,1]
## NVDA 0.0000
## MSFT 0.4441
## JPM 0.0000
## BAC 0.0219
## AMT 0.5340
## PLD 0.0000
## attr(,"names")
## [1] "NVDA" "MSFT" "JPM" "BAC" "AMT" "PLD"
cat("Annualized Return:", round(tan.ret * 100, 2), "%\n")
## Annualized Return: 31.63 %
cat("Annualized Std Dev:", round(as.numeric(tan.sd) * 100, 2), "%\n")
## Annualized Std Dev: 16.74 %
cat("Sharpe Ratio:", round(tan.sr, 4), "\n")
## Sharpe Ratio: 1.77
# Part A - Section 1d: Efficient Frontier
n.points <- 100
ret.range <- seq(min(mu.vec) * 252, max(mu.vec) * 252, length.out = n.points)
front.sd <- rep(NA, n.points)
for (i in seq_along(ret.range)) {
target <- ret.range[i] / 252
Amat.f <- cbind(rep(1, n), mu.vec, diag(n))
bvec.f <- c(1, target, rep(0, n))
tryCatch({
sol <- solve.QP(Dmat, dvec, Amat.f, bvec.f, meq = 2)
w.f <- sol$solution
front.sd[i] <- sqrt(t(w.f) %*% vcov.mat %*% w.f) * sqrt(252)
}, error = function(e) NULL)
}
ind.ret <- mu.vec * 252
ind.sd <- sqrt(diag(vcov.mat)) * sqrt(252)
plot(front.sd * 100, ret.range * 100,
type = "l", col = "steelblue", lwd = 2.5,
xlim = c(0, max(ind.sd * 100) * 1.15),
ylim = c(min(c(ret.range, ind.ret)) * 100 - 5,
max(c(ret.range, ind.ret)) * 100 + 5),
xlab = "Annualized Std Dev (%)",
ylab = "Annualized Return (%)",
main = "Mean-Variance Efficient Frontier\nFifth Year (April 2018 - March 2019)",
cex.lab = 1.0)
points(ind.sd * 100, ind.ret * 100, pch = 19, col = "gray40", cex = 1.3)
text(ind.sd * 100, ind.ret * 100, labels = tickers,
pos = 4, cex = 0.75, col = "gray30")
points(as.numeric(mvp.sd) * 100, mvp.ret * 100,
pch = 18, col = "darkgreen", cex = 2.5)
text(as.numeric(mvp.sd) * 100, mvp.ret * 100,
labels = "MVP", pos = 4, cex = 0.85, col = "darkgreen", font = 2)
points(as.numeric(tan.sd) * 100, tan.ret * 100,
pch = 17, col = "firebrick", cex = 2.5)
text(as.numeric(tan.sd) * 100, tan.ret * 100,
labels = "TP", pos = 4, cex = 0.85, col = "firebrick", font = 2)
cml.sd <- seq(0, max(ind.sd * 100) * 1.15, length.out = 100)
cml.ret <- 2 + tan.sr * cml.sd
lines(cml.sd, cml.ret, lty = 2, col = "firebrick", lwd = 1.5)
legend("topleft",
legend = c("Efficient Frontier", "Capital Market Line",
"Individual Stocks", "Min Variance Portfolio",
"Tangency Portfolio"),
col = c("steelblue","firebrick","gray40","darkgreen","firebrick"),
lty = c(1, 2, NA, NA, NA),
pch = c(NA, NA, 19, 18, 17),
pt.cex = c(NA, NA, 1.3, 2, 2),
lwd = c(2.5, 1.5, NA, NA, NA),
bty = "n",
cex = 0.8)
# Part B - Section 1: Global Minimum Variance Portfolio
gmv.sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
gmv.w <- gmv.sol$solution
gmv.w <- gmv.w / sum(gmv.w)
names(gmv.w) <- tickers
gmv.ret <- sum(gmv.w * mu.vec) * 252
gmv.sd <- sqrt(t(gmv.w) %*% vcov.mat %*% gmv.w) * sqrt(252)
cat("=== Global Minimum Variance Portfolio ===\n")
## === Global Minimum Variance Portfolio ===
cat("Weights:\n"); print(round(gmv.w, 4))
## Weights:
## NVDA MSFT JPM BAC AMT PLD
## 0.0000 0.0232 0.3826 0.0000 0.4283 0.1658
cat("Annualized Return:", round(gmv.ret * 100, 2), "%\n")
## Annualized Return: 16.93 %
cat("Annualized Std Dev:", round(as.numeric(gmv.sd) * 100, 2), "%\n")
## Annualized Std Dev: 13.78 %
# Part B - Section 2: Most Diversified Portfolio
sigma.vec <- sqrt(diag(vcov.mat))
neg.dr <- function(w) {
w <- abs(w) / sum(abs(w))
port.vol <- sqrt(t(w) %*% vcov.mat %*% w)
dr <- sum(w * sigma.vec) / port.vol
return(-as.numeric(dr))
}
lower <- rep(0, n)
upper <- rep(1, n)
set.seed(123)
mdp.sol <- DEoptim(neg.dr, lower, upper,
control = DEoptim.control(itermax = 5000,
NP = 200,
trace = FALSE,
CR = 0.9,
F = 0.8))
mdp.w <- mdp.sol$optim$bestmem
mdp.w <- mdp.w / sum(mdp.w)
names(mdp.w) <- tickers
if (max(mdp.w) > 0.95) warning("MDP may be degenerate")
mdp.ret <- sum(mdp.w * mu.vec) * 252
mdp.sd <- sqrt(t(mdp.w) %*% vcov.mat %*% mdp.w) * sqrt(252)
cat("=== Most Diversified Portfolio ===\n")
## === Most Diversified Portfolio ===
cat("Weights:\n"); print(round(mdp.w, 4))
## Weights:
## NVDA MSFT JPM BAC AMT PLD
## 0.1241 0.0037 0.0691 0.2272 0.3941 0.1818
cat("Diversification Ratio:", round(-neg.dr(mdp.w), 4), "\n")
## Diversification Ratio: 1.5549
cat("Annualized Return:", round(mdp.ret * 100, 2), "%\n")
## Annualized Return: 13.17 %
cat("Annualized Std Dev:", round(as.numeric(mdp.sd) * 100, 2), "%\n")
## Annualized Std Dev: 15.55 %
# Part B - Section 3: Equal Risk Contribution Portfolio
erc.obj <- function(w) {
w <- abs(w) / sum(abs(w))
port.var <- as.numeric(t(w) %*% vcov.mat %*% w)
mrc <- as.numeric(vcov.mat %*% w)
rc <- w * mrc / port.var
return(sum((rc - mean(rc))^2))
}
set.seed(456)
erc.sol <- DEoptim(erc.obj, lower, upper,
control = DEoptim.control(itermax = 5000,
NP = 200,
trace = FALSE,
CR = 0.9,
F = 0.8))
erc.w <- erc.sol$optim$bestmem
erc.w <- erc.w / sum(erc.w)
names(erc.w) <- tickers
if (max(erc.w) > 0.95) warning("ERC may be degenerate")
erc.ret <- sum(erc.w * mu.vec) * 252
erc.sd <- sqrt(t(erc.w) %*% vcov.mat %*% erc.w) * sqrt(252)
port.var <- as.numeric(t(erc.w) %*% vcov.mat %*% erc.w)
mrc <- as.numeric(vcov.mat %*% erc.w)
rc <- erc.w * mrc / port.var
cat("=== Equal Risk Contribution Portfolio ===\n")
## === Equal Risk Contribution Portfolio ===
cat("Weights:\n"); print(round(erc.w, 4))
## Weights:
## NVDA MSFT JPM BAC AMT PLD
## 0.0800 0.1308 0.1763 0.1525 0.2602 0.2002
cat("Risk Contributions (should be ~equal):\n"); print(round(rc, 4))
## Risk Contributions (should be ~equal):
## NVDA MSFT JPM BAC AMT PLD
## 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667
cat("Annualized Return:", round(erc.ret * 100, 2), "%\n")
## Annualized Return: 13.67 %
cat("Annualized Std Dev:", round(as.numeric(erc.sd) * 100, 2), "%\n")
## Annualized Std Dev: 15.48 %
# Part B - Section 4: Minimum Tail-Dependent Portfolio
ret.mat.scaled <- apply(ret.mat, 2, rank) / (nrow(ret.mat) + 1)
tdc.matrix <- matrix(0, n, n)
threshold <- 0.10
for (i in 1:n) {
for (j in 1:n) {
if (i != j) {
u1 <- ret.mat.scaled[, i]
u2 <- ret.mat.scaled[, j]
joint.tail <- sum(u1 <= threshold & u2 <= threshold)
margin.tail <- sum(u2 <= threshold)
tdc.matrix[i, j] <- joint.tail / margin.tail
}
}
}
rownames(tdc.matrix) <- tickers
colnames(tdc.matrix) <- tickers
cat("=== Pairwise Lower Tail Dependence Coefficients ===\n")
## === Pairwise Lower Tail Dependence Coefficients ===
print(round(tdc.matrix, 4))
## NVDA MSFT JPM BAC AMT PLD
## NVDA 0.00 0.48 0.36 0.36 0.20 0.32
## MSFT 0.48 0.00 0.44 0.40 0.28 0.36
## JPM 0.36 0.44 0.00 0.80 0.24 0.28
## BAC 0.36 0.40 0.80 0.00 0.24 0.24
## AMT 0.20 0.28 0.24 0.24 0.00 0.44
## PLD 0.32 0.36 0.28 0.24 0.44 0.00
mtd.obj <- function(w) {
w <- abs(w) / sum(abs(w))
tdc.val <- as.numeric(t(w) %*% tdc.matrix %*% w)
herfindahl <- sum(w^2)
return(tdc.val + 0.5 * herfindahl)
}
upper.mtd <- rep(0.6, n)
set.seed(789)
mtd.sol <- DEoptim(mtd.obj, lower, upper.mtd,
control = DEoptim.control(itermax = 5000,
NP = 200,
trace = FALSE,
CR = 0.9,
F = 0.8))
mtd.w <- mtd.sol$optim$bestmem
mtd.w <- mtd.w / sum(mtd.w)
names(mtd.w) <- tickers
if (max(mtd.w) > 0.95) warning("MTD may be degenerate")
mtd.ret <- sum(mtd.w * mu.vec) * 252
mtd.sd <- sqrt(t(mtd.w) %*% vcov.mat %*% mtd.w) * sqrt(252)
cat("=== Minimum Tail-Dependent Portfolio ===\n")
## === Minimum Tail-Dependent Portfolio ===
cat("Weights:\n"); print(round(mtd.w, 4))
## Weights:
## NVDA MSFT JPM BAC AMT PLD
## 0.3451 0.0000 0.2212 0.0000 0.4336 0.0000
cat("Weighted Avg Tail Dependence:", round(mtd.obj(mtd.w), 4), "\n")
## Weighted Avg Tail Dependence: 0.3389
cat("Annualized Return:", round(mtd.ret * 100, 2), "%\n")
## Annualized Return: 6.75 %
cat("Annualized Std Dev:", round(as.numeric(mtd.sd) * 100, 2), "%\n")
## Annualized Std Dev: 22.04 %
# Part B - Graphical Comparison: Weight Bar Chart
all.weights <- rbind(gmv.w, mdp.w, erc.w, mtd.w)
rownames(all.weights) <- c("GMV","MDP","ERC","MTD")
port.names <- c("GMV","MDP","ERC","MTD")
port.colors <- c("steelblue","firebrick","darkgreen","darkorange")
stock.colors <- c("#4E79A7","#F28E2B","#E15759","#76B7B2","#59A14F","#EDC948")
par(mfrow = c(1,1), mar = c(5, 5, 4, 10), xpd = TRUE)
barplot(t(all.weights),
beside = TRUE,
col = stock.colors,
names.arg = port.names,
xlab = "Portfolio",
ylab = "Weight",
main = "Portfolio Weights by Stock\nFifth Year (April 2018 - March 2019)",
ylim = c(0, max(all.weights) * 1.3),
cex.names = 1.0,
cex.axis = 0.9,
border = "white")
legend("topright",
inset = c(-0.32, 0),
legend = tickers,
fill = stock.colors,
border = "white",
bty = "n",
cex = 0.85,
title = "Stock",
title.adj = 0)
par(xpd = FALSE)
# Part B - Graphical Comparison: Risk-Return Scatter
port.ret.vec <- c(gmv.ret, mdp.ret, erc.ret, mtd.ret) * 100
port.sd.vec <- c(as.numeric(gmv.sd), as.numeric(mdp.sd),
as.numeric(erc.sd), as.numeric(mtd.sd)) * 100
port.shapes <- c(16, 17, 15, 18)
ord <- order(port.sd.vec)
par(mar = c(5, 5, 4, 2))
plot(port.sd.vec, port.ret.vec,
type = "n",
xlim = c(min(port.sd.vec) - 3, max(port.sd.vec) + 3),
ylim = c(min(port.ret.vec) - 4, max(port.ret.vec) + 4),
xlab = "Annualized Std Dev (%)",
ylab = "Annualized Return (%)",
main = "Risk-Return Comparison\nFifth Year (April 2018 - March 2019)",
cex.lab = 1.0, cex.axis = 0.9, las = 1)
abline(h = seq(0, 25, by = 5), lty = 3, col = "gray88")
abline(v = seq(10, 30, by = 5), lty = 3, col = "gray88")
lines(port.sd.vec[ord], port.ret.vec[ord],
col = "gray60", lty = 2, lwd = 1.5)
points(port.sd.vec, port.ret.vec,
pch = port.shapes, col = port.colors, cex = 3)
offsets.x <- c(-0.5, 1.8, 1.8, 0.5)
offsets.y <- c( 1.5, 1.2, -1.8, 1.5)
for (i in 1:4) {
text(port.sd.vec[i] + offsets.x[i],
port.ret.vec[i] + offsets.y[i],
labels = port.names[i],
col = port.colors[i], cex = 0.95, font = 2)
}
legend("topright",
legend = paste0(port.names, " (",
round(port.ret.vec, 1), "% / ",
round(port.sd.vec, 1), "%)"),
col = port.colors,
pch = port.shapes,
pt.cex = 1.5,
bty = "n",
cex = 0.82,
title = "Portfolio (Return / Std Dev)",
title.adj = 0)
# Part B - Graphical Comparison: Pie Charts
par(mfrow = c(2, 2), mar = c(1, 1, 3, 1))
for (i in 1:4) {
w <- all.weights[i, ]
lbl <- ifelse(w > 0.02,
paste0(tickers, "\n", round(w * 100, 1), "%"), "")
pie(w, labels = lbl, col = stock.colors,
main = paste(port.names[i], "Portfolio"),
cex = 0.85, border = "white")
}
par(mfrow = c(1, 1), mar = c(5, 4, 4, 2))