# Load Required Libraries
library(xts)
library(quantmod)
library(zoo)
library(ghyp)
# Create the Data Files
start.date <- "2014-04-01"
end.date <- "2019-04-01"
getSymbols("NVDA", from = start.date, to = end.date, src = "yahoo", auto.assign = TRUE)
## [1] "NVDA"
getSymbols("BAC", from = start.date, to = end.date, src = "yahoo", auto.assign = TRUE)
## [1] "BAC"
data.NVDA <- NVDA
data.BAC <- BAC
names(data.NVDA) <- paste0("NVDA.", c("Open","High","Low","Close","Volume","Adjusted"))
names(data.BAC) <- paste0("BAC.", c("Open","High","Low","Close","Volume","Adjusted"))
# Subset to Fourth Year
hw3.start <- "2017-04-01"
hw3.end <- "2018-03-31"
NVDA.hw3 <- window(data.NVDA, start = hw3.start, end = hw3.end)
BAC.hw3 <- window(data.BAC, start = hw3.start, end = hw3.end)
# Part A - Variance-Covariance Matrix
# Daily log returns
NVDA.ret <- na.omit(diff(log(NVDA.hw3$NVDA.Adjusted)))
BAC.ret <- na.omit(diff(log(BAC.hw3$BAC.Adjusted)))
# Merge on common dates
ret.merged <- merge(NVDA.ret, BAC.ret, join = "inner")
colnames(ret.merged) <- c("NVDA", "BAC")
ret.matrix <- as.matrix(ret.merged)
# Variance-Covariance Matrix
vcov.matrix <- cov(ret.matrix)
cat("=== Variance-Covariance Matrix (Daily Log Returns, Year 4) ===\n")
## === Variance-Covariance Matrix (Daily Log Returns, Year 4) ===
print(round(vcov.matrix, 8))
## NVDA BAC
## NVDA 0.00071228 0.00013230
## BAC 0.00013230 0.00019919
cat("\nVariance of NVDA:", round(vcov.matrix["NVDA","NVDA"], 8), "\n")
##
## Variance of NVDA: 0.00071228
cat("Variance of BAC: ", round(vcov.matrix["BAC","BAC"], 8), "\n")
## Variance of BAC: 0.00019919
cat("Covariance: ", round(vcov.matrix["NVDA","BAC"], 8), "\n")
## Covariance: 0.0001323
cat("Correlation: ", round(vcov.matrix["NVDA","BAC"] /
sqrt(vcov.matrix["NVDA","NVDA"] * vcov.matrix["BAC","BAC"]), 4), "\n")
## Correlation: 0.3512
cat("\n=== Annualized Variance-Covariance Matrix ===\n")
##
## === Annualized Variance-Covariance Matrix ===
print(round(vcov.matrix * 252, 6))
## NVDA BAC
## NVDA 0.179495 0.033339
## BAC 0.033339 0.050196
cat("\nAnnualized Std Dev NVDA:", round(sqrt(vcov.matrix["NVDA","NVDA"] * 252) * 100, 2), "%\n")
##
## Annualized Std Dev NVDA: 42.37 %
cat("Annualized Std Dev BAC: ", round(sqrt(vcov.matrix["BAC","BAC"] * 252) * 100, 2), "%\n")
## Annualized Std Dev BAC: 22.4 %
# Part B - Section 1: Fit GHD, HYP, and NIG
nvda.vec <- as.numeric(ret.merged[, "NVDA"])
fit.ghd <- fit.ghypuv(nvda.vec, symmetric = FALSE, silent = TRUE)
fit.hyp <- fit.hypuv(nvda.vec, symmetric = FALSE, silent = TRUE)
fit.nig <- fit.NIGuv(nvda.vec, symmetric = FALSE, silent = TRUE)
cat("--- GHD Parameters ---\n"); print(fit.ghd)
## --- GHD Parameters ---
## Asymmetric Generalized Hyperbolic Distribution:
##
## Parameters:
## lambda alpha.bar mu sigma gamma
## -0.1925642112 0.4663137226 0.0034417053 0.0265777097 -0.0003814184
##
## log-likelihood:
## 577.2492
##
##
## Call:
## fit.ghypuv(data = nvda.vec, symmetric = FALSE, silent = TRUE)
cat("\n--- HYP Parameters ---\n"); print(fit.hyp)
##
## --- HYP Parameters ---
## Asymmetric Hyperbolic Distribution:
##
## Parameters:
## alpha.bar mu sigma gamma
## 0.0001313448 0.0024083411 0.0256551425 0.0006571446
##
## log-likelihood:
## 576.6797
##
##
## Call:
## fit.hypuv(data = nvda.vec, symmetric = FALSE, silent = TRUE)
cat("\n--- NIG Parameters ---\n"); print(fit.nig)
##
## --- NIG Parameters ---
## Asymmetric Normal Inverse Gaussian Distribution:
##
## Parameters:
## alpha.bar mu sigma gamma
## 0.4618002457 0.0035731862 0.0267362868 -0.0005072733
##
## log-likelihood:
## 577.2038
##
##
## Call:
## fit.NIGuv(data = nvda.vec, symmetric = FALSE, silent = TRUE)
# Part B - Section 2: Combined Density Plot
nvda.vec <- as.numeric(ret.merged[, "NVDA"])
x.seq <- seq(min(nvda.vec) - 0.01, max(nvda.vec) + 0.01, length.out = 500)
dens.ghd <- dghyp(x.seq, fit.ghd)
dens.hyp <- dghyp(x.seq, fit.hyp)
dens.nig <- dghyp(x.seq, fit.nig)
# Remove NAs and Inf before calculating y.max
all.dens <- c(dens.ghd, dens.hyp, dens.nig, density(nvda.vec)$y)
all.dens <- all.dens[is.finite(all.dens)]
y.max <- max(all.dens, na.rm = TRUE) * 1.05
hist(nvda.vec, breaks = 50, freq = FALSE,
col = "gray90",
border = "white",
xlab = "Daily Log Return",
main = "NVDA Daily Log Returns: GHD, HYP, and NIG Density Fits\nYear 4 (April 2017 - March 2018)",
ylim = c(0, y.max))
lines(x.seq, dens.ghd, col = "steelblue", lwd = 2)
lines(x.seq, dens.hyp, col = "firebrick", lwd = 2, lty = 2)
lines(x.seq, dens.nig, col = "darkgreen", lwd = 2, lty = 3)
lines(density(nvda.vec), col = "black", lwd = 1.5, lty = 4)
legend("topright",
legend = c("GHD", "HYP", "NIG", "Empirical Density"),
col = c("steelblue","firebrick","darkgreen","black"),
lwd = 2, lty = c(1,2,3,4), cex = 0.85)
# Part B - Section 3: Q-Q Plots
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
qqghyp(fit.ghd, plot.it = TRUE,
main = "Q-Q Plot: GHD\nNVDA Year 4",
ghyp.col = "steelblue", line = TRUE)
qqghyp(fit.hyp, plot.it = TRUE,
main = "Q-Q Plot: HYP\nNVDA Year 4",
ghyp.col = "firebrick", line = TRUE)
qqghyp(fit.nig, plot.it = TRUE,
main = "Q-Q Plot: NIG\nNVDA Year 4",
ghyp.col = "darkgreen", line = TRUE)
par(mfrow = c(1, 1))
# Part B - Section 4: Likelihood Ratio Test
lrt.ghd.hyp <- lik.ratio.test(fit.ghd, fit.hyp)
lrt.ghd.nig <- lik.ratio.test(fit.ghd, fit.nig)
cat("GHD vs HYP:\n"); print(lrt.ghd.hyp)
## GHD vs HYP:
## $statistic
## L
## 0.5657949
##
## $p.value
## [1] 0.2858537
##
## $df
## [1] 1
##
## $H0
## [1] TRUE
cat("\nGHD vs NIG:\n"); print(lrt.ghd.nig)
##
## GHD vs NIG:
## $statistic
## L
## 0.9555958
##
## $p.value
## [1] 0.7631114
##
## $df
## [1] 1
##
## $H0
## [1] TRUE
cat("\n=== Log-Likelihoods ===\n")
##
## === Log-Likelihoods ===
cat("GHD:", round(logLik(fit.ghd), 4), "\n")
## GHD: 577.2492
cat("HYP:", round(logLik(fit.hyp), 4), "\n")
## HYP: 576.6797
cat("NIG:", round(logLik(fit.nig), 4), "\n")
## NIG: 577.2038
cat("\n=== AIC (lower is better) ===\n")
##
## === AIC (lower is better) ===
cat("GHD:", round(AIC(fit.ghd), 4), "\n")
## GHD: -1144.498
cat("HYP:", round(AIC(fit.hyp), 4), "\n")
## HYP: -1145.359
cat("NIG:", round(AIC(fit.nig), 4), "\n")
## NIG: -1146.408
# Part B - Section 5: VaR Calculation and Plot
var.ghd.95 <- qghyp(0.05, fit.ghd)
var.ghd.99 <- qghyp(0.01, fit.ghd)
var.hyp.95 <- qghyp(0.05, fit.hyp)
var.hyp.99 <- qghyp(0.01, fit.hyp)
var.nig.95 <- qghyp(0.05, fit.nig)
var.nig.99 <- qghyp(0.01, fit.nig)
var.hist.95 <- as.numeric(quantile(nvda.vec, 0.05))
var.hist.99 <- as.numeric(quantile(nvda.vec, 0.01))
cat("=== Value at Risk ===\n")
## === Value at Risk ===
cat(sprintf("%-12s %10s %10s\n", "Model", "VaR 95%", "VaR 99%"))
## Model VaR 95% VaR 99%
cat(sprintf("%-12s %10.4f %10.4f\n", "GHD", var.ghd.95, var.ghd.99))
## GHD -0.0386 -0.0742
cat(sprintf("%-12s %10.4f %10.4f\n", "HYP", var.hyp.95, var.hyp.99))
## HYP -0.0383 -0.0670
cat(sprintf("%-12s %10.4f %10.4f\n", "NIG", var.nig.95, var.nig.99))
## NIG -0.0383 -0.0752
cat(sprintf("%-12s %10.4f %10.4f\n", "Historical", var.hist.95, var.hist.99))
## Historical -0.0381 -0.0715
var.data <- matrix(c(abs(var.ghd.95), abs(var.ghd.99),
abs(var.hyp.95), abs(var.hyp.99),
abs(var.nig.95), abs(var.nig.99),
abs(var.hist.95), abs(var.hist.99)),
nrow = 2, byrow = FALSE)
colnames(var.data) <- c("GHD", "HYP", "NIG", "Historical")
rownames(var.data) <- c("95%", "99%")
bp <- barplot(var.data, beside = TRUE,
col = c("steelblue","firebrick"),
ylim = c(0, max(var.data) * 1.4),
main = "VaR Comparison by Model\nNVDA Year 4 (April 2017 - March 2018)",
ylab = "VaR (Absolute Daily Log Return)",
xlab = "Model")
legend("topright",
legend = c("95% VaR", "99% VaR"),
fill = c("steelblue","firebrick"),
cex = 0.85)
text(bp, var.data + 0.001,
labels = round(var.data, 4),
cex = 0.7, pos = 3)
# Part B - Section 6: ES Calculation and Plot
es.ghd.95 <- ESghyp(0.05, fit.ghd)
es.ghd.99 <- ESghyp(0.01, fit.ghd)
es.hyp.95 <- ESghyp(0.05, fit.hyp)
es.hyp.99 <- ESghyp(0.01, fit.hyp)
es.nig.95 <- ESghyp(0.05, fit.nig)
es.nig.99 <- ESghyp(0.01, fit.nig)
es.hist.95 <- mean(nvda.vec[nvda.vec <= var.hist.95])
es.hist.99 <- mean(nvda.vec[nvda.vec <= var.hist.99])
cat("=== Expected Shortfall ===\n")
## === Expected Shortfall ===
cat(sprintf("%-12s %10s %10s\n", "Model", "ES 95%", "ES 99%"))
## Model ES 95% ES 99%
cat(sprintf("%-12s %10.4f %10.4f\n", "GHD", es.ghd.95, es.ghd.99))
## GHD -0.0609 -0.0991
cat(sprintf("%-12s %10.4f %10.4f\n", "HYP", es.hyp.95, es.hyp.99))
## HYP -0.0561 -0.0848
cat(sprintf("%-12s %10.4f %10.4f\n", "NIG", es.nig.95, es.nig.99))
## NIG -0.0616 -0.1024
cat(sprintf("%-12s %10.4f %10.4f\n", "Historical", es.hist.95, es.hist.99))
## Historical -0.0600 -0.0807
es.data <- matrix(c(abs(es.ghd.95), abs(es.ghd.99),
abs(es.hyp.95), abs(es.hyp.99),
abs(es.nig.95), abs(es.nig.99),
abs(es.hist.95), abs(es.hist.99)),
nrow = 2, byrow = FALSE)
colnames(es.data) <- c("GHD", "HYP", "NIG", "Historical")
rownames(es.data) <- c("95%", "99%")
bp2 <- barplot(es.data, beside = TRUE,
col = c("steelblue","firebrick"),
ylim = c(0, max(es.data) * 1.4),
main = "Expected Shortfall (ES) Comparison by Model\nNVDA Year 4 (April 2017 - March 2018)",
ylab = "ES (Absolute Daily Log Return)",
xlab = "Model")
legend("topright",
legend = c("95% ES", "99% ES"),
fill = c("steelblue","firebrick"),
cex = 0.85)
text(bp2, es.data + 0.001,
labels = round(es.data, 4),
cex = 0.7, pos = 3)
