# Pricing a Vanilla European Call Option


```python
#Import Libraries
import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import norm
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.ticker as ticker
from matplotlib.axis import Axis
from scipy.stats import ncx2
```

<center><h4>Abstract</h4></center>
In this paper, we price a vanilla European call option under the Heston model and then simulate the monthly share price over a year using the Constant Elasticity of Variance (CEV) model, with the assumption of constant volatility each month. Monte Carlo simulations with varying sample sizes are run and the results are plotted against the closed form value for com- parison.
<br>
<br>

#### Keyboard:
European Options, Heston model, Constant Elasticity of Variance (CEV) model and Stochastic Volatility.

##  Introduction

In this paper we go beyond the constant volatility assumption in the Black-Scholes model. Most of code implemented in this submission is derived from Module 5 [7] of the course.
We initialize most variables as given by the question.

- Maturity = 1 year
- Strike price as 100\$
- $S_0=100$ share price is \$100
- $r=8%$  (compound interest rate risk free)
- $\sigma=30%$ as volatility


```python
#Share Specific Information
r = 0.08 #risk-free continuously compounded interest rate 
S0 = 100 #current share price 
sigma = 0.3

#Call Specific Information
K = 100
T = 1
k_log = np.log(K)
```

## Part 1: Fourier pricing


```python

```

Black-Scholes option pricing models assumed volatility of the underlying assets to be constant or a predetermined function of time, we will now implement a model which features instantaneous variance of asset price using volatility that evolves stochastically in time.  Although there are several models incorporating stochastic volatility including introduced by Hull and White [1], Stein and Stein [2] and  Heston [3] [4],  we will implement Heston’s constant interest rate model.

With the assumption that the underlying share follows the Heston model dynamics, the additional parameters required are specified as follows:


```python
#Share Specific Information
v0 = 0.06
kappa = 9
theta = 0.06
sigma = 0.3
rho = -0.4

#Approximation Information
t_max = 30
N = 100
```

The Characteristic function is implemented using a function presented by Albrecher et al [5].  The function is written as:
$$\phi_{S_T} = exp(C(\tau;u)+D(\tau;u)v_t + iu \log(S_t) )$$

Where,

$$C(\tau;u) = ri\tau u + \theta \kappa [\tau x\_ - \frac{1}{a}\log(\frac{1-ge^{d\tau}}{1-g})],$$

$$D(\tau;u) = (\frac{1-e^{d\tau}}{1-ge^{d\tau}})x\_,$$


\begin{equation}
  \begin{split}
    \tau = T -t \\
    g = \frac{x\_}{x_+}\\
    x_\pm = \frac{b\pm d}{2a}
  \end{split}
\quad\quad
  \begin{split}
    d = \sqrt{b^2 - 4ac}\\
    c = - \frac{u^2 + Ui}{2} \\
    b = \kappa - \rho \sigma iu \\
    a = \frac {\sigma^2}{2} \\
  \end{split}
\end{equation}


```python
#Characteristic function code
a = sigma**2/2

def c(u):
    return -(u**2+1j*u)/2

def d(u):
    return np.sqrt(b(u)**2-4*a*c(u))

def b(u):
    return kappa - rho*sigma*1j*u

def g(u):
    return xminus(u)/xplus(u)

def C(u):
    val1 = T*xminus(u)-np.log((1-g(u)*np.exp(-T*d(u)))/(1-g(u)))/a
    return r*T*1j*u + theta*kappa*val1

def D(u):
    val1 = 1-np.exp(-T*d(u))
    val2 = 1-g(u)*np.exp(-T*d(u))
    return (val1/val2)*xminus(u)

def log_char(u):
    return np.exp(C(u) + D(u)*v0 + 1j*u*np.log(S0))

def xminus(u):
    return (b(u)-d(u))/(2*a)

def xplus(u):
    return (b(u)+d(u))/(2*a)

def adj_char(u):
    return log_char(u-1j)/log_char(-1j)
```

Now we vectorize the code, calculate an estimate for integrals and calculate the Fourier estimate of our call price.


```python
delta_t = t_max/N
from_1_to_N = np.linspace(1,N,N)
t_n = (from_1_to_N-1/2)*delta_t

first_integral = sum((((np.exp(-1j*t_n*k_log)*adj_char(t_n)).imag)/t_n)*delta_t)
second_integral = sum((((np.exp(-1j*t_n*k_log)*log_char(t_n)).imag)/t_n)*delta_t)

fourier_call_val = S0*(1/2 + first_integral/np.pi)-np.exp(-r*T)*K*(1/2 + second_integral/np.pi)
fourier_call_val
```




    13.734895692109049



