Herron Topic 5 - Simulations

FINA 6333 for Spring 2024

Author

Richard Herron

In finance, we typically perform simulations with Monte Carlo methods. Monte Carlo methods are used throughout science, engineering, and mathematics, and they are especially useful in finance due to the randomness of asset prices. This lecture notebook provides an introduction to Monte Carlo methods in finance.

The basic idea behind Monte Carlo methods in finance is to create many possible paths of financial variables (e.g., stock prices, interest rates, exchange rates) based on their historical behavior, statistical properties, and any other relevant factors. We use these paths to simulate the performance of financial instruments or portfolios, and then we aggregate these results to estimate metrics of interest. This approach helps us model and analyze complex financial problems that may be difficult or impossible to solve analytically.

We could spend a semester (or more) on simulations and Monte Carlo methods. In this lecture notebook, we will limit our focus to:

  1. Option pricing: Monte Carlo methods can be used to estimate the fair value of financial derivatives, such as options and warrants. By simulating the potential future paths of the underlying asset and calculating the payoffs of the derivative at each path, the expected payoff can be computed and discounted to present value to determine the option’s price.
  2. Value at Risk (VaR): Monte Carlo simulations can be used to compute VaR, a widely used risk management metric that estimates the potential loss in the value of a portfolio over a specific time horizon, given a certain level of confidence (e.g., 95% or 99%). By simulating numerous scenarios and observing the distribution of potential losses, the VaR can be calculated as the threshold below which a certain percentage of losses fall.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pandas_datareader as pdr
from scipy.stats import norm
import warnings
import yfinance as yf
%precision 2
pd.options.display.float_format = '{:.2f}'.format
%config InlineBackend.figure_format = 'retina'

1 Option Pricing

We can use Monte Carlo methods to value stock options. First, we simulate several hundred or thousand possible (but random) price paths for the underlying stock. Then, we calculate the payoff for each path. On some paths, the option will expire “in the money” with \(S_T > K\) and pay \(S_T - K\). On other paths, the option will expire “out of the money” with \(S_T < K\) and pay 0. We average these payoffs and discount them to today. The present value of these payoffs is the option price. This is an illustrative example, and there is a lot of depth to Monte Carlo methods for option pricing.

1.1 Simulating Stock Prices

We can simulate stock prices with the following stochastic differential equation (SDE) for Geometric Brownian Motion (GBM): \(dS = \mu S dt + \sigma S dW_t\). GBM does not account for mean-reversion and time-dependent volatility. So GBM is often used for stocks and not for bond prices, which tend to display long-term reversion to the face value. In the GBM SDE:

  1. \(S\) is the stock price
  2. \(\mu\) is the drift coefficient (i.e., instantaneous expected return)
  3. \(\sigma\) is the diffusion coefficient (i.e., volatility of the drift)
  4. \(W_t\) is the Wiener Process or Brownian Motion

The GBM SDE has a closed-form solution: \(S(t) = S_0 \exp\left(\left(\mu - \frac{1}{2}\sigma^2\right)t + \sigma W_t\right)\). We can apply this closed form solution recursively: \(S(t_{i+1}) = S(t_i) \exp\left(\left(\mu - \frac{1}{2}\sigma^2\right)(t_{i+1} - t_i) + \sigma \sqrt{t_{i+1} - t_i} Z_{i+1}\right)\). Here, \(Z_i\) is a Standard Normal random variable (\(dW_t\) are independent and normally distributed) and \(i = 0, \ldots, T-1\) is the time index.

We can use this closed form solution to simulate stock prices for AAPL.

