# 9. Fourier-based Option Pricing¶

For several reasons, it is beneficial to have available alternative valuation and pricing approaches to the Monte Carlo simulation approach. One application area is to benchmark Monte Carlo-based valuation results against other (potentially more accurate) results. Another area is model calibration to liquidly traded vanilla instruments where generally faster numerial methods can be applied.

This part introduces Fouried-based valuation functions and benchmarks valuation results from the “standard”, simulation-based DX Analytics modeling approach to output of those functions.

[1]:

import dx
import datetime as dt


## 9.1. Risk Factors¶

The examples and benchmarks to follow rely on four different models:

• geometric Brownian motion (Black-Scholes-Merton 1973)
• jump diffusion (Merton 1976)
• stochastic volatility (Heston 1993)
• stochastic volatility jump diffusion (Bates 1996)

For details on these models and the Fourier-based option pricing approach refer to Hilpisch (2015) (cf. http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119037999.html).

We first define the single market and valuation environments.

[2]:

# constant short rate
r = dx.constant_short_rate('r', 0.01)

[3]:

# geometric Brownian motion
me = dx.market_environment('me', dt.datetime(2015, 1, 1))

[4]:

# jump component

[5]:

# stochastic volatiltiy component

[6]:

# valuation environment
val_env = dx.market_environment('val_env', dt.datetime(2015, 1, 1))
# 25,000 paths
# weekly frequency

[7]:

# add valuation environment to market environment


Equipped with the single market environments and the valuation environment, we can instantiate the simulation model objects.

[8]:

gbm = dx.geometric_brownian_motion('gbm', me)

[9]:

jd = dx.jump_diffusion('jd', me)

[10]:

sv = dx.stochastic_volatility('sv', me)

[11]:

svjd = dx.stoch_vol_jump_diffusion('svjd', me)


## 9.2. Plain Vanilla Put and Call Options¶

Based on the just defined risk factors, we define 8 diffent options—a European put and call option per risk factor, respectively.

[12]:

# market environment for the options
me_option = dx.market_environment('option', dt.datetime(2015, 1, 1))

[13]:

euro_put_gbm = dx.valuation_mcs_european_single('euro_put', gbm, me_option,
'np.maximum(strike - maturity_value, 0)')
euro_call_gbm = dx.valuation_mcs_european_single('euro_call', gbm, me_option,
'np.maximum(maturity_value - strike, 0)')

[14]:

euro_put_jd = dx.valuation_mcs_european_single('euro_put', jd, me_option,
'np.maximum(strike - maturity_value, 0)')
euro_call_jd = dx.valuation_mcs_european_single('euro_call', jd, me_option,
'np.maximum(maturity_value - strike, 0)')

[15]:

euro_put_sv = dx.valuation_mcs_european_single('euro_put', sv, me_option,
'np.maximum(strike - maturity_value, 0)')
euro_call_sv = dx.valuation_mcs_european_single('euro_call', sv, me_option,
'np.maximum(maturity_value - strike, 0)')

[16]:

euro_put_svjd = dx.valuation_mcs_european_single('euro_put', svjd, me_option,
'np.maximum(strike - maturity_value, 0)')
euro_call_svjd = dx.valuation_mcs_european_single('euro_call', svjd, me_option,
'np.maximum(maturity_value - strike, 0)')


## 9.3. Valuation Benchmarking¶

In this sub-section, we benchmark the Monte Carlo value estimates against the Fourier-based pricing results.

[17]:

import numpy as np
import pandas as pd


We first define some parameters used throughout.

[18]:

freq = '2m'  # used for maturity definitions
periods = 3  # number of intervals for maturity grid
strikes = 5  # number of strikes per maturity
initial_value = 100  # initial value for all risk factors
start = 0.8  # lowest strike in percent of spot
end = 1.2  # highest strike in percent of spot
start_date = '2015/3/1'  # start date for simulation/pricing


### 9.3.1. Geometric Brownian Motion¶

We need to initialize the valuation object first.

[19]:

euro_put_gbm.present_value()
# method call needed for initialization

[19]:

7.401144


There is a valuation class for European put and call options in the Black-Scholes-Merton model available called BSM_european_option. It is based on the analytical pricing formula for that model and is instantiated as follows:

[20]:

bsm_option = dx.BSM_european_option('bsm_opt', me_option)