To see the effectiveness of the pricing option under Heston dynamics we will also price the call option under Black-Scholes assumption.


```python
d_1_stock = (np.log(S0/K)+(r + sigma**2/2)*(T))/(sigma*np.sqrt(T))
d_2_stock = d_1_stock - sigma*np.sqrt(T)

analytic_callprice = S0*norm.cdf(d_1_stock)-K*np.exp(-r*(T))*norm.cdf(d_2_stock)
analytic_callprice
```




    15.711312547892973



## Part 2: CEV Simulations

We use the following formula to simulate share price walk :
$$S_{t_{i+1}} = S_{t_i}e^{(r - \frac{\sigma^2 (t_i,t_{i+1})}{2})(t_{i+1}-t_i)+\sigma(t_i, t_{i+1})\sqrt{t_{i+1}-t_i}Z}$$
where $S_{ti}$ is the share price at time $t_i$, $\sigma(t_i, t_{i+1})$ is the volatility for the period $[t_i,t_{i+1}]$, r is the risk-free interest rate, and $Z \sim N(0,1)$

For estimations we will use sample size from 1000 upto 50,000 for monthly simulation over 1 year. More over
we define $\sigma(t_i,t_{i+1}) = \sigma(S_{ti})^{\gamma-1}$, 
with $\sigma = 0.3$ and $\gamma = 0.75$.


We are attempting to run simulations using the CEV model. However, while the CEV model assumes that volatility is a continuous function of time and share price, we are making a simplifying assumption that volatility is constant over each simulation period.

Cox [6] developed the constant elasticity of variance (CEV) option pricing model, it attempts to capture stochastic volatility and is given by:
<br>
<br>
<center>$dS_t = \mu S_t dt + \sigma {S_t}^{\gamma} dW_t$</center>
<br>
If $\gamma = 1$ this model return the same value as Black-Scholes model, however if the value of $\gamma < 1$ we experience an effect called leverage effect where the volatility increases as the price decreases over subsequent time periods.


```python
#from question
sigma_const = 0.30
gamma = 0.75

#assumption
r = 0.08
S0 = 100
T = 1
sample_sizes = range(1000, 50001, 1000)

timesteps = 12
```

First, we define our helper functions. The next_share_price function is used to calculate the evolution of the share price at $t+1$, from the share price at t. We generate the random variable $Z$ from within this function. The effective sigma is also written as a function of the share price at $t$. 

Just for exploration purposes, we also added a varying_vol flag, to allow us to switch between a constant volatility and varying volatility. We use the initial stock price instead of the previous price when assuming a contant volatility.

The other function is the generate_share_price_path function. We first create an empty numpy array of shape (sample_size x timesteps+1). Note the addition of one element to the timestep, as the first element is equal to the initial stock price. We then iterate through each path, and each timestep, applying the next_share_price function against the previous share price. We then convert the result into a pandas DataFrame and return the results.


```python
def next_share_price(prev_price, r, dT, sigma_const, gamma, sample_size, varying_vol = True):
    Z = stats.norm.rvs(size=sample_size)
    if varying_vol:
        sigma = sigma_const*(prev_price)**(gamma-1)
    else:
        sigma = sigma_const*(S0)**(gamma-1)
        
    return prev_price*np.exp((r-(sigma**2)/2)*(dT)+(sigma)*(np.sqrt(dT))*Z)

def generate_share_price_path(S0, r, T, sigma_const, gamma, sample_size, timesteps, varying_vol = True):
    df = pd.DataFrame([S0]*sample_size)
    for t in range(1, timesteps+1):
        df[t] = next_share_price(df[t-1], r, 1/timesteps, sigma_const, gamma, sample_size, varying_vol)
    return df.T
```


```python
T = 10
sample_size = 20

share_price_path_cev = generate_share_price_path(S0, r, T, sigma_const, gamma, sample_size, timesteps)
share_price_path_black_scholes = \
  generate_share_price_path(S0, r, T, sigma_const, 1.0, sample_size, timesteps, varying_vol=False)

red_patch = mpatches.Patch(color='red', label='CEV')
blue_patch = mpatches.Patch(color='blue', label='BSM')

fig, ax = plt.subplots(figsize=(12,4))
ax.plot(share_price_path_cev, color='red')
ax.plot(share_price_path_black_scholes, color='blue')
plt.xlabel("Timestep")
plt.ylabel("Share price")
plt.legend(handles=[red_patch, blue_patch], loc='upper left')
plt.title("CEV vs BSM Stock Price Paths")
plt.show()
```

