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. .. code:: python import dx import datetime as dt 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**. .. code:: python # constant short rate r = dx.constant_short_rate('r', 0.01) .. code:: python # 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') .. code:: python # jump component me.add_constant('lambda', 0.4) me.add_constant('mu', -0.6) me.add_constant('delta', 0.2) .. code:: python # 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) .. code:: python # 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)) .. code:: python # 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**. .. code:: python gbm = dx.geometric_brownian_motion('gbm', me) .. code:: python jd = dx.jump_diffusion('jd', me) .. code:: python sv = dx.stochastic_volatility('sv', me) .. code:: python svjd = dx.stoch_vol_jump_diffusion('svjd', me) 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. .. code:: python # 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) .. code:: python 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)') .. code:: python 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)') .. code:: python 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)') .. code:: python 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)') Valuation Benchmarking ---------------------- In this sub-section, we benchmark the **Monte Carlo value estimates** against the **Fourier-based pricing results**. .. code:: python import numpy as np import pandas as pd We first define some parameters used throughout. .. code:: python 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 Geometric Brownian Motion ~~~~~~~~~~~~~~~~~~~~~~~~~ We need to initialize the valuation object first. .. code:: python euro_put_gbm.present_value() # method call needed for initialization .. parsed-literal:: 7.417061 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: .. code:: python 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. .. code:: python %%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) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 0.0357 | 0.0338 | 0.0019 | 5.51 0.244 | 90.000 | 0.6502 | 0.6524 | -0.0023 | -0.35 0.244 | 100.000 | 3.8025 | 3.8130 | -0.0104 | -0.27 0.244 | 110.000 | 10.6245 | 10.6957 | -0.0711 | -0.67 0.244 | 120.000 | 19.7058 | 19.8537 | -0.1479 | -0.74 0.411 | 80.000 | 0.1741 | 0.1748 | -0.0007 | -0.40 0.411 | 90.000 | 1.3236 | 1.3241 | -0.0004 | -0.03 0.411 | 100.000 | 4.8395 | 4.8985 | -0.0590 | -1.20 0.411 | 110.000 | 11.3407 | 11.4275 | -0.0869 | -0.76 0.411 | 120.000 | 19.9276 | 20.0325 | -0.1049 | -0.52 0.578 | 80.000 | 0.3882 | 0.3917 | -0.0035 | -0.89 0.578 | 90.000 | 1.9368 | 1.9466 | -0.0098 | -0.50 0.578 | 100.000 | 5.7718 | 5.7593 | 0.0125 | 0.22 0.578 | 110.000 | 12.0367 | 12.0934 | -0.0567 | -0.47 0.578 | 120.000 | 20.2347 | 20.3153 | -0.0806 | -0.40 CPU times: user 14.6 s, sys: 898 ms, total: 15.5 s Wall time: 15.5 s The same now for the **European call option**. .. code:: python euro_call_gbm.present_value() # method call needed for initialization .. parsed-literal:: 8.369015 .. code:: python %%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) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 20.0767 | 20.2286 | -0.1520 | -0.75 0.244 | 90.000 | 10.7843 | 10.8716 | -0.0873 | -0.80 0.244 | 100.000 | 4.0583 | 4.0565 | 0.0018 | 0.04 0.244 | 110.000 | 0.9560 | 0.9636 | -0.0076 | -0.78 0.244 | 120.000 | 0.1395 | 0.1460 | -0.0065 | -4.46 0.411 | 80.000 | 20.3924 | 20.5029 | -0.1105 | -0.54 0.411 | 90.000 | 11.6414 | 11.6932 | -0.0518 | -0.44 0.411 | 100.000 | 5.2774 | 5.3086 | -0.0312 | -0.59 0.411 | 110.000 | 1.8614 | 1.8787 | -0.0173 | -0.92 0.411 | 120.000 | 0.5270 | 0.5246 | 0.0024 | 0.46 0.578 | 80.000 | 20.7512 | 20.8528 | -0.1016 | -0.49 0.578 | 90.000 | 12.4161 | 12.4654 | -0.0493 | -0.40 0.578 | 100.000 | 6.3267 | 6.3357 | -0.0090 | -0.14 0.578 | 110.000 | 2.7082 | 2.7274 | -0.0192 | -0.71 0.578 | 120.000 | 1.0061 | 1.0070 | -0.0009 | -0.09 CPU times: user 14.1 s, sys: 783 ms, total: 14.9 s Wall time: 14.9 s 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: .. code:: python 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) Jump Diffusion ~~~~~~~~~~~~~~ The next model is the jump diffusion as proposed by **Merton (1976)**. .. code:: python euro_put_jd.present_value() # method call needed for initialization .. parsed-literal:: 14.272832 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. .. code:: python %time valuation_benchmarking(euro_put_jd, dx.M76_put_value) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 2.1959 | 2.1594 | 0.0364 | 1.687 0.244 | 90.000 | 3.3174 | 3.2826 | 0.0348 | 1.060 0.244 | 100.000 | 5.6949 | 5.8842 | -0.1893 | -3.217 0.244 | 110.000 | 11.5185 | 11.6115 | -0.0929 | -0.800 0.244 | 120.000 | 19.9645 | 20.0857 | -0.1213 | -0.604 0.411 | 80.000 | 3.3886 | 3.4505 | -0.0619 | -1.795 0.411 | 90.000 | 5.2458 | 5.2162 | 0.0296 | 0.568 0.411 | 100.000 | 8.1833 | 8.2266 | -0.0433 | -0.527 0.411 | 110.000 | 13.3401 | 13.4430 | -0.1029 | -0.765 0.411 | 120.000 | 20.6134 | 20.9238 | -0.3104 | -1.484 0.578 | 80.000 | 4.6632 | 4.6090 | 0.0542 | 1.177 0.578 | 90.000 | 6.8978 | 6.8782 | 0.0196 | 0.285 0.578 | 100.000 | 10.0860 | 10.2077 | -0.1217 | -1.192 0.578 | 110.000 | 15.2581 | 15.2251 | 0.0330 | 0.217 0.578 | 120.000 | 21.9876 | 22.0886 | -0.1010 | -0.457 CPU times: user 36.7 s, sys: 1.5 s, total: 38.2 s Wall time: 38.2 s Accordingly, the benchmarking for the **European call options** based on the Fourier-based ``M76_call_value`` function. .. code:: python euro_call_jd.present_value() # method call needed for initialization .. parsed-literal:: 15.141689 .. code:: python %time valuation_benchmarking(euro_call_jd, dx.M76_call_value) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 22.1912 | 22.3543 | -0.1631 | -0.730 0.244 | 90.000 | 13.3980 | 13.5018 | -0.1038 | -0.769 0.244 | 100.000 | 6.0899 | 6.1277 | -0.0378 | -0.618 0.244 | 110.000 | 1.8246 | 1.8794 | -0.0548 | -2.917 0.244 | 120.000 | 0.3882 | 0.3780 | 0.0102 | 2.700 0.411 | 80.000 | 23.6387 | 23.7786 | -0.1399 | -0.589 0.411 | 90.000 | 15.4714 | 15.5853 | -0.1139 | -0.731 0.411 | 100.000 | 8.5522 | 8.6367 | -0.0846 | -0.979 0.411 | 110.000 | 3.8717 | 3.8941 | -0.0224 | -0.574 0.411 | 120.000 | 1.4045 | 1.4160 | -0.0114 | -0.807 0.578 | 80.000 | 25.0104 | 25.0701 | -0.0597 | -0.238 0.578 | 90.000 | 17.3636 | 17.3970 | -0.0334 | -0.192 0.578 | 100.000 | 10.7576 | 10.7841 | -0.0265 | -0.246 0.578 | 110.000 | 5.8877 | 5.8591 | 0.0285 | 0.487 0.578 | 120.000 | 2.8004 | 2.7803 | 0.0201 | 0.723 CPU times: user 37.3 s, sys: 1.52 s, total: 38.8 s Wall time: 38.8 s 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. .. code:: python euro_put_sv.present_value() # method call needed for initialization .. parsed-literal:: 5.317924 .. code:: python %time valuation_benchmarking(euro_put_sv, dx.H93_put_value) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 0.0446 | 0.0504 | -0.0058 | -11.556 0.244 | 90.000 | 0.5315 | 0.5741 | -0.0426 | -7.418 0.244 | 100.000 | 3.2501 | 3.3204 | -0.0703 | -2.118 0.244 | 110.000 | 10.1363 | 10.2404 | -0.1042 | -1.017 0.244 | 120.000 | 19.5830 | 19.7354 | -0.1524 | -0.772 0.411 | 80.000 | 0.1385 | 0.1603 | -0.0218 | -13.620 0.411 | 90.000 | 0.9712 | 1.0063 | -0.0350 | -3.483 0.411 | 100.000 | 3.9225 | 4.0234 | -0.1009 | -2.508 0.411 | 110.000 | 10.4090 | 10.5487 | -0.1397 | -1.324 0.411 | 120.000 | 19.5018 | 19.6293 | -0.1275 | -0.649 0.578 | 80.000 | 0.2417 | 0.2838 | -0.0421 | -14.841 0.578 | 90.000 | 1.2879 | 1.3631 | -0.0751 | -5.512 0.578 | 100.000 | 4.4287 | 4.5467 | -0.1180 | -2.595 0.578 | 110.000 | 10.7018 | 10.8389 | -0.1370 | -1.264 0.578 | 120.000 | 19.4616 | 19.5767 | -0.1151 | -0.588 CPU times: user 33.3 s, sys: 3.95 s, total: 37.2 s Wall time: 37.2 s Second, the benchmarking for the **European call options** based on the Fourier-based ``H93_call_value`` function. .. code:: python euro_call_sv.present_value() # method call needed for initialization .. parsed-literal:: 6.248828 .. code:: python %time valuation_benchmarking(euro_call_sv, dx.H93_call_value) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 20.0898 | 20.2453 | -0.1555 | -0.768 0.244 | 90.000 | 10.6877 | 10.7933 | -0.1056 | -0.978 0.244 | 100.000 | 3.5002 | 3.5639 | -0.0637 | -1.787 0.244 | 110.000 | 0.4871 | 0.5083 | -0.0212 | -4.170 0.244 | 120.000 | 0.0231 | 0.0276 | -0.0045 | -16.147 0.411 | 80.000 | 20.3443 | 20.4884 | -0.1441 | -0.703 0.411 | 90.000 | 11.2545 | 11.3754 | -0.1208 | -1.062 0.411 | 100.000 | 4.3477 | 4.4335 | -0.0858 | -1.935 0.411 | 110.000 | 0.9492 | 0.9998 | -0.0506 | -5.061 0.411 | 120.000 | 0.1152 | 0.1214 | -0.0062 | -5.107 0.578 | 80.000 | 20.6306 | 20.7450 | -0.1144 | -0.551 0.578 | 90.000 | 11.7448 | 11.8818 | -0.1370 | -1.153 0.578 | 100.000 | 4.9837 | 5.1231 | -0.1394 | -2.721 0.578 | 110.000 | 1.3770 | 1.4729 | -0.0959 | -6.512 0.578 | 120.000 | 0.2390 | 0.2684 | -0.0294 | -10.951 CPU times: user 35.6 s, sys: 4.3 s, total: 39.9 s Wall time: 39.9 s 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``. .. code:: python euro_put_svjd.present_value() # method call needed for initialization .. parsed-literal:: 13.046603 .. code:: python %time valuation_benchmarking(euro_put_svjd, dx.B96_put_value) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 2.0128 | 2.1638 | -0.1510 | -6.979 0.244 | 90.000 | 3.1453 | 3.2761 | -0.1308 | -3.991 0.244 | 100.000 | 5.5157 | 5.5889 | -0.0732 | -1.309 0.244 | 110.000 | 10.8980 | 11.0733 | -0.1753 | -1.583 0.244 | 120.000 | 19.6077 | 19.8344 | -0.2266 | -1.143 0.411 | 80.000 | 3.2972 | 3.4366 | -0.1394 | -4.055 0.411 | 90.000 | 4.9726 | 5.1339 | -0.1613 | -3.142 0.411 | 100.000 | 7.4453 | 7.7747 | -0.3293 | -4.236 0.411 | 110.000 | 12.4031 | 12.5744 | -0.1713 | -1.362 0.411 | 120.000 | 19.8908 | 20.1827 | -0.2919 | -1.446 0.578 | 80.000 | 4.3281 | 4.5480 | -0.2199 | -4.834 0.578 | 90.000 | 6.5681 | 6.7163 | -0.1482 | -2.206 0.578 | 100.000 | 9.4546 | 9.6585 | -0.2039 | -2.111 0.578 | 110.000 | 13.9146 | 14.1776 | -0.2630 | -1.855 0.578 | 120.000 | 20.5290 | 20.9423 | -0.4132 | -1.973 CPU times: user 55.8 s, sys: 2.6 s, total: 58.4 s Wall time: 58.4 s The Fourier-based counterpart function for **European call options** is called ``B96_call_value``. .. code:: python euro_call_svjd.present_value() # method call needed for initialization .. parsed-literal:: 13.96036 .. code:: python %time valuation_benchmarking(euro_call_svjd, dx.B96_call_value) .. parsed-literal:: T | strike | mcs | fou | dif | rel 0.244 | 80.000 | 22.1387 | 22.3587 | -0.2200 | -0.984 0.244 | 90.000 | 13.3566 | 13.4953 | -0.1386 | -1.027 0.244 | 100.000 | 5.7380 | 5.8325 | -0.0945 | -1.620 0.244 | 110.000 | 1.2884 | 1.3411 | -0.0527 | -3.930 0.244 | 120.000 | 0.1113 | 0.1266 | -0.0153 | -12.051 0.411 | 80.000 | 23.5057 | 23.7647 | -0.2590 | -1.090 0.411 | 90.000 | 15.3243 | 15.5030 | -0.1787 | -1.153 0.411 | 100.000 | 7.9859 | 8.1848 | -0.1989 | -2.430 0.411 | 110.000 | 2.8734 | 3.0255 | -0.1521 | -5.029 0.411 | 120.000 | 0.6190 | 0.6749 | -0.0558 | -8.273 0.578 | 80.000 | 24.7489 | 25.0091 | -0.2603 | -1.041 0.578 | 90.000 | 16.9918 | 17.2351 | -0.2432 | -1.411 0.578 | 100.000 | 9.9490 | 10.2349 | -0.2859 | -2.793 0.578 | 110.000 | 4.5873 | 4.8117 | -0.2244 | -4.663 0.578 | 120.000 | 1.4985 | 1.6340 | -0.1354 | -8.288 CPU times: user 54.6 s, sys: 2.53 s, total: 57.1 s Wall time: 57.1 s 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") 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 \| team@tpq.io \| http://twitter.com/dyjh **Quant Platform** \| http://quant-platform.com **Derivatives Analytics with Python (Wiley Finance)** \| http://derivatives-analytics-with-python.com **Python for Finance (O'Reilly)** \| http://python-for-finance.com