# Financial Engineering and Risk Management tutorials: Option Pricing

This tutorial is a guide for pricing European put option using either numerical integration or FFT. The second approach is useful since one only needs to specify the model’s characteristic function, instead of the density which may note be readily available.

C’est un support des travaux dirigés pour le cours de M2 “Financial Engineering and Risk Management” que j’ai enseigné aux étudiants de l’Université Paris-Sud en 2017 et 2018 (devenue désormais Paris-Saclay).

Part one

Consider a speculator who believes that the stock of Amazon will increase in value in the next 6 months. The current stock price is $2000. Suppose that the speculator buys 100 stocks today. Consider the following scenarios for the portfolio’s P&L (profit and loss): Scenario Today price Future price (6m) P&L 1$2,000 $1,900 2$2,000 $1,950 3$2,000 $2,000 4$2,000 $2,050 5$2,000 $2,100 Q1. What is the P&L under Scenario 2? (select with mouse to reveal the answer) -5000 Q2. What is the P&L under Scenario 5? (select with mouse to reveal the answer) 10000 Instead of buying stocks, suppose that the speculator decides to buy 100 6-month call options on the Amazon stock with strike K=$2025. Assume that the option price is $10. Consider the following scenarios for the portfolio’s P&L (profit and loss): Scenario Today price Future price (6m) Option payoff P&L 1$2,000 $1,900 2$2,000 $1,950 3$2,000 $2,000 4$2,000 $2,050 5$2,000 $2,100 Q3. What is the P&L under Scenario 2? Answer: -1000 Q4. What is the P&L under Scenario 5? Answer: 6500 Suppose that your investment portfolio consists of 100 stocks of Amazon and assume that you do not buy any options to hedge your position. Consider the following scenarios for the portfolio’s value in 6 months: Scenario Today price Future price (6m) Future portfolio value (6m) 1$2,000 $1,600 2$2,000 $1,800 3$2,000 $2,000 4$2,000 $2,200 5$2,000 $2,400 Q5. What will be your future portfolio value under Scenario 1? Answer: 160000 Q6. What will be your future portfolio value under Scenario 4? Answer: 220000 Instead, suppose that you decide to buy 100 6-month put options on the Amazon stock at strike K=1900 (you still own the stocks). Assume that the put price is$12.

Consider the following scenarios for the portfolio’s value in 6 months:

Scenario Today price Future price (6m) Option payoff Future portfolio value (6m)
1 $2,000$1,600
2 $2,000$1,800
3 $2,000$2,000
4 $2,000$2,200
5 $2,000$2,400

Q7. What will be your future portfolio value under Scenario 1? Note: please take the cost of options into account when calculating portfolio value. Answer: 188800

Q8. What will be your future portfolio value under Scenario 4? Note: take the cost of options into account when calculating portfolio value. Answer: 218800

To answer the following questions, use a Jupyter notebook.

Assume that the stock price follows a lognormal distribution. Using the numerical integration method given during the lecture, estimate the price of a European put with the following parameters: spot price $$S_0​=100$$, strike $$K=90$$, risk-free rate $$r=4%$$, dividend rate $$q=2%$$, volatility $$σ=25%$$, and maturity $$T=1$$ year.

Do the numerical integration with the number of grid points equal to $$N=2n, n=1,2,...,15$$. For each $$N$$ your code should report $$η$$ and $$P_0$$ ​the price of the put:

$$N$$ $$η$$ $$P_0$$​
21
22
215

Q9. What is the value of $$η$$ when $$n=3$$? Answer: 11.25

Q10 What is the value of $$η$$ when $$n=10$$? Answer: 0.088

Q11. What is the value of $$P_0​$$ when $$n=4$$? Answer: 4.482

Q12. What is the value of $$P_0$$​ when $$n=15$$? Answer: 4.524

Prepare the code to price put options using FFT for different models.

For the following questions, assume that the spot price is $$S_0​=100$$, the strike is $$K=80$$, the risk-free rate is $$r=5%$$, the dividend rate is $$1%$$, and the maturity is $$T=1$$.

