Market Risk Analytics: Multi-Method VaR & Backtesting Framework

Portfolio: SPY, TLT, GLD, HYG

# 1. LOAD LIBRARIES
library(quantmod)   # Financial data
library(MASS)       # Monte Carlo simulation
library(ggplot2)    # Visualisation

# 2. DATA COLLECTION
# Multi-asset portfolio
assets = c("SPY", "TLT", "GLD", "HYG")

# Pull historical data
getSymbols(assets, from = "2020-01-01")

## [1] "SPY" "TLT" "GLD" "HYG"

# Extract closing prices and merge
prices = na.omit(merge(
  Cl(SPY), Cl(TLT), Cl(GLD), Cl(HYG)
))

# 3. RETURN CALCULATION
# Compute log returns (standard in finance)
returns = na.omit(diff(log(prices)))

# 4. PORTFOLIO CONSTRUCTION
# Equal weights
weights = c(0.25, 0.25, 0.25, 0.25)

# Portfolio returns (matrix multiplication)
portfolio_returns = returns %*% weights

# 5. PARAMETRIC VaR FUNCTION
parametric_var = function(returns, weights, alpha = 0.99, V = 1e6) {
  
  mu = colMeans(returns)              # Mean returns per asset
  cov_mat = cov(returns)              # Covariance matrix
  
  port_mean = sum(weights * mu)       # Portfolio mean
  port_sd = sqrt(t(weights) %*% cov_mat %*% weights)  # Portfolio standard deviation
  
  z = qnorm(1 - alpha)                # Normal quantile
  
  VaR = V * (port_mean + z * port_sd) # Convert to monetary VaR
  
  return(abs(VaR))
}

# 6. HISTORICAL VaR FUNCTION
historical_var = function(portfolio_returns, alpha = 0.99, V = 1e6) {
  
  VaR = quantile(portfolio_returns, probs = 1 - alpha)
  
  return(abs(V * VaR))
}

# 7. MONTE CARLO VaR FUNCTION
monte_carlo_var = function(returns, weights, alpha = 0.99, 
                            V = 1e6, n_sim = 10000) {
  
  mu = colMeans(returns)
  cov_mat = cov(returns)
  
  # Simulate correlated returns
  sim_returns = mvrnorm(n_sim, mu = mu, Sigma = cov_mat)
  
  # Portfolio simulation
  port_sim = sim_returns %*% weights
  
  VaR = quantile(port_sim, probs = 1 - alpha)
  
  return(abs(V * VaR))
}

# 8. COMPUTE VaR VALUES
V = 1e6   # Portfolio value in dollars

var_param = as.numeric(parametric_var(returns, weights, 0.99, V))
var_hist  = historical_var(portfolio_returns, 0.99, V)
var_mc    = monte_carlo_var(returns, weights, 0.99, V)

cat("Parametric VaR:", var_param, "\n")

## Parametric VaR: 14812.46

cat("Historical VaR:", var_hist, "\n")

## Historical VaR: 16184.21

cat("Monte Carlo VaR:", var_mc, "\n")

## Monte Carlo VaR: 14836

# 9. BACKTESTING FUNCTION
backtest_var = function(returns, VaR, V = 1e6) {
  
  returns = as.numeric(returns)
  VaR_return = VaR / V
  
  violations = returns < -VaR_return
  x = sum(violations)
  n = length(returns)
  violation_rate = x / n
  
  return(list(
    violations = x,
    total_obs = n,
    violation_rate = violation_rate
  ))
}

# 10. KUPIEC TEST FUNCTION
kupiec_test = function(returns, VaR, alpha = 0.99, V = 1e6) {
  
  returns = as.numeric(returns)
  VaR_return = VaR / V
  
  violations = returns < -VaR_return
  x = sum(violations)
  n = length(returns)
  
  p_hat = x / n
  p = 1 - alpha
  
  LR = -2 * log(((1 - p)^(n - x) * p^x) /
                   ((1 - p_hat)^(n - x) * p_hat^x))
  
  p_value = 1 - pchisq(LR, df = 1)
  
  return(list(
    violations = x,
    violation_rate = p_hat,
    expected_rate = p,
    LR_stat = LR,
    p_value = p_value
  ))
}

# 11. RUN BACKTESTING
bt_param = kupiec_test(portfolio_returns, var_param, 0.99, V)
bt_hist  = kupiec_test(portfolio_returns, var_hist, 0.99, V)
bt_mc    = kupiec_test(portfolio_returns, var_mc, 0.99, V)

print(bt_param)

## $violations
## [1] 23
## 
## $violation_rate
## [1] 0.01471529
## 
## $expected_rate
## [1] 0.01
## 
## $LR_stat
## [1] 3.065054
## 
## $p_value
## [1] 0.07999251

print(bt_hist)

## $violations
## [1] 16
## 
## $violation_rate
## [1] 0.01023672
## 
## $expected_rate
## [1] 0.01
## 
## $LR_stat
## [1] 0.00877897
## 
## $p_value
## [1] 0.9253505

print(bt_mc)

## $violations
## [1] 23
## 
## $violation_rate
## [1] 0.01471529
## 
## $expected_rate
## [1] 0.01
## 
## $LR_stat
## [1] 3.065054
## 
## $p_value
## [1] 0.07999251

# 12. EXPECTED SHORTFALL (CVaR)
expected_shortfall = function(returns, alpha = 0.99) {
  
  returns = as.numeric(returns)
  VaR = quantile(returns, probs = 1 - alpha)
  ES = mean(returns[returns <= VaR])
  
  return(abs(ES))
}

es_value = expected_shortfall(portfolio_returns)
cat("Expected Shortfall (ES):", es_value, "\n")

## Expected Shortfall (ES): 0.02550301

# 13. VISUALISATION

df = data.frame(returns = portfolio_returns)

ggplot(df, aes(x = returns)) +
  geom_histogram(bins = 50, fill="steelblue", color="black") +
  ggtitle("Portfolio Return Distribution") +
  xlab("Daily Log Returns") +
  ylab("Frequency")

# Save plot
ggsave("var_distribution_new_data.png")

## Saving 7 x 5 in image