Next, we create a python dictionary called share_price_paths, to hold our results for part 2. The key of this dictionary would be the number of sample price paths, while the values would be the dataframes containing the price paths. We also track the rough processing time required at each step, by printing the time when each iteration completes.


```python
import time

T = 1
sample_sizes = range(1000, 50001, 1000)

share_price_paths = {}

print("Start generating share price path")
start = time.time()
for sample_size in sample_sizes:
    share_val = generate_share_price_path(S0, r, T, sigma_const, gamma, sample_size, timesteps, varying_vol=False) 
     
    share_price_paths[sample_size] = share_val
    #print("Updated for sample size {} at {}".format(sample_size, datetime.datetime.now().strftime('%H:%M')))
end = time.time()
print(f"Generating all samples paths takes {(end - start):.2f}s")
```

    Start generating share price path
    Generating all samples paths takes 1.59s


To display our output, we show the first 10 price paths generated by our iteration with 1000 samples:


```python
share_price_paths[1000].iloc[:, 0:10]
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
      <td>100.000000</td>
    </tr>
    <tr>
      <th>1</th>
      <td>100.857562</td>
      <td>100.106398</td>
      <td>103.296220</td>
      <td>103.237244</td>
      <td>102.477151</td>
      <td>100.141983</td>
      <td>99.860355</td>
      <td>99.265029</td>
      <td>100.465829</td>
      <td>96.664721</td>
    </tr>
    <tr>
      <th>2</th>
      <td>102.299627</td>
      <td>102.170988</td>
      <td>103.124669</td>
      <td>106.038341</td>
      <td>102.533543</td>
      <td>101.845595</td>
      <td>100.491371</td>
      <td>99.395467</td>
      <td>105.083317</td>
      <td>92.029417</td>
    </tr>
    <tr>
      <th>3</th>
      <td>101.867161</td>
      <td>93.865529</td>
      <td>101.634586</td>
      <td>104.976242</td>
      <td>108.473032</td>
      <td>103.778385</td>
      <td>101.742397</td>
      <td>94.643315</td>
      <td>102.962963</td>
      <td>94.429574</td>
    </tr>
    <tr>
      <th>4</th>
      <td>102.813675</td>
      <td>91.982826</td>
      <td>105.212492</td>
      <td>102.155602</td>
      <td>106.666274</td>
      <td>101.342648</td>
      <td>102.277561</td>
      <td>97.341252</td>
      <td>100.677285</td>
      <td>100.997206</td>
    </tr>
    <tr>
      <th>5</th>
      <td>100.933753</td>
      <td>96.627653</td>
      <td>106.951672</td>
      <td>104.840510</td>
      <td>107.463153</td>
      <td>101.973944</td>
      <td>98.958936</td>
      <td>98.323329</td>
      <td>101.680525</td>
      <td>102.179669</td>
    </tr>
    <tr>
      <th>6</th>
      <td>101.957430</td>
      <td>95.208273</td>
      <td>104.464336</td>
      <td>106.214328</td>
      <td>107.054801</td>
      <td>106.243952</td>
      <td>100.637030</td>
      <td>103.183938</td>
      <td>104.453154</td>
      <td>107.263225</td>
    </tr>
    <tr>
      <th>7</th>
      <td>104.571916</td>
      <td>92.148076</td>
      <td>102.967291</td>
      <td>106.638492</td>
      <td>104.328094</td>
      <td>110.008954</td>
      <td>97.368919</td>
      <td>107.018006</td>
      <td>102.202276</td>
      <td>102.605055</td>
    </tr>
    <tr>
      <th>8</th>
      <td>106.576893</td>
      <td>90.127079</td>
      <td>104.435163</td>
      <td>106.043822</td>
      <td>103.826505</td>
      <td>111.997720</td>
      <td>100.257464</td>
      <td>106.680848</td>
      <td>98.383468</td>
      <td>103.449951</td>
    </tr>
    <tr>
      <th>9</th>
      <td>107.267472</td>
      <td>95.636057</td>
      <td>108.110705</td>
      <td>108.198575</td>
      <td>104.420906</td>
      <td>117.704747</td>
      <td>99.341943</td>
      <td>108.905238</td>
      <td>102.333938</td>
      <td>103.678585</td>
    </tr>
    <tr>
      <th>10</th>
      <td>115.118830</td>
      <td>97.666874</td>
      <td>111.960859</td>
      <td>110.175120</td>
      <td>103.001066</td>
      <td>114.141261</td>
      <td>105.439317</td>
      <td>110.180501</td>
      <td>102.920079</td>
      <td>104.825189</td>
    </tr>
    <tr>
      <th>11</th>
      <td>117.368677</td>
      <td>97.603712</td>
      <td>115.867144</td>
      <td>106.104115</td>
      <td>105.544988</td>
      <td>112.539562</td>
      <td>107.405076</td>
      <td>107.156807</td>
      <td>100.694101</td>
      <td>102.665990</td>
    </tr>
    <tr>
      <th>12</th>
      <td>121.017598</td>
      <td>97.368054</td>
      <td>113.049343</td>
      <td>103.781586</td>
      <td>100.745486</td>
      <td>112.449897</td>
      <td>109.702714</td>
      <td>109.048143</td>
      <td>105.944799</td>
      <td>100.248490</td>
    </tr>
  </tbody>
</table>
</div>



Next, we plot the price paths for the iteration with 1000 samples, 


```python
fig, ax = plt.subplots(figsize=(12,4))

