library(readxl)
library(ggplot2)
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(zoo)
## Warning: package 'zoo' was built under R version 4.3.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(TTR)
## Warning: package 'TTR' was built under R version 4.3.3
library(e1071)
## Warning: package 'e1071' was built under R version 4.3.3
library(sn)
## Warning: package 'sn' was built under R version 4.3.3
## Loading required package: stats4
##
## Attaching package: 'sn'
## The following object is masked from 'package:stats':
##
## sd
library(moments)
##
## Attaching package: 'moments'
## The following objects are masked from 'package:e1071':
##
## kurtosis, moment, skewness
library(PearsonDS)
## Warning: package 'PearsonDS' was built under R version 4.3.3
library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.3.3
FINM_8100_DATA_1_ <- read_excel("C:/Users/thech/Downloads/FINM 8100 DATA (1).xlsx")
## New names:
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...11`
## • `` -> `...13`
## • `` -> `...15`
## • `` -> `...17`
## • `` -> `...19`
## • `` -> `...21`
## • `` -> `...23`
## • `` -> `...25`
## • `` -> `...27`
## • `` -> `...29`
## • `` -> `...31`
## • `` -> `...32`
#################Commonwealth Bank#########################
commonwealth<-FINM_8100_DATA_1_[,2:3]
commonwealth[1,1]<-"year"
colnames(commonwealth) <- as.character(unlist(commonwealth[1, ]))
cm<-commonwealth[-1, ]
cmdotcom<-cm%>%filter(year>=20000301 & year <=20030331)
cmm <- mean(as.numeric(cmdotcom$Change.cm), na.rm = TRUE)
cms <- sd(as.numeric(cmdotcom$Change.cm), na.rm = TRUE)
hist(as.numeric(cmdotcom$Change.cm), breaks = 100, freq = FALSE, main = " Distribution of Commonwealth in Dotcom", xlab = " Daily return")
# Add normal distribution curve
curve(
dnorm(x, mean = cmm, sd = cms),
col = "lightblue", lwd = 2,
add = TRUE
)
# Add true distribution curve
dens <- density(as.numeric(cmdotcom$Change.cm), bw = 0.1, na.rm = TRUE)
lines(dens, col = "blue", lwd = 2)
##################Variance-Covariance Approach
cmchange <- as.numeric(cmdotcom$Change.cm)
##50 days
#----------EQUALLY WEIGHTED
win50 <- 50
alpha <- 0.01
z <- qnorm(alpha)
cmmu_roll50 <- rollapply(cmchange, width = win50, FUN = mean, align = "right", fill = NA, na.rm = TRUE)
cmsd_roll50 <- rollapply(cmchange, width = win50, FUN = sd, align = "right", fill = NA, na.rm = TRUE)
cmVaR50d_99 <- (cmmu_roll50 + z * cmsd_roll50)
##120 days
win120 <- 120
alpha <- 0.01
z <- qnorm(alpha)
cmmu_roll120 <- rollapply(cmchange, width = win120, FUN = mean, align = "right", fill = NA, na.rm = TRUE)
cmsd_roll120 <- rollapply(cmchange, width = win120, FUN = sd, align = "right", fill = NA, na.rm = TRUE)
cmVaR120d_99 <- (cmmu_roll120 + z * cmsd_roll120)
##250 days
win250 <- 250
alpha <- 0.01
z <- qnorm(alpha)
cmmu_roll250 <- rollapply(cmchange, width = win250, FUN = mean, align = "right", fill = NA, na.rm = TRUE)
cmsd_roll250 <- rollapply(cmchange, width = win250, FUN = sd, align = "right", fill = NA, na.rm = TRUE)
cmVaR250d_99 <- (cmmu_roll250 + z * cmsd_roll250)
#plot cm
# Make sure you have a Date column
dates <- as.Date(as.character(cmdotcom$year), format = "%Y%m%d")
# Base plot with dates on x-axis
plot(dates, cmchange, type = "l", col = "grey",
main = "Commonwealth Bank Daily Returns with 99% VaR, equally weighted",
xlab = "Date", ylab = "Return & VaR")
# Overlay the VaR lines
lines(dates, cmVaR50d_99, col = "black")
lines(dates, cmVaR120d_99, col = "blue")
lines(dates, cmVaR250d_99, col = "purple")
# Format x-axis: tick every 6 months
ticks <- as.Date(c("2001-01-01", "2002-01-01", "2003-01-01"))
# Draw axis
axis.Date(1, at = ticks, format = "%Y")
#--------------exponentially weighted
library(TTR)
lam50 <- 49/51
lam120 <- 119/121
lam250 <- 249/251
# Returns
rcm <- as.numeric(cmchange)
alpha <- 0.01
z <- qnorm(alpha)
# choose lambda (decay) directly, e.g. 0.94 is common
lambda <- 0.94
# Step 1: de-mean returns (you can also just use r if mean≈0)
eps <- rcm - mean(rcm, na.rm = TRUE)
# Step 2: compute EWMA variance with EMA() on squared eps
ewma_var50 <- EMA(eps^2, ratio = 1 - lam50)
ewma_var120 <- EMA(eps^2, ratio = 1 - lam120)
ewma_var250 <- EMA(eps^2, ratio = 1 - lam250)
cmmu_roll50 <- rollapply(cmchange, width = win50, FUN = mean, align = "right", fill = NA, na.rm = TRUE)
cmmu_roll120 <- rollapply(cmchange, width = win120, FUN = mean, align = "right", fill = NA, na.rm = TRUE)
cmmu_roll250 <- rollapply(cmchange, width = win250, FUN = mean, align = "right", fill = NA, na.rm = TRUE)
# Step 3: EWMA volatility
ewma_sd50 <- sqrt(ewma_var50)
ewma_sd120 <- sqrt(ewma_var120)
ewma_sd250 <- sqrt(ewma_var250)
# Step 4: VaR series
VaR_ewma50 <- cmmu_roll50 + z * ewma_sd50
VaR_ewma120 <- cmmu_roll120 + z * ewma_sd120
VaR_ewma250 <- cmmu_roll250 + z * ewma_sd250
# Plot
dates <- as.Date(as.character(cmdotcom$year), format = "%Y%m%d")
# Base plot with dates on x-axis
plot(dates, cmchange, type = "l", col = "grey",
main = "Commonwealth Bank Daily Returns with 99% VaR,exponentially weighted",
xlab = "Date", ylab = "Return & VaR")
lines(dates,VaR_ewma50, col = "black")
lines(dates, VaR_ewma120, col = "blue")
lines(dates, VaR_ewma250, col = "purple")
# Format x-axis: tick every 6 months
#plot cm
# Make sure you have a Date column
ticks <- as.Date(c("2001-01-01", "2002-01-01", "2003-01-01"))
axis.Date(1, at = ticks, format = "%Y")
################## Historical Approach
#50days
win50 <- 50
alpha <- 0.01
cmVaR50h_99 <- rollapply(cmchange, width = win50, FUN = function(cmchange) quantile(cmchange, probs = alpha, na.rm = TRUE), align = "right", fill = NA)
#120 days
win120 <- 120
alpha <- 0.01
cmVaR120h_99 <- rollapply(cmchange, width = win120, FUN = function(cmchange) quantile(cmchange, probs = alpha, na.rm = TRUE), align = "right", fill = NA)
#250 days
win250 <- 250
alpha <- 0.01
cmVaR250h_99 <- rollapply(cmchange, width = win250, FUN = function(cmchange) quantile(cmchange, probs = alpha, na.rm = TRUE), align = "right", fill = NA)
dates <- as.Date(as.character(cmdotcom$year), format = "%Y%m%d")
# Base plot with dates on x-axis
plot(dates, cmchange, type = "l", col = "grey",
main = "Commonwealth Bank Daily Returns with 99% VaR, historical",
xlab = "Date", ylab = "Return & VaR")
lines(dates, cmVaR50h_99, col = "black")
lines(dates, cmVaR120h_99,col = "blue")
lines(dates, cmVaR250h_99, col = "purple")
# Format x-axis: tick every 6 months
ticks <- as.Date(c("2001-01-01", "2002-01-01", "2003-01-01"))
# Draw axis
axis.Date(1, at = ticks, format = "%Y")
####################Monte Carlo Approach
stackr<-FINM_8100_DATA_1_[,2:9]
stackr[1,1]<-"year"
colnames(stackr) <- as.character(unlist(stackr[1, ]))
stackr<-stackr[-1, ]
dotcom<-stackr%>%filter(year>=20000301 & year <=20030331)
coldotcom <- as.numeric(c(as.matrix(dotcom[2:8])))
skewness(coldotcom, na.rm = TRUE)
## [1] -1.379403
stackr <- as.numeric(stackr$Change.cm)
stackr <- stackr[is.finite(stackr)]
coldotcom <- as.numeric(coldotcom)
coldotcom <- coldotcom[is.finite(coldotcom)]
cmm<- mean(rcm,na.rm = TRUE)
cms<- sd(rcm, na.rm = TRUE)
cmsk<- skewness(coldotcom, na.rm = TRUE)
cmk<- kurtosis(coldotcom, na.rm = TRUE)
# install.packages("PearsonDS") # run once
## --- your inputs (clean vectors) ---
# bank returns (for mean & sd)
bank <- as.numeric(stackr) # e.g., stackr <- as.numeric(stackr$Change.cm)
bank <- bank[is.finite(bank)]
# industry pooled returns (for skew & kurt)
ind <- as.numeric(coldotcom)
ind <- ind[is.finite(ind)]
## --- target moments -----------------
mu <- mean(cmchange, na.rm = TRUE) # mean
sdv <- stats::sd(rcm, na.rm = TRUE) # sd (use stats::sd to avoid masking)
sk <- skewness(ind, na.rm = TRUE) # skewness
ku <- kurtosis(ind, na.rm = TRUE) # **regular** kurtosis (not excess)
# (optional) tiny nudge to satisfy Pearson existence condition: beta2 > beta1 + 1
if (ku <= sk^2 + 1) ku <- sk^2 + 1 + 1e-6
## --- fit Pearson distribution from moments ---
fitPRS<- pearsonFitM(moments = list(mean = mu, variance = sdv^2,
skewness = sk, kurtosis = ku))
## --- Monte Carlo sampling (choose N) ---
set.seed(123)
N <- nrow(dotcom)
x_sim <- rpearson(N, params = fitPRS)
## quick check
c(mean = mean(x_sim), sd = sd(x_sim),
skew = skewness(x_sim), kurt = kurtosis(x_sim))
## mean sd skew kurt
## -0.0008058715 0.1103552042 -1.5876367207 11.7074746594
## example plot
hist(x_sim, breaks = 80, freq = FALSE, col = "lightblue",
main = "PearsonDS Monte Carlo")
lines(density(x_sim, bw = 1, na.rm = TRUE), lwd = 2)
#50days
win50 <- 50
alpha <- 0.01
cmVaR50m_99 <- rollapply(x_sim, width = win50, FUN = function(x_sim) quantile(x_sim, probs = alpha, na.rm = TRUE), align = "right", fill = NA)
#120 days
win120 <- 120
alpha <- 0.01
cmVaR120m_99 <- rollapply(x_sim, width = win120, FUN = function(x_sim) quantile(x_sim, probs = alpha, na.rm = TRUE), align = "right", fill = NA)
#250 days
win250 <- 250
alpha <- 0.01
cmVaR250m_99 <- rollapply(x_sim, width = win250, FUN = function(x_sim) quantile(x_sim, probs = alpha, na.rm = TRUE), align = "right", fill = NA)
dates <- as.Date(as.character(cmdotcom$year), format = "%Y%m%d")
# Base plot with dates on x-axis
plot(dates, x_sim, type = "l", col = "grey",
main = "Commonwealth Bank Daily Returns with 99% VaR, monte-carlo",
xlab = "Date", ylab = "Return & VaR")
lines(dates, cmVaR50m_99, col = "black")
lines(dates, cmVaR120m_99,col = "blue")
lines(dates, cmVaR250m_99, col = "purple")
ticks <- as.Date(c("2001-01-01", "2002-01-01", "2003-01-01"))
# Draw axis
axis.Date(1, at = ticks, format = "%Y")
# install.packages("writexl") # run once
library(writexl)
## Warning: package 'writexl' was built under R version 4.3.3
# 1) Dates (from your yyyymmdd column)
dates <- as.Date(as.character(cmdotcom$year), format = "%Y%m%d")
# 2) Map "pretty column name" -> "object name in R"
var_map <- c(
VarCov_50 = "cmVaR50d_99",
VarCov_120 = "cmVaR120d_99",
VarCov_250 = "cmVaR250d_99",
EWMA_50 = "VaR_ewma50",
EWMA_120 = "VaR_ewma120",
EWMA_250 = "VaR_ewma250",
Hist_50 = "cmVaR50h_99",
Hist_120 = "cmVaR120h_99",
Hist_250 = "cmVaR250h_99",
MC_50 = "cmVaR50m_99",
MC_120 = "cmVaR120m_99",
MC_250 = "cmVaR250m_99"
)
# 3) Start with Date + Return
cols <- list(
Date = dates,
Return = as.numeric(cmchange)
)
# 4) Add each VaR column if it exists
for (pretty_name in names(var_map)) {
obj_name <- var_map[[pretty_name]]
if (exists(obj_name, inherits = FALSE)) {
cols[[pretty_name]] <- get(obj_name)
}
}
# 5) Build data frame & (optionally) drop all-NA columns
VaR_export <- as.data.frame(cols)
VaR_export <- VaR_export[ , !sapply(VaR_export, function(v) all(!is.finite(v) | is.na(v))), drop = FALSE]
# 6) Write to Excel (single sheet)
#write_xlsx(VaR_export, path = "Commonwealth_VaR_Export.xlsx")
# Where is it saved?
#getwd() # prints the folder; that's where the file is
#--------------------------------------Test----------------------------------------------#
#####count
valid <- is.finite(cmVaR50d_99) & is.finite(cmchange)
exceptionse50 <- (cmchange < cmVaR50d_99) & valid
ce50<-sum(exceptionse50, na.rm = TRUE)
pe50<-ce50/sum(valid)
# z-Test at 99%
thres99<- 0.01
z<- (pe50-thres99)/sqrt(thres99*(1-thres99)/sum(valid))
cn<-qnorm(1-thres99)
z
## [1] 2.359221
cn
## [1] 2.326348
#--> reject null: more than 0.01 -> more bleached than expected or underestimate risks but not much
#####magnitude
exceptionse50 <-cmchange < cmVaR50d_99
loss_exc <- cmchange[exceptionse50]
var_exc <- cmVaR50d_99[exceptionse50]
ok <- is.finite(loss_exc) & is.finite(var_exc)
d <- var_exc[ok] - loss_exc[ok] # differences
n <- length(d)
mean_d <- mean(d, na.rm = TRUE)
sd_d <- sd(d, na.rm = TRUE)
t <- mean_d / (sd_d / sqrt(n)) # test statistic
ct<-qt(1-thres99,n-1)
t
## [1] 3.17601
ct
## [1] 2.650309
#--> rejest null: the difference is not 0 -> exceeded by large
#####loss function
y <- as.numeric(cmchange)
q <- cmVaR50d_99
thres99 <- 0.01
# Loss function
s <- ((1-thres99) - (y < q)) * (y - q)
loss<-sum(s,na.rm = TRUE)
loss
## [1] 184.8002
alpha <- 0.01
# 1) Backtest helper for ONE VaR series
backtest_var <- function(y, q, alpha = 0.01) {
valid <- is.finite(y) & is.finite(q)
yv <- y[valid]; qv <- q[valid]
n <- length(yv)
if (n == 0) return(data.frame(n=0, exceptions=NA, p_exc=NA, z=NA, zcrit=NA,
t_short=NA, tcrit=NA,
sum_sf=NA, total_loss=NA, mean_loss=NA))
# exceptions (Kupiec-style z)
exc <- (yv < qv)
ce <- sum(exc)
pe <- ce / n
z <- (pe - alpha) / sqrt(alpha*(1-alpha)/n)
zc <- qnorm(1 - alpha)
# shortfall (how far below VaR)
sf <- qv[exc] - yv[exc] # positive = exceedance magnitude
if (length(sf) >= 3) {
tstat <- mean(sf) / (sd(sf) / sqrt(length(sf)))
tcrit <- qt(1 - alpha, df = length(sf) - 1)
sumsf <- sum(sf, na.rm = TRUE)
} else {
tstat <- NA; tcrit <- NA; sumsf <- NA
}
# quantile loss
loss_vec <- (alpha - (yv < qv)) * (yv - qv)
total_loss <- sum(loss_vec, na.rm = TRUE)
mean_loss <- mean(loss_vec, na.rm = TRUE)
data.frame(
n = n,
exceptions = ce,
p_exc = pe,
z = z, zcrit = zc,
t_short = tstat, tcrit = tcrit,
sum_sf = sumsf,
total_loss = total_loss,
mean_loss = mean_loss
)
}
# 2) Collect VaR series you have
var_names <- c(
"VarCov_50" = "cmVaR50d_99",
"VarCov_120" = "cmVaR120d_99",
"VarCov_250" = "cmVaR250d_99",
"EWMA_50" = "VaR_ewma50",
"EWMA_120" = "VaR_ewma120",
"EWMA_250" = "VaR_ewma250",
"Hist_50" = "cmVaR50h_99",
"Hist_120" = "cmVaR120h_99",
"Hist_250" = "cmVaR250h_99",
"MC_50" = "cmVaR50m_99",
"MC_120" = "cmVaR120m_99",
"MC_250" = "cmVaR250m_99"
)
var_list <- lapply(var_names, function(obj) if (exists(obj, inherits = TRUE)) get(obj) else NULL)
var_list <- var_list[!sapply(var_list, is.null)] # drop missing
# 3) Run backtests for ALL VaR series
results <- do.call(
rbind,
lapply(names(var_list), function(nm) {
cbind(model = nm, backtest_var(y = as.numeric(cmchange), q = var_list[[nm]], alpha = alpha))
})
)
row.names(results) <- NULL
print(results)
## model n exceptions p_exc z zcrit t_short tcrit
## 1 VarCov_50 755 14 0.018543046 2.3592211 2.326348 3.176010 2.650309
## 2 VarCov_120 685 11 0.016058394 1.5936216 2.326348 2.658118 2.763769
## 3 VarCov_250 555 7 0.012612613 0.6185915 2.326348 2.594303 3.142668
## 4 EWMA_50 755 8 0.010596026 0.1645968 2.326348 2.652299 2.997952
## 5 EWMA_120 685 8 0.011678832 0.4416060 2.326348 2.484191 2.997952
## 6 EWMA_250 555 10 0.018018018 1.8984358 2.326348 2.178316 2.821438
## 7 Hist_50 755 18 0.023841060 3.8223039 2.326348 3.730660 2.566934
## 8 Hist_120 685 13 0.018978102 2.3616320 2.326348 2.330366 2.680998
## 9 Hist_250 555 6 0.010810811 0.1919767 2.326348 2.164940 3.364930
## 10 MC_50 755 12 0.015894040 1.6276797 2.326348 3.305582 2.718079
## 11 MC_120 685 6 0.008759124 -0.3264044 2.326348 1.988431 3.364930
## 12 MC_250 555 4 0.007207207 -0.6612529 2.326348 1.872184 4.540703
## sum_sf total_loss mean_loss
## 1 0.7123264 2.571800 0.003406357
## 2 0.6267987 2.347286 0.003426696
## 3 0.6338765 2.045312 0.003685247
## 4 0.3849883 2.250383 0.002980640
## 5 0.4945944 2.230696 0.003256490
## 6 0.5373925 1.960110 0.003531729
## 7 0.5591988 2.329849 0.003085893
## 8 0.6492843 2.313772 0.003377769
## 9 0.5353620 1.990325 0.003586172
## 10 0.6791893 3.020245 0.004000324
## 11 0.3732375 2.582347 0.003769850
## 12 0.2118250 2.099235 0.003782405
# 4) Export to Excel
#write_xlsx(as.data.frame(results), "VaR_backtests.xlsx")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
library(readxl)
library(dplyr)
library(tidyr)
library(zoo)
library(TTR)
library(PearsonDS)
library(moments)
library(purrr)
## Warning: package 'purrr' was built under R version 4.3.3
#----------------------- load & tidy -----------------------
raw <- read_excel("C:/Users/thech/Downloads/FINM 8100 DATA (1).xlsx")
## New names:
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...11`
## • `` -> `...13`
## • `` -> `...15`
## • `` -> `...17`
## • `` -> `...19`
## • `` -> `...21`
## • `` -> `...23`
## • `` -> `...25`
## • `` -> `...27`
## • `` -> `...29`
## • `` -> `...31`
## • `` -> `...32`
# Make first row the header, like you did, but for the whole sheet
names(raw) <- as.character(unlist(raw[1, ]))
## Warning: The `value` argument of `names<-()` can't be empty as of tibble 3.0.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
df <- raw[-1, ]
# Make first row the header, then drop it
names(raw) <- as.character(unlist(raw[1, ]))
df <- raw[-1, ]
# ---- FIX BAD COLUMN NAMES ----
nm <- names(df)
nm <- trimws(nm)
nm[is.na(nm) | nm == ""] <- paste0("col_", which(is.na(nm) | nm == ""))
nm <- make.names(nm, unique = TRUE) # also fixes duplicates, spaces, etc.
names(df) <- nm
# --------------------------------
# Ensure the yyyymmdd column is called 'year' (your sheet has it in the 2nd col)
if (!"year" %in% names(df)) names(df)[2] <- "year"
# Now this will work
df <- df %>% mutate(year = as.numeric(year))
# Ensure the yyyymmdd column is named 'year' and is numeric
if (!"year" %in% names(df)) {
# If your second column is the yyyymmdd, rename it
names(df)[2] <- "year"
}
df <- df %>% mutate(year = as.numeric(year))
# Grab all bank columns automatically
bank_cols <- grep("^Change\\.", names(df), value = TRUE)
# Numeric returns
df <- df %>%
mutate(across(all_of(bank_cols), ~ suppressWarnings(as.numeric(.x))))
# Focus on dot-com window
dot <- df %>% filter(year >= 20000301 & year <= 20030331)
# For MC skew/kurt we’ll use the pooled industry (all banks stacked)
industry_pool <- dot %>% select(all_of(bank_cols)) %>% as.matrix() %>% as.numeric()
industry_pool <- industry_pool[is.finite(industry_pool)]
sk_ind <- skewness(industry_pool, na.rm = TRUE)
ku_ind <- kurtosis(industry_pool, na.rm = TRUE)
if (ku_ind <= sk_ind^2 + 1) ku_ind <- sk_ind^2 + 1 + 1e-6 # Pearson existence tweak
alpha <- 0.05
z_q <- qnorm(alpha)
wins <- c(50, 120, 250)
lam_from_win <- function(w) (w - 1) / (w + 1) # your choice earlier
names(wins) <- paste0(wins, "d")
dates <- as.Date(as.character(dot$year), format = "%Y%m%d")
#---------------------- helpers ---------------------------
roll_mean <- function(x, w) rollapply(x, w, mean, align = "right", fill = NA, na.rm = TRUE)
roll_sd <- function(x, w) rollapply(x, w, sd, align = "right", fill = NA, na.rm = TRUE)
roll_qtl <- function(x, w, p) rollapply(x, w, function(y) quantile(y, probs = p, na.rm = TRUE),
align = "right", fill = NA)
backtest_var <- function(y, q, alpha = 0.05) {
valid <- is.finite(y) & is.finite(q)
yv <- y[valid]; qv <- q[valid]; n <- length(yv)
if (n == 0) return(tibble(n=0, exceptions=NA, p_exc=NA, z=NA, zcrit=NA,
t_short=NA, tcrit=NA, sum_sf=NA, total_loss=NA, mean_loss=NA))
exc <- (yv < qv)
ce <- sum(exc); pe <- ce/n
z <- (pe - alpha)/sqrt(alpha*(1-alpha)/n); zc <- qnorm(1 - alpha)
sf <- qv[exc] - yv[exc]
if (length(sf) >= 3) {
tstat <- mean(sf) / (sd(sf)/sqrt(length(sf)))
tcrit <- qt(1 - alpha, df = length(sf) - 1)
sumsf <- sum(sf)
} else {
tstat <- NA; tcrit <- NA; sumsf <- NA
}
loss_vec <- (alpha - (yv < qv)) * (yv - qv)
tibble(
n = n, exceptions = ce, p_exc = pe, z = z, zcrit = zc,
t_short = tstat, tcrit = tcrit, sum_sf = sumsf,
total_loss = sum(loss_vec, na.rm = TRUE),
mean_loss = mean(loss_vec, na.rm = TRUE)
)
}
compute_bank <- function(ret_vec) {
r <- as.numeric(ret_vec)
# Var-Cov (equal weights)
vc <- map(wins, ~ roll_mean(r, .x) + z_q * roll_sd(r, .x))
names(vc) <- paste0("VarCov_", names(vc))
# EWMA (lambda from window, like your code)
ewma_list <- imap(wins, function(w, nm) {
lam <- lam_from_win(w)
eps <- r - mean(r, na.rm = TRUE)
sd_ewma <- sqrt(EMA(eps^2, ratio = 1 - lam))
mu_roll <- roll_mean(r, w)
mu_roll + z_q * sd_ewma
})
names(ewma_list) <- paste0("EWMA_", names(ewma_list))
# Historical
hist_list <- map(wins, ~ roll_qtl(r, .x, alpha))
names(hist_list) <- paste0("Hist_", names(hist_list))
# PearsonDS Monte Carlo (match sample length N, then roll quantiles)
N <- length(r)
mu <- mean(r, na.rm = TRUE)
sdv <- sd(r, na.rm = TRUE)
fit <- pearsonFitM(moments = list(mean = mu, variance = sdv^2,
skewness = sk_ind, kurtosis = ku_ind))
set.seed(123)
x_sim <- rpearson(N, params = fit)
mc_list <- map(wins, ~ roll_qtl(x_sim, .x, alpha))
names(mc_list) <- paste0("MC_", names(mc_list))
# Assemble tidy export table
export_df <- tibble(
Date = dates,
Return = r
) %>% bind_cols(as_tibble(vc),
as_tibble(ewma_list),
as_tibble(hist_list),
as_tibble(mc_list))
# Backtests for all series
bt <- bind_rows(lapply(names(export_df)[-(1:2)], function(colnm) {
res <- backtest_var(y = export_df$Return, q = export_df[[colnm]], alpha = alpha)
mutate(res, model = colnm)
}))
list(export = export_df, backtests = bt)
}
#------------------- run for ALL banks --------------------
# long form for clarity, then split by bank and map
long <- dot %>%
select(year, all_of(bank_cols)) %>%
pivot_longer(-year, names_to = "bank", values_to = "ret")
by_bank <- split(long, long$bank)
results_per_bank <- imap(by_bank, function(df_bank, bank_name) {
out <- compute_bank(df_bank$ret)
out$backtests <- mutate(out$backtests, bank = bank_name, .before = 1)
out
})
# 1) One big backtest table across banks & models
backtests_all <- bind_rows(lapply(results_per_bank, `[[`, "backtests"))
# 2) (Optional) write one Excel per bank
# library(writexl)
# iwalk(results_per_bank, function(x, nm) {
# write_xlsx(x$export, paste0(gsub("^Change\\.", "", nm), "_VaR_Export.xlsx"))
# })
# 3) (Optional) plot function for any bank + model set
plot_bank <- function(bank = "Change.cm") {
ex <- results_per_bank[[bank]]$export
matplot(ex$Date, ex[, -c(1,2)], type = "l", lty = 1,
main = paste0(bank, " — 95% VaR (all methods)"),
xlab = "Date", ylab = "Return / VaR")
lines(ex$Date, ex$Return, col = "grey")
}
# Quick peek:
# head(backtests_all)
# plot_bank("Change.cm")
library(writexl)
# Combine all VaR export tables (add bank column)
VaR_all <- dplyr::bind_rows(
lapply(names(results_per_bank), function(bk) {
df <- results_per_bank[[bk]]$export
df$bank <- bk
df
}),
.id = NULL
)
# Backtests already combined earlier
Backtests_all <- backtests_all
# Save as one Excel with 2 sheets
write_xlsx(
list(
"VaR_Exports" = VaR_all,
"Backtests" = Backtests_all
),
path = "AllBanks_VaR_and_Backtestsdot95.xlsx"
)
getwd()
## [1] "C:/Users/thech/OneDrive/Documents/ANU/SEM 4/applied project/VaR"