Portfolio Risk Management: A Practical Guide to VaR and Monte Carlo Simulation with Python

• 1. Introduction: Speaking Risk in Numbers
• 2. VaR Fundamentals and Three Methodologies
  • What Is VaR?
  • Three Key VaR Parameters
  • Overview of Three Calculation Methodologies
• 3. Parametric VaR Implementation
  • Multi-Asset Extension of Parametric VaR: Using the Covariance Matrix
• 4. Historical Simulation VaR
• 5. Monte Carlo Simulation VaR
  • Geometric Brownian Motion (GBM) Based Simulation
  • Determining Simulation Count and Convergence
• 6. CVaR (Expected Shortfall) Extension
• 7. Portfolio Optimization and Risk Integration
• 8. VaR Methodology Comparison Table
  • Methodology Selection Guidelines
• 9. Backtesting and Model Validation
  • Kupiec Test (POF Test)
  • Christoffersen Test (Conditional Coverage Test)
• 10. Operational Considerations
  • Data Quality Issues
  • Model Limitations and Assumption Violations
  • Regulatory Considerations
  • Production Operations Checkpoint
• 11. Failure Cases and Recovery Procedures
  • Case 1: The 2008 Financial Crisis and VaR's Limitations
  • Case 2: Risk Underestimation Due to Implementation Errors
  • Case 3: Failure to Detect Volatility Regime Transition
• 12. Checklist
  • Model Construction Phase
  • Backtesting Phase
  • Production Operations Phase
  • Governance Phase
• 13. References

1. Introduction: Speaking Risk in Numbers

"How much risk does this portfolio carry?" Answering "it's somewhat risky" or "it has high volatility" is meaningless in practice. Risk managers, regulators, and investors all want concrete numbers. "At a 95% confidence level, the maximum one-day loss is 230 million won." This is why Value at Risk (VaR) exists.

Since VaR was popularized in the early 1990s through JP Morgan's RiskMetrics system, it has become a core metric in the Basel regulatory framework (Basel II/III/IV). However, the 2008 financial crisis starkly exposed VaR's limitations, making it essential to perform integrated risk frameworks that combine CVaR (Expected Shortfall), stress testing, and backtesting rather than relying on a single VaR number.

This article implements each of the three major VaR calculation methodologies -- Parametric (Variance-Covariance), Historical Simulation, and Monte Carlo Simulation -- in Python code, validates models through backtesting, and covers in detail the precautions and failure cases you must know when operating risk systems in production environments.

The target audience is developers or junior quants familiar with Python basics and pandas/numpy. Financial domain knowledge is explained from scratch wherever possible, so you can follow along even without background knowledge.

2. VaR Fundamentals and Three Methodologies

What Is VaR?

VaR (Value at Risk) represents the maximum expected loss a portfolio can experience under a given confidence level and holding period. Mathematically expressed:

VaR_alpha = -inf{x : P(L > x) <= 1 - alpha}

Here, alpha is the confidence level (typically 95% or 99%), L is the loss distribution, and x is the loss amount. Intuitively, "95% VaR = 100 million won" means "there is a 95% probability that losses will not exceed 100 million won over the next day." Conversely, it also means there is a 5% probability of losses exceeding 100 million won.

Three Key VaR Parameters

• Confidence Level: 95%, 99%, 99.5%, etc. Basel regulations use 99% as the default.
• Holding Period: 1 day, 10 days, 1 month, etc. Trading books typically use 1 day, while investment portfolios use 10 days to 1 month.
• Observation Window: The historical data period used for VaR calculation. Typically 250 days (1 year) to 500 days (2 years).

Overview of Three Calculation Methodologies

VaR calculation methods are broadly classified into three categories.

Parametric (Variance-Covariance) Method: Assumes returns follow a normal distribution and calculates VaR using only the mean and standard deviation (or covariance matrix). Fast but fails to capture fat-tail phenomena.

Historical Simulation Method: Directly uses past return data to construct the loss distribution. No distribution assumption is needed, but it cannot reflect scenarios that have not occurred in the past.

Monte Carlo Simulation Method: Generates tens of thousands to hundreds of thousands of random scenarios from a stochastic model (e.g., Geometric Brownian Motion) to estimate the loss distribution. Most flexible but computationally expensive.

3. Parametric VaR Implementation

Parametric VaR is the simplest and fastest method. Assuming portfolio returns follow a normal distribution N(mu, sigma^2), VaR is calculated as follows:

VaR = -(mu + z*alpha * sigma) _ Portfolio_Value

Here, z_alpha is the quantile of the standard normal distribution (e.g., z = -1.645 for 95% confidence level).

import numpy as np
import pandas as pd
from scipy import stats
import yfinance as yf
from datetime import datetime, timedelta