ax.plot(share_price_paths[1000])
plt.xlabel('Timestep')
plt.ylabel('Share price')
plt.title('Sample of 1000 Share Price Paths')
plt.show()
```

## Part 3: Monte-Carlo pricing

Using Monte Carlo, we calculate the price of the vanilla call option as follows:


```python
price_estimate = []
price_std = []


for size in sample_sizes:
    S_Ts = share_price_paths[size].iloc[12, :]
    payoff = np.maximum(S_Ts - K, 0)
    discounted_price = np.exp(-r*T)*payoff
    price_estimate.append(discounted_price.mean())
    price_std.append(discounted_price.std()/np.sqrt(size))
```


```python
print("The price estimated by Monte Carlo when using sample size of 50,000 is : {:.3f}".format(price_estimate[-1])) 
```

    The price estimated by Monte Carlo when using sample size of 50,000 is : 8.713


### Comparison with calculation of CEV model using noncentral chi-squared distribution 


```python
S0 = 100
sigma = 0.3
gamma = 0.75
r = 0.08
T = 1
```


```python
z = 2 + 1/(1-gamma)
def C(t,K):
    kappa = 2*r/(sigma**2*(1-gamma)*(np.exp(2*r*(1-gamma)*t)-1))
    x = kappa*S0**(2*(1-gamma))*np.exp(2*r*(1-gamma)*t)
    y = kappa*K**(2*(1-gamma))
    return S0*(1-ncx2.cdf(y,z,x))-K*np.exp(-r*t)*ncx2.cdf(x,z-2,y)
```


```python
cev_call_price = C(T, 100)
```


```python
print("The price calculated via CEV model using noncentral chi-squared distribution is : {:.3f}".format(cev_call_price))
```

    The price calculated via CEV model using noncentral chi-squared distribution is : 8.702


## Part 4:  Monte Carlo estimates

We plot the prices of our Monte Carlo estimates against the CEV noncentral chi-squared distribution prices calculated from Part 1 above:


```python
fig, ax = plt.subplots(figsize=(12,4))
ticks_x = ticker.FuncFormatter(lambda x, pos: '{0:g}'.format(x * 1000))

ax.plot([cev_call_price]*50, label='CEV Call Price')
ax.plot(price_estimate, '.', label='Monte Carlo Estimates')
ax.plot(cev_call_price + 3*np.array(price_std), 'r')
ax.plot(cev_call_price - 3*np.array(price_std), 'r')
ax.xaxis.set_major_formatter(ticks_x)
plt.legend()
plt.show()
```

We also plot the prices of our Monte Carlo estimates against the Fourier Call prices calculated from Part 1 above:


```python
fig, ax = plt.subplots(figsize=(12,4))

ax.plot([fourier_call_val]*50, label='Fourier Call Price')
ax.plot(price_estimate, '.', label='Monte Carlo Estimates')
ax.plot(fourier_call_val + 3*np.array(price_std), 'r')
ax.plot(fourier_call_val - 3*np.array(price_std), 'r')
ax.xaxis.set_major_formatter(ticks_x)

