2
$\begingroup$

I have this exercise. $\\\\$ Look for realistic values ​​of the parameters and calculate the price of a European Call with maturity $T = 0.5$ and $S_0 = 1$ for the strike values $​​K = 0.5,0.6, ......, 1.5$. Display the volatility smile.

I poste my code in python, but I don't understand why the volatility smiley is not displayed.

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from py_vollib_vectorized import vectorized_implied_volatility as implied_vol


# Parameters
# simulation dependent
S0 = 1           # asset price
T = 0.5              # time in years
r = 0.00003            # risk-free rate
N = 125             # number of time steps in simulation
M = 10000         # number of simulations

# Heston dependent parameters
kappa = 1.1              # rate of mean reversion of variance under risk-neutral dynamics
theta = 0.02        # long-term mean of variance under risk-neutral dynamics
v0 = 0.025          # initial variance under risk-neutral dynamics
rho = -0.7              # correlation between returns and variances under risk-neutral dynamics
sigma = 0.1          # volatility of volatility


def heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M):
    """
    Inputs:
     - S0, v0: initial parameters for asset and variance
     - rho   : correlation between asset returns and variance
     - kappa : rate of mean reversion in variance process
     - theta : long-term mean of variance process
     - sigma : vol of vol / volatility of variance process
     - T     : time of simulation
     - N     : number of time steps
     - M     : number of scenarios / simulations
    
    Outputs:
    - asset prices over time (numpy array)
    - variance over time (numpy array)
    """
    # initialise other parameters
    dt = T/N
    mu = np.array([0,0])
    cov = np.array([[1,rho],
                    [rho,1]])

    # arrays for storing prices and variances
    S = np.full(shape=(N+1,M), fill_value=S0)
    v = np.full(shape=(N+1,M), fill_value=v0)

    # sampling correlated brownian motions under risk-neutral measure
    Z = np.random.multivariate_normal(mu, cov, (N,M))

    for i in range(1,N+1):
        S[i] = S[i-1] * np.exp( (r - 0.5*v[i-1])*dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0] )
        v[i] = np.maximum(v[i-1] + kappa*(theta-v[i-1])*dt + sigma*np.sqrt(v[i-1]*dt)*Z[i-1,:,1],0)
    
    return S, v


rho_p = 0.98
rho_n = -0.98

S_p,v_p = heston_model_sim(S0, v0, rho_p, kappa, theta, sigma,T, N, M)
S_n,v_n = heston_model_sim(S0, v0, rho_n, kappa, theta, sigma,T, N, M)

fig, (ax1, ax2)  = plt.subplots(1, 2, figsize=(12,5))
time = np.linspace(0,T,N+1)
ax1.plot(time,S_p)
ax1.set_title('Heston Model Asset Prices')
ax1.set_xlabel('Time')
ax1.set_ylabel('Asset Prices')

ax2.plot(time,v_p)
ax2.set_title('Volatility Monte Carlo Simulation')
ax2.set_xlabel('Time')
ax2.set_ylabel('Variance')

plt.show()
# simulate gbm process at time T
gbm = S0*np.exp( (r - theta**2/2)*T + np.sqrt(theta)*np.sqrt(T)*np.random.normal(0,1,M) )

fig, ax = plt.subplots()

ax = sns.kdeplot(S_p[-1], label=r"$\rho= 0.98$", ax=ax)
ax = sns.kdeplot(S_n[-1], label=r"$\rho= -0.98$", ax=ax)
ax = sns.kdeplot(gbm, label="GBM", ax=ax)

plt.title(r'Asset Price Density under Heston Model')
plt.xlim([0.5, 1.5])
plt.xlabel('$S_T$')
plt.ylabel('Density')
plt.legend()
plt.show()


rho = -0.7
S,v = heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M)

# Set strikes and complete MC option price for different strikes
K = np.arange(0.5,1.5)



calls = np.array([np.exp(-r*T)*np.mean(np.maximum(S-k,0)) for k in K])


call_ivs = implied_vol(calls, S0, K, T, r, flag='c', q=0, return_as='numpy', on_error='ignore')

plt.plot(K, call_ivs, label=r'Call')


plt.ylabel('Implied Volatility')
plt.xlabel('Strike')

plt.title('Implied Volatility Smile from Heston Model')
plt.legend()
plt.show()

Thanks for the support.

$\endgroup$

1 Answer 1

0
$\begingroup$

Spotted issues:

  • The main issue: You define S0 = 1 which in Python is of type Int. Your stock price trajectories are built upon S0 and thus S will be of type Int as-well (more specifically, Int32). This implies that any decimals in your estimates will be truncated. To fix this simply define S0 = 1.0, which is now of type float.

  • Minor issue: When defining your set of strikes $K$ you need to define the step-size (automatically set to 1) in the numpy arange function, K = np.arange(0.1,1.6,0.1). This will give you the specified set of strikes.

  • While testing the code, I have spotted that implied_vol function might run into some optimization issues, that generates a weird smile. I got a decent smile (shown below) after a couple of re-runs of the function.