# Portfolio tickers and weights
tickers = ['AAPL', 'MSFT', 'GOOGL', 'AMZN', 'NVDA']
weights = np.array([0.25, 0.20, 0.20, 0.20, 0.15])
portfolio_value = 1_000_000_000  # 1 billion won

# Collect 2 years of price data
end_date = datetime(2026, 3, 1)
start_date = end_date - timedelta(days=730)
prices = yf.download(tickers, start=start_date, end=end_date)['Adj Close']

# Calculate daily log returns
log_returns = np.log(prices / prices.shift(1)).dropna()

# Portfolio returns
portfolio_returns = log_returns.dot(weights)

# Parametric VaR calculation
confidence_levels = [0.90, 0.95, 0.99]
mean_return = portfolio_returns.mean()
std_return = portfolio_returns.std()

print("=" * 60)
print("Parametric VaR (Variance-Covariance Method)")
print("=" * 60)
print(f"Daily mean return:    {mean_return:.6f}")
print(f"Daily std deviation:  {std_return:.6f}")
print(f"Portfolio value:      {portfolio_value:,.0f} won")
print("-" * 60)

for cl in confidence_levels:
    z_score = stats.norm.ppf(1 - cl)
    var_pct = -(mean_return + z_score * std_return)
    var_amount = var_pct * portfolio_value
    print(f"{cl*100:.0f}% VaR: {var_pct:.4%} = {var_amount:,.0f} won")

# 10-day VaR (Basel standard): sqrt(10) scaling
print("\n10-day VaR (sqrt(T) scaling):")
for cl in confidence_levels:
    z_score = stats.norm.ppf(1 - cl)
    var_1d = -(mean_return + z_score * std_return)
    var_10d = var_1d * np.sqrt(10)
    var_10d_amount = var_10d * portfolio_value
    print(f"{cl*100:.0f}% 10-day VaR: {var_10d:.4%} = {var_10d_amount:,.0f} won")

Let us highlight a few key points from the code above. First, we use log returns. Log returns are additive over time, allowing multi-period return aggregation, and better fit the normal distribution assumption. Second, the sqrt(10) rule is applied for scaling to 10-day VaR. This is only accurate under the assumption that returns are independent and identically distributed (i.i.d.). Since real financial time series exhibit autocorrelation and volatility clustering, this scaling is merely an approximation.

Multi-Asset Extension of Parametric VaR: Using the Covariance Matrix

Instead of directly computing the standard deviation of the single portfolio return, using the covariance matrix of individual assets enables risk contribution analysis.

# Covariance matrix-based portfolio VaR
cov_matrix = log_returns.cov()

# Portfolio variance = w^T * Sigma * w
portfolio_variance = weights @ cov_matrix.values @ weights
portfolio_std = np.sqrt(portfolio_variance)

# Marginal VaR contribution per asset (Component VaR)
marginal_var = cov_matrix.values @ weights / portfolio_std
component_var = weights * marginal_var

print("\n" + "=" * 60)
print("Component VaR Analysis (99% Confidence Level)")
print("=" * 60)
z_99 = stats.norm.ppf(0.99)

for i, ticker in enumerate(tickers):
    cvar_pct = component_var[i] * z_99
    cvar_amount = cvar_pct * portfolio_value
    contribution_pct = component_var[i] / portfolio_std * 100
    print(f"{ticker}: Weight {weights[i]*100:.0f}% | "
          f"Risk contribution {contribution_pct:.1f}% | "
          f"Component VaR {cvar_amount:,.0f} won")

# Verify diversification benefit
individual_vars = np.array([
    weights[i] * np.sqrt(cov_matrix.values[i, i]) * z_99
    for i in range(len(tickers))
])
undiversified_var = individual_vars.sum() * portfolio_value
diversified_var = portfolio_std * z_99 * portfolio_value
diversification_benefit = undiversified_var - diversified_var

print(f"\nUndiversified VaR (sum):    {undiversified_var:,.0f} won")
print(f"Diversified VaR (portfolio): {diversified_var:,.0f} won")
print(f"Diversification benefit:     {diversification_benefit:,.0f} won "
      f"({diversification_benefit/undiversified_var*100:.1f}% reduction)")

The key insight from this analysis is the diversification benefit. The portfolio's overall VaR is smaller than the simple sum of individual stock VaRs, and this difference represents the risk reduction effect of diversification. The lower the correlations, the greater the diversification benefit, and as correlations approach 1 (e.g., correlation convergence during crises), diversification benefit shrinks rapidly.

4. Historical Simulation VaR

Historical simulation is more realistic than the parametric method in that it makes no distribution assumptions. Since past return data is directly used as the loss distribution, fat tails and skewness are naturally reflected.