The following routine benchmarks the Monte Carlo value estimates for the European put option against the output from the valuation object based on the analytical pricing formula. The results are quite good since this model is quite easy to discretize exactly and therefore generally shows good convergence of the Monte Carlo estimates.

[21]:

%%time
# European put
print('%4s  | %7s | %7s | %7s | %7s | %7s' % ('T', 'strike', 'mcs', 'fou', 'dif', 'rel'))
for maturity in pd.date_range(start=start_date, freq=freq, periods=periods):
bsm_option.maturity = maturity
euro_put_gbm.update(maturity=maturity)
for strike in np.linspace(start, end, strikes) * initial_value:
T = (maturity - me_option.pricing_date).days / 365.
euro_put_gbm.update(strike=strike)
mcs = euro_put_gbm.present_value()
bsm_option.strike = strike
ana = bsm_option.put_value()
print('%4.3f | %7.3f | %7.4f | %7.4f | %7.4f | %7.2f '
% (T, strike, mcs, ana, mcs - ana, (mcs - ana) / ana * 100))

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 |  0.0327 |  0.0338 | -0.0011 |   -3.22
0.244 |  90.000 |  0.6532 |  0.6524 |  0.0008 |    0.12
0.244 | 100.000 |  3.7811 |  3.8130 | -0.0318 |   -0.84
0.244 | 110.000 | 10.6205 | 10.6957 | -0.0752 |   -0.70
0.244 | 120.000 | 19.7026 | 19.8537 | -0.1512 |   -0.76
0.411 |  80.000 |  0.1755 |  0.1748 |  0.0007 |    0.42
0.411 |  90.000 |  1.3393 |  1.3241 |  0.0152 |    1.15
0.411 | 100.000 |  4.8421 |  4.8985 | -0.0564 |   -1.15
0.411 | 110.000 | 11.3570 | 11.4275 | -0.0705 |   -0.62
0.411 | 120.000 | 19.9190 | 20.0325 | -0.1135 |   -0.57
0.578 |  80.000 |  0.3820 |  0.3917 | -0.0097 |   -2.47
0.578 |  90.000 |  1.9342 |  1.9466 | -0.0124 |   -0.64
0.578 | 100.000 |  5.7446 |  5.7593 | -0.0147 |   -0.25
0.578 | 110.000 | 12.0049 | 12.0934 | -0.0885 |   -0.73
0.578 | 120.000 | 20.2267 | 20.3153 | -0.0886 |   -0.44
CPU times: user 38.5 s, sys: 6.24 s, total: 44.7 s
Wall time: 9.89 s


The same now for the European call option.

[22]:

euro_call_gbm.present_value()
# method call needed for initialization

[22]:

8.426388

[23]:

%%time
# European calls
print('%4s  | %7s | %7s | %7s | %7s | %7s' % ('T', 'strike', 'mcs', 'fou', 'dif', 'rel'))
for maturity in pd.date_range(start=start_date, freq=freq, periods=periods):
euro_call_gbm.update(maturity=maturity)
for strike in np.linspace(start, end, strikes) * initial_value:
T = (maturity - me_option.pricing_date).days / 365.
euro_call_gbm.update(strike=strike)
mcs = euro_call_gbm.present_value()
bsm_option.strike = strike
bsm_option.maturity = maturity
ana = bsm_option.call_value()
print('%4.3f | %7.3f | %7.4f | %7.4f | %7.4f | %7.2f ' \
% (T, strike, mcs, ana, mcs - ana, (mcs - ana) / ana * 100))

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 | 20.0731 | 20.2286 | -0.1556 |   -0.77
0.244 |  90.000 | 10.7908 | 10.8716 | -0.0808 |   -0.74
0.244 | 100.000 |  4.0200 |  4.0565 | -0.0365 |   -0.90
0.244 | 110.000 |  0.9631 |  0.9636 | -0.0005 |   -0.05
0.244 | 120.000 |  0.1356 |  0.1460 | -0.0104 |   -7.14
0.411 |  80.000 | 20.3867 | 20.5029 | -0.1162 |   -0.57
0.411 |  90.000 | 11.6556 | 11.6932 | -0.0376 |   -0.32
0.411 | 100.000 |  5.3106 |  5.3086 |  0.0020 |    0.04
0.411 | 110.000 |  1.8540 |  1.8787 | -0.0247 |   -1.31
0.411 | 120.000 |  0.5041 |  0.5246 | -0.0205 |   -3.91
0.578 |  80.000 | 20.7445 | 20.8528 | -0.1083 |   -0.52
0.578 |  90.000 | 12.3760 | 12.4654 | -0.0894 |   -0.72
0.578 | 100.000 |  6.3499 |  6.3357 |  0.0142 |    0.22
0.578 | 110.000 |  2.7110 |  2.7274 | -0.0165 |   -0.60
0.578 | 120.000 |  1.0002 |  1.0070 | -0.0068 |   -0.67
CPU times: user 36.9 s, sys: 5.8 s, total: 42.7 s
Wall time: 9.39 s


