Derivative Security Pricing Project
Exotic Underperformance Option Pricing
Technical stack: Python, NumPy, pandas, Matplotlib, SciPy
1. Data and Environment Setup
Import libraries and data of:
- End-of-day bid and ask prices of options on the S&P 500 index, traded on the CBOE on 9 March 2022
- Daily closing prices of S&P 500 index from 9 March 2021 to 8 March 2023
- End-of-day bid and ask prices of options on Moderna Inc., traded on CBOE on 9 March 2022
- Daily closing prices on the Moderna Inc. stock from 9 March 2021 to 9 March 2023
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.optimize import fsolve
from scipy.stats import norm
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
mrna = pd.read_excel("MRNA.xlsx",)
mrna["Date"] = pd.to_datetime(mrna["Date"], format="%d/%m/%Y")
mrna['Price'] = mrna['Close']
mrna = mrna[['Date', 'Price']]
mrna = mrna.sort_values(by = 'Date')
mrna = mrna.reset_index(drop = True)
print("The data of MRNA stocks:")
display(mrna)
The data of MRNA stocks:
| Date | Price | |
|---|---|---|
| 0 | 2021-03-09 | 130.87 |
| 1 | 2021-03-10 | 129.75 |
| 2 | 2021-03-11 | 140.47 |
| 3 | 2021-03-12 | 136.99 |
| 4 | 2021-03-15 | 143.66 |
| ... | ... | ... |
| 500 | 2023-03-03 | 143.20 |
| 501 | 2023-03-06 | 144.03 |
| 502 | 2023-03-07 | 141.05 |
| 503 | 2023-03-08 | 142.08 |
| 504 | 2023-03-09 | 137.36 |
505 rows × 2 columns
sp500 = pd.read_excel("SP500.xlsx", skiprows=2)
sp500["Date"] = pd.to_datetime(sp500["Date"], format="%b %d, %Y")
sp500['Price'] = sp500['Close*'] # Do not remove the "*"
sp500 = sp500[['Date', 'Price']]
sp500 = sp500.sort_values(by = 'Date')
sp500 = sp500.reset_index(drop = True)
print("The data of SP500 Index:")
display(sp500)
The data of SP500 Index:
| Date | Price | |
|---|---|---|
| 0 | 2021-03-09 | 3875.44 |
| 1 | 2021-03-10 | 3898.81 |
| 2 | 2021-03-11 | 3939.34 |
| 3 | 2021-03-12 | 3943.34 |
| 4 | 2021-03-15 | 3968.94 |
| ... | ... | ... |
| 499 | 2023-03-02 | 3981.35 |
| 500 | 2023-03-03 | 4045.64 |
| 501 | 2023-03-06 | 4048.42 |
| 502 | 2023-03-07 | 3986.37 |
| 503 | 2023-03-08 | 3992.01 |
504 rows × 2 columns
mrna_quote = pd.read_excel("mrna_quotedata20220309.xlsx")
mrna_quote["Expiration Date"] = pd.to_datetime(mrna_quote["Expiration Date"], format = "%a %b %d %Y")
list_date = mrna_quote['Expiration Date'].unique()
delete_index_list = []
for t in list_date:
subset = mrna_quote[(mrna_quote['Expiration Date'] == t)]
if len(subset) > 1:
list_strike = subset['Strike'].unique()
for s in list_strike:
subset2 = subset[(subset['Strike'] == s)]
if len(subset2) > 1:
max_bid_c = subset2['Bid'].max()
max_bid_p = subset2['Bid.1'].max()
min_ask_c = subset2['Ask'].min()
min_ask_p = subset2['Ask.1'].min()
list_index = subset2.index
mrna_quote.loc[list_index[0], 'Bid'] = max_bid_c
mrna_quote.loc[list_index[0], 'Bid.1'] = max_bid_p
mrna_quote.loc[list_index[0], 'Ask'] = min_ask_c
mrna_quote.loc[list_index[0], 'Ask.1'] = min_ask_p
delete_index_list.extend(list_index[1:])
mrna_quote = mrna_quote.drop(delete_index_list)
mrna_quote['Call Price'] = (mrna_quote['Ask'] + mrna_quote['Bid']) / 2
mrna_quote['Put Price'] = (mrna_quote['Ask.1'] + mrna_quote['Bid.1']) / 2
mrna_quote['Call Volume'] = mrna_quote['Volume']
mrna_quote['Put Volume'] = mrna_quote['Volume.1']
mrna_quote = mrna_quote[['Expiration Date', 'Call Price', 'Put Price', 'Strike', 'Call Volume', 'Put Volume']]
mrna_quote = mrna_quote.sort_values(by = 'Expiration Date')
mrna_quote = mrna_quote.reset_index(drop = True)
print("The data of MRNA options:")
display(mrna_quote)
The data of MRNA options:
| Expiration Date | Call Price | Put Price | Strike | Call Volume | Put Volume | |
|---|---|---|---|---|---|---|
| 0 | 2022-03-11 | 67.700 | 0.005 | 75.0 | 0 | 54 |
| 1 | 2022-03-11 | 0.710 | 12.500 | 155.0 | 712 | 7 |
| 2 | 2022-03-11 | 1.070 | 11.125 | 152.5 | 1783 | 1 |
| 3 | 2022-03-11 | 1.310 | 10.850 | 150.0 | 5541 | 250 |
| 4 | 2022-03-11 | 2.340 | 8.500 | 149.0 | 295 | 11 |
| ... | ... | ... | ... | ... | ... | ... |
| 640 | 2024-01-19 | 52.900 | 38.400 | 135.0 | 8 | 0 |
| 641 | 2024-01-19 | 48.575 | 42.775 | 140.0 | 40 | 7 |
| 642 | 2024-01-19 | 47.000 | 42.550 | 145.0 | 12 | 1 |
| 643 | 2024-01-19 | 43.000 | 50.675 | 155.0 | 2 | 2 |
| 644 | 2024-01-19 | 4.180 | 498.000 | 640.0 | 3 | 0 |
645 rows × 6 columns
sp500_quote = pd.read_excel("spx_quotedata20220309_all.xlsx")
sp500_quote["Expiration Date"] = pd.to_datetime(sp500_quote["Expiration Date"], format = "%a %b %d %Y")
list_date = sp500_quote['Expiration Date'].unique()
delete_index_list = []
for t in list_date:
subset = sp500_quote[(sp500_quote['Expiration Date'] == t)]
if len(subset) > 1:
list_strike = subset['Strike'].unique()
for s in list_strike:
subset2 = subset[(subset['Strike'] == s)]
if len(subset2) > 1:
max_bid_c = subset2['Bid'].max()
max_bid_p = subset2['Bid.1'].max()
min_ask_c = subset2['Ask'].min()
min_ask_p = subset2['Ask.1'].min()
list_index = subset2.index
sp500_quote.loc[list_index[0], 'Bid'] = max_bid_c
sp500_quote.loc[list_index[0], 'Bid.1'] = max_bid_p
sp500_quote.loc[list_index[0], 'Ask'] = min_ask_c
sp500_quote.loc[list_index[0], 'Ask.1'] = min_ask_p
delete_index_list.extend(list_index[1:])
sp500_quote = sp500_quote.drop(delete_index_list)
sp500_quote['Call Price'] = (sp500_quote['Ask'] + sp500_quote['Bid']) / 2
sp500_quote['Put Price'] = (sp500_quote['Ask.1'] + sp500_quote['Bid.1']) / 2
sp500_quote['Call Volume'] = sp500_quote['Volume']
sp500_quote['Put Volume'] = sp500_quote['Volume.1']
sp500_quote = sp500_quote[['Expiration Date', 'Call Price', 'Put Price', 'Strike', 'Call Volume', 'Put Volume']]
sp500_quote = sp500_quote.sort_values(by = 'Expiration Date')
sp500_quote = sp500_quote.reset_index(drop = True)
print("The data of SP500 options:")
display(sp500_quote)
The data of SP500 options:
| Expiration Date | Call Price | Put Price | Strike | Call Volume | Put Volume | |
|---|---|---|---|---|---|---|
| 0 | 2022-03-09 | 875.800 | 0.025 | 3400 | 0 | 1 |
| 1 | 2022-03-09 | 0.025 | 117.100 | 4395 | 332 | 2 |
| 2 | 2022-03-09 | 0.025 | 124.600 | 4400 | 2209 | 56 |
| 3 | 2022-03-09 | 0.025 | 129.600 | 4405 | 369 | 17 |
| 4 | 2022-03-09 | 0.025 | 134.600 | 4410 | 1287 | 15 |
| ... | ... | ... | ... | ... | ... | ... |
| 4503 | 2025-12-19 | 3913.400 | 3.400 | 200 | 0 | 1 |
| 4504 | 2025-12-19 | 2341.400 | 95.300 | 2000 | 0 | 10 |
| 4505 | 2025-12-19 | 2178.400 | 117.550 | 2200 | 0 | 10 |
| 4506 | 2025-12-19 | 801.450 | 595.000 | 4200 | 0 | 50 |
| 4507 | 2026-12-18 | 239.450 | 1631.000 | 6000 | 0 | 10 |
4508 rows × 6 columns
2. Assumptions
The following assumptions have been made:
The daily observations in the data are equally spaced, even though due to weekends and public holidays, this is just an approximation.
Risk-free interest rates are deterministic, but time-dependent. Thus, the continuously compounded short rate $r(t)$ is a deterministic function of $t$.
The underlying S&P 500 index pays dividends, which can be adequately approximated by a continuous dividend rate $q(t)$, which is also deterministic and time-dependent.
Where necessary, we use loglinear interpolation of zero coupon bond prices and "dividend discount factors" (the latter are defined below). I.e., given $B(0,T_1)$ and $B(0,T_2)$,for $T_1 < T < T_2$ loglinear interpolation for $B(0,T)$,
$$B(0,T) = B(0,T_1)^{\frac{T_2 - T}{T_2 - T_1}} \cdot B(0,T_2)^{\frac{T - T_1}{T_2 - T_1}}$$
- The dynamics of the S&P 500 index and the Modenra Inc. stock price are represented by the stochastic processes $X$ and $Y$, respectively, where $X$ and $Y$ follow the dynamics (under the USD risk-neutral measure) of:
$$dX(t) = X(t)((r(t) - q(t)) dt + \sigma (t) dW(t))$$ $$dY(t) = Y(t)(r(t) dt + \eta (t) dW(t))$$
Here, $W(t)$ is a two-dimensional vector-valued Brownian motion with independent components, and $\sigma(t)$ and $\eta(t)$ are two-dimensional (time-dependent) vectors. For the specific numerical calculations, we assume that
$$\sigma(t) = \nu(t) \begin{bmatrix} \sigma_1 \\ 0 \end{bmatrix} \quad \qquad \eta(t) = \xi(t) \begin{bmatrix} \eta_1 \\ \eta_2 \end{bmatrix} $$
Where $\sigma^2_1 = 1$ and $\eta_1^2 + \eta_2^2 = 1$.
3. Stage 1: Extracting Term Structure Using Put/Call Parity
Determine discount factors and dividend yields for the S&P 500 index and Moderna Inc. stock using put-call parity, and plot the term structure.
For this stage, we have removed all dates with missing observations - the entire row of data is removed. This is because when minimizer is used, it requires at least two observations to estimate results.
today = pd.to_datetime("2022-03-09")
def Stage_1():
list_date = sp500_quote[sp500_quote['Expiration Date'] > today + pd.DateOffset(days=6)]["Expiration Date"].unique()
list_valid_date = []
for date in list_date:
option_temp = sp500_quote[(sp500_quote["Expiration Date"] == date) & (sp500_quote["Call Volume"] >0) & (sp500_quote["Put Volume"] > 0)]
if len(option_temp) > 1:
list_valid_date.append(date)
D = []
B = []
S_0 = sp500[sp500['Date'] == today]['Price'].values[0]
def target(x, call, put, K):
D, B = x
error = D*S_0 - K*B - (call - put)
return np.sum(error**2)
for date in list_valid_date:
option_temp = sp500_quote[(sp500_quote["Expiration Date"] == date) & (sp500_quote["Call Volume"] >0) & (sp500_quote["Put Volume"] > 0)]
K = option_temp["Strike"].to_numpy()
call = option_temp["Call Price"].to_numpy()
put = option_temp["Put Price"].to_numpy()
D_result, B_result = minimize(target, [0.9, 0.9], args=(call, put, K), bounds=((0, 2), (0, 2))).x
D.append(D_result)
B.append(B_result)
df = pd.DataFrame({
'Maturity': list_valid_date,
'Zero Coupon Prices': B,
'Dividend Discount Factors': D})
display(df)
def year_transform(x):
return (x - today).days / 365
list_date_value = [year_transform(x) for x in list_valid_date]
list_date_value = [0] + list_date_value
B = [1] + B
D = [1] + D
log_B_interpolate = interp1d(list_date_value, np.log(B), kind='linear')
log_D_interpolate = interp1d(list_date_value, np.log(D), kind='linear')
def B_cal (T):
return np.exp(log_B_interpolate(T))
def D_cal (T):
return np.exp(log_D_interpolate(T))
T = np.linspace(0, list_date_value[-1], 100)
B_plot = B_cal(T)
D_plot = D_cal(T)
plt.figure(figsize=(12, 5))
plt.plot(T, B_plot, label='B(0,T)')
plt.plot(T, D_plot, label='D(0,T)')
plt.xlabel('In Years')
plt.ylabel('Value')
plt.title('Zero Coupon Prices B(0,T) and Dividend Discount Factors D(0,T) Over Time')
plt.legend()
plt.grid(True)
plt.show()
return B_cal, D_cal, year_transform, list_valid_date
B_cal, D_cal, year_transform, term_structure_day = Stage_1()
| Maturity | Zero Coupon Prices | Dividend Discount Factors | |
|---|---|---|---|
| 0 | 2022-03-16 | 0.999686 | 0.998793 |
| 1 | 2022-03-18 | 0.999688 | 0.998771 |
| 2 | 2022-03-21 | 0.999703 | 0.998702 |
| 3 | 2022-03-23 | 1.000410 | 0.999398 |
| 4 | 2022-03-25 | 0.999732 | 0.998665 |
| 5 | 2022-03-28 | 0.999388 | 0.998318 |
| 6 | 2022-03-30 | 0.998308 | 0.997145 |
| 7 | 2022-03-31 | 0.999415 | 0.998180 |
| 8 | 2022-04-01 | 0.999491 | 0.998292 |
| 9 | 2022-04-04 | 0.999160 | 0.997932 |
| 10 | 2022-04-06 | 0.998730 | 0.997408 |
| 11 | 2022-04-08 | 0.998834 | 0.997373 |
| 12 | 2022-04-14 | 0.999441 | 0.997949 |
| 13 | 2022-04-22 | 0.998503 | 0.996984 |
| 14 | 2022-04-29 | 0.998756 | 0.997205 |
| 15 | 2022-05-20 | 0.999111 | 0.996814 |
| 16 | 2022-05-31 | 0.996427 | 0.993924 |
| 17 | 2022-06-17 | 0.998533 | 0.995587 |
| 18 | 2022-06-30 | 0.998167 | 0.995262 |
| 19 | 2022-07-15 | 0.993915 | 0.991039 |
| 20 | 2022-08-19 | 0.995478 | 0.992735 |
| 21 | 2022-09-16 | 0.995024 | 0.992178 |
| 22 | 2022-09-30 | 0.992713 | 0.990207 |
| 23 | 2022-10-21 | 0.993521 | 0.991494 |
| 24 | 2022-12-16 | 0.990846 | 0.989205 |
| 25 | 2022-12-30 | 0.988554 | 0.987751 |
| 26 | 2023-01-20 | 0.989062 | 0.989030 |
| 27 | 2023-02-17 | 0.984697 | 0.985646 |
| 28 | 2023-03-17 | 0.985942 | 0.987298 |
| 29 | 2023-06-16 | 0.977896 | 0.982033 |
| 30 | 2023-12-15 | 0.969584 | 0.979622 |
4. Stage 2: Piecewise Constant Volatility
Determine an implied term structure of volatility such that the Black/Scholes prices resulting from this volatility equal the observed prices.
For each maturity greater than six days in the index option and Moderna Inc. op-tion data sets, choose strikes closest to the forward price for that maturity. Usingcall option “mid” prices for these strikes, determine an implied term structure ofvolatility, i.e., determine a piecewise constant functionsν(t) andξ(t) such thatthe Black/Scholes prices resulting from this volatility equal the observed prices.
Forward Price
The forward price $F(0,T)$ for the option maturity $T$, with the current price $S(0)$ and continuously compounded domestic interest rate $r(t)$, is given by: $$\boxed{ F(0,T) = \frac{S(0) D(0,T)}{B(0,T)} = S(0) e^{ \int_0^T(r(s) - q(s))\,ds} }$$
For SP500
Using Ito-Lemma, we have
$$ \ln X(T) = \ln X(0) + \int_0^T \big(r(s) - q(s) - \frac{1}{2} \sigma(s)^\top \sigma(s)\big)\, ds + \int_0^T \sigma (t)^\top\, dW(s) $$
Or
$$X(T) = \frac{X(0) D(0,T)}{B(0,T)} \exp\!\left\{-\frac{1}{2} \int_0^T \sigma(s)^\top \sigma(s)\, ds + \int_0^T \sigma(s)^\top\, dW(s)\right\} $$
We can define an at-the-money SP500 call option formula as follows: $$ C(0, X(0)) = B(0,T)\,\mathbb{E}_{\mathbb{Q}}\!\left[ (X(T) - K)^+ \right] = B(0,T)\,\mathbb{E}_{\mathbb{Q}}\!\left[ X(T)\,\mathbf{1}_{\{X(T) > K\}} \right] - B(0,T)\,K\,\mathbb{E}_{\mathbb{Q}}\!\left[\,\mathbf{1}_{\{X(T) > K\}}\right] $$
Using Radon-Nikodym measure change technique to change measure from $\mathbb{Q}$ to $\mathbb{\tilde{Q}}$: $$ \frac{d\tilde{\mathbb{Q}}}{d\mathbb{Q}} = \frac{X(T) B(0,T)}{X(0)}\,\mathbf{C} = D(0,T)\, \exp\!\left\{ -\frac{1}{2}\!\int_0^T \sigma(s)^\top \sigma(s)\, ds + \int_0^T \sigma(s)^\top\,dW(s) \right\} \mathbf{C_1} $$
is an exponential martingale when $\mathbf{C_1} = \frac{1}{D(0,T)}$
Using Girsanov's theorem, we have $ d\tilde{W}(t) = dW(t) - \sigma(t)\,dt$, we have the following result:
$$ \ln X(T) \;\overset{\mathbb{\tilde{Q}}}{\sim}\; \mathcal{N}\!\left( \ln\!\left(\frac{X(0)D(0,T)}{B(0,T)}\right) +\frac{1}{2}\!\int_0^T \sigma(s)^\top \sigma(s)\,ds,\; \int_0^T \sigma(s)^\top \sigma(s)\,ds \right) $$
It is given that the piecewise constant function of $\sigma(t) = \sigma(t) = \nu(t) \begin{pmatrix} \sigma_1 \\ 0\end{pmatrix} $
$$ \int_0^T \sigma(s)^\top \sigma(s)\,ds = \int_0^T \nu(s)^2\,ds = \int_0^{T_1} \nu(s)^2\,ds + \int_{T_1}^{T_2} \nu(s)^2\,ds + ... + \int_{T_i}^{T} \nu(s)^2\,ds $$
Assuming that $\nu(t) = \mathbb{c}_i$ when $T_i \le t < T_{i+1}$, we have:
$$ \int_0^T \sigma(s)^\top \sigma(s)\,ds = \sum_{j=0}^{i-1}\int_{T_j}^{T_{j+1}} \mathbb{c_j}^2\,ds + \int_{T_i}^T \mathbb{c_i}\,ds $$
However since we do not let $T$ reach $T_{i+1}$, we simply drop the second integral, and thus:
$$ \int_0^T \sigma(s)^\top \sigma(s)\,ds = \sum_{j=0}^{i-1} \mathbb{c_j}^2 (T_{j+1} - T_j) $$
The call option formula is as follows:
$$ C(0, X(0)) = B(0,T_i)\,\mathbb{E}_{\mathbb{\tilde{Q}}}\!\left[ X(T_i)\,\mathbf{1}_{\{X(T_i) > K_i \}} \frac{X(0)}{X(T_i) B(0,T_i)} \frac{1}{\mathbf{C_1}} \right] - B(0,T_i)\,K\,\mathbb{E}_{\mathbb{Q}}\!\left[\,\mathbf{1}_{\{X(T_i) > K_i\}}\right] $$
$$\boxed{ C(0, X(0)) = D(0,T_i) X(0) \Phi(d_{1,i}) - B(0,T_i) K \Phi(d_{2,i})} $$
$$\boxed{ d_{1,i} = \frac{ \ln\!\left(\dfrac{X(0)\,D(0,T_i)}{B(0,T_i)\,K_i}\right) + \frac{1}{2}\displaystyle\sum_{j=0}^{i-1} \mathbb{c_j}^{2}\,\big(T_{j+1}-T_j\big) }{ \sqrt{\displaystyle\sum_{j=0}^{i-1} \mathbb{c_j}^{2}\,\big(T_{j+1}-T_j\big)} }} $$
$$\boxed{ d_{2,i} = \frac{ \ln\!\left(\dfrac{X(0)\,D(0,T_i)}{B(0,T_i)\,K_i}\right) - \frac{1}{2}\displaystyle\sum_{j=0}^{i-1} \mathbb{c_j}^{2}\,\big(T_{j+1}-T_j\big) }{ \sqrt{\displaystyle\sum_{j=0}^{i-1} \mathbb{c_j}^{\,2}\,\big(T_{j+1}-T_j\big)} }} $$
$\quad$
For Moderna stock
Using Ito-Lemma, we have: $$ \ln Y(T) = \ln X(0) + \int_0^T \big(r(s) - \frac{1}{2} \eta(s)^\top \eta(s)\big)\, ds + \int_0^T \eta (t)^\top\, dW(s) $$
Or
$$Y(T) = \frac{Y(0)}{B(0,T)} \exp\!\left\{-\frac{1}{2} \int_0^T \eta(s)^\top \eta(s)\, ds + \int_0^T \eta(s)^\top\, dW(s)\right\} $$
We can define an at-the-money Moderna call option formula as follows:
$$ C(0, Y(0)) = B(0,T)\,\mathbb{E}_{\mathbb{Q}}\!\left[ (Y(T) - K)^+ \right] = B(0,T)\,\mathbb{E}_{\mathbb{Q}}\!\left[ Y(T)\,\mathbf{1}_{\{Y(T) > K\}} \right] - B(0,T)\,K\,\mathbb{E}_{\mathbb{Q}}\!\left[\,\mathbf{1}_{\{Y(T) > K\}}\right] $$
Using Radon-Nikodym measure change technique to change measure from $\mathbb{Q}$ to $\hat{\mathbb{Q}}$: $$ \frac{d\hat{\mathbb{Q}}}{d\mathbb{Q}} = \frac{Y(T) B(0,T)}{Y(0)}\,\mathbf{C} = \exp\!\left\{ -\frac{1}{2}\!\int_0^T \eta(s)^\top \eta(s)\, ds + \int_0^T \eta(s)^\top\,dW(s) \right\} \mathbf{C_2} $$
is an exponential martingale when $\mathbf{C_2} = 1$
Using Girsanov's theorem, we have $ d\hat{W}(t) = dW(t) - \eta(t)\,dt$, we have the following result:
$$ \ln Y(T) \;\overset{\mathbb{\hat{Q}}}{\sim}\; \mathcal{N}\!\left( \ln\!\left(\frac{Y(0)}{B(0,T)}\right) +\frac{1}{2}\!\int_0^T \eta(s)^\top \eta(s)\,ds,\; \int_0^T \eta(s)^\top \eta(s)\,ds \right) $$
The call option formula is as follows (note that since Moderna stock does not pay dividends, its American call option is typically the same as an European call option):
$$ C(0, Y(0)) = B(0,T_i)\,\mathbb{E}_{\mathbb{\hat{Q}}}\!\left[ Y(T_i)\,\mathbf{1}_{\{Y(T_i) > K_i \}} \frac{Y(0)}{Y(T_i) B(0,T_i)} \frac{1}{\mathbf{C_2}} \right] - B(0,T_i)\,K\,\mathbb{E}_{\mathbb{Q}}\!\left[\,\mathbf{1}_{\{Y(T_i) > K_i\}}\right] $$
$$\boxed{C(0, Y(0)) = Y(0) \Phi(d_{1,i}) - B(0,T_i) K \Phi(d_{2,i})} $$
$$\boxed{ d_{1,i} = \frac{ \ln\!\left(\dfrac{Y(0)}{B(0,T_i)\,K_i}\right) + \frac{1}{2}\displaystyle\sum_{j=0}^{i-1} \mathbb{p}_j^{\,2}\,\big(T_{j+1}-T_j\big) }{ \sqrt{\displaystyle\sum_{j=0}^{i-1} \mathbb{p}_j^{\,2}\,\big(T_{j+1}-T_j\big)} }} $$
$$\boxed{ d_{2,i} = \frac{ \ln\!\left(\dfrac{Y(0)}{B(0,T_i)\,K_i}\right) - \frac{1}{2}\displaystyle\sum_{j=0}^{i-1} \mathbb{p}_j^{\,2}\,\big(T_{j+1}-T_j\big) }{ \sqrt{\displaystyle\sum_{j=0}^{i-1} \mathbb{p}_j^{\,2}\,\big(T_{j+1}-T_j\big)} }} $$
Upon reviewing the initial results of Stage 2, several assumptions were made and corresponding code adjustments were implemented. These assumptions are:
All Moderna options with nonzero volumes in both puts and calls were removed.
When the forward price matches two strikes on the same maturity date, the mid-price is computed as the average of the highest bid and lowest ask for that date.
In cases where the implied volatility for a later period is lower than the prior period, contrary to the monotonicity condition implied by the derived formula, the volatility for that period is set equal to the previous period.
We also filter the list of maturities based on the
term_structure_daydefined in Stage 1. This means that we only calculate implied volatilities for maturities that exist withinterm_structure_day.After generating the implied volatilities for both the S&P 500 and Moderna, we merged their maturities using
sorted()andset(). We then completed the piecewise constant function table usingbfill(), which fills missing values by carrying forward the next available data within the period.
def Stage_2():
def iv_compute(price, df, if_index):
list_date_dat1 = sp500_quote[sp500_quote['Expiration Date'] > today + pd.DateOffset(days=6)]["Expiration Date"].unique()
list_date_dat2 = mrna_quote[mrna_quote['Expiration Date'] > today + pd.DateOffset(days=6)]["Expiration Date"].unique()
current_price = (price[price['Date'] == today]['Price'].values[0])
if if_index == 'True':
list_date = list_date_dat1
else:
list_date = list_date_dat2[list_date_dat2.isin(list_date_dat1)].unique()
list_valid_date = []
info_df = pd.DataFrame({'Call':[], 'Strike':[], 'Discount Factor':[], 'Dividend Factor':[]})
for date in list_date:
option_temp = df[(df["Expiration Date"] == date) & (df["Call Volume"] > 0) & (df["Put Volume"] > 0)]
if (len(option_temp) > 0) and (date <= term_structure_day[-1]):
list_valid_date.append(date)
for date in list_valid_date:
option_temp = df[(df["Expiration Date"] == date) & (df["Call Volume"] > 0) & (df["Put Volume"] > 0)].copy()
T = year_transform(date)
if if_index == 'True':
forward_price = current_price * D_cal(T) / B_cal(T)
else:
forward_price = current_price / B_cal(T)
option_temp['Moneyness'] = np.abs(option_temp['Strike'] - forward_price)
atm_option = option_temp.loc[option_temp['Moneyness'].idxmin()]
info_df = pd.concat([info_df, pd.DataFrame({
'Call': [atm_option['Call Price']],
'Strike': [atm_option['Strike']],
'Discount Factor': [B_cal(T)],
'Dividend Factor': [D_cal(T)]
})], ignore_index=True)
def target(iv):
B = info_df["Discount Factor"].to_numpy(dtype=float)
D = info_df["Dividend Factor"].to_numpy(dtype=float) if if_index == 'True' else np.ones(len(info_df))
K = info_df["Strike"].to_numpy(dtype=float)
call = info_df["Call"].to_numpy(dtype=float)
logarg = (current_price * D) / (K * B)
d_1 = (np.log(logarg) + 0.5 * iv) / np.sqrt(iv)
d_2 = d_1 - np.sqrt(iv)
price = D * current_price * norm.cdf(d_1) - K * B * norm.cdf(d_2)
return call - price
x0 = 0.002*np.ones(len(info_df))
iv_cal = fsolve(target, x0=x0)
for i in range(1, len(iv_cal)):
if iv_cal[i] < iv_cal[i-1]:
iv_cal[i] = iv_cal[i-1]
component_list = np.diff(iv_cal)
component_list = np.insert(component_list, 0, iv_cal[0])
t = [year_transform(x) for x in list_valid_date]
t = [0] + t
component = np.sqrt(component_list / np.diff(t))
for j in range(len(component)):
if component[j] == 0:
component[j] = component[j-1]
else:
component[j] = component[j]
return component, list_valid_date, info_df
sp500_iv, sp500_date, sp500_info_df = iv_compute(sp500, sp500_quote, "True")
mrna_iv, mrna_date, mrna_info_df = iv_compute(mrna, mrna_quote, "False")
merged_date = sorted(set(sp500_date + mrna_date))
df1 = pd.DataFrame({
'Maturity': sp500_date,
'Implied Volatility': sp500_iv
})
df2 = pd.DataFrame({
'Maturity': mrna_date,
'Implied Volatility': mrna_iv
})
table = pd.DataFrame({
'Maturity': merged_date
})
table = table.merge(df2[['Maturity', 'Implied Volatility']], on='Maturity', how='left')
table.rename(columns={'Implied Volatility': 'Moderna IV'}, inplace=True)
table = table.merge(df1[['Maturity', 'Implied Volatility']], on='Maturity', how='left')
table.rename(columns={'Implied Volatility': 'SP500 IV'}, inplace=True)
table['Moderna IV'] = table['Moderna IV'].bfill()
table = table.dropna(subset=['Moderna IV'])
table['Maturity'] = [year_transform(x) for x in table['Maturity']]
def variance_covariance(MaturityT, type="covariance"):
if MaturityT > table['Maturity'].max():
raise ValueError("Maturity must be before ", table['Maturity'].max())
elif MaturityT == 0:
return 0
elif type == "covariance":
v = table['Moderna IV'].to_numpy(float)
xi = table['SP500 IV'].to_numpy(float)
T = table['Maturity'].to_numpy(float)
if MaturityT < T[0]:
return v[0] * xi[0] * MaturityT
closest_T = T[T <= MaturityT].max()
idx_closest_T = np.where(table['Maturity'] == closest_T)[0][0]
result = np.array([])
for i in range(len(T[T <= closest_T])):
if i == 0:
vxi = v[i] * xi[i] * T[i]
result = np.append(result, vxi)
else:
vxi = v[i] * xi[i] * (T[i] - T[i-1])
result = np.append(result, vxi)
if MaturityT != T[-1]:
result = np.append(
result,
table.loc[idx_closest_T + 1, 'Moderna IV']
* table.loc[idx_closest_T +1 , 'SP500 IV']
* (MaturityT - closest_T))
else:
result = np.append(
result,
table.loc[idx_closest_T, 'Moderna IV']
* table.loc[idx_closest_T , 'SP500 IV']
* (MaturityT - closest_T))
result_ = np.sum(result)
return result_
elif type == "variance_index":
v = table['Moderna IV'].to_numpy(float)
T = table['Maturity'].to_numpy(float)
if MaturityT < T[0]:
return v[0] ** 2 * MaturityT
closest_T = T[T <= MaturityT].max()
idx_closest_T = np.where(table['Maturity'] == closest_T)[0][0]
result = np.array([])
for i in range(len(T[T <= closest_T])):
if i == 0:
v_temp = (v[i] ** 2) * T[i]
result = np.append(result, v_temp)
else:
v_temp = (v[i] ** 2) * (T[i] - T[i-1])
result = np.append(result, v_temp)
if MaturityT != T[-1]:
result = np.append(
result,
(table.loc[idx_closest_T + 1, 'Moderna IV'] ** 2)
* (MaturityT - closest_T))
else:
result = np.append(
result,
(table.loc[idx_closest_T, 'Moderna IV'] ** 2)
* (MaturityT - closest_T))
result_ = np.sum(result)
return result_
elif type == "variance_moderna":
xi = table['SP500 IV'].to_numpy(float)
T = table['Maturity'].to_numpy(float)
if MaturityT < T[0]:
return xi[0] ** 2 * MaturityT
closest_T = T[T <= MaturityT].max()
idx_closest_T = np.where(table['Maturity'] == closest_T)[0][0]
result = np.array([])
for i in range(len(T[T <= closest_T])):
if i == 0:
xi_temp = (xi[i] ** 2) * T[i]
result = np.append(result, xi_temp)
else:
xi_temp = (xi[i] ** 2) * (T[i] - T[i-1])
result = np.append(result, xi_temp)
if MaturityT != T[-1]:
result = np.append(
result,
(table.loc[idx_closest_T + 1, 'SP500 IV'] ** 2)
* (MaturityT - closest_T))
else:
result = np.append(
result,
(table.loc[idx_closest_T, 'SP500 IV'] ** 2)
* (MaturityT - closest_T))
result_ = np.sum(result)
return result_
display(table)
plt.figure(figsize = (12, 5))
for i in range(len(table)):
if i == 0:
plt.plot([0, table.loc[i, 'Maturity']], [table['Moderna IV'].iloc[i], table['Moderna IV'].iloc[i]], color = 'blue', label = 'Moderna')
plt.plot([0, table.loc[i, 'Maturity']], [table['SP500 IV'].iloc[i], table['SP500 IV'].iloc[i]], color = 'orange', label = 'SP500')
else:
plt.plot([table.loc[i - 1, 'Maturity'], table.loc[i, 'Maturity']], [table['Moderna IV'].iloc[i], table['Moderna IV'].iloc[i]], color = 'blue')
plt.plot([table.loc[i - 1, 'Maturity'], table.loc[i, 'Maturity']], [table['SP500 IV'].iloc[i], table['SP500 IV'].iloc[i]], color = 'orange')
if i != len(table) - 1:
plt.plot([table.loc[i, 'Maturity'], table.loc[i, 'Maturity']], [table['Moderna IV'].iloc[i], table['Moderna IV'].iloc[i+1]], '--', color = 'blue')
plt.plot([table.loc[i, 'Maturity'], table.loc[i, 'Maturity']], [table['SP500 IV'].iloc[i], table['SP500 IV'].iloc[i+1]], '--', color = 'orange')
plt.xlabel('Time to Maturity')
plt.ylabel('Volatilities')
plt.title('Implied Volatilities')
plt.grid(True)
plt.legend()
plt.show()
return table, variance_covariance
volatility_structure, cov_cal = Stage_2()
| Maturity | Moderna IV | SP500 IV | |
|---|---|---|---|
| 0 | 0.019178 | 0.767189 | 0.304656 |
| 1 | 0.024658 | 0.767189 | 0.307581 |
| 2 | 0.032877 | 0.900254 | 0.216518 |
| 3 | 0.038356 | 0.900254 | 0.304644 |
| 4 | 0.043836 | 0.900254 | 0.297142 |
| 5 | 0.052055 | 0.746459 | 0.199395 |
| 6 | 0.057534 | 0.746459 | 0.317569 |
| 7 | 0.060274 | 0.746459 | 0.278991 |
| 8 | 0.063014 | 0.746459 | 0.300477 |
| 9 | 0.071233 | 0.737439 | 0.238322 |
| 10 | 0.076712 | 0.737439 | 0.220320 |
| 11 | 0.082192 | 0.737439 | 0.326301 |
| 12 | 0.098630 | 0.741217 | 0.247077 |
| 13 | 0.120548 | 0.661346 | 0.232064 |
| 14 | 0.139726 | 0.723257 | 0.273073 |
| 15 | 0.197260 | 0.723257 | 0.256016 |
| 16 | 0.227397 | 0.681508 | 0.222435 |
| 17 | 0.273973 | 0.681508 | 0.256595 |
| 18 | 0.309589 | 0.640984 | 0.248618 |
| 19 | 0.350685 | 0.640984 | 0.241108 |
| 20 | 0.389041 | 0.764470 | 0.263931 |
| 21 | 0.446575 | 0.764470 | 0.198135 |
| 22 | 0.523288 | 0.764470 | 0.268144 |
| 23 | 0.561644 | 0.366353 | 0.258924 |
| 24 | 0.619178 | 0.366353 | 0.258924 |
| 25 | 0.695890 | 0.581112 | 0.514537 |
| 26 | 0.772603 | 0.581112 | 0.514537 |
| 27 | 0.810959 | 0.581112 | 0.514537 |
| 28 | 0.868493 | 0.581112 | 0.514537 |
5. Stage 3: Building Blocks for Option Pricing
Derive analytical expressions for standard deviation of the daily logarithmic returns of the S&P 500 index and Moderna Inc. stock, and the correlation between them, under the risk-neutral measure, and determine the parameters of the implied volatility of Moderna Inc. stock.
We first find the daily logarithmic returns of the S&P 500 index and Moderna Inc. stock, which are defined as follows:
$$ r_{SP500}(t) = \ln\!\left(\frac{S_{SP500}(t)}{S_{SP500}(t-1)}\right) \quad \text{and} \quad r_{MRNA}(t) = \ln\!\left(\frac{S_{MRNA}(t)}{S_{MRNA}(t-1)}\right) $$
where $S_{SP500}(t)$ and $S_{MRNA}(t)$ are the prices of the S&P 500 index and Moderna Inc. stock at time $t$, respectively.
Then we can determine their standard deviations and correlation under the risk-neutral measure using the following formulas:
$$\sigma_{SP500} = \sqrt{\frac{1}{N-1} \sum_{t=1}^{N} (r_{SP500}(t) - \bar{r}_{SP500})^2}$$ $$\sigma_{MRNA} = \sqrt{\frac{1}{N-1} \sum_{t=1}^{N} (r_{MRNA}(t) - \bar{r}_{MRNA})^2}$$ $$\rho = \frac{\sum_{t=1}^{N} (r_{SP500}(t) - \bar{r}_{SP500})(r_{MRNA}(t) - \bar{r}_{MRNA})}{(N-1) \sigma_{SP500} \sigma_{MRNA}}$$
where $\bar{r}_{SP500}$ and $\bar{r}_{MRNA}$ are the mean daily logarithmic returns of the S&P 500 index and Moderna Inc. stock, respectively, and $N$ is the total number of observations.
def Stage_3a():
def year_transform(x):
return (today - x).days / 365
mrna_history = mrna[mrna['Date'] < today].copy()
mrna_history.sort_values(by='Date', inplace=True)
mrna_history = mrna_history.reset_index(drop=True)
daily_logarithmic_returns_mrna = np.log(mrna_history['Price'] / mrna_history['Price'].shift(1))
sp500_history = sp500[sp500['Date'] < today].copy()
sp500_history.sort_values(by='Date', inplace=True)
sp500_history = sp500_history.reset_index(drop=True)
daily_logarithmic_returns_sp500 = np.log(sp500_history['Price'] / sp500_history['Price'].shift(1))
correlation = daily_logarithmic_returns_mrna.cov(daily_logarithmic_returns_sp500)/(daily_logarithmic_returns_mrna.std()*daily_logarithmic_returns_sp500.std())
mrna_history['Date'] = mrna_history['Date'].apply(lambda x: year_transform(x))
sp500_history['Date'] = sp500_history['Date'].apply(lambda x: year_transform(x))
print(f"Annualized Standard Deviation of Daily Logarithmic Returns for MRNA: {daily_logarithmic_returns_mrna.std()*np.sqrt(365)}")
print(f"Annualized Standard Deviation of Daily Logarithmic Returns for SP500: {daily_logarithmic_returns_sp500.std()*np.sqrt(365)}")
print(f"Annualized Covariance between Daily Logarithmic Returns of MRNA and SP500: {daily_logarithmic_returns_mrna.cov(daily_logarithmic_returns_sp500)*365}")
print(f"Annualized Variance of Daily Logarithmic Returns for MRNA: {(daily_logarithmic_returns_mrna.std()*np.sqrt(365))**2}")
print(f"Annualized Variance of Daily Logarithmic Returns for SP500: {(daily_logarithmic_returns_sp500.std()*np.sqrt(365))**2}")
print(f"Correlation between Daily Logarithmic Returns of MRNA and SP500: {correlation}")
return correlation
correlation_log_ret = Stage_3a()
Annualized Standard Deviation of Daily Logarithmic Returns for MRNA: 0.9817400532068851 Annualized Standard Deviation of Daily Logarithmic Returns for SP500: 0.17158858830252613 Annualized Covariance between Daily Logarithmic Returns of MRNA and SP500: 0.03366555162428605 Annualized Variance of Daily Logarithmic Returns for MRNA: 0.9638135320706576 Annualized Variance of Daily Logarithmic Returns for SP500: 0.02944264363565381 Correlation between Daily Logarithmic Returns of MRNA and SP500: 0.19984846826387667
Using Ito-Lemma, we have for Y:
$$ \ln Y(t+\Delta t )-\ln Y(t)= \int_{t}^{t+\Delta t} \big(r(s) - \frac{1}{2} \eta(s)^\top\eta(s)\big)\,ds + \int_{t}^{t+\Delta t} \eta(s)^\top\,dW(s) $$
$$ \ln Y(t+\Delta t )-\ln Y(t) \;\overset{\mathcal{F}_t}{\sim}\; \mathcal{N}\! \left( \int_t^{t+\Delta t}\big(r(s)-\frac{1}{2}\eta(s)^\top\eta(s)\big)\,ds, \int_t^{t+\Delta t}\eta(s)^\top\eta(s)\,ds) \right) $$
Assuming that $T_i \le t + \Delta t < T_{i+1}$. Then, we have:
$$ \sqrt{\ln Y (t+\Delta t) -\ln Y(t)|\mathcal{F_t}} = \sqrt{\int_t^{t+\Delta t}\eta(s)^\top\eta(s)\,ds} = \sqrt{\int_t^{t+\Delta t} \xi(s)^2\,ds} = \sqrt{\mathbb{p_i^2}(\Delta t)} $$
Assuming that $T_i \le t < T_{i+1} \le t+\Delta t$. Then, we have: $$ \sqrt{\text{Var}\big(\ln Y(t+\Delta t) - \ln Y(t)\big)|\mathcal{F_t}} = \sqrt{\int_t^{t+\Delta t} \eta(s)^\top\eta(s)\,ds} = \sqrt{\int_t^{T_{i+1}}\xi(s)^2\,ds + \int_{T_{i+1}}^{t+\Delta t}\xi(s)^2\,ds } = \sqrt{\mathbb{p_i^2}(T_{i+1} - t) + \mathbb{p_{i+1}^2}(t+\Delta t - T_{i+1})} $$
However, since we are only interested in daily changes, $t+\Delta t$ can't exceed $T_{i+1}$, so we can simply drop this assumption.
$\quad$
Using Ito-Lemma, we have for X:
$$ \ln X(t+\Delta t )-\ln X(t)= \int_{t}^{t+\Delta t} \big(r(s) - q(s) - \frac{1}{2} \sigma(s)^\top\sigma(s)\big)\,ds + \int_{t}^{t+\Delta t} \sigma(s)^\top\,dW(s) $$
$$ \ln X(t+\Delta t )-\ln X(t) \;\overset{\mathcal{F}_t}{\sim}\; \mathcal{N}\! \left( \int_t^{t+\Delta t}\big(r(s)-q(s)-\frac{1}{2}\sigma(s)^\top\sigma(s)\big)\,ds, \int_t^{t+\Delta t}\sigma(s)^\top\sigma(s)\,ds) \right) $$
Assuming that $T_i \le t + \Delta t < T_{i+1}$, we have:
$$ \sqrt{\ln X (t+\Delta t) -\ln X(t)|\mathcal{F_t}} = \sqrt{\int_t^{t+\Delta t}\sigma(s)^\top\sigma(s)\,ds} = \sqrt{\int_t^{t+\Delta t} \nu(s)^2\,ds} = \sqrt{\mathbb{c_i^2}(\Delta t)} $$
Assuming that $T_i \le t < T_{i+1} \le t+\Delta t$, we have: $$ \sqrt{\text{Var}\big(\ln X(t+\Delta t) - \ln X(t)\big)|\mathcal{F_t}} = \sqrt{\int_t^{t+\Delta t} \sigma(s)^\top\sigma(s)\,ds} = \sqrt{\int_t^{T_{i+1}}\nu(s)^2\,ds + \int_{T_{i+1}}^{t+\Delta t}\nu(s)^2\,ds } = \sqrt{\mathbb{c_i^2}(T_{i+1} - t) + \mathbb{c_{i+1}^2}(t+\Delta t - T_{i+1})} $$
Similarly, we can drop this assumption.
Thus,
$$\boxed{ \sqrt{\ln Y (t+\Delta t) -\ln Y(t)|\mathcal{F_t}} = \sqrt{\mathbb{p_i^2}(\Delta t)}} $$
And
$$\boxed{ \sqrt{\ln X (t+\Delta t) -\ln X(t)|\mathcal{F_t}} = \sqrt{\mathbb{c_i^2}(\Delta t)}} $$
We have: $ \text{Var}\big[\ln X(t+\Delta t)Y(t+\Delta t)- \ln X(t)Y(t) | \mathcal{F_t}\big] = \text{Var}\big[\ln X(t+\Delta t) - ln X(t) + \ln Y(t+\Delta t)- \ln Y(t)\big] | \mathcal{F_t}\big] $
Since X and Y are NOT independent (Moderna is in SP500) and there is no information that X and Y are indentically distributed, we can conclude that:
$ \text{Var}\big[\ln X(t+\Delta t) - ln X(t) + \ln Y(t+\Delta t)- \ln Y(t)\big] | \mathcal{F_t}\big] $
$$=\text{Var}\big[\ln X(t+\Delta t) -\ln X(t)\big] + \text{Var}\big[\ln Y(t+\Delta t) -\ln Y(t)\big] + 2 \text{Cov}\big[\ln X(t+\Delta t) -\ln X(t), \ln Y(t+\Delta t) -\ln Y(t)\big] $$
Where the covariance is:
$ \text{Cov}\big[\ln X(t+\Delta t) -\ln X(t), \ln Y(t+\Delta t) -\ln Y(t)\big] $
$ = \text{Cov} \big[\int_t^{T_{i+1}}\sigma(s)^\top\,dW(s) , \int_t^{T_{i+1}}\eta(s)^\top\,dW(s) \big] $
$ = \int_t^{T_{i+1}} \sigma(s)^\top \eta(s)\,ds $
$ = \int_t^{T_{i+1}} \nu(s) \xi(s) \begin{pmatrix} \sigma_1 & 0 \end{pmatrix}\ \begin{pmatrix} \eta_1 \\ \eta_2 \end{pmatrix}\,ds $
$ = \int_t^{T_{i+1}} \nu(s) \xi(s)\ \sigma_1 \eta_1\,ds $
Assuming that $T_i \le t + \Delta t < T_{i+1}$, we have:
$$\boxed{ \text{Var}\big[\ln X(t+\Delta t) - \ln X(t) + \ln Y(t+\Delta t)- \ln Y(t)\big] | \mathcal{F_t}\big] = (\mathbb{c_i^2} + \mathbb{p_i^2} + 2 \mathbb{c_i} \mathbb{p_i} \sigma_1 \eta_1)\Delta t } $$
We do not consider where $T_i \le t < T_{i+1} \le t+\Delta t$. However, just in case, the result for the assumption is as follows:
$$ \text{Var}\big[\ln X(t+\Delta t) - \ln X(t) + \ln Y(t+\Delta t)- \ln Y(t)\big] | \mathcal{F_t}\big] = (\mathbb{c_i^2} + \mathbb{p_i^2} + 2 \mathbb{c_i} \mathbb{p_i} \sigma_1 \eta_1) (T_{i+1}-t) + (\mathbb{c_{i+1}^2} + \mathbb{p_{i+1}^2} + 2 \mathbb{c_{i+1}} \mathbb{p_{i+1}} \sigma_1 \eta_1)(t+\Delta t - T_{i+1}) $$
Assuming that volatilities are constant in time such at $\nu(t)\equiv\ \hat{\nu}$ and $\xi(t) \equiv\ \hat{\xi}$.
From above, we have calculated the correlation $\rho_{X,Y}$ of the daily logarithmic returns from the datasets. We will use this correlation $\rho_{X,Y}$ to determine the appropriate choices for the model parameters $\eta_1$ and $\eta_2$.
$$ \text{Corr} \Big[\ln \Big(\frac{\ln X(t+\Delta t)}{\ln X(t)}\Big) , \ln \Big(\frac{\ln Y(t+\Delta t)}{\ln Y(t)}\Big)\Big] = \frac {\text{Cov} \Big[\ln \Big(\frac{\ln X(t+\Delta t)}{\ln X(t)}\Big) , \ln \Big(\frac{\ln Y(t+\Delta t)}{\ln Y(t)}\Big)\Big]} {{\sqrt{\text{Var}(\ln X(t+\Delta t) - \ln X(t))} \sqrt{\text{Var}(\ln X(t+\Delta t) - \ln X(t))}}} $$
$$\Leftrightarrow \text{Corr} \Big[\ln \Big(\frac{\ln X(t+\Delta t)}{\ln X(t)}\Big) , \ln \Big(\frac{\ln Y(t+\Delta t)}{\ln Y(t)}\Big)\Big] = \frac {\int_t^{t+\Delta t} \sigma(s)^\top \eta(s)\,ds} {{\sqrt{\int_t^{t+\Delta t}\sigma(s)^\top \sigma(s)\,ds} \sqrt{\int_t^{t+\Delta t}\eta(s)^\top \eta(s)\,ds}}} $$
$$\Leftrightarrow \text{Corr} \Big[\ln \Big(\frac{\ln X(t+\Delta t)}{\ln X(t)}\Big) , \ln \Big(\frac{\ln Y(t+\Delta t)}{\ln Y(t)}\Big)\Big] = \frac {\hat{\nu}\sigma_1 \hat{\xi} \eta_1 \Delta t} {{\sqrt{ \hat{\nu}^2 \Delta t} \sqrt{\hat{\xi}^2 \Delta t}}} $$
$$\Leftrightarrow \boxed{ \text{Corr} \Big[\ln \Big(\frac{\ln X(t+\Delta t)}{\ln X(t)}\Big) , \ln \Big(\frac{\ln Y(t+\Delta t)}{\ln Y(t)}\Big)\Big] = \sigma_1 \eta_1} $$
Given that $\sigma_1 = 1$ and $\eta_1^2 + \eta_2^2 = 1$, we can see that the choice of $\eta_1$ depends on values of $\rho_{X,Y}$ and $\sigma_1$. Assuming that $\rho_{X,Y}$ returns a postive value, we will have $\eta_1 = \rho_{X,Y}$ when $\sigma_1 = 1$, and $\eta_1 = -\rho_{X,Y}$ when $\sigma_1 = -1$.
Alternatively, when $\rho_{X,Y}$ returns a negative value, we will have $\eta_1 = -\rho_{X,Y}$ when $\sigma_1 = 1$, and $\eta_1 = \rho_{X,Y}$ when $\sigma_1 = -1$.
For each of the scenarios above, there are two choices of $\eta_2$ such that $\eta_2 = |\sqrt{1-\eta_1^2}|$.
def Stage_3b(sgm_1):
correlation = correlation_log_ret
eta_1 = correlation / sgm_1
abs_eta_2 = np.abs(np.sqrt(1 - eta_1**2))
return eta_1, abs_eta_2
choice_eta1, choice_eta2 = Stage_3b(1)
alt_choice_eta1, alt_choice_eta2 = Stage_3b(-1)
# Combination of
print("The choices of model parameters eta_1 and eta_2:")
print(f"eta_1 is {choice_eta1} and eta_2 is {choice_eta2}")
print(f"eta_1 is {-choice_eta1} and eta_2 is {-choice_eta2}")
print(f"eta_1 is {alt_choice_eta1} and eta_2 is {alt_choice_eta2}")
print(f"eta_1 is {-alt_choice_eta1} and eta_2 is {-alt_choice_eta2}")
The choices of model parameters eta_1 and eta_2: eta_1 is 0.19984846826387667 and eta_2 is 0.9798268161938528 eta_1 is -0.19984846826387667 and eta_2 is -0.9798268161938528 eta_1 is -0.19984846826387667 and eta_2 is 0.9798268161938528 eta_1 is 0.19984846826387667 and eta_2 is -0.9798268161938528
6. Stage 4: Pricing An Underperformance Option On Moderna Inc.
Attempt to determine the "fair" price of this option on 9 March 2022 using the model parameters obtained in the previous stages, where the option expires on 20 January 2023, and the "strike factor" is chosen such that on 9 March 2022, the option is at-the-money based on forward prices.
An underperformance option on Moderna Inc. has payoff at time $T$: $$ \text{Payoff} = \max(0,\, K X(T) - Y(T)) $$
The "fair" price at time $T_0$ (9 March 2022) of this option, which expires on 20 January 2023, is given by:
$$ C_{\text{moderna}}(0, X(0), Y(0)) = B(0,T)\,\mathbb{E}_{\mathbb{Q}}\!\left[(K X(T) - Y(T))^+\right] = B(0,T)\,\mathbb{E}_{\mathbb{Q}}\!\left[(K X(T) - Y(T))\,\mathbf{1}_{\{Y(T)/X(T)<K\}}\right] $$
$$ $$
Using Radon-Nikodym Derivative method for $X(T)$ to change measure from $\mathbb{Q}$ to $\tilde{\mathbb{Q}}$, we have:
$$ \frac{d\tilde{\mathbb{Q}}}{d\mathbb{Q}} = \frac{X(T)B(0,T)}{X(0)B(T,T)}\,\mathbf{C_1} = D(0,T)\, \exp\!\left\{ -\frac{1}{2}\!\int_0^T \sigma(s)^\top \sigma(s)\,ds + \int_0^T \sigma(s)^\top dW(s) \right\} \mathbf{C_1} $$
This is an exponential martingale when $\mathbf{C_1} = \frac{1}{D(0,T)}$.
Using Girsanov, we have the Brownian motion under $\tilde{\mathbb{Q}}$: $\quad$ $ d\tilde{W}(t) = dW(t) - \sigma(t)\,dt $
Under $\tilde{\mathbb{Q}}$:
$$ \ln X(T) = \ln\!\left(\frac{X(0)D(0,T)}{B(0,T)}\right) + \int_0^T \frac{1}{2} \sigma(s)^\top \sigma(s)\,ds + \int_0^T \sigma(s)^\top d\tilde{W}(s) $$
$$ \ln Y(T) = \ln\!\left(\frac{Y(0)}{B(0,T)}\right) + \int_0^T \left[\sigma(s)^\top \eta(s) - \frac{1}{2}\eta(s)^\top \eta(s)\right]ds + \int_0^T \eta(s)^\top d\tilde{W}(s) $$
Hence,
$$ \ln\!\frac{Y(T)}{X(T)} = \ln\!\left(\frac{Y(0)}{X(0)D(0,T)}\right) + \int_0^T \!\Big[\sigma(s)^\top \eta(s) - \frac{1}{2}\eta(s)^\top \eta(s) - \frac{1}{2}\sigma(s)^\top \sigma(s)\Big]\,ds + \int_0^T (\eta(s)-\sigma(s))^\top d\tilde{W}(s) $$
Thus,
$$ \ln\!\frac{Y(T)}{X(T)} \;\overset{\tilde{\mathbb{Q}}|\mathbf{F_t}}{\sim}\; \mathcal{N}\!\left( \ln\!\left(\frac{Y(0)}{X(0)D(0,T)}\right) + \int_0^T [\sigma(s)^\top \eta(s) - \frac{1}{2}\eta(s)^\top \eta(s) - \frac{1}{2}\sigma(s)^\top \sigma(s)]\,ds,\; \int_0^T \big[(\eta(s)-\sigma(s))^\top (\eta(s)-\sigma(s)\big]ds \right) $$
For $Y(T)$, we use R-N to change measure from $\mathbb{Q}$ to $\hat{\mathbb{Q}}$:
$$ \frac{d\hat{\mathbb{Q}}}{d\mathbb{Q}} = \frac{Y(T)B(0,T)}{Y(0)} = \exp\left\{ -\frac{1}{2} \int_0^T \eta(s)^T \eta(s) \, ds + \int_0^T \eta(s)^\top \, dW(s) \right\} $$
This is an exponential martingale.
Using Girsanov's theorem, we have $d\hat{W}(t) = dW(t)-\eta(t)\,dt$
$$ \ln X(T) = \ln\left( \frac{X(0)D(0,T)}{B(0,T)} \right) + \int_0^T \left[ \sigma(s)^\top \eta(t) - \frac{1}{2} \sigma(s)^\top \sigma(s) \right] ds + \int_0^T \sigma(s)^\top \, d\hat{W}(s) $$
$$ \ln Y(T) = \ln\left( \frac{Y(0)}{B(0,T)} \right) + \int_0^T \frac{1}{2} \eta(s)^\top \eta(s)\,ds + \int_0^T \eta(s)^\top\,d\hat{W}(s) $$
$$ \ln \frac{Y(T)}{X(T)} = \ln\left( \frac{Y(0)}{X(0)D(0,T)} \right) + \int_0^T \left[ \frac{1}{2} \eta(s)^\top \eta(s) - \sigma(s)^\top \eta(t) + \frac{1}{2} \sigma(s)^\top \sigma(s) \right]\,ds + \int_0^T (\eta(s) - \sigma(s))^\top \, d\hat{W}(s) $$
$$ \ln \left( \frac{Y(T)}{X(T)} \right)\;\overset{\tilde{\mathbb{Q}}|\mathbf{F_t}}{\sim}\; \mathcal{N}\left( \ln\left( \frac{Y(0)}{X(0)D(0,T)} \right) + \int_0^T \left[ \frac{1}{2} \eta(s)^\top \eta(s) - \sigma(s)^\top \eta(s) + \frac{1}{2} \sigma(s)^\top \sigma(s) \right] ds,\ \int_0^T (\eta(s) - \sigma(s))^\top (\eta(s) - \sigma(s)) \, ds \right) $$
Changing measure via the numeraire technique, we have the price of the "under-performing" Moderna option as:
$$ C_{\text{moderna}}(0, X(0), Y(0)) = D(0,T)\,K\,X(0)\,\mathbb{E}_{\tilde{\mathbb{Q}}}\!\left[\mathbf{1}_{\left\{\frac{Y(T)}{X(T)} < K\right\}}\right] - Y(0)\,\mathbb{E}_{\hat{\mathbb{Q}}}\!\left[\mathbf{1}_{\left\{\frac{Y(T)}{X(T)} <K\right\}}\right] $$
Finally,
$$\boxed{ C_{\text{moderna}}(0, X(0), Y(0)) = D(0,T)\,K X(0)\,\Phi(d_1) - Y(0)\,\Phi(d_2)} $$
$$\boxed{ d_1 = \frac{ \ln\left( \frac{K X(0) D(0,T)}{Y(0)} \right) + \frac{1}{2}\int_0^T \left[\eta(t)^\top \eta(t) -2 \sigma(s)^\top \eta(t) + \sigma(s)^\top \sigma(s) \right] ds }{ \sqrt{ \int_0^T \left( \eta(s)^\top \eta(s) - 2 \sigma(s)^\top \eta(s) + \sigma(s)^\top \sigma(s) \right) ds }}} $$
$$\boxed{ d_2 = \frac{ \ln\left( \frac{K X(0) D(0,T)}{Y(0)} \right) - \frac{1}{2}\int_0^T \left[\eta(t)^\top \eta(t) -2 \sigma(s)^\top \eta(t) + \sigma(s)^\top \sigma(s) \right] ds }{ \sqrt{ \int_0^T \left( \eta(s)^\top \eta(s) - 2 \sigma(s)^\top \eta(s) + \sigma(s)^\top \sigma(s) \right) ds }}} $$
with $ \boxed{K = \frac{Y(0)}{X(0) D(0,T)}} $ , given that the "strike factor" $K$ is chosen such that at time $T_0$ (9 March 2022), the option is at-the-money based on forward prices.
def Stage_4():
X = sp500[sp500['Date'] == today]['Price'].values[0]
Y = mrna[mrna['Date'] == today]['Price'].values[0]
maturity = pd.to_datetime("2023-01-20")
delta_t = year_transform(maturity)
K = Y / (X * D_cal(delta_t))
vol_T = (cov_cal(delta_t, type = "variance_index")
+ cov_cal(delta_t, type = "variance_moderna")
- 2 * cov_cal(delta_t, type = "covariance")
* correlation_log_ret)
d1 = (np.log((D_cal(delta_t) * K * X) / Y) + 0.5 * vol_T) / np.sqrt(vol_T)
d2 = d1 - np.sqrt(vol_T)
price = D_cal(delta_t) * K * X * norm.cdf(d1) - Y * norm.cdf(d2)
return price
print("The price of the under-performing option on Moderna Inc is: ", Stage_4())
The price of the under-performing option on Moderna Inc is: 35.75113444654894
7. Stage 5: Dynamic Hedging.
Conduct a dynamic hedging strategy using self-financing replicating portfolios for the underperformance option on Moderna Inc. stock, and evaluate the performance of this strategy over the life of the option.
We suppose that an investment bank sells this option on 9 March 2022 and conducts a dynamic hedging strategy based on the assumptions behind the price calculated in Stage 4, with daily rebalancing of the hedge.
A hypothetical scenario is considered: What is the bank’s profit or loss on 20 January 2023?
OPTION PRICING Form above we have already derived the price of the under-performing Moderna option at time $t < T$ as:
We have the under-performing Moderna option price at time $t < T$ is:
$$\boxed{ C_{\text{moderna}}(t, X(t), Y(t)) = D(t,T)\,K X(t)\,\Phi(d_1) - Y(t)\,\Phi(d_2)} $$
With
$$ \boxed{ d_1 = \frac{ \ln\left( \frac{K X(t) D(t,T)}{Y(t)} \right) + \frac{1}{2} \int_t^T \left[ \sigma(s)^\top \sigma(s) - 2 \sigma(s)^\top \eta(t) + \eta(t)^\top \eta(t) \right] ds} {\sqrt{ \int_t^T \left( \eta(s)^\top \eta(s) - 2\sigma(s)^\top \eta(s) + \sigma(s)^\top \sigma(s) \right) ds }} } $$
$$ \boxed{ d_2 = \frac{ \ln\left( \frac{K X(t) D(t,T)}{Y(t)} \right) - \frac{1}{2} \int_t^T \left[ \sigma(s)^\top \sigma(s) - 2 \sigma(s)^\top \eta(t) + \eta(t)^\top \eta(t) \right] ds}{ \sqrt{ \int_0^T \left( \eta(s)^\top \eta(s) - 2\sigma(s)^\top \eta(s) + \sigma(s)^\top \sigma(s) \right) ds } } } $$
By Ito-Lemma, we have the dynamics of the option price as:
$$ d C_{\text{moderna}} (t, X(t), Y(t)) = D(t,T)\,K\,X(t)\Phi(d_1)\,dX(t) - \Phi(d_2)\,dY(t) + \big[...\big]\,dt $$
We will need to use log-linear interpolation of zero coupon prices and dividend discount factors obtained in Stage 1, such that for $T_i \le T < T_{i+1}$, we have at time T (and analogously at time t)
$$ B(t,T) = \frac{B(0,T)}{B(0,t)} = B(0,T_i)^{\frac{T-t}{T_{i+1}-T_i}} \cdot B(0,T_{i+1})^{\frac{T-t}{T_{i+1}-T_i}} $$
$\quad$
SELF-FINANCING REPLICATING PORTFOLIO
The function of the self-financing replicating portfolio with daily rebalancing of the hedge can be defined as:
$$ V(t,X(t),Y(t)) = Q_1(t) X(t) + Q_2(t) Y(t) + Q_B(t) B(t,T) $$
with $ Q_1(t) , Q_2(t)$ and $Q_B(t)$ are quantities of SPX500 index, Morderna stock and the Zero Coupon Bond, respectively.
Where the dynamics of self-financing replicating portfolio is:
$$ dV(t,X(t),Y(t)) = Q_1(t) dX(t) + Q_2(t) dY(t) + [...]dt $$
We only consider Delta-hedging for the replicating self-financing portfolio as there is insufficient data to on the prices of S&P 500 and Moderna Inc. options beyond 9 March 2022 (our dataset contains prices observed on 9 March 2022 only). Also as we are working on an exotic option, there are no other options traded in the market, on the same underlying assets, but with different characteristics to support gamma- or vega-based hedges. Accordingly, Delta/gamma-hedging and Delta/vega-hedging are not considered.
Under Delta-hedging, the quantities of X and Y in the portfolio are as follows:
$$\boxed { Q_1(t) = D(t,T) K \Phi(d_1) } $$
$$\boxed { Q_2(t) = -\Phi(d_2) } $$
$$\boxed { Q_B(t) = \frac {C_{\text{moderna}} (t, X(t), Y(t)) - Q_1(t) X(t) - Q_2(t) Y(t)}{B(t,T)} = 0 } $$
The self-financing replicating portfolio that adjusts the hedge once per day is thus reduced to:
$$\boxed{ V(t,X(t),Y(t)) = Q_1(t) X(t) + Q_2(t) Y(t)} $$
with $ \boxed{K = \frac{Y(0)}{X(0) D(0,T)}} $ , given that the "strike factor" $K$ is chosen such that at time $T_0$ (9 March 2022), the option is at-the-money based on forward prices.
$\quad$
DYNAMIC HEDGING STRATEGY
Given the positions of S&P 500 index and Moderna Inc. stock above, we will conduct the dynamic hedging strategy daily, from 9 March 2022 to the day before the option maturity on 20 January 2023.
To do this, we will conduct daily hedging where we will buy the positions in the portfolio and sell it on the following day and then buy the portfolio again based on the new positions. The process is repeated until we are at the date just before the maturity date.
Each day's profit & loss is immediately reinvested in zero-coupond bonds. Moreover, the daily profit & loss is divided by the discount factor to get forward the value on 20 January 2023. We will then sum on the profit & loss on 20 January 2023 to get the overall profit & loss.
def Stage_5():
X = sp500[sp500['Date'] == today]['Price'].values[0]
Y = mrna[mrna['Date'] == today]['Price'].values[0]
maturity = pd.to_datetime("2023-01-20")
T = year_transform(maturity)
K = Y / (X * D_cal(T))
def delta_cal(x):
t = year_transform(x)
X_t = sp500[sp500['Date'] == x]['Price']
Y_t = mrna[mrna['Date'] == x]['Price']
vol_T = cov_cal(T, type = "variance_index") + cov_cal(T, type = "variance_moderna") - 2 * cov_cal(T, type = "covariance") * correlation_log_ret
vol_t = cov_cal(t, type = "variance_index") + cov_cal(t, type = "variance_moderna") - 2 * cov_cal(t, type = "covariance") * correlation_log_ret
vol = vol_T - vol_t
D = D_cal(T) / D_cal(t)
d1 = (np.log((D * K * X_t) / Y_t) + 0.5 * vol) / np.sqrt(vol)
d2 = d1 - np.sqrt(vol)
q1 = D * K * norm.cdf(d1)
q2 = -norm.cdf(d2)
return q1, q2
list_date = mrna[(mrna['Date'] <= maturity) & (mrna['Date'] > today)]['Date'].values
current_pos = delta_cal(today)
daily_pnl = []
T_pnl = []
for t in list_date:
X_t = sp500[sp500['Date'] == t]['Price']
Y_t = mrna[mrna['Date'] == t]['Price']
sell_price = current_pos[0] * X_t + current_pos[1] * Y_t
if t < maturity:
new_pos = delta_cal(t)
buy_price = new_pos[0] * X_t + new_pos[1] * Y_t
daily_pnl.append(sell_price - buy_price)
B = B_cal(T) / B_cal(year_transform(t))
T_pnl.append(daily_pnl[-1] / B)
current_pos = new_pos
else:
buy_price = np.maximum(K * X_t - Y_t, 0)
daily_pnl.append(sell_price - buy_price)
T_pnl.append(daily_pnl[-1])
print('Profit and Loss from Hedging is ', np.sum(T_pnl))
pnl_cum_sum = np.cumsum(T_pnl)
plt.figure(figsize = (12, 5))
plt.plot(list_date, daily_pnl, label = "Daily Profit & Loss")
plt.plot(list_date, pnl_cum_sum, label = "Accumulated Profit & Loss over time")
plt.xlabel('Date')
plt.ylabel('Profit and Loss')
plt.legend()
plt.title('Profit and Loss from Hedging')
plt.show()
Stage_5()
Profit and Loss from Hedging is 4.472563229909699
8. Stage 6: Plotting The "Fair" Price.
Holding all other conditions fixed, we will plot the “fair” price of this option on 9 March2022 as a function of $\eta_1$, for $\eta_1 \in [-1, 1]$, as obtained in Stage 3.
Stage_6_corr_values = np.linspace(-1, 1, 200)
Stage_6_option_prices = []
def Stage_6(corr):
X = sp500[sp500['Date'] == today]['Price'].values[0]
Y = mrna[mrna['Date'] == today]['Price'].values[0]
maturity = pd.to_datetime("2023-01-20")
delta_t = year_transform(maturity)
K = Y / (X * D_cal(delta_t))
vol_T = (cov_cal(delta_t, type = "variance_index")
+ cov_cal(delta_t, type = "variance_moderna")
- 2 * cov_cal(delta_t, type = "covariance")
* corr)
d1 = (np.log((D_cal(delta_t) * K * X) / Y) + 0.5 * vol_T) / np.sqrt(vol_T)
d2 = d1 - np.sqrt(vol_T)
price = D_cal(delta_t) * K * X * norm.cdf(d1) - Y * norm.cdf(d2)
return price
for corr in Stage_6_corr_values:
Stage_6_option_prices.append(Stage_6(corr))
plt.figure(figsize=(12, 5))
plt.plot(Stage_6_corr_values, Stage_6_option_prices)
plt.xlabel('Correlation between MRNA and SP500')
plt.ylabel('Option Price')
plt.title('Option Price vs Correlation')
plt.grid(True)
plt.show()
Project Summary
The results indicate that the option value decreases as the correlation increases. When MRNA is negatively correlated with the S&P 500, the option benefits from diversification because the long S&P 500 position and the short Moderna position offset each other more effectively. As a result, the payoff is higher and the option is worth more (approximately 50 when $\rho_{X,Y} = -1$).
In contrast, a positive correlation results in a lower option price. When MRNA and the S&P 500 move in the same direction (as $\rho_{X,Y}$ tends to 1), the diversification benefit is reduced. The payoff is then lower and the option is worth less (approximately 20 when $\rho_{X,Y} = 1$).