def historical_var(returns, confidence_level=0.95, portfolio_value=1_000_000_000):
    """
    Calculate VaR using Historical Simulation method.

    Parameters
    ----------
    returns : pd.Series
        Portfolio daily returns time series
    confidence_level : float
        Confidence level (e.g., 0.95, 0.99)
    portfolio_value : float
        Total portfolio value (won)

    Returns
    -------
    dict
        VaR calculation results
    """
    sorted_returns = returns.sort_values()
    n = len(sorted_returns)

    # Lower (1-confidence_level) percentile
    percentile_index = int(np.floor((1 - confidence_level) * n))
    var_pct = -sorted_returns.iloc[percentile_index]

    # Interpolated version (more accurate)
    var_pct_interp = -np.percentile(returns, (1 - confidence_level) * 100)

    return {
        'var_pct': var_pct,
        'var_pct_interp': var_pct_interp,
        'var_amount': var_pct * portfolio_value,
        'var_amount_interp': var_pct_interp * portfolio_value,
        'n_observations': n,
        'worst_loss': -sorted_returns.iloc[0],
        'worst_loss_date': sorted_returns.index[0],
    }

# Calculate Historical VaR
results = {}
for cl in [0.90, 0.95, 0.99]:
    result = historical_var(portfolio_returns, cl)
    results[cl] = result
    print(f"\n{cl*100:.0f}% Historical VaR:")
    print(f"  Percentile method: {result['var_pct']:.4%} = {result['var_amount']:,.0f} won")
    print(f"  Interpolation:     {result['var_pct_interp']:.4%} = {result['var_amount_interp']:,.0f} won")

print(f"\nMaximum loss in observation period: {results[0.95]['worst_loss']:.4%} "
      f"({results[0.95]['worst_loss_date'].strftime('%Y-%m-%d')})")

# Rolling analysis of Historical VaR
rolling_window = 250  # 1-year window
rolling_var_95 = portfolio_returns.rolling(rolling_window).apply(
    lambda x: -np.percentile(x, 5), raw=True
)

print(f"\nRolling VaR (250-day window) statistics:")
print(f"  Mean 95% VaR: {rolling_var_95.mean():.4%}")
print(f"  Max 95% VaR:  {rolling_var_95.max():.4%} (highest risk period)")
print(f"  Min 95% VaR:  {rolling_var_95.min():.4%} (lowest risk period)")

The biggest weakness of historical simulation is the "ghost effect." When an extreme event is included in the observation window, VaR stays elevated, then drops abruptly on the day that event falls out of the window. For example, when using a 250-day window, if there was a large crash exactly 250 days ago, VaR suddenly decreases the next day when that crash exits the window. This is an artifact where the measured value changes abruptly even though actual risk has not changed.

A solution to this problem is the exponentially weighted moving average (EWMA) based weighted historical simulation, which mitigates the ghost effect by assigning higher weights to more recent observations.

5. Monte Carlo Simulation VaR

Monte Carlo simulation is the most flexible and powerful method for VaR calculation. The basic idea is to randomly generate thousands to tens of thousands of price paths from a stochastic model (typically Geometric Brownian Motion, GBM), then calculate portfolio P&L from each path to construct the loss distribution.

Geometric Brownian Motion (GBM) Based Simulation

In Geometric Brownian Motion, the change in stock price S follows this stochastic differential equation:

dS = mu _ S _ dt + sigma _ S _ dW

Here, mu is the drift (expected return), sigma is volatility, and dW is the increment of the Wiener process. Discretized:

S(t+dt) = S(t) _ exp((mu - sigma^2/2) _ dt + sigma _ sqrt(dt) _ Z)

Where Z is a random number following the standard normal distribution.

def monte_carlo_var(
    returns_df,
    weights,
    portfolio_value,
    n_simulations=50_000,
    n_days=1,
    confidence_level=0.95,
    seed=42
):
    """
    Calculate portfolio VaR using Monte Carlo simulation.

    Generates scenarios from a multivariate normal distribution
    with correlation structure based on GBM.

    Parameters
    ----------
    returns_df : pd.DataFrame
        Individual asset daily log returns DataFrame
    weights : np.ndarray
        Asset weight vector
    portfolio_value : float
        Total portfolio value (won)
    n_simulations : int
        Number of simulations
    n_days : int
        Holding period (days)
    confidence_level : float
        Confidence level
    seed : int
        Random seed (for reproducibility)

    Returns
    -------
    dict
        VaR, CVaR, simulation results, etc.
    """
    np.random.seed(seed)

    # Individual asset mean returns and covariance matrix
    mean_returns = returns_df.mean().values
    cov_matrix = returns_df.cov().values
    n_assets = len(weights)

    # Cholesky decomposition for correlated random number generation
    L = np.linalg.cholesky(cov_matrix)

    # Run simulation
    portfolio_sim_returns = np.zeros(n_simulations)

    for sim in range(n_simulations):
        cumulative_return = 0
        for day in range(n_days):
            # Generate independent standard normal random numbers
            Z = np.random.standard_normal(n_assets)
            # Apply correlation
            correlated_returns = mean_returns + L @ Z
            # Portfolio return
            daily_port_return = weights @ correlated_returns
            cumulative_return += daily_port_return

        portfolio_sim_returns[sim] = cumulative_return

    # P&L distribution
    sim_pnl = portfolio_sim_returns * portfolio_value

    # VaR calculation
    var_pct = -np.percentile(portfolio_sim_returns, (1 - confidence_level) * 100)
    var_amount = var_pct * portfolio_value

    # CVaR (Expected Shortfall) calculation
    var_threshold = np.percentile(portfolio_sim_returns, (1 - confidence_level) * 100)
    tail_losses = portfolio_sim_returns[portfolio_sim_returns <= var_threshold]
    cvar_pct = -tail_losses.mean()
    cvar_amount = cvar_pct * portfolio_value

    return {
        'var_pct': var_pct,
        'var_amount': var_amount,
        'cvar_pct': cvar_pct,
        'cvar_amount': cvar_amount,
        'sim_returns': portfolio_sim_returns,
        'sim_pnl': sim_pnl,
        'n_simulations': n_simulations,
        'n_days': n_days,
        'confidence_level': confidence_level,
    }