There might be some other small issues, but fixing the above problems gives me a working solution.

Verification:

I took the liberty of fixing the above issues and ran the code. However, I used the original Heston parameters which are found in the Youtube Tutorial that you clearly followed to produce the simulations, density and smile. This results in the following graphs:

plot123

And the following smile:

plot4

I will leave interpretation of the graphs up to you. I hope this helps.


Appendix: Python Code

# Parameters
# simulation dependent
S0 = 1.0          # asset price
T = 0.5              # time in years
r = 0.00003            # risk-free rate
#r=0.02
N = 252             # number of time steps in simulation
M = 10000         # number of simulations

# Heston dependent parameters
#kappa = 1.1              # rate of mean reversion of variance under risk-neutral dynamics
#theta = 0.02        # long-term mean of variance under risk-neutral dynamics
#v0 = 0.025          # initial variance under risk-neutral dynamics
#rho = -0.7              # correlation between returns and variances under risk-neutral dynamics
#sigma = 0.1          # volatility of volatility

kappa = 3
theta = 0.20**2
v0=0.25**2
#rho = 0.7
sigma = 0.6 



def heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M):
    """
    Inputs:
     - S0, v0: initial parameters for asset and variance
     - rho   : correlation between asset returns and variance
     - kappa : rate of mean reversion in variance process
     - theta : long-term mean of variance process
     - sigma : vol of vol / volatility of variance process
     - T     : time of simulation
     - N     : number of time steps
     - M     : number of scenarios / simulations
    
    Outputs:
    - asset prices over time (numpy array)
    - variance over time (numpy array)
    """
    # initialise other parameters
    dt = T/N
    mu = np.array([0,0])
    cov = np.array([[1,rho],
                    [rho,1]])

    # arrays for storing prices and variances
    S = np.full(shape=(N+1,M), fill_value=S0)
    v = np.full(shape=(N+1,M), fill_value=v0)

    # sampling correlated brownian motions under risk-neutral measure
    Z = np.random.multivariate_normal(mu, cov, (N,M))

    for i in range(1,N+1):
        S[i] = S[i-1] *  np.exp((r - 0.5*v[i-1])*dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0])
        v[i] = np.maximum(v[i-1] + kappa*(theta-v[i-1])*dt + sigma*np.sqrt(v[i-1]*dt)*Z[i-1,:,1],0)
    
    return S, v


rho_p = 0.98
rho_n = -0.98

S_p,v_p = heston_model_sim(S0, v0, rho_p, kappa, theta, sigma,T, N, M)
S_n,v_n = heston_model_sim(S0, v0, rho_n, kappa, theta, sigma,T, N, M)

fig, (ax1, ax2)  = plt.subplots(1, 2, figsize=(12,5))
time = np.linspace(0,T,N+1)
ax1.plot(time,S_p)
ax1.set_title('Heston Model Asset Prices')
ax1.set_xlabel('Time')
ax1.set_ylabel('Asset Prices')

ax2.plot(time,v_p)
ax2.set_title('Volatility Monte Carlo Simulation')
ax2.set_xlabel('Time')
ax2.set_ylabel('Variance')

plt.show()
# simulate gbm process at time T
gbm = S0*np.exp( (r - theta**2/2)*T + np.sqrt(theta)*np.sqrt(T)*np.random.normal(0,1,M) )

fig, ax = plt.subplots()

ax = sns.kdeplot(S_p[-1], label=r"$\rho= 0.98$", ax=ax)
ax = sns.kdeplot(S_n[-1], label=r"$\rho= -0.98$", ax=ax)
ax = sns.kdeplot(gbm, label="GBM", ax=ax)

plt.title(r'Asset Price Density under Heston Model')
plt.xlim([0.5, 1.5])
plt.xlabel('$S_T$')
plt.ylabel('Density')
plt.legend()
plt.show()

rho = -0.7
S,v = heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M)

# Set strikes and complete MC option price for different strikes

K = np.arange(0.1,1.6,0.1)

calls = np.array([np.exp(-r*T)*np.mean(np.maximum(S-k,0)) for k in K])


call_ivs = implied_vol(calls, S0, K, T, r, flag='c', q=0, return_as='numpy', on_error='ignore')

plt.plot(K, call_ivs, label=r'Call')


plt.ylabel('Implied Volatility')
plt.xlabel('Strike')

plt.title('Implied Volatility Smile from Heston Model')
plt.legend()
plt.show()
$\endgroup$
2
  • 1
    $\begingroup$ Thanks very very much Pleb, I resolve my problem! Can I ask you to take a look at the Euler discretization code too? i just didn't understand how to set it up in python. thank you so much! $\endgroup$
    – GloBag578
    May 31 at 19:55
  • $\begingroup$ quant.stackexchange.com/questions/71024/… $\endgroup$
    – GloBag578
    May 31 at 20:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.