First, we will price put options using the BMS (Black-Merton-Scholes) model. For this model, assume that the parameters are $$Θ={σ}={0.3}$$. Compute the put prices according to the following table:

$$α$$ $$η=0.1$$
$$N=26$$
$$η=0.1$$
$$N=210$$
$$η=0.25$$
$$N=26$$
$$η=0.25$$
$$N=210$$
-1.01
-1.25
-1.5
-1.75
-2
-5

Q13. Given your numerical results what is the value of the put option under the BMS model? Answer: 2.708

Next, we will price put options using the Heston model. For this model, assume that the parameters are $$Θ={κ,θ,λ,ρ,v_0​}={2,0.05,0.3,-0.7,0.04}$$. Compute the put prices according to the following table:

$$α$$ $$η=0.1$$
$$N=26$$
$$η=0.1$$
$$N=210$$
$$η=0.25$$
$$N=26$$
$$η=0.25$$
$$N=210$$
-1.01
-1.25
-1.5
-1.75
-2
-5

Q14. Given your numerical results what is the value of the put option under the Heston model? Answer: 1.341

Finally, we will price put options using the VG (Variance Gamma) model. For this model, assume that the parameters are $$Θ={σ,ν,θ}={0.3,0.5,-0.4}$$. Compute the put prices according to the following table:

$$α$$ $$η=0.1$$
$$N=26$$
$$η=0.1$$
$$N=210$$
$$η=0.25$$
$$N=26$$
$$η=0.25$$
$$N=210$$
-1.01
-1.25
-1.5
-1.75
-2
-5

Q15. Given your numerical results what is the value of the put option under the VG model? Answer: 5.3137

### Solutions

A1: $$P\&L = 100 · (S_6-S_0)$$
A2: $$P\&L = 100 · (S_6-S_0)$$
A3: $$P\&L = 100 · ((S_6-K)^{+}-C)$$
A4: $$P\&L = 100 · ((S_6-K)^{+}-C)$$
A5: $$V=100 · S_6$$
A6: $$V=100 · (S_6 + (K-S_6)^{+}-P)$$
A7: $$V=100 · (S_6 + (K-S_6)^{+}-P)$$
A8: $$\eta = K/N = K/2^n$$ for puts.
A9: $$\eta = K/N = K/2^n$$ for puts.

Modeling:

import numpy as np # for fast vector computations
from time import time # to obtain the running time of some functions


We assume that the stock price follows a lognormal distribution. The following is a function that computes the lognormal density.

def logNormal(S, r, q, sig, S0, T):
''' Computes the log normal density function for the stock price. '''
f = np.exp(-0.5*((np.log(S/S0)-(r-q-sig**2/2)*T)/(sig*np.sqrt(T)))**2) / (sig*S*np.sqrt(2*np.pi*T))
return f


Warning: be careful when evaluating this function at $$S=0.0$$. The theoretical value should be $$f=0.0$$, but the function above will give you a warning and a NaN value. A quick fix is to evaluate it at $$S=\epsilon$$, where $$\epsilon$$ is a small positive number such as 1e-8.

We want to estimate the price of a European put with the following parameters:

# Parameters for put
S0 = 100
K = 90
r = 0.04
q = 0.02
sig = 0.25
T = 1.0


Now, we will use the numerical integration method to price this put. For this, we create the following function that takes as arguments the above parameters and a positive integer $$N$$ which is the number of grid points.

def numerical_integral_put(r, q, S0, K, sig, T, N):
''' Numerical integration for puts. '''

# temporary value for output
eta = 0.0 # spacing of the grid
priceP = 0.0 # price of the put

#discount factor
df = np.exp(-r * T)
# step size
eta = 1. * K / N
# vector of stock prices
S = np.arange(0, N) * eta
# avoid numerical issues
S[0] = 1e-8
# vector of weights
w = np.ones(N) * eta
w[0] = eta / 2
# lognormal densities
logN = np.zeros(N)
logN = logNormal(S, r, q, sig, S0, T)
# numerical integral
sumP = np.sum((K - S) * logN * w)
# price put
priceP = df * sumP

return eta, priceP