plt.legend()
plt.show()
```

We notice that there is a discrepency. Upon further investigation, for the Heston model, ```v_0 = 0.06``` ```(v_0 = stock volatility ^ 0.5)```, so ```stock volatility = .06**0.5= 0.2449```. Note that sigma under Heston model refers to volatility of stock volatility

Under our stock price Monte Carlo calculation, the default ```stock volatility is 0.0948``` (given by  ```0.3*(100)*-.25``` ). Thus its much less than Heston. To have the same initial stock volatility, we find a new value for $\sigma$ by equating $\sigma(S_{ti})^{\gamma-1} = \sqrt{0.06}$, giving us $\sigma = 0.775$. 

We then find that our newly calculated Monte Carlo calculated call prices are aligned with the Fourier prices.


```python
sigma_const = 0.775

T = 1
sample_sizes = range(1000, 50001, 1000)

share_price_paths = {}

print("Start generating share price path")
start = time.time()
for sample_size in sample_sizes:
    share_val = generate_share_price_path(S0, r, T, sigma_const, gamma, sample_size, timesteps, varying_vol=False) 
     
    share_price_paths[sample_size] = share_val
    #print("Updated for sample size {} at {}".format(sample_size, datetime.datetime.now().strftime('%H:%M')))
end = time.time()
print(f"Generating all samples paths takes {(end - start):.2f}s")
```

    Start generating share price path
    Generating all samples paths takes 1.63s



```python
price_estimate = []
price_std = []


for size in sample_sizes:
    S_Ts = share_price_paths[size].iloc[12, :]
    payoff = np.maximum(S_Ts - K, 0)
    discounted_price = np.exp(-r*T)*payoff
    price_estimate.append(discounted_price.mean())
    price_std.append(discounted_price.std()/np.sqrt(size))
```


```python
fig, ax = plt.subplots(figsize=(12,4))

ax.plot([fourier_call_val]*50, label='Fourier Call Price')
ax.plot(price_estimate, '.', label='Monte Carlo Estimates')
ax.plot(fourier_call_val + 3*np.array(price_std), 'r')
ax.plot(fourier_call_val - 3*np.array(price_std), 'r')
ax.xaxis.set_major_formatter(ticks_x)
plt.legend()
plt.show()
```

#### Conclusion
A vanilla European call option was priced for a fluctuating volatility condition, using the Heston model. The call price was evaluated as USD 13.73 which was verified to be similar to - but lower than - a constant volatility estimate using the Black-Scholes model at USD 15.71.


Initially the values calculated for the underlying share price did not agree as expected, with the estimates well outside the expected error range. This was because the volatility term was not equivalent in the two different calculation methods. After resolving this, the CEV model calcu- lated a price of USD 8.70 for the underlying share, which agreed closely with the Monte Carlo estimate of USD 8.69. This represents an absolute error of 1c and a relative discrepancy of approximately 0.1


The agreement between the results from various methods, although expected, provides confidence to choose the most suitable method for a situation with the knowledge that the results are accurate within a small margin of error.

#### References

[1] Hull, J. and White, A. (1987). The pricing of options on assets with stochastic volatilities. The journal of finance, 42(2):281–300.

[2] Stein, E. M. and Stein, J. C. (1991). Stock price distributions with stochastic volatility: an analytic approach. Review of financial Studies, 4(4):727–752.

[3] Heston, S. L. (1993). A closed-form solution for options with stochastic volatility with applications to bond and currency options. Review of financial studies, 6(2):327–343. 

[4] Heston, S. L. (1997). A simple new formula for options with stochastic volatility. 

[5] Albrecher, H., Mayer, P., Schoutens, W. and Tistaert, J. (2007). “The Little Heston Trap”, Wilmott (1): 83–92.




3d50fb0b6fb4f2a0fe58ea17a4876f67b358d8dd
ipynb
Jupyter Notebook
Pricing a Vanilla European Call Option.ipynb
Raanaji/ComputationalProjects
8c90747cfde55ccb664579628c26dc62fcfa4c70
MIT
Pricing a Vanilla European Call Option.ipynb
Raanaji/ComputationalProjects
8c90747cfde55ccb664579628c26dc62fcfa4c70
MIT
Pricing a Vanilla European Call Option.ipynb
Raanaji/ComputationalProjects
8c90747cfde55ccb664579628c26dc62fcfa4c70
MIT
Qwen/Qwen-72B
1. YES
2. YES