【Quant】Options Pricing with Monte Carlo Simulation

Find theoretical price for options via Monte Carlo simulation

Photo by Rishi Jhajharia on Unsplash

Highlight

• Difficulty：★★★★☆
• Price options via Monte Carlo simulation
• Introduce variance reduction technique for enhancing princing result

Preface

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.

Programming environment and Module required

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_base="http://10.10.10.66"
tejapi.ApiConfig.api_key = api_key
tejapi.ApiConfig.ignoretz = True

Database

Stock trading database: Unadjusted daily stock price, database code is (TWN/APRCD).
Derivatives database: Options daily transaction information, database code is (TWN/AOPTION).

Import data

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']
    }
)

# set index
stocks = stocks.set_index('mdate')
puts = puts.set_index('mdate')
calls = calls.set_index('mdate')

Data processing

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()

Stock price prediction

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:

• S0: Asset price at initial date
• r: Asset`s historical return
• sigma: Asset`s standard deviation of historical return
• T: Time to maturity
• Nsteps: Numbers of time interval
• Nrep: Numbers of stock path

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 SPATH

After 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()

figure 1: Monte Carlo simulation on stock price

Options pricing

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 【Quant】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.

figure 2: Comparison between the price of MC and BS

Antithetic variate method

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 P0

Next, 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))

figure 3: Comparison between AV and MC

Control variate method

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:

• CallOrPut: Choose between call and put
• K: Strike price
• S0: Asset price at initial time point
• r: Asset`s historical return
• sigma: Asset`s standard deviation of historical return
• T: Time to maturity
• Nsteps: Numbers of time interval
• Nrep: Numbers of stock path
• Npilot: Time period for obtaining optimal c*

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 CVP

By 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.

figure 4: Comparison between AV and CV

Practical example

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.

figure 5: Comparison for all method and TEJ BS price

Conclusion

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 E Shop to build their own 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.

Source code

Extended reading

• 【Quant】CRR Model
• 【Quant】Black Scholes model and Greeks

Related link

• TEJ API
• TEJ E-Shop

You could give us encouragement by …
We will share financial database applications every week.
If you think today’s article is good, you can click on the applause icon once.
If you think it is awesome, you can hold the applause icon until 50 times.
Any feedback is welcome, please feel free to leave a comment below.