We use the previous function to price the put for different values of $$N=2^n$$, $$n=1,2,\ldots,15$$.

# vector with all values of n
n_min = 1
n_max = 15
n_vec = np.arange(n_min, n_max + 1, dtype=int)

# compute the numerical integration for each value of n
start_time = time()
# will store the results in vectors
eta_vec = np.zeros(n_max)
put_vec = np.zeros(n_max)
for i in range(n_max):
N = 2** n_vec[i]

[eta_vec[i], put_vec[i]] = numerical_integral_put(r, q, S0, K, sig, T, N)

end_time = time()
print('Elapsed time (seconds): ' + str(end_time - start_time))

Elapsed time (seconds): 0.8806734085083008


We print a table with the results.

# print a table with the numerical integration for each value of n
print('N\teta\tP_0')
for i in range(n_max):
print('2^%i\t%.3f\t%.4f' % (n_vec[i], eta_vec[i], put_vec[i]))

N	eta	P_0
2^1	45.000	0.4847
2^2	22.500	3.8251
2^3	11.250	4.3546
2^4	5.625	4.4822
2^5	2.812	4.5137
2^6	1.406	4.5216
2^7	0.703	4.5236
2^8	0.352	4.5241
2^9	0.176	4.5242
2^10	0.088	4.5242
2^11	0.044	4.5242
2^12	0.022	4.5242
2^13	0.011	4.5242
2^14	0.005	4.5242
2^15	0.003	4.5242


A9:

# eta, when n = 3
n_question = 3
# remember that Python indices start at 0
print(eta_vec[n_question-1])

11.25


A10:

# eta, when n = 10
n_question = 10
# remember that Python indices start at 0
print(eta_vec[n_question-1])

0.087890625


A11:

# put, when n = 4
n_question = 4
# remember that Python indices start at 0
print(put_vec[n_question-1])

4.48223458094208


A12:

# put, when n = 15
n_question = 15
# remember that Python indices start at 0
print(put_vec[n_question-1])

4.524213707706823


During the lectures, we saw how to price call options using FFT (the extension to put options is straighforward).

We first need a function that computes the characteristic function of the log-stock price for 3 different models: BMS, Heston, and VG (each one with different sets of parameters). The function below will do the job! The part of Heston and VG is already complete, you just need to implement the BMS part.

def generic_CF(u, params, S0, r, q, T, model):
''' Computes the characteristic function for different models (BMS, Heston, VG). '''

if (model == 'BMS'):
# unpack parameters
sig = params[0]

mu = np.log(S0) + (r-q-sig**2/2)*T
a = sig*np.sqrt(T)
phi = np.exp(1j*mu*u-(a*u)**2/2)

elif(model == 'Heston'):
# unpack parameters
kappa  = params[0]
theta  = params[1]
sigma  = params[2]
rho    = params[3]
v0     = params[4]
# cf
tmp = (kappa-1j*rho*sigma*u)
g = np.sqrt((sigma**2)*(u**2+1j*u)+tmp**2)
pow1 = 2*kappa*theta/(sigma**2)
numer1 = (kappa*theta*T*tmp)/(sigma**2) + 1j*u*T*r + 1j*u*np.log(S0)
log_denum1 = pow1 * np.log(np.cosh(g*T/2)+(tmp/g)*np.sinh(g*T/2))
tmp2 = ((u*u+1j*u)*v0)/(g/np.tanh(g*T/2)+tmp)
log_phi = numer1 - log_denum1 - tmp2
phi = np.exp(log_phi)

elif (model == 'VG'):
# unpack parameters
sigma  = params[0];
nu     = params[1];
theta  = params[2];
# cf
if (nu == 0):
mu = np.log(S0) + (r-q - theta -0.5*sigma**2)*T
phi  = np.exp(1j*u*mu) * np.exp((1j*theta*u-0.5*sigma**2*u**2)*T)
else:
mu  = np.log(S0) + (r-q + np.log(1-theta*nu-0.5*sigma**2*nu)/nu)*T
phi = np.exp(1j*u*mu)*((1-1j*nu*theta*u+0.5*nu*sigma**2*u**2)**(-T/nu))

