Monte Carlo Method: Value at Risk (VaR) In Python

In this blog post, we will demonstrate how to perform a Value at Risk (VaR) simulation based on the Monte Carlo Method for a portfolio of stocks using Python and the yfinance library. We will use historical stock prices to calculate the expected return and standard deviation of the portfolio, and then run a Monte Carlo simulation to estimate the VaR.

Importing Libraries

import numpy as np
import pandas as pd
import datetime as dt
import yfinance as yf
import matplotlib.pyplot as plt
from scipy.stats import norm

We begin by importing the necessary libraries: numpy, pandas, datetime, yfinance, matplotlib, and scipy.stats.

Setting Time Range and Ticker List

### Set time from to a certain number of years
years = 15

endDate = dt.datetime.now()
startDate = endDate - dt.timedelta(days = 365*years)

### Create a list of tickers
tickers = ['SPY','BND','GLD','QQQ','VTI']

In this section, we set the time range for our historical data to 15 years and create a list of stock tickers for our portfolio.

Downloading Adjusted Close Prices

The code explicitly requests auto_adjust=False to retain Adj Close and multi_level_index=False for a flat single-ticker result. Use a current yfinance version supporting those arguments. Check that downloads contain the expected tickers and enough common trading dates; a failed or missing download should be corrected before interpreting VaR.

### Download the daily adjusted close prices for the tickers
adj_close_df = pd.DataFrame()

for ticker in tickers:
    data = yf.download(ticker, start=startDate, end=endDate,
                       auto_adjust=False, multi_level_index=False, progress=False)
    adj_close_df[ticker] = data['Adj Close']

print(adj_close_df)

Next, we download the daily adjusted close prices for each stock ticker using the yfinance library and store the data in a DataFrame.

Calculating Daily Log Returns

### Calculate the daily log returns and drop any NAs
log_returns = np.log(adj_close_df/adj_close_df.shift(1))
log_returns  = log_returns.dropna()

print(log_returns)

We calculate the daily log returns for each stock in our portfolio and drop any missing values.

Defining Functions to Calculate Portfolio Expected Return and Standard Deviation

### Create a function that will be used to calculate portfolio expected return
# Historical parameters are estimates, not guarantees of future returns.
def expected_return(weights, log_returns):
    portfolio_log_returns = np.log1p(np.expm1(log_returns).dot(weights))
    return portfolio_log_returns.mean()

### Create a function that will be used to calculate portfolio standard deviation
def standard_deviation(weights, log_returns):
    portfolio_log_returns = np.log1p(np.expm1(log_returns).dot(weights))
    return portfolio_log_returns.std(ddof=1)

These functions estimate daily mean and standard deviation of the portfolio log return. We first weight asset simple returns and then use log1p, assuming costless daily rebalancing to fixed weights. A weighted average of asset log returns is not exactly the portfolio log return. Historical estimates and the normal-return model can miss future regime changes and tail events.

Creating a Covariance Matrix

### Create a covariance matrix for all the securities
cov_matrix = np.expm1(log_returns).cov()
print(cov_matrix)

This simple-return covariance matrix shows the cross-asset dependence in the data. The functions above measure the actual aggregated portfolio log-return volatility used in the simulation; that is not exactly a weighted asset-log-return covariance calculation.

Calculating Portfolio Expected Return and Standard Deviation

### Create an equally weighted portfolio and find total portfolio expected return and standard deviation
portfolio_value = 1000000
weights = np.array([1/len(tickers)]*len(tickers))
portfolio_expected_return = expected_return(weights, log_returns)
portfolio_std_dev = standard_deviation(weights, log_returns)

We create an equally weighted portfolio and calculate the portfolio’s expected return and standard deviation using the functions we defined earlier.

Defining Functions for Monte Carlo Simulation

def random_z_score():
    return np.random.normal(0, 1)

### Create a function to calculate scenarioGainLoss
days = 20

def scenario_gain_loss(portfolio_value, portfolio_std_dev, z_score, days):
    scenario_log_return = (portfolio_expected_return * days
                           + portfolio_std_dev * z_score * np.sqrt(days))
    return portfolio_value * np.expm1(scenario_log_return)

The first function draws a standard-normal shock. The second models the horizon portfolio log return with mean proportional to days and volatility proportional to the square root of days, then converts it to dollar P&L using expm1. This assumes independent normal daily log returns with stable parameters. It is the Monte Carlo version of a parametric lognormal model, so with many draws it approaches that model’s analytical quantile. It does not demonstrate nonlinear option repricing or a nonnormal simulation.

Running Monte Carlo Simulation

### Run 10000 simulations
simulations = 10000
scenarioReturn = []

for i in range(simulations):
    z_score = random_z_score()
    scenarioReturn.append(scenario_gain_loss(portfolio_value, portfolio_std_dev, z_score, days))

We run 10,000 Monte Carlo simulations, calculating the scenario gain/loss for each simulation and storing the results in a list.

Calculating Value at Risk (VaR)

### Specify a confidence interval and calculate the Value at Risk (VaR)
confidence_interval = 0.99
VaR = -np.percentile(scenarioReturn, 100 * (1 - confidence_interval))
print(VaR)

We specify a confidence interval of 99% and calculate Value at Risk (VaR) using the results of our simulations.

Plotting the Results

### Plot the results of all 10000 scenarios
plt.hist(scenarioReturn, bins=50, density=True)
plt.xlabel('Scenario Gain/Loss ($)')
plt.ylabel('Frequency')
plt.title(f'Distribution of Portfolio Gain/Loss Over {days} Days')
plt.axvline(-VaR, color='r', linestyle='dashed', linewidth=2, label=f'VaR at {confidence_interval:.0%} confidence level')
plt.legend()
plt.show()

Finally, we create a histogram to visualize the distribution of portfolio gain/loss over the specified number of days. We also add a vertical dashed line to indicate the VaR at our chosen confidence level.

Conclusion