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))
me.add_constant('initial_value', 100.)
me.add_constant('volatility', 0.2)
me.add_constant('final_date', dt.datetime(2015, 12, 31))
me.add_constant('currency', 'EUR')

[4]:

# jump component
me.add_constant('lambda', 0.4)
me.add_constant('mu', -0.6)
me.add_constant('delta', 0.2)

[5]:

# stochastic volatiltiy component
me.add_constant('rho', -.5)
me.add_constant('kappa', 5.0)
me.add_constant('theta', 0.02)
me.add_constant('vol_vol', 0.3)

[6]:

# valuation environment
val_env = dx.market_environment('val_env', dt.datetime(2015, 1, 1))
val_env.add_constant('paths', 55000)
    # 25,000 paths
val_env.add_constant('frequency', 'D')
    # weekly frequency
val_env.add_curve('discount_curve', r)
val_env.add_constant('starting_date', dt.datetime(2015, 1, 1))
val_env.add_constant('final_date', dt.datetime(2015, 12, 31))

[7]:

# add valuation environment to market environment
me.add_environment(val_env)

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))
me_option.add_constant('maturity', dt.datetime(2015, 12, 31))
me_option.add_constant('strike', 100.)
me_option.add_constant('currency', 'EUR')
me_option.add_environment(me)
me_option.add_environment(val_env)

[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)
        me_option.add_constant('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()
            me_option.add_constant('strike', strike)
            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

Copyright, License & Disclaimer

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

DX Analytics (the “dx library” or “dx package”) is licensed under the GNU Affero General Public License version 3 or later (see http://www.gnu.org/licenses/).

DX Analytics comes with no representations or warranties, to the extent permitted by applicable law.

http://tpq.io | dx@tpq.io | http://twitter.com/dyjh

Quant Platform | http://pqp.io

Python for Finance Training | http://training.tpq.io

Certificate in Computational Finance | http://compfinance.tpq.io

Derivatives Analytics with Python (Wiley Finance) | http://dawp.tpq.io

Python for Finance (2nd ed., O’Reilly) | http://py4fi.tpq.io