return phi


We now provide you with a function that uses FFT to price European option for any of the 3 models we have discussed. The same function works for both calls and puts (you should ask yourself: how and why?). The function is complete, you don’t need to add any code. It return two vectors, one with the log-strikes and one with option prices.

def genericFFT(params, S0, K, r, q, T, alpha, eta, n, model):
''' Option pricing using FFT (model-free). '''

N = 2**n

# step-size in log strike space
lda = (2 * np.pi / N) / eta

# choice of beta
#beta = np.log(S0)-N*lda/2 # the log strike we want is in the middle of the array
beta = np.log(K) # the log strike we want is the first element of the array

# forming vector x and strikes km for m=1,...,N
km = np.zeros(N)
xX = np.zeros(N)

# discount factor
df = np.exp(-r*T)

nuJ = np.arange(N) * eta
psi_nuJ = generic_CF(nuJ - (alpha + 1) * 1j, params, S0, r, q, T, model) / ((alpha + 1j*nuJ)*(alpha+1+1j*nuJ))

km = beta + lda * np.arange(N)
w = eta * np.ones(N)
w[0] = eta / 2
xX = np.exp(-1j * beta * nuJ) * df * psi_nuJ * w

yY = np.fft.fft(xX)
cT_km = np.zeros(N)
multiplier = np.exp(-alpha * km) / np.pi
cT_km = multiplier * np.real(yY)

return km, cT_km


We want to estimate a European put with the following parameters:

# parameters
S0 = 100
K = 80
r = 0.05
q = 0.01
T = 1.0


For the 3 models, we consider the same values of $$\alpha$$, $$\eta$$, and $$n$$.

# parameters for fft
eta_vec = np.array([0.1, 0.25])
n_vec = np.array([6, 10])
alpha_vec = np.array([-1.01, -1.25, -1.5, -1.75, -2., -5.])
num_prices = len(eta_vec) * len(n_vec) * len(alpha_vec)


We write a function that given a model and the above parameters will compute the put price.

def price_all_puts(params, S0, K, r, q, T, model, alpha_vec, eta_vec, n_vec):
num_prices = len(eta_vec) * len(n_vec) * len(alpha_vec)
# output is a matrix, the columns correspond to eta, n, alpha, and put price
put_matrix = np.zeros([num_prices, 4])
i = 0
for eta in eta_vec:
for n in n_vec:
for alpha in alpha_vec:
# pricing via fft
put = 0 # store here the put price
###############################################################
###############################################################
###############################################################
k_vec, option_vec = genericFFT(params, S0, K, r, q, T, alpha, eta, n, model)
put = option_vec[0] # only interested in one strike
###############################################################
###############################################################
###############################################################
put_matrix[i] = np.array([eta, n, alpha, put])
i += 1
return put_matrix


#### BMS

These are the parameters for the model.

# model-specific parameters
mod = 'BMS'
sig = 0.3
params = [sig]

start_time = time()
bms_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)
end_time = time()
print('Elapsed time (seconds): ' + str(end_time - start_time))

Elapsed time (seconds): 0.11588907241821289

# print results in table form
print('Model = ' + mod)
print('eta\tN\talpha\tput')
for i in range(num_prices):
print('%.2f\t2^%i\t%.2f\t%.4f' % (bms_matrix[i,0], bms_matrix[i,1], bms_matrix[i,2], bms_matrix[i,3]))

Model = BMS
eta	N	alpha	put
0.10	2^6	-1.01	89.7431
0.10	2^6	-1.25	2.7396
0.10	2^6	-1.50	2.7569
0.10	2^6	-1.75	2.7701
0.10	2^6	-2.00	2.7789
0.10	2^6	-5.00	2.6727
0.10	2^10	-1.01	89.7316
0.10	2^10	-1.25	2.7080
0.10	2^10	-1.50	2.7080
0.10	2^10	-1.75	2.7080
0.10	2^10	-2.00	2.7080
0.10	2^10	-5.00	2.7080
0.25	2^6	-1.01	269.0367
0.25	2^6	-1.25	2.8504
0.25	2^6	-1.50	2.7083
0.25	2^6	-1.75	2.7080
0.25	2^6	-2.00	2.7080
0.25	2^6	-5.00	2.7080
0.25	2^10	-1.01	269.0367
0.25	2^10	-1.25	2.8504
0.25	2^10	-1.50	2.7083
0.25	2^10	-1.75	2.7080
0.25	2^10	-2.00	2.7080
0.25	2^10	-5.00	2.7080


