Quickly compute Value at Risk with Monte Carlo

Value at risk (VaR) is a tool professional traders use to manage risk. It estimates how much a portfolio might lose, given normal market conditions, over a set time period.

There are three ways to compute VaR: the parametric method, the historical method, and the Monte Carlo method.

In contrast to the parametric and historical methods which are backward looking, Monte Carlo is forward looking.

In today’s newsletter, you’ll be able to simulate the equity curve of a portfolio of ETFs and compute the VaR using the Monte Carlo Method.

If you’re ready, let’s go!

Quickly compute Value at Risk with Monte Carlo

The Monte Carlo method for computing VaR involves generating a large number (sometimes millions) of hypothetical simulations of the evolution of a portfolio.

The scenarios correspond to a potential future value of a portfolio considering asset prices, weights, and covariances.

By simulating these scenarios, we can estimate the distribution of the portfolio’s future value.

From there, we can figure out VaR by identifying the worst losses at a specific confidence level, usually the 95th percentile.

The VaR represents the maximum amount of money we can expect to lose over a set time horizon.

Luckily for us, Python makes it easy to do.

Imports and set up

We only need three libraries for today’s issue. NumPy for the linear algebra, pandas for building DataFrames, and OpenBB for data. We’ll start by downloading data for a mock portfolio of 25 sector ETFs.

import numpy as np
import pandas as pd
from openbb import obb

sectors = [
    "XLE", 
    "XLF", 
    "XLU", 
    "XLI", 
    "GDX", 
    "XLK", 
    "XLV", 
    "XLY", 
    "XLP", 
    "XLB", 
    "XOP", 
    "IYR", 
    "XHB", 
    "ITB", 
    "VNQ", 
    "GDXJ", 
    "IYE", 
    "OIH", 
    "XME", 
    "XRT", 
    "SMH", 
    "IBB", 
    "KBE", 
    "KRE", 
    "XTL", 
]

data = obb.equity.price.historical(
    sectors, 
    start_date="2022-01-01", 
    provider="yfinance"
).to_df()

Next, we’ll compute the historic mean returns, weights, and covariance matrix.

data["returns"] = data.groupby("symbol").close.pct_change()

portfolio_stats = data.groupby("symbol").agg(
    daily_returns=("returns", "mean"),
)

portfolio_stats["weights"] = 1 / len(sectors)

covariance_matrix = (
    data
    .pivot(
        columns="symbol", 
        values="returns"
    )
    .dropna()
    .cov()
)

In today’s example, we use a simple average of past returns to represent future expected returns. We also use a static, equal weighted portfolio. From our returns, we can use pandas to compute the covariance matrix between the sector historic returns. We’ll use the covariance matrix in the Monte Carlo to simulate the covariance between assets in the future price paths.

Set up the Monte Carlo simulation

Now that we have our historical returns and covariance matrix, we can generate the simulated price paths.

simulations = 1_000
days = len(data.index.unique())
initial_capital = 100_000

portfolio = np.zeros((days, simulations))

historical_returns = np.full(
    shape=(days, len(sectors)), 
    fill_value=portfolio_stats.daily_returns
)

We start by defining the number of simulations, the number of days, and the initial portfolio capital. We set up arrays to store the portfolio values and historical returns. Next, we run the simulation.

L = np.linalg.cholesky(covariance_matrix)

for i in range(0, simulations):
    Z = np.random.normal(size=(days, len(sectors)))
    daily_returns = historical_returns + np.dot(L, Z.T).T
    portfolio[:, i] = (
        np.cumprod(np.dot(daily_returns, portfolio_stats.weights) + 1) * initial_capital
    )

simulated_portfolio = pd.DataFrame(portfolio)

First, we calculate the Cholesky decomposition of the covariance matrix. Sounds complicated but all it does is generate correlated random variables. In the loop, we generate a normally distributed random variable and adjust it by the Cholesky factor to simulate daily returns.

These returns are then used to calculate our cumulative portfolio returns over the time period for each simulation. Finally, we use the NumPy array output to construct a pandas DataFrame.

Analyze the results

Now that we have our simulate results, we can compute VaR and conditional VaR.

alpha = 5

def montecarlo_var(alpha):
    sim_val = simulated_portfolio.iloc[-1, :]
    return np.percentile(sim_val, alpha)

def conditional_var(alpha):
    sim_val = simulated_portfolio.iloc[-1, :]
    return sim_val[sim_val <= montecarlo_var(alpha)].mean()

mc_var = montecarlo_var(alpha)
cond_var = conditional_var(alpha)

ax = simulated_portfolio.plot(lw=0.25, legend=False)
ax.axhline(mc_var, lw=0.5, c="r")
ax.axhline(cond_var, lw=0.5, c="g")

VaR is the value at the lower 5th percentile of the final portfolio values. It represents the dollar amount the portfolio will lose with 95% confidence.

The Conditional Value at Risk (CVaR) is also known as Expected Shortfall. CVaR is computed as the mean of the simulated values at the last time step that are less than or equal to VaR. Since CVaR captures the worst case losses beyond that of the point estimate provided by VaR, it’s often seen as a superior to VaR.

The plot visualizes all price paths, VaR, and CVaR. We can see that CVaR (the green line) is less than VaR which demonstrates the more conservative calculation.

Next steps

This simulation assumes static expected returns, no drift, no volatility, and a static covariance matrix. As a next step, use a look back window for the covariance matrix. You can also try introducing drift and volatility into the price paths.