### 9.3.2. Benchmarking Function¶

All other valuation benchmarks are generated with Fourier-based pricing functions for which the handling is identical. We therefore use the following function for the benchmarks from now on:

[24]:

def valuation_benchmarking(valuation_object, fourier_function):
print('%4s  | %7s | %7s | %7s | %7s | %7s' % ('T', 'strike', 'mcs', 'fou', 'dif', 'rel'))
for maturity in pd.date_range(start=start_date, freq=freq, periods=periods):
valuation_object.update(maturity=maturity)
for strike in np.linspace(start, end, strikes) * initial_value:
T = (maturity - me_option.pricing_date).days / 365.
valuation_object.update(strike=strike)
mcs = valuation_object.present_value()
fou = fourier_function(me_option)
print('%4.3f | %7.3f | %7.4f | %7.4f | %7.4f | %7.3f '
% (T, strike, mcs, fou, mcs - fou, (mcs - fou) / fou * 100))


### 9.3.3. Jump Diffusion¶

The next model is the jump diffusion as proposed by Merton (1976).

[25]:

euro_put_jd.present_value()
# method call needed for initialization

[25]:

14.106036


There is a Fourier-based pricing function available which is called M76_put_value and which is used for the benchmarking for the European put options that follows.

[26]:

%time valuation_benchmarking(euro_put_jd, dx.M76_put_value)

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 |  2.1771 |  2.1594 |  0.0177 |   0.819
0.244 |  90.000 |  3.2142 |  3.2826 | -0.0685 |  -2.086
0.244 | 100.000 |  5.8588 |  5.8842 | -0.0254 |  -0.431
0.244 | 110.000 | 11.4707 | 11.6115 | -0.1408 |  -1.213
0.244 | 120.000 | 19.9428 | 20.0857 | -0.1429 |  -0.712
0.411 |  80.000 |  3.4273 |  3.4505 | -0.0233 |  -0.674
0.411 |  90.000 |  5.1224 |  5.2162 | -0.0938 |  -1.798
0.411 | 100.000 |  8.2061 |  8.2266 | -0.0205 |  -0.249
0.411 | 110.000 | 13.3971 | 13.4430 | -0.0458 |  -0.341
0.411 | 120.000 | 20.8313 | 20.9238 | -0.0925 |  -0.442
0.578 |  80.000 |  4.5777 |  4.6090 | -0.0313 |  -0.679
0.578 |  90.000 |  6.8477 |  6.8782 | -0.0304 |  -0.443
0.578 | 100.000 | 10.2062 | 10.2077 | -0.0015 |  -0.014
0.578 | 110.000 | 15.0388 | 15.2251 | -0.1862 |  -1.223
0.578 | 120.000 | 21.8634 | 22.0886 | -0.2252 |  -1.020
CPU times: user 1min 29s, sys: 11.5 s, total: 1min 40s
Wall time: 22.6 s


Accordingly, the benchmarking for the European call options based on the Fourier-based M76_call_value function.

[27]:

euro_call_jd.present_value()
# method call needed for initialization

[27]:

15.076098

[28]:

%time valuation_benchmarking(euro_call_jd, dx.M76_call_value)

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 | 22.1538 | 22.3543 | -0.2005 |  -0.897
0.244 |  90.000 | 13.3560 | 13.5018 | -0.1458 |  -1.080
0.244 | 100.000 |  6.0947 |  6.1277 | -0.0330 |  -0.539
0.244 | 110.000 |  1.8695 |  1.8794 | -0.0098 |  -0.524
0.244 | 120.000 |  0.3701 |  0.3780 | -0.0079 |  -2.086
0.411 |  80.000 | 23.6404 | 23.7786 | -0.1383 |  -0.582
0.411 |  90.000 | 15.4757 | 15.5853 | -0.1096 |  -0.703
0.411 | 100.000 |  8.6258 |  8.6367 | -0.0109 |  -0.126
0.411 | 110.000 |  3.8863 |  3.8941 | -0.0078 |  -0.200
0.411 | 120.000 |  1.4004 |  1.4160 | -0.0156 |  -1.101
0.578 |  80.000 | 25.0494 | 25.0701 | -0.0207 |  -0.083
0.578 |  90.000 | 17.3956 | 17.3970 | -0.0014 |  -0.008
0.578 | 100.000 | 10.7217 | 10.7841 | -0.0624 |  -0.579
0.578 | 110.000 |  5.8156 |  5.8591 | -0.0435 |  -0.743
0.578 | 120.000 |  2.7383 |  2.7803 | -0.0420 |  -1.511
CPU times: user 1min 30s, sys: 11.9 s, total: 1min 42s
Wall time: 23 s