aapl = (
    yf.download(tickers='AAPL')
    .assign(Return=lambda x: x['Adj Close'].pct_change())
    .rename_axis(columns=['Variable'])
)
[*********************100%%**********************]  1 of 1 completed
aapl.describe()
Variable Open High Low Close Adj Close Volume Return
count 10926.00 10926.00 10926.00 10926.00 10926.00 10926.00 10925.00
mean 21.13 21.35 20.91 21.14 20.37 319740189.00 0.00
std 43.79 44.26 43.35 43.82 43.42 335918295.20 0.03
min 0.05 0.05 0.05 0.05 0.04 0.00 -0.52
25% 0.30 0.30 0.29 0.30 0.24 114434300.00 -0.01
50% 0.52 0.53 0.51 0.52 0.42 207588200.00 0.00
75% 19.45 19.60 19.23 19.40 16.89 399854000.00 0.01
max 198.02 199.62 197.00 198.11 197.86 7421640800.00 0.33

We will use returns from 2021 to predict prices in 2022 (i.e., train in 2021, and test in 2022).

train = aapl.loc['2021']
test = aapl.loc['2022']

We will use the following function to simulate price paths. Throughout this lecture notebook, we will use one-trading-day steps (i.e., dt=1).

def simulate_gbm(S_0, mu, sigma, n_steps, dt=1, seed=42):
    '''
    Function to simulate stock prices following Geometric Brownian Motion (GBM).
    
    Parameters
    ------------
    S_0 : float
        Initial stock price
    mu : float
        Drift coefficient
    sigma : float
        Diffusion coefficient
    n_steps : int
        Length of the forecast horizon in time increments, so T = n_steps * dt
    dt : int
        Time increment, typically one day
    seed : int
        Random seed for reproducibility

    Returns
    -----------
    S_t : np.ndarray
        Array (length: n_steps + 1) of simulated prices
    '''

    np.random.seed(seed)
    dW = np.random.normal(scale=np.sqrt(dt), size=n_steps)
    W = dW.cumsum()
    
    t = np.linspace(dt, n_steps * dt, n_steps)
    
    S_t = S_0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * W)
    S_t = np.insert(S_t, 0, S_0)
    
    return S_t

Now we will simulate price paths. Here is one simulated price path:

simulate_gbm(
    S_0=train['Adj Close'].iloc[-1],
    mu=train['Return'].pipe(np.log1p).mean(),
    sigma=train['Return'].pipe(np.log1p).std(),
    n_steps=test.shape[0]
)
array([175.34, 176.91, 176.71, 178.72, 183.26, 182.78, 182.3 , 187.1 ,
       189.59, 188.38, 190.21, 189.02, 187.83, 188.75, 183.32, 178.58,
       177.18, 174.55, 175.61, 173.29, 169.64, 173.8 , 173.36, 173.73,
       170.04, 168.76, 169.24, 166.36, 167.53, 166.12, 165.53, 164.14,
       169.2 , 169.34, 166.71, 169.07, 166.01, 166.74, 161.82, 158.63,
       159.29, 161.33, 161.94, 161.81, 161.21, 157.66, 156.04, 155.07,
       157.85, 158.88, 154.67, 155.63, 154.85, 153.36, 155.02, 157.73,
       160.24, 158.29, 157.69, 158.68, 161.32, 160.27, 159.97, 157.36,
       154.58, 156.74, 160.31, 160.3 , 163.03, 164.14, 162.64, 163.75,
       167.96, 168.04, 172.43, 165.61, 167.95, 168.36, 167.74, 168.16,
       163.14, 162.74, 163.84, 167.89, 166.69, 164.75, 163.62, 166.18,
       167.22, 166.  , 167.53, 167.97, 170.74, 169.03, 168.34, 167.47,
       163.82, 164.76, 165.61, 165.8 , 165.36, 161.88, 160.97, 160.27,
       158.42, 158.18, 159.37, 164.36, 164.99, 165.84, 165.82, 161.03,
       161.14, 161.46, 168.05, 167.71, 168.69, 168.78, 165.87, 169.07,
       171.27, 173.61, 171.31, 175.34, 171.67, 173.46, 179.76, 177.15,
       175.76, 176.22, 175.01, 170.95, 171.32, 168.64, 170.09, 167.81,
       172.16, 170.22, 169.53, 171.91, 168.78, 169.56, 173.28, 169.12,
       169.79, 170.67, 172.97, 169.8 , 166.47, 168.03, 169.  , 169.85,
       170.96, 169.31, 170.11, 171.08, 169.34, 174.59, 176.09, 172.99,
       174.98, 172.48, 174.83, 178.25, 176.14, 179.03, 180.39, 182.94,
       188.71, 188.18, 186.15, 183.74, 181.58, 181.55, 182.72, 183.72,
       186.33, 186.57, 191.1 , 190.51, 199.09, 201.28, 198.78, 195.65,
       197.36, 196.87, 199.31, 201.02, 201.  , 198.54, 194.05, 192.88,
       195.72, 196.59, 192.96, 193.69, 195.08, 192.58, 193.25, 193.63,
       190.37, 191.65, 193.56, 197.11, 200.63, 196.52, 193.83, 195.62,
       197.43, 199.25, 211.99, 214.13, 218.24, 221.79, 224.32, 223.44,
       226.38, 223.87, 223.27, 221.79, 222.32, 230.84, 224.37, 227.05,
       221.57, 220.16, 224.22, 224.68, 221.12, 218.86, 221.46, 219.15,
       220.13, 220.53, 218.5 , 226.27, 228.79, 221.81, 222.7 , 220.62,
       223.85, 221.29, 221.13, 223.13, 226.45, 222.43, 221.49, 220.06,
       218.03, 224.44, 226.12, 221.89])