# Run Monte Carlo VaR
mc_results = {}
for cl in [0.90, 0.95, 0.99]:
    result = monte_carlo_var(
        log_returns, weights, portfolio_value,
        n_simulations=100_000,
        n_days=1,
        confidence_level=cl,
        seed=42
    )
    mc_results[cl] = result
    print(f"\n{cl*100:.0f}% Monte Carlo VaR (100,000 simulations):")
    print(f"  VaR:  {result['var_pct']:.4%} = {result['var_amount']:,.0f} won")
    print(f"  CVaR: {result['cvar_pct']:.4%} = {result['cvar_amount']:,.0f} won")

# Check VaR convergence by number of simulations
print("\n95% VaR convergence by number of simulations:")
for n_sim in [1_000, 5_000, 10_000, 50_000, 100_000, 500_000]:
    result = monte_carlo_var(
        log_returns, weights, portfolio_value,
        n_simulations=n_sim, confidence_level=0.95, seed=42
    )
    print(f"  N={n_sim:>7,}: VaR = {result['var_pct']:.6%}")

Determining Simulation Count and Convergence

The number of simulations (N) in Monte Carlo VaR directly affects accuracy. The standard error of the VaR estimate is roughly proportional to 1/sqrt(N), so quadrupling the number of simulations halves the standard error. In practice, a minimum of 10,000 is recommended, and 50,000 to 100,000 for production. When using extreme confidence levels like 99% VaR, more simulations are needed since there are fewer samples in the tail region.

Cholesky decomposition decomposes the covariance matrix into a lower triangular matrix L such that L * L^T = Sigma. Multiplying an independent standard normal random vector Z by L yields a correlated random vector with the original covariance structure. If the covariance matrix is not positive definite, Cholesky decomposition fails, which occurs when the data period is shorter than the number of assets or when there are linear dependencies. In such cases, eigendecomposition should be used, or the Higham algorithm should be applied to approximate the nearest positive definite matrix.

6. CVaR (Expected Shortfall) Extension

The biggest limitation of VaR is that it provides no information about "how large losses beyond VaR could be." Even if 95% VaR is 2%, losses in the remaining 5% could be 3% or 30%. CVaR (Conditional VaR), also known as Expected Shortfall (ES), is defined as the average of losses exceeding VaR and captures tail risk much more effectively.

The mathematical definition is:

CVaR_alpha = E[L | L > VaR_alpha]

Unlike VaR, CVaR is a coherent risk measure that satisfies subadditivity. Subadditivity means that the risk of combining two portfolios is less than or equal to the sum of their individual risks. VaR does not satisfy this property, so optimizing based on VaR may underestimate the benefits of diversification.

Since Basel III, the regulatory framework has adopted Expected Shortfall as the primary market risk measure instead of VaR (FRTB: Fundamental Review of the Trading Book).