### 9.3.4. Stochastic Volatility¶

Stochastic volatility models like the one of Heston (1993) are popular to reproduce implied volatility smiles observed in markets. First, the benchmarking for the European put options using the Fourier-based H93_put_value function.

[29]:

euro_put_sv.present_value()
# method call needed for initialization

[29]:

5.3329

[30]:

%time valuation_benchmarking(euro_put_sv, dx.H93_put_value)

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 |  0.0504 |  0.0504 | -0.0000 |  -0.011
0.244 |  90.000 |  0.5512 |  0.5741 | -0.0229 |  -3.984
0.244 | 100.000 |  3.2609 |  3.3204 | -0.0595 |  -1.791
0.244 | 110.000 | 10.1399 | 10.2404 | -0.1005 |  -0.981
0.244 | 120.000 | 19.5842 | 19.7354 | -0.1512 |  -0.766
0.411 |  80.000 |  0.1444 |  0.1603 | -0.0159 |  -9.910
0.411 |  90.000 |  0.9420 |  1.0063 | -0.0643 |  -6.386
0.411 | 100.000 |  3.9185 |  4.0234 | -0.1049 |  -2.607
0.411 | 110.000 | 10.4422 | 10.5487 | -0.1065 |  -1.010
0.411 | 120.000 | 19.5004 | 19.6293 | -0.1288 |  -0.656
0.578 |  80.000 |  0.2502 |  0.2838 | -0.0336 | -11.847
0.578 |  90.000 |  1.2868 |  1.3631 | -0.0763 |  -5.596
0.578 | 100.000 |  4.4023 |  4.5467 | -0.1444 |  -3.176
0.578 | 110.000 | 10.6897 | 10.8389 | -0.1492 |  -1.376
0.578 | 120.000 | 19.4625 | 19.5767 | -0.1142 |  -0.583
CPU times: user 1min 43s, sys: 16.2 s, total: 1min 59s
Wall time: 26.4 s


Second, the benchmarking for the European call options based on the Fourier-based H93_call_value function.

[31]:

euro_call_sv.present_value()
# method call needed for initialization

[31]:

6.307286

[32]:

%time valuation_benchmarking(euro_call_sv, dx.H93_call_value)

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 | 20.0909 | 20.2453 | -0.1544 |  -0.763
0.244 |  90.000 | 10.6943 | 10.7933 | -0.0989 |  -0.917
0.244 | 100.000 |  3.5110 |  3.5639 | -0.0529 |  -1.483
0.244 | 110.000 |  0.4705 |  0.5083 | -0.0378 |  -7.441
0.244 | 120.000 |  0.0274 |  0.0276 | -0.0002 |  -0.549
0.411 |  80.000 | 20.3465 | 20.4884 | -0.1419 |  -0.693
0.411 |  90.000 | 11.2741 | 11.3754 | -0.1012 |  -0.890
0.411 | 100.000 |  4.3208 |  4.4335 | -0.1127 |  -2.543
0.411 | 110.000 |  0.9430 |  0.9998 | -0.0568 |  -5.685
0.411 | 120.000 |  0.1072 |  0.1214 | -0.0142 | -11.702
0.578 |  80.000 | 20.6329 | 20.7450 | -0.1121 |  -0.540
0.578 |  90.000 | 11.7564 | 11.8818 | -0.1254 |  -1.056
0.578 | 100.000 |  5.0226 |  5.1231 | -0.1005 |  -1.961
0.578 | 110.000 |  1.3453 |  1.4729 | -0.1276 |  -8.665
0.578 | 120.000 |  0.2373 |  0.2684 | -0.0311 | -11.598
CPU times: user 1min 44s, sys: 16.5 s, total: 2min
Wall time: 26.6 s