A13: 2.7080

#### Heston

# model-specific parameters
mod = 'Heston'
kappa = 2.
theta = 0.05
lda = 0.3
rho = -0.7
v0 = 0.04
params = [kappa, theta, lda, rho, v0]

start_time = time()
heston_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)
end_time = time()
print('Elapsed time (seconds): ' + str(end_time - start_time))

Elapsed time (seconds): 0.04112815856933594

# print results in table form
print('Model = ' + mod)
print('eta\tN\talpha\tput')
for i in range(num_prices):
print('%.2f\t2^%i\t%.2f\t%.4f' % (heston_matrix[i,0], heston_matrix[i,1], heston_matrix[i,2], heston_matrix[i,3]))

Model = Heston
eta	N	alpha	put
0.10	2^6	-1.01	88.0744
0.10	2^6	-1.25	1.1102
0.10	2^6	-1.50	1.1666
0.10	2^6	-1.75	1.2167
0.10	2^6	-2.00	1.2605
0.10	2^6	-5.00	1.3979
0.10	2^10	-1.01	88.3648
0.10	2^10	-1.25	1.3412
0.10	2^10	-1.50	1.3412
0.10	2^10	-1.75	1.3412
0.10	2^10	-2.00	1.3412
0.10	2^10	-5.00	1.3412
0.25	2^6	-1.01	267.6738
0.25	2^6	-1.25	1.4873
0.25	2^6	-1.50	1.3450
0.25	2^6	-1.75	1.3445
0.25	2^6	-2.00	1.3442
0.25	2^6	-5.00	1.3415
0.25	2^10	-1.01	267.6698
0.25	2^10	-1.25	1.4835
0.25	2^10	-1.50	1.3414
0.25	2^10	-1.75	1.3412
0.25	2^10	-2.00	1.3412
0.25	2^10	-5.00	1.3412


A14: 1.3412

#### VG

# model-specific parameters
mod = 'VG'
sigma = 0.3
nu = 0.5
theta = -0.4
params = [sigma, nu, theta]

start_time = time()
vg_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)
end_time = time()
print('Elapsed time (seconds): ' + str(end_time - start_time))

Elapsed time (seconds): 0.020387887954711914

# print results in table form
print('Model = ' + mod)
print('eta\tN\talpha\tput')
for i in range(num_prices):
print('%.2f\t2^%i\t%.2f\t%.4f' % (vg_matrix[i,0], vg_matrix[i,1], vg_matrix[i,2], vg_matrix[i,3]))

Model = VG
eta	N	alpha	put
0.10	2^6	-1.01	92.2103
0.10	2^6	-1.25	5.2047
0.10	2^6	-1.50	5.2230
0.10	2^6	-1.75	5.2401
0.10	2^6	-2.00	5.2556
0.10	2^6	-5.00	0.0030
0.10	2^10	-1.01	92.3373
0.10	2^10	-1.25	5.3137
0.10	2^10	-1.50	5.3137
0.10	2^10	-1.75	5.3137
0.10	2^10	-2.00	5.3137
0.10	2^10	-5.00	0.0000
0.25	2^6	-1.01	271.6399
0.25	2^6	-1.25	5.4539
0.25	2^6	-1.50	5.3121
0.25	2^6	-1.75	5.3122
0.25	2^6	-2.00	5.3124
0.25	2^6	-5.00	0.0019
0.25	2^10	-1.01	271.6423
0.25	2^10	-1.25	5.4560
0.25	2^10	-1.50	5.3139
0.25	2^10	-1.75	5.3137
0.25	2^10	-2.00	5.3137
0.25	2^10	-5.00	0.0020


A15: 5.3137