def comprehensive_risk_metrics(returns, portfolio_value, confidence_level=0.95):
    """
    Calculate comprehensive risk metrics for a portfolio.

    Computes VaR, CVaR, plus maximum drawdown (MDD), volatility,
    skewness, kurtosis, etc. to provide a complete risk picture.
    """
    sorted_returns = returns.sort_values()
    n = len(sorted_returns)

    # VaR
    var_pct = -np.percentile(returns, (1 - confidence_level) * 100)

    # CVaR (Expected Shortfall)
    var_threshold = np.percentile(returns, (1 - confidence_level) * 100)
    tail_returns = returns[returns <= var_threshold]
    cvar_pct = -tail_returns.mean()

    # Volatility (annualized)
    daily_vol = returns.std()
    annual_vol = daily_vol * np.sqrt(252)

    # Skewness and Kurtosis
    skewness = returns.skew()
    kurtosis = returns.kurtosis()  # Excess kurtosis (normal distribution = 0)

    # Maximum Drawdown
    cumulative = (1 + returns).cumprod()
    running_max = cumulative.cummax()
    drawdown = (cumulative - running_max) / running_max
    max_drawdown = drawdown.min()

    # VaR ratio vs normal distribution (tail risk indicator)
    parametric_var = -(returns.mean() + stats.norm.ppf(1 - confidence_level) * returns.std())
    tail_ratio = var_pct / parametric_var if parametric_var > 0 else float('inf')

    results = {
        'VaR': var_pct,
        'CVaR': cvar_pct,
        'VaR (amount)': var_pct * portfolio_value,
        'CVaR (amount)': cvar_pct * portfolio_value,
        'CVaR/VaR ratio': cvar_pct / var_pct if var_pct > 0 else 0,
        'Daily volatility': daily_vol,
        'Annual volatility': annual_vol,
        'Skewness': skewness,
        'Excess kurtosis': kurtosis,
        'Maximum drawdown': max_drawdown,
        'Tail ratio (actual/parametric)': tail_ratio,
    }

    print("=" * 60)
    print(f"Comprehensive Risk Analysis ({confidence_level*100:.0f}% Confidence Level)")
    print("=" * 60)
    for key, value in results.items():
        if 'amount' in key:
            print(f"  {key}: {value:,.0f} won")
        elif 'ratio' in key or 'volatility' in key or 'drawdown' in key:
            print(f"  {key}: {value:.4%}")
        else:
            print(f"  {key}: {value:.6f}")

    # Interpretation guide
    print("\n[Interpretation Guide]")
    if kurtosis > 1:
        print(f"  - Excess kurtosis {kurtosis:.2f}: Tails are heavier than normal "
              f"distribution (fat tails). Parametric VaR likely underestimates risk.")
    if skewness < -0.5:
        print(f"  - Skewness {skewness:.2f}: Negative asymmetry. "
              f"Large losses tend to occur more frequently/severely than large gains.")
    if tail_ratio > 1.2:
        print(f"  - Tail ratio {tail_ratio:.2f}: Actual tail risk is "
              f"{(tail_ratio-1)*100:.0f}% greater than the normal distribution assumption.")

    return results

# Execute
risk_metrics = comprehensive_risk_metrics(portfolio_returns, portfolio_value, 0.95)

The CVaR/VaR ratio is a useful indicator of tail risk severity. Under a normal distribution at a 95% confidence level, the CVaR/VaR ratio is approximately 1.25. If this ratio exceeds 1.5, it signals the presence of extreme tail risk, and you should strengthen stress testing and consider portfolio adjustments.

7. Portfolio Optimization and Risk Integration

Moving beyond simply measuring VaR and CVaR, incorporating them into portfolio construction enables risk-based optimization. The PyPortfolioOpt library provides various methods including mean-variance optimization (Markowitz), minimum CVaR optimization, and Hierarchical Risk Parity (HRP).

from pypfopt import EfficientFrontier, risk_models, expected_returns
from pypfopt import HRPOpt
from pypfopt.efficient_frontier import EfficientCVaR

# Estimate expected returns and covariance matrix
mu = expected_returns.mean_historical_return(prices)
S = risk_models.sample_cov(prices)

# 1. Mean-variance optimization (minimum volatility portfolio)
ef_minvol = EfficientFrontier(mu, S)
ef_minvol.min_volatility()
w_minvol = ef_minvol.clean_weights()
perf_minvol = ef_minvol.portfolio_performance(verbose=False)

# 2. Maximum Sharpe ratio portfolio
ef_sharpe = EfficientFrontier(mu, S)
ef_sharpe.max_sharpe(risk_free_rate=0.035)
w_sharpe = ef_sharpe.clean_weights()
perf_sharpe = ef_sharpe.portfolio_performance(verbose=False)

# 3. Hierarchical Risk Parity (HRP)
hrp = HRPOpt(log_returns)
hrp.optimize()
w_hrp = hrp.clean_weights()
perf_hrp = hrp.portfolio_performance(verbose=False)

# 4. Minimum CVaR portfolio
ef_cvar = EfficientCVaR(mu, log_returns)
ef_cvar.min_cvar()
w_cvar = ef_cvar.clean_weights()

print("=" * 70)
print("Portfolio Optimization Results Comparison")
print("=" * 70)
print(f"{'Strategy':<20} {'Exp Return':>10} {'Volatility':>10} {'Sharpe':>8}")
print("-" * 70)
print(f"{'Min Volatility':<20} {perf_minvol[0]:>10.2%} {perf_minvol[1]:>10.2%} "
      f"{perf_minvol[2]:>8.3f}")
print(f"{'Max Sharpe':<20} {perf_sharpe[0]:>10.2%} {perf_sharpe[1]:>10.2%} "
      f"{perf_sharpe[2]:>8.3f}")