We will combine simulate_gbm() with a list comprehension and pd.concat() to simulate and collect many price paths. To simplify this combination, we will write a helper function simulate_gbm_series() that:

  1. Returns a series
  2. Helps us vary the seed argument (otherwise seed will remain the default value of 42, and every returned price path will be the same)
def simulate_gbm_series(seed, train=train, test=test):
    S_t = simulate_gbm(
        S_0=train['Adj Close'].iloc[-1],
        mu=train['Return'].pipe(np.log1p).mean(),
        sigma=train['Return'].pipe(np.log1p).std(),
        n_steps=test.shape[0],
        seed=seed
    )
    return pd.Series(data=S_t, index=test.index.insert(0, train.index[-1]))
n = 100
S_t = pd.concat(
    objs=[simulate_gbm_series(seed=seed) for seed in range(n)],
    axis=1,
    keys=range(n),
    names=['Simulation']
)
S_t
Simulation 0 1 2 3 4 5 6 7 8 9 ... 90 91 92 93 94 95 96 97 98 99
Date
2021-12-31 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 ... 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34
2022-01-03 180.49 180.09 174.37 180.56 175.67 176.75 174.66 180.28 175.78 175.53 ... 174.90 173.85 176.17 180.11 177.60 173.69 175.72 170.73 176.84 175.13
2022-01-04 181.83 178.55 174.40 182.00 177.25 176.02 176.87 179.15 179.03 174.91 ... 174.46 171.14 177.21 183.46 182.11 172.28 175.77 169.76 181.67 181.11
2022-01-05 184.86 177.25 168.79 182.47 174.66 183.11 177.67 179.43 173.78 172.03 ... 173.02 169.84 175.35 189.74 180.28 173.39 174.13 167.94 184.76 182.12
2022-01-06 191.72 174.45 173.41 177.36 176.77 182.57 175.35 180.78 170.20 172.18 ... 171.03 169.28 168.40 191.34 185.54 175.23 175.30 170.82 182.90 186.18
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2022-12-23 262.97 307.38 208.18 242.58 279.49 226.84 237.21 207.60 215.97 252.47 ... 238.62 348.72 236.58 238.55 210.66 239.32 221.16 182.12 212.46 250.85
2022-12-27 267.77 305.56 209.65 250.91 278.46 226.17 234.01 205.90 212.35 258.89 ... 232.39 348.33 238.03 240.24 211.97 239.87 225.12 183.42 212.66 250.37
2022-12-28 264.63 314.76 211.58 255.46 275.50 223.94 238.34 205.77 211.07 255.26 ... 236.72 346.96 241.94 237.04 211.12 232.35 234.97 185.88 210.80 248.57
2022-12-29 258.84 308.66 210.69 256.40 274.30 225.66 237.77 201.47 210.91 258.48 ... 241.33 361.64 238.74 238.67 210.80 227.92 237.16 180.01 207.84 250.33
2022-12-30 261.26 307.30 209.97 263.75 275.18 225.39 242.19 201.15 212.93 261.03 ... 236.42 370.19 241.31 238.46 204.33 228.72 236.61 182.26 212.64 248.27

