Table of Contents
Monte Carlo simulation has been widely adopted in the field of financial research. In 【Quant(19)】Prediction of Portfolio Performance, we have introduced how to use Monte Carlo simulation for stock price prediction. In today`s article, we will extend the application to more complex options pricing. In 【Quant】CRR Model and 【Quant】Black Scholes model and Greeks, we respectively used the principles of binomial trees and the Black-Scholes formula to calculate theoretical prices of options. These articles also explained many fundamental concepts related to options. For those who have a limited understanding of options, it is recommended to read these two articles first before continuing with the current one. In the subsequent sections, we will demonstrate how to use Monte Carlo simulation to predict stock prices. We will then extend the discussion to pricing European options and finally introduce some variance reduction techniques to assist us in using Monte Carlo simulation.
Windows 11 and Jupyter Notebook is used as editor
# Import required packages
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
import time
import tejapi
plt.style.use('bmh')
# Log in TEJ API
api_key = 'YOUR_KEY'
tejapi.ApiConfig.api_key = api_key
tejapi.ApiConfig.ignoretz = TrueStock trading database: Unadjusted daily stock price, database code is (TWN/APRCD).
Derivatives database: Options daily transaction information, database code is (TWN/AOPTION).
Using the unadjusted closing prices of the Taiwan Weighted Stock Index (Y9999) within the time period from January 31, 2021, to April 19, 2023. We will also load the Taiwan Weighted Index call and put options (TXO202304C15500, TXO202304P15500). These options are European-style options, with a start trading date of January 31 and an expiration date of April 19. The strike price is set at 15,500. Furthermore, we will set the “mdate” (date) column as the index.
# Import required data
gte, lte = '2021-03-16', '2023-04-20'
stocks = tejapi.get('TWN/APRCD', # stock price
paginate = True,
coid = 'Y9999',
mdate = {'gte':gte, 'lte':lte},
opts = {
'columns':[ 'mdate','close_d']
}
)
# Get options price
puts = tejapi.get( # puts price
'TWN/AOPTION',
paginate = True,
coid = 'TXO202304P15500',
mdate = {'gte':gte, 'lte':lte},
opts = {
'columns':['mdate', 'coid','settle', 'kk', 'theoremp', 'acls', 'ex_price', 'td1y', 'avolt', 'rtime']
}
)
calls = tejapi.get( # calls price
'TWN/AOPTION',
paginate = True,
coid = 'TXO202304C15500',
mdate = {'gte':gte, 'lte':lte},
opts = {
'columns':['mdate', 'coid','settle', 'kk', 'theoremp', 'acls', 'ex_price', 'td1y', 'avolt', 'rtime']
}
)
Calculating daily return and moving volatility, set 252 days as the window.
# Calculate daily return
stocks['daily return'] = np.log(stocks['close_d']) - np.log(stocks['close_d'].shift(1))
stocks['moving volatility'] = stocks['daily return'].rolling(252).std()We utilize Monte Carlo simulation to simulate stock price paths. The concept of Monte Carlo simulation is quite simple. It involves obtaining the return process of the asset and discretizing it, then using small time intervals to calculate the changes in asset prices. For example, considering stock prices, their returns follow a Geometric Brownian motion. Thus, we can obtain a discretized stochastic differential equation (Equation 1), where Wt represents a Wiener process. After applying Itô’s formula, we obtain Equation 2 as the main equation for Monte Carlo simulation to predict stock prices, where Zt follows a standard normal distribution.
Next, we can use Python to program the above equation. The essence of the Monte Carlo method lies in simultaneously estimating multiple stock price paths using Equation 2. Finally, by summing and averaging the last stock price of each path, we obtain the predicted stock price. Here, we define the following variables:
We codify the above equation into “mc_asset” function for executing Monte Carlo simulation.
def mc_asset(S0, r, sigma, T, Nsteps, Nrep):
SPATH = np.zeros((Nrep, 1 + Nsteps))
SPATH[:, 0] = S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)
for i in range(0,Nrep):
for j in range(0,Nsteps):
SPATH[i,j+1] = SPATH[i,j] * np.exp(nudt + sidt * np.random.normal())
return SPATHAfter setting up the function, we can set the arguments to testify whether the function can work or not. The result can be visualized as figure 1. Each line in figure 1 represents a stock simulation path.
S0 = 100
K = 110
CallOrPut = 'call'
r = 0.03
sigma = 0.25
T = 0.5
Nsteps = 10000
Nrep = 1000
SPATH = mc_asset(S0, r, sigma, T, Nsteps, Nrep)
plt.figure(figsize = (10,8))
for i in range(len(SPATH)):
plt.plot(SPATH[i])
plt.xlabel('Numbers of steps')
plt.ylabel('Stock price')
plt.title('Monte Carlo Simulation for Stock Price')
plt.show()
We can use the aforementioned method to predict the stock price at options expiration. Then, we calculate the intrinsic value of the option at expiration for each path. Finally, by discounting the intrinsic value back to present value and taking the average, we can obtain the theoretical option price. Please refer to the specific code below:
def mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep):
SPATH = mc_asset(S0, r, sigma, T, Nsteps, Nrep)
if CallOrPut == 'call':
payoffs = np.maximum(SPATH[:,-1] - K, 0)
return np.mean(payoffs)*np.exp(-r*T)
else:
payoffs = np.maximum(K - SPATH[:,-1], 0)
return np.mean(payoffs)np.exp(-rT)We can choose to calculate the price of call or put by simply modifying the argument of “CallOrPut”. Further, we provide all arguments appropriate value for verification. Moreover, we compare the result with the theoretical price calculated from Black Scholes model and Greeks.
S0 = 100
K = 110
CallOrPut = 'put'
r = 0.03
sigma = 0.25
T = 0.5
Nsteps = 10000
Nrep = 1000
p_ = mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep)
mybs = BS_formula(S0, K, r, sigma, T)
c, p = mybs.BS_price()
print(f'Monte Carlo price: {c_}')
print(f'Black Scholes price: {p}')Results are displayed down below in figure 2. As you can see, there is little difference between the price of Monte Carlo simulation and Black-Scholes formula.
Start Predict Market Pricing with High-Quality Investment Database by TEJ!
Due to the nature of Monte Carlo simulation in option pricing, which involves generating a large number of theoretical prices and taking the average, there is a potential issue with higher volatility in the simulated prices. This means that the simulated prices may exhibit extreme values, leading to a deviation in the simulation results. To address this problem, we can employ the Antithetic Variate method to reduce volatility.
The concept behind this method is to generate a price path for the underlying asset (Equation 3) and simultaneously generate a path with opposite returns (Equation 4). In this case, the correlation between the two paths is -1, resulting in the minimum covariance when combining the two paths. This, in turn, reduces the volatility in estimating the option price.
The specific code implementation is as follows. We generate two matrices for calculation: SPATH1 represents the forward path, while SPATH2 represents the path in the opposite direction. It can be observed that in each iteration when generating random numbers (epsilon), both matrices share the same random numbers. As a result, the computational effort is reduced, and the volatility in predicting option prices is decreased.
def mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep):
SPATH1 = np.zeros((int(Nrep/2), 1 + Nsteps))
SPATH2 = np.zeros((int(Nrep/2), 1 + Nsteps))
SPATH1[:, 0], SPATH2[:, 0] = S0, S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)
for i in range(0,int(Nrep/2)):
for j in range(0,Nsteps):
epsilon = np.random.normal()
SPATH1[i,j+1] = SPATH1[i,j] * np.exp(nudt + sidt * epsilon)
SPATH2[i,j+1] = SPATH2[i,j] * np.exp(nudt - sidt * epsilon)
if CallOrPut == 'call':
C1 = np.maximum(SPATH1[:, -1] - K, 0)
C2 = np.maximum(SPATH2[:, -1] - K, 0)
C = np.mean(0.5 * (C1 + C2))
C0 = np.exp(-r*T) * C
return C0
else:
P1 = np.maximum(K - SPATH1[:, -1], 0)
P2 = np.maximum(K - SPATH2[:, -1], 0)
P = np.mean(0.5 * (P1 + P2))
P0 = np.exp(-r*T) * P
return P0Next, we input the values and verify whether the option prices obtained using the Antithetic Variate method are consistent with the results from the regular Monte Carlo simulation. The results can be seen in Figure 3, and it can be observed that they are indeed very close to each other.
CallOrPut = 'put'
K = 110
S0 = 100
r = 0.03
sigma = 0.25
T = 0.5
Nrep = 10000
Nsteps = 1000
print('Price under AV: ', mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Price under MC: ', mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
In addition to the Antithetic Variate method, we can also reduce the volatility of option theoretical prices using the Control Variate method. Suppose we have two random variables, X and Y, where calculating the mean and variance of variable Y is straightforward. Let’s assume that these two variables can be combined to form a new variable, Z (Equation 5). In this case, the expected value of Z is the same as the expected value of X (as shown in Equation 6), while the variance is determined by the parameter c. Therefore, we can find an optimal value for c* that minimizes the variance of Z, as shown in Equation 7. We consider X and Y as the option and underlying stock prices for each path. By using the covariance between historical option and stock prices and the variance of stock prices, we can calculate the optimal c* and then compute the option theoretical price (E(Z)) using Equation 5. This allows us to achieve the goal of reducing the volatility in Monte Carlo simulations.
The programming result is shown as follows:
def mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot):
# Calculate covariance between stock and options price
SPATH = np.zeros((Npilot, 1 + Nsteps))
SPATH[:, 0] = S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)
for i in range(0,Npilot):
for j in range(0,Nsteps):
SPATH[i,j+1] = SPATH[i,j] * np.exp(nudt + sidt * np.random.normal())
Sn = SPATH[:, -1]
if CallOrPut == 'call':
Cn = np.maximum(SPATH[:,-1] - K, 0) * np.exp(-r*T)
MatCov = np.cov(Sn, Cn)[0,1]
VarY = S0 ** 2 * np.exp(2 * r * T) * (np.exp(T * sigma ** 2) - 1)
c = -MatCov / VarY
ExpY = S0 * np.exp(r*T)
else:
Pn = np.maximum(K - SPATH[:,-1], 0) * np.exp(-r*T)
MatCov = np.cov(Sn, Pn)[0,1]
VarY = S0 ** 2 * np.exp(2 * r * T) * (np.exp(T * sigma ** 2) - 1)
c = -MatCov / VarY
ExpY = S0 * np.exp(r*T)
# Applying control variate function with optimal c*
SPATH2 = np.zeros((Nrep, 1 + Nsteps))
SPATH2[:, 0] =S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)
for i in range(0,Nrep):
for j in range(0,Nsteps):
SPATH2[i,j+1] = SPATH2[i,j] * np.exp(nudt + sidt * np.random.normal())
S = SPATH2[:, -1]
if CallOrPut == 'call':
C = np.maximum(SPATH2[:,-1] - K, 0) * np.exp(-r*T)
CVC = np.mean(C + c * (S - ExpY))
return CVC
else:
P = np.maximum(K - SPATH2[:,-1], 0) * np.exp(-r*T)
CVP = np.mean(P + c * (S - ExpY))
return CVPBy inputting numbers into arguments, we testify the difference between options price from Antithetic Variate and Control Variate method. The result is shown in figure 4. We can say they are identical.
CallOrPut = 'put'
K = 110
S0 = 100
r = 0.03
sigma = 0.25
T = 0.5
Nrep = 5000
Nsteps = 1000
Npilot= 5000
print('Price under AV: ', mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Price under CV: ', mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot))
In summary, we have learned three methods for calculating option theoretical prices using Monte Carlo simulation: the conventional method, the Antithetic Variate method, and the Control Variate method. Next, we will introduce real-life examples to compare the theoretical prices obtained from each method with the Black-Scholes prices calculated by TEJ, to see if there are significant differences between them.
S0 = stocks.loc['2023-01-31']['close_d']
K = 15500
r = stocks['daily return'].rolling(252).mean().loc['2023-01-31'] # average return of stock
T = 51 / 252
sigma = stocks.loc['2023-01-31']['moving volatility'] * np.sqrt(252)
Nstep = 50000
Nrep = 50000
Npilot = 5000
CallOrPut = 'put'
print('Monte Carlo theoretical price: ', mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Monte Carlo with AV theoretical price: ', mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Monte Carlo with CV theoretical price: ', mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot))
print('TEJ Black Scholes price: ', puts.loc['2023-01-31']['theoremp'])
print('Real price: ', puts.loc['2023-01-31']['settle'])We use a Taiwan Stock Exchange put with a strike price of 15,500 and a time period from January 1, 2023, to April 19, 2023. Sigma is calculated as the standard deviation of Taiwan Stock Exchange returns with a window of 252 days, and r is the average return over the past 252 days. We assume today is January 31, 2023. The results are shown in figure 5, where we can observe that the prices obtained using the three methods are closer to the actual prices compared to TEJ’s calculated prices.
Monte Carlo pricing method relies more on the law of large numbers compared to the CRR model and Black-Scholes model. It gradually approximates reasonable theoretical prices by simulating a large number of stock price paths. In today’s era of significant advancements in computer performance, data-driven algorithms and pricing models are expected to become increasingly prevalent. When engaging in options trading, investors may also consider incorporating the Monte Carlo method into their considerations. In the future, we will provide more knowledge related to options and derivative products. Please continue to follow us. Additionally, readers and investors are welcome to purchase solutions from TEJ Quantitative Solutions to obtain a comprehensive option pricing programs. Please note that the pricing formulas and options products mentioned above are for demonstration purposes only and do not constitute any investment or recommendation by us.