print(f"{'HRP':<20} {perf_hrp[0]:>10.2%} {perf_hrp[1]:>10.2%} "
      f"{perf_hrp[2]:>8.3f}")

print("\nWeight comparison:")
for ticker in tickers:
    print(f"  {ticker}: MinVol {w_minvol.get(ticker,0):.1%} | "
          f"MaxSharpe {w_sharpe.get(ticker,0):.1%} | "
          f"HRP {w_hrp.get(ticker,0):.1%} | "
          f"MinCVaR {w_cvar.get(ticker,0):.1%}")

The minimum volatility portfolio minimizes overall risk but may also have lower expected returns. The maximum Sharpe portfolio has the best risk-adjusted return but tends to concentrate in specific stocks. HRP (Hierarchical Risk Parity) does not require inversion of the covariance matrix, making it less sensitive to estimation errors and providing robust results when the number of assets is large or the covariance matrix is unstable.

In practice, rather than relying on a single optimization method, it is common to ensemble results from multiple methods or combine investor views using the Black-Litterman model.

8. VaR Methodology Comparison Table

A systematic comparison of the three VaR methodologies is as follows:

Methodology Selection Guidelines

• Simple stock/bond portfolio: Start with Historical VaR and compare Parametric VaR as a benchmark
• Including options/derivatives: Monte Carlo VaR is mandatory (must reflect non-linear payoffs)
• Real-time risk monitoring: Parametric VaR (speed priority)
• Regulatory reporting purposes: Historical or Monte Carlo VaR + backtesting is mandatory
• Risk limit setting: CVaR-based is recommended (VaR alone does not reflect tail risk)

9. Backtesting and Model Validation

Once you have built a VaR model, you must validate model accuracy through backtesting. Backtesting is the procedure of calculating VaR on historical data, then statistically testing whether the number of times actual losses exceed VaR (VaR breaches) matches theoretical expectations.

For example, backtesting 95% VaR over 250 days yields an expected number of breaches of approximately 12.5 (250 x 5%). If the number of breaches is significantly more than this, the model is underestimating risk; if fewer, it is overestimating.

Kupiec Test (POF Test)

The Kupiec test uses a likelihood ratio test to check whether the proportion of failures matches the theoretical ratio.

Christoffersen Test (Conditional Coverage Test)

The Christoffersen test also tests independence of breaches. If VaR breaches occur consecutively (clustering), it indicates the model is not properly reflecting volatility clustering.

def backtest_var(returns, var_series, confidence_level=0.95, method_name=""):
    """
    VaR Backtesting: Kupiec Test

    Parameters
    ----------
    returns : pd.Series
        Actual portfolio returns
    var_series : pd.Series
        VaR estimates for the same period (positive values)
    confidence_level : float
        VaR confidence level
    method_name : str
        Method name (for output)
    """
    # Determine VaR breaches
    breaches = returns < -var_series
    n_breaches = breaches.sum()
    n_total = len(returns)
    breach_rate = n_breaches / n_total
    expected_rate = 1 - confidence_level
    expected_breaches = expected_rate * n_total

    print(f"\n{'=' * 60}")
    print(f"VaR Backtesting Results: {method_name}")
    print(f"{'=' * 60}")
    print(f"  Observation period:    {n_total} trading days")
    print(f"  VaR breaches:          {n_breaches}")
    print(f"  Expected breaches:     {expected_breaches:.1f}")
    print(f"  Actual breach rate:    {breach_rate:.2%}")
    print(f"  Expected breach rate:  {expected_rate:.2%}")

    # Kupiec likelihood ratio test
    p = expected_rate
    p_hat = breach_rate if breach_rate > 0 else 1e-10

    # LR statistic calculation
    if n_breaches == 0:
        lr_stat = -2 * (n_total * np.log(1 - p) - n_total * np.log(1 - p_hat))
    elif n_breaches == n_total:
        lr_stat = -2 * (n_total * np.log(p) - n_total * np.log(p_hat))
    else:
        lr_stat = -2 * (
            n_breaches * np.log(p) +
            (n_total - n_breaches) * np.log(1 - p) -
            n_breaches * np.log(p_hat) -
            (n_total - n_breaches) * np.log(1 - p_hat)
        )

    # Compare with chi-squared distribution (1 degree of freedom)
    p_value = 1 - stats.chi2.cdf(lr_stat, df=1)

    print(f"\n  Kupiec LR statistic: {lr_stat:.4f}")
    print(f"  p-value:             {p_value:.4f}")

    if p_value > 0.05:
        verdict = "PASS (model adequate)"
    else:
        verdict = "FAIL (model inadequate - review needed)"
    print(f"  Verdict (5% significance): {verdict}")

    # Basel Traffic Light System
    # Based on 250 days: Green (0-4), Yellow (5-9), Red (10+)
    if n_total >= 250:
        scale_factor = n_total / 250
        scaled_breaches = n_breaches / scale_factor
        if scaled_breaches <= 4:
            zone = "Green Zone - Good"
        elif scaled_breaches <= 9:
            zone = "Yellow Zone - Caution needed"
        else:
            zone = "Red Zone - Model revision required"
        print(f"  Basel Traffic Light: {zone}")

    # Consecutive breach analysis (Christoffersen independence preliminary analysis)
    if n_breaches > 0:
        consecutive = 0
        max_consecutive = 0
        for b in breaches:
            if b:
                consecutive += 1
                max_consecutive = max(max_consecutive, consecutive)
            else:
                consecutive = 0
        print(f"\n  Max consecutive breaches: {max_consecutive} trading days")
        if max_consecutive >= 3:
            print(f"  Warning: {max_consecutive} consecutive breaches. "
                  f"Possible failure to capture volatility clustering.")

    return {
        'n_breaches': n_breaches,
        'breach_rate': breach_rate,
        'lr_stat': lr_stat,
        'p_value': p_value,
        'pass': p_value > 0.05,
    }