### 9.3.5. Stochastic Volatility Jump-Diffusion¶

Finally, we consider the combination of the stochastic volatility and jump diffusion models from before as proposed by Bates (1996). The Fourier-based pricing function for European put options is called B96_put_value.

[33]:

euro_put_svjd.present_value()
# method call needed for initialization

[33]:

12.938869

[34]:

%time valuation_benchmarking(euro_put_svjd, dx.B96_put_value)

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 |  2.1374 |  2.1638 | -0.0265 |  -1.224
0.244 |  90.000 |  3.1814 |  3.2761 | -0.0947 |  -2.891
0.244 | 100.000 |  5.4111 |  5.5889 | -0.1778 |  -3.182
0.244 | 110.000 | 11.0267 | 11.0733 | -0.0466 |  -0.421
0.244 | 120.000 | 19.5912 | 19.8344 | -0.2431 |  -1.226
0.411 |  80.000 |  3.3949 |  3.4366 | -0.0417 |  -1.213
0.411 |  90.000 |  5.0453 |  5.1339 | -0.0887 |  -1.727
0.411 | 100.000 |  7.5677 |  7.7747 | -0.2069 |  -2.662
0.411 | 110.000 | 12.5075 | 12.5744 | -0.0669 |  -0.532
0.411 | 120.000 | 19.8595 | 20.1827 | -0.3233 |  -1.602
0.578 |  80.000 |  4.3390 |  4.5480 | -0.2090 |  -4.595
0.578 |  90.000 |  6.4859 |  6.7163 | -0.2304 |  -3.431
0.578 | 100.000 |  9.2934 |  9.6585 | -0.3651 |  -3.780
0.578 | 110.000 | 13.9231 | 14.1776 | -0.2545 |  -1.795
0.578 | 120.000 | 20.6726 | 20.9423 | -0.2697 |  -1.288
CPU times: user 2min 42s, sys: 22.1 s, total: 3min 4s
Wall time: 39.8 s


The Fourier-based counterpart function for European call options is called B96_call_value.

[35]:

euro_call_svjd.present_value()
# method call needed for initialization

[35]:

13.836211

[36]:

%time valuation_benchmarking(euro_call_svjd, dx.B96_call_value)

   T  |  strike |     mcs |     fou |     dif |     rel
0.244 |  80.000 | 22.1335 | 22.3587 | -0.2252 |  -1.007
0.244 |  90.000 | 13.3312 | 13.4953 | -0.1640 |  -1.216
0.244 | 100.000 |  5.6897 |  5.8325 | -0.1428 |  -2.448
0.244 | 110.000 |  1.2988 |  1.3411 | -0.0423 |  -3.156
0.244 | 120.000 |  0.1163 |  0.1266 | -0.0103 |  -8.110
0.411 |  80.000 | 23.5008 | 23.7647 | -0.2639 |  -1.110
0.411 |  90.000 | 15.3156 | 15.5030 | -0.1875 |  -1.209
0.411 | 100.000 |  7.9903 |  8.1848 | -0.1944 |  -2.376
0.411 | 110.000 |  2.8889 |  3.0255 | -0.1366 |  -4.514
0.411 | 120.000 |  0.6102 |  0.6749 | -0.0647 |  -9.583
0.578 |  80.000 | 24.7946 | 25.0091 | -0.2145 |  -0.858
0.578 |  90.000 | 16.9820 | 17.2351 | -0.2531 |  -1.469
0.578 | 100.000 |  9.9994 | 10.2349 | -0.2355 |  -2.300
0.578 | 110.000 |  4.5780 |  4.8117 | -0.2337 |  -4.857
0.578 | 120.000 |  1.4891 |  1.6340 | -0.1448 |  -8.864
CPU times: user 2min 46s, sys: 22.9 s, total: 3min 8s
Wall time: 40.9 s


## 9.4. Sources of Errors¶

Numerical methods like Monte Carlo simulation might suffer from different sources of errors, like for example:

• discretization error: every discretization of a continuous time interval—or a continuous state space to this end—leads to a so-called discretization error
• approximation errors: DX Analytics uses in several places approximative, Euler-based discretization schemes (e.g. for performance reasons and to allow for consistent correlation modeling) which are known to be biased
• numerical errors: the approximation of a continuous probability distribution by a finite, discrete set of (pseudo-) random numbers introduces also errors

© Dr. Yves J. Hilpisch | The Python Quants GmbH