252 rows × 100 columns

Below, we prefix the simulated price path column names with _ to hide them from the legend. However, this feature triggers a warning 100 times! We will suppress these 100 warnings with the warnings package. Generally, we should avoid suppressing warnings, however it is the easiest option here.

fig, ax = plt.subplots(1,1)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    S_t.add_prefix('_').plot(alpha=0.1, ax=ax)
S_t.mean(axis=1).plot(label='Mean', ax=ax)
aapl.loc[S_t.index, ['Adj Close']].plot(label='Observed', ax=ax)
plt.legend()
plt.ylabel('Price ($)')
plt.title(
    'Apple Simulated and Observed Prices' + 
    f'\nTrained from {train.index[0]:%Y-%m-%d} to {train.index[-1]:%Y-%m-%d}'
)
plt.show()

1.2 Pricing Stock Options

We can use simulated price paths to price options! We will use the Black and Scholes (1973) formula as a benchmark. Black and Scholes (1973) provide a closed form (i.e., analytic) solution to price European options.

def price_bs(S_0, K, T, r, sigma, type='call'):
    '''
    Function used for calculating the price of European options using the analytical form of the Black-Scholes model.
    
    Parameters
    ------------
    S_0 : float
        Initial stock price
    K : float
        Strike price
    T : float
        Time to expiration in days
    r : float
        Daily risk-free rate
    sigma : float
        Standard deviation of daily stock returns
    type : str
        Type of the option. Allowable: ['call', 'put']
    
    Returns
    -----------
    option_premium : float
        The premium on the option calculated using the Black-Scholes model
    '''
    
    d1 = (np.log(S_0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if type == 'call':
        val = (norm.cdf(d1, 0, 1) * S_0) - (norm.cdf(d2, 0, 1) * K * np.exp(-r * T))
    elif type == 'put':
        val = (norm.cdf(-d2, 0, 1) * K * np.exp(-r * T)) - (norm.cdf(-d1, 0, 1) * S_0)
    else:
        raise ValueError('Wrong input for type!')
        
    return val

We can use the AAPL parameters above to price a European call option on AAPL stock. We will calculate its price at the end of 2021 with an expiration at the end of 2022, assuming a 5% risk-free rate.

S_0 = train['Adj Close'].iloc[-1]
K = 100
T = test.shape[0]
r = 0.05/252
sigma = train['Return'].pipe(np.log1p).std()
S_0
175.34
_ = price_bs(
    S_0=S_0,
    K=K,
    T=T,
    r=r,
    sigma=sigma
)

print(f'Black and Scholes (1973) option price: {_:0.2f}')
Black and Scholes (1973) option price: 80.28

To simulate the Black and Scholes (1973) option price, we must simulate AAPL prices with the same 5% risk-free rate as the drift. We will write a new helper function to use the same inputs as above.

def simulate_gbm_series(seed, S_0=S_0, T=T, r=r, sigma=sigma, train=train, test=test):
    S_t = simulate_gbm(
        S_0=S_0,
        mu=r,
        sigma=sigma,
        n_steps=T,
        seed=seed
    )
    return pd.Series(data=S_t, index=test.index.insert(0, train.index[-1]))

We will simulate more price paths to increase the precision of out option price.

n = 10_000
S_t = pd.concat(
    objs=[simulate_gbm_series(seed=seed) for seed in range(n)],
    axis=1,
    keys=range(n),
    names=['Simulation']
)
S_t
Simulation 0 1 2 3 4 5 6 7 8 9 ... 9990 9991 9992 9993 9994 9995 9996 9997 9998 9999
Date
2021-12-31 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 ... 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34 175.34
2022-01-03 180.31 179.91 174.20 180.38 175.49 176.58 174.49 180.10 175.61 175.36 ... 177.19 175.01 173.26 173.03 172.07 171.35 172.64 175.64 176.66 173.93
2022-01-04 181.47 178.20 174.06 181.64 176.90 175.67 176.53 178.79 178.67 174.57 ... 176.51 175.87 174.64 175.25 170.52 170.45 177.24 174.25 175.47 175.58
2022-01-05 184.31 176.73 168.29 181.93 174.15 182.57 177.15 178.90 173.27 171.53 ... 179.24 174.44 175.30 173.07 166.81 168.76 179.92 174.76 175.79 175.73
2022-01-06 190.97 173.77 172.73 176.67 176.08 181.85 174.66 180.07 169.53 171.51 ... 182.15 175.25 174.19 175.99 167.06 171.90 181.76 175.84 171.63 174.39
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2022-12-23 206.32 241.17 163.34 190.32 219.29 177.98 186.11 162.88 169.45 198.09 ... 155.35 165.02 136.69 250.88 143.47 128.72 241.58 203.06 193.55 185.52
2022-12-27 209.89 239.51 164.33 196.67 218.27 177.28 183.42 161.39 166.44 202.92 ... 159.47 163.95 137.86 254.08 144.77 133.00 239.79 203.96 186.88 191.02
2022-12-28 207.22 246.48 165.68 200.04 215.73 175.36 186.64 161.13 165.28 199.88 ... 157.21 163.16 139.14 258.26 145.83 130.34 239.68 212.17 190.77 191.76
2022-12-29 202.49 241.46 164.82 200.58 214.58 176.54 186.01 157.61 164.99 202.20 ... 155.79 164.04 138.34 260.85 143.42 128.22 241.18 216.13 186.57 188.86
2022-12-30 204.18 240.17 164.10 206.12 215.06 176.14 189.28 157.20 166.41 204.00 ... 155.60 163.55 134.08 267.53 140.05 130.25 243.21 222.57 186.76 191.31

252 rows × 10000 columns

We can compare this price to a simulated price. The payoff on the call option is \(S_T - K\) or \(0\), whichever is higher. The price of the option is the present value of the mean payoff, discounted at the risk-free rate.

payoff = np.maximum(S_t.iloc[-1] - K, 0)
print(f'Simulated option price: {payoff.mean() * np.exp(-r * T):0.2f}')
Simulated option price: 80.98

The option prices do not match exactly. However, we can simulate more price paths to bring our simulated option price closer to the analytic solution.

2 Estimating Value-at-Risk using Monte Carlo

Value-at-Risk (VaR) measures the risk associated with a portfolio VaR reports the worst expected loss, at a given level of confidence, over a certain horizon under normal market conditions. For example, say the 1-day 95% VaR of our portfolio is 100 dollars. This implies that that 95% of the time (under normal market conditions), we should not lose more than 100 dollars over one day. We typically present VaR as a positive value, so a VaR of 100 dollars implies a loss of less than 100 dollars.

We can calculate VaR several ways, including:

  • Parametric Approach (Variance-Covariance)
  • Historical Simulation Approach
  • Monte Carlo simulations

We only consider the last method to calculate the 1-day VaR of an portfolio of 20 shares each of META and GOOG.

tickers = ['GOOG', 'META']
shares = np.array([20, 20])
T = 1
n_sims = 10_000

However, we will download all data from Yahoo! Finance and subset our data later.

df = (
    yf.download(tickers=tickers)
    .rename_axis(columns=['Variable', 'Ticker'])
)
[*********************100%%**********************]  2 of 2 completed

Next, we calculate daily returns during 2022. Choosing the window to define “normal market conditions” is part art, part science, and beyod the scope of this lecture notebook.

returns = df['Adj Close'].pct_change().loc['2022']
returns
Ticker GOOG META
Date
2022-01-03 0.00 0.01
2022-01-04 -0.00 -0.01
2022-01-05 -0.05 -0.04
2022-01-06 -0.00 0.03
2022-01-07 -0.00 -0.00
... ... ...
2022-12-23 0.02 0.01
2022-12-27 -0.02 -0.01
2022-12-28 -0.02 -0.01
2022-12-29 0.03 0.04
2022-12-30 -0.00 0.00

251 rows × 2 columns

We will need the variance-covariance matrix.

cov_mat = returns.cov()
cov_mat * 100**2
Ticker GOOG META
Ticker
GOOG 5.96 6.74
META 6.74 16.38

We will use the variance-covariance matrix to calculate the Cholesky decomposition.

chol_mat = np.linalg.cholesky(cov_mat)
chol_mat
array([[0.02, 0.  ],
       [0.03, 0.03]])

The Cholesky decomposition helps us generate random variables with the same variance and covariance as the observed data.

rv = np.random.normal(size=(n_sims, len(tickers)))
correlated_rv = (chol_mat @ rv.T).T
correlated_rv
array([[ 0.04,  0.05],
       [ 0.02,  0.02],
       [-0.02, -0.04],
       ...,
       [-0.04,  0.01],
       [-0.  , -0.09],
       [-0.02,  0.  ]])

These random variables have a variance-covariance matrix similar to the real data.

np.cov(correlated_rv.T)  * 100**2
array([[ 5.97,  6.85],
       [ 6.85, 16.47]])
np.allclose(cov_mat, np.cov(correlated_rv.T), rtol=0.05)
True
np.mean(correlated_rv, axis=0) * 100
array([0.02, 0.02])
returns.mean().values * 100
array([-0.16, -0.32])

Here are the parameters for the simulated price paths:

mu = returns.mean().values
sigma = returns.std().values
S_0 = df.loc['2021', 'Adj Close'].iloc[-1].values
P_0 = S_0.dot(shares)

Calculate terminal prices using the GBM formula above:

S_T = S_0 * np.exp((r - 0.5 * sigma ** 2) * T + sigma * np.sqrt(T) * correlated_rv)

S_T
array([[144.79, 336.46],
       [144.74, 336.03],
       [144.61, 335.27],
       ...,
       [144.53, 335.94],
       [144.65, 334.6 ],
       [144.6 , 335.84]])

Calculate terminal portfolio values and returns. Note that these are dollar values, since VaR is typically expressed in dollar values.

P_T = S_T.dot(shares)
P_T
array([9625.09, 9615.38, 9597.54, ..., 9609.35, 9584.95, 9608.85])
P_diff = P_T - P_0
P_diff
array([ 11.63,   1.92, -15.92, ...,  -4.11, -28.51,  -4.61])
P_diff.mean()
-4.37

Next, we calculate VaR.

percentiles = [0.01, 0.05, 0.1]
var = np.percentile(P_diff, percentiles)

for x, y in zip(percentiles, var):
    print(f'1-day VaR with {100-x}% confidence: ${-y:.2f}')
1-day VaR with 99.99% confidence: $47.97
1-day VaR with 99.95% confidence: $44.61
1-day VaR with 99.9% confidence: $42.59

Finally, we will plot VaR:

fig, ax = plt.subplots()
ax.hist(P_diff, bins=100, density=True)
ax.set_title(f'Distribution of 1-Day Changes in Portfolio Value\n from {n_sims} Simulations')
ax.axvline(x=var[2], color='red', ls='--')
ax.text(x=var[2], y=1, s='99% 1-Day VaR', color='red', ha='right', va='top', rotation=90, transform=ax.get_xaxis_transform())
ax.set_ylabel('Density')
ax.set_xlabel('1-Day Change in Portfolio Value ($)')
plt.show()