# Perform backtesting with rolling VaR
window = 250
test_returns = portfolio_returns.iloc[window:]

# Rolling parametric VaR calculation
rolling_mean = portfolio_returns.rolling(window).mean().iloc[window:]
rolling_std = portfolio_returns.rolling(window).std().iloc[window:]
z_95 = stats.norm.ppf(0.95)
parametric_var_series = -(rolling_mean - z_95 * rolling_std)

# Rolling historical VaR calculation
historical_var_series = portfolio_returns.rolling(window).apply(
    lambda x: -np.percentile(x, 5), raw=True
).iloc[window:]

# Run backtesting
bt_param = backtest_var(test_returns, parametric_var_series, 0.95, "Parametric VaR")
bt_hist = backtest_var(test_returns, historical_var_series, 0.95, "Historical VaR")

When backtesting reveals model failure, you must immediately analyze the cause. Common causes and responses are as follows:

• Too many breaches (risk underestimation): Extend the observation window or expand data to include stress periods. Apply GARCH models to reflect time-varying volatility. Or operate conservatively by raising the confidence level.
• Too few breaches (risk overestimation): Capital efficiency decreases, so refine model parameters or improve volatility estimation methods.
• Breaches occurring in clusters: Integrate volatility clustering models such as GARCH/EGARCH into the VaR framework.

10. Operational Considerations

Data Quality Issues

VaR model accuracy is entirely dependent on input data quality. Common data issues encountered in practice include:

• Missing price data: If prices are missing on certain trading days, return calculations become distorted. Simple forward fill understates volatility. If the missing ratio exceeds 5%, the VaR reliability for that stock should be reassessed.
• Unadjusted dividends/splits: Not using adjusted prices causes artificial crashes on ex-dividend or split dates, overestimating VaR.
• Survivorship bias: If delisted stocks are excluded from the data, risk is underestimated. This is particularly serious in historical VaR.
• Stale prices: For illiquid assets, the last traded price may not reflect current value, distorting both volatility and correlations.

Model Limitations and Assumption Violations

Regulatory Considerations

• Basel III/IV: Market risk capital calculated based on 97.5% ES (Expected Shortfall). Independent backtesting must pass to use internal models.
• FRTB (Fundamental Review of the Trading Book): Transition from VaR-based to ES-based. Requires use of ES from stress periods.
• Korea FSS: Financial institutions with approved internal models can use their own VaR models but must submit regular backtesting reports.

Production Operations Checkpoint

• Automatically record daily VaR calculation results and trigger alerts when there are sharp changes compared to the previous day (e.g., over 30% change).
• Perform backtesting on a monthly basis to continuously monitor model accuracy.
• Review model parameters (observation window, volatility model, etc.) quarterly and adjust to market conditions.
• Update stress test scenarios at least annually. Add new crisis cases (e.g., COVID pandemic, SVB incident) to scenarios.

11. Failure Cases and Recovery Procedures

Case 1: The 2008 Financial Crisis and VaR's Limitations

The 2008 global financial crisis was a historic event that demonstrated the fundamental limitations of VaR-based risk management. Before the crisis, most investment banks used 95% or 99% VaR as risk limits, but actual losses exceeded VaR by orders of magnitude.

Root Cause Analysis:

• Assumption violations: Correlations of underlying assets in subprime mortgage-related structured products (CDOs, etc.) rapidly converged to 1 during the crisis. VaR models were based on peacetime low correlations and thus grossly overestimated diversification benefits.
• Liquidity risk not reflected: VaR assumes assets can be liquidated normally, but market liquidity completely evaporated during the crisis. Theoretical VaR was meaningless in a market with no bid prices.
• Model risk: The Gaussian copula model used for CDO pricing underestimated tail dependence.

Lessons:

• VaR is not an upper bound on risk, but merely "a loss estimate under normal market conditions."
• Using VaR alone is dangerous; CVaR + stress testing + liquidity risk assessment must always be performed together.
• The robustness of correlation assumptions must be tested under stress scenarios.

Case 2: Risk Underestimation Due to Implementation Errors

In practice, VaR model failures occur more frequently from implementation bugs than theoretical limitations. The following are error patterns that easily occur in real-world implementations.

Common bugs:

• Return calculation direction error: Computing (P_{t-1} - P_t) / P_t instead of (P_t - P_{t-1}) / P_{t-1} produces errors in both sign and magnitude.
• Annualization mistake: Multiplying daily volatility by 252 instead of sqrt(252) overestimates volatility by approximately 16 times.
• Array index error: Confusing ascending/descending order in VaR percentile calculations can cause the minimum loss to be reported as VaR -- a serious error.
• Timezone mismatch: Ignoring different closing time points for US and European stocks distorts correlations.

Recovery procedure:

1. Immediately inspect the data pipeline when VaR changes abruptly (missing data, outliers, duplicates)
2. Cross-verify manual calculations against code results for a single asset
3. Unit test code with synthetic data generated from known distributions (e.g., standard normal)
4. Reproduce VaR for specific past dates to verify temporal consistency
5. Cross-validate with an independent second implementation (or vendor system)

Case 3: Failure to Detect Volatility Regime Transition

In early March 2020 during the COVID pandemic onset, the VIX index exceeded 80, rendering VaR models trained on 2018-2019 data powerless overnight. Such extreme volatility was not included in the 250-day (1-year) observation window.

Response measures:

• Apply GARCH(1,1) model to immediately reflect volatility spikes from recent days into VaR.
• Use longer windows (500 days, 750 days) in parallel to include past crisis data such as 2008.
• Add regime-switching (Markov regime-switching) capability to the VaR model to automatically detect whether the current state is a high-volatility or low-volatility regime.
• Monitor external volatility indicators like VIX and proactively adjust risk limits when they spike.

12. Checklist

Check the following items when building a VaR-based risk management system.

Model Construction Phase

• Confirm return calculation method (log returns vs arithmetic returns, choose appropriately for the use case)
• Validate data quality (missing values, outliers, dividend/split adjustment)
• Calculate at least 2 VaR methodologies in parallel (Parametric + Historical, or Historical + Monte Carlo)
• Calculate CVaR (Expected Shortfall) alongside to understand tail risk
• Perform sensitivity analysis on observation window length (compare 250, 500, 750 days)
• Verify covariance matrix is positive definite (test Cholesky decomposition feasibility)

Backtesting Phase

• Confirm Kupiec test (POF test) passes
• Analyze breach pattern clustering (Christoffersen independence test)
• Verify Green Zone maintenance under Basel Traffic Light system
• Ensure at least 250 trading days backtesting period

Production Operations Phase

• Build daily VaR automated calculation and recording pipeline
• Set up automatic alert mechanism for sudden VaR changes (threshold exceedance vs previous day)
• Auto-generate monthly backtesting reports
• Establish quarterly model parameter review schedule
• Annual stress test scenario updates
• Unit tests on code changes (synthetic data-based VaR reproduction tests)
• Regular cross-validation with a second system or vendor

Governance Phase

• Document risk limit framework and VaR linkage methodology
• Establish escalation procedures for VaR exceedances
• Establish model change management process
• Register VaR model in model risk inventory and conduct periodic reviews

13. References

• Investopedia: Value at Risk (VaR) Explained - Suitable introductory material for understanding VaR fundamentals. Provides intuitive explanations of key parameters such as confidence level and holding period.
• Risk Engineering: Value at Risk - Academic-level material on the mathematical definition and three methodologies of VaR. Includes discussion of limitations and alternatives.
• Interactive Brokers: Risk Metrics in Python - VaR and CVaR Guide - Practical guide for implementing VaR and CVaR in Python. Provides code examples from data collection using yfinance to backtesting.
• PyPortfolioOpt GitHub Repository - Python implementations of various portfolio optimization algorithms including mean-variance, HRP, and Black-Litterman. Can be used for risk-integrated optimization.
• PyQuant News: Quickly Compute Value at Risk with Monte Carlo - A concise newsletter summarizing the essentials of Monte Carlo simulation-based VaR calculation. Includes practical advice on Cholesky decomposition and simulation convergence.
• Hull, J.C. (2022). Options, Futures, and Other Derivatives. 11th Edition. Pearson. - The standard financial engineering textbook, rigorously covering the mathematical foundations of three methodologies in the VaR chapter.
• Jorion, P. (2006). Value at Risk: The New Benchmark for Managing Financial Risk. 3rd Edition. McGraw-Hill. - The most comprehensive monograph on VaR, covering the full range from methodology to backtesting to regulatory application.