diff --git a/tf_quant_finance/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..5f62270a2 Binary files /dev/null and b/tf_quant_finance/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..de4a2d147 Binary files /dev/null and b/tf_quant_finance/math/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/__pycache__/piecewise.cpython-37.pyc b/tf_quant_finance/math/__pycache__/piecewise.cpython-37.pyc new file mode 100644 index 000000000..f801b8315 Binary files /dev/null and b/tf_quant_finance/math/__pycache__/piecewise.cpython-37.pyc differ diff --git a/tf_quant_finance/math/interpolation/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/interpolation/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..1ad46e498 Binary files /dev/null and b/tf_quant_finance/math/interpolation/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/interpolation/cubic/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/interpolation/cubic/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..ba78ce5db Binary files /dev/null and b/tf_quant_finance/math/interpolation/cubic/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/interpolation/cubic/__pycache__/cubic_interpolation.cpython-37.pyc b/tf_quant_finance/math/interpolation/cubic/__pycache__/cubic_interpolation.cpython-37.pyc new file mode 100644 index 000000000..48782b6ce Binary files /dev/null and b/tf_quant_finance/math/interpolation/cubic/__pycache__/cubic_interpolation.cpython-37.pyc differ diff --git a/tf_quant_finance/math/interpolation/linear/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/interpolation/linear/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..d3a269127 Binary files /dev/null and b/tf_quant_finance/math/interpolation/linear/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/interpolation/linear/__pycache__/linear_interpolation.cpython-37.pyc b/tf_quant_finance/math/interpolation/linear/__pycache__/linear_interpolation.cpython-37.pyc new file mode 100644 index 000000000..f55f28b93 Binary files /dev/null and b/tf_quant_finance/math/interpolation/linear/__pycache__/linear_interpolation.cpython-37.pyc differ diff --git a/tf_quant_finance/math/optimizer/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/optimizer/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..2423a3f3e Binary files /dev/null and b/tf_quant_finance/math/optimizer/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/optimizer/__pycache__/conjugate_gradient.cpython-37.pyc b/tf_quant_finance/math/optimizer/__pycache__/conjugate_gradient.cpython-37.pyc new file mode 100644 index 000000000..3de122d0b Binary files /dev/null and b/tf_quant_finance/math/optimizer/__pycache__/conjugate_gradient.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/pde/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..1eb6c8eb5 Binary files /dev/null and b/tf_quant_finance/math/pde/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/__pycache__/grid_stepper.cpython-37.pyc b/tf_quant_finance/math/pde/__pycache__/grid_stepper.cpython-37.pyc new file mode 100644 index 000000000..c572e7b78 Binary files /dev/null and b/tf_quant_finance/math/pde/__pycache__/grid_stepper.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/__pycache__/pde_kernels.cpython-37.pyc b/tf_quant_finance/math/pde/__pycache__/pde_kernels.cpython-37.pyc new file mode 100644 index 000000000..144dadbf4 Binary files /dev/null and b/tf_quant_finance/math/pde/__pycache__/pde_kernels.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/grids/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/pde/grids/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..3c9b1accf Binary files /dev/null and b/tf_quant_finance/math/pde/grids/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/grids/__pycache__/grids_impl.cpython-37.pyc b/tf_quant_finance/math/pde/grids/__pycache__/grids_impl.cpython-37.pyc new file mode 100644 index 000000000..f610b55eb Binary files /dev/null and b/tf_quant_finance/math/pde/grids/__pycache__/grids_impl.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/internal/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/pde/internal/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..6dc08485a Binary files /dev/null and b/tf_quant_finance/math/pde/internal/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/internal/__pycache__/pde_time_marching_schemes_internal.cpython-37.pyc b/tf_quant_finance/math/pde/internal/__pycache__/pde_time_marching_schemes_internal.cpython-37.pyc new file mode 100644 index 000000000..81951a5b0 Binary files /dev/null and b/tf_quant_finance/math/pde/internal/__pycache__/pde_time_marching_schemes_internal.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..68a5cc6d3 Binary files /dev/null and b/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/pde_time_marching_scheme.cpython-37.pyc b/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/pde_time_marching_scheme.cpython-37.pyc new file mode 100644 index 000000000..c1ea09263 Binary files /dev/null and b/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/pde_time_marching_scheme.cpython-37.pyc differ diff --git a/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/pde_time_marching_schemes.cpython-37.pyc b/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/pde_time_marching_schemes.cpython-37.pyc new file mode 100644 index 000000000..534b0f97d Binary files /dev/null and b/tf_quant_finance/math/pde/time_marching_schemes/__pycache__/pde_time_marching_schemes.cpython-37.pyc differ diff --git a/tf_quant_finance/math/random_ops/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/random_ops/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..76bbadbd0 Binary files /dev/null and b/tf_quant_finance/math/random_ops/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/random_ops/__pycache__/stateless.cpython-37.pyc b/tf_quant_finance/math/random_ops/__pycache__/stateless.cpython-37.pyc new file mode 100644 index 000000000..dda429086 Binary files /dev/null and b/tf_quant_finance/math/random_ops/__pycache__/stateless.cpython-37.pyc differ diff --git a/tf_quant_finance/math/random_ops/halton/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/random_ops/halton/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..5841059dc Binary files /dev/null and b/tf_quant_finance/math/random_ops/halton/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/random_ops/halton/__pycache__/halton_impl.cpython-37.pyc b/tf_quant_finance/math/random_ops/halton/__pycache__/halton_impl.cpython-37.pyc new file mode 100644 index 000000000..96382d72c Binary files /dev/null and b/tf_quant_finance/math/random_ops/halton/__pycache__/halton_impl.cpython-37.pyc differ diff --git a/tf_quant_finance/math/random_ops/sobol/__pycache__/__init__.cpython-37.pyc b/tf_quant_finance/math/random_ops/sobol/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 000000000..9fa55d3b4 Binary files /dev/null and b/tf_quant_finance/math/random_ops/sobol/__pycache__/__init__.cpython-37.pyc differ diff --git a/tf_quant_finance/math/random_ops/sobol/__pycache__/sobol_impl.cpython-37.pyc b/tf_quant_finance/math/random_ops/sobol/__pycache__/sobol_impl.cpython-37.pyc new file mode 100644 index 000000000..f8d6d9eb2 Binary files /dev/null and b/tf_quant_finance/math/random_ops/sobol/__pycache__/sobol_impl.cpython-37.pyc differ diff --git a/tf_quant_finance/volatility/bachelier_tf.py b/tf_quant_finance/volatility/bachelier_tf.py new file mode 100644 index 000000000..e7a9d7af7 --- /dev/null +++ b/tf_quant_finance/volatility/bachelier_tf.py @@ -0,0 +1,252 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import tensorflow.compat.v2 as tf +import numpy as np + + +def _ncdf(x): + return (tf.math.erf(x / _SQRT_2) + 1.) / 2.0 + +def _npdf(x): + return tf.math.exp(-x**2/2)/_SQRT_2/_SQRT_pi + +_SQRT_2 = tf.math.sqrt(tf.constant(2.0,dtype=tf.float64)) #1.4142135623730951 +_SQRT_pi = tf.math.sqrt(tf.constant(np.pi,dtype=tf.float64)) + +# straight fwd implementation of the Bachelier pricing +# there is a version with just one call to exp !! + + +def bachelier_option_price(spots, + strikes, + volatilities, + expiries, + discount_rates = None, + discount_factors = None, + is_call_options=None, + dtype = None, + name = None): + """ computes the Bachelier price for a batch of European options. + We assume a standard Brownian motion of the form + dS = r dt + sigma dW + for the underlying + + ## References: + [1] Kienitz, J. "Interest Rate Derivatives Explained I", Palgrave McMillan (2014) p.119 + Link: https://www.palgrave.com/gp/book/9781137360069 + [2] Terakado, Satoshi: On the Option Pricing Formula Based on the Bachelier Model + Link: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3428994 + + #### Examples + + spots = np.array([0.03, 0.02]) + strikes = np.array([.02, .02]) + volatilities = np.array([.004, .005]) + expiries = 2.0 + discount_rates = [0.02, 0.01] + computed_prices = bachelier_option_price( + spots, + strikes, + volatilities, + expiries, + dtype=tf.float64) + # Expected print output of computed prices: + # + + Args: + + spots : A real `Tensor` of any shape. The current spot prices to + expiry. + strikes : A real `Tensor` of the same shape and dtype as `spots`. The + strikes of the options to be priced. + volatilities : A real `Tensor` of same shape and dtype as `spots`. The + volatility to expiry. + expiries : A real `Tensor` of same shape and dtype as `spots`. + discount_rates : rates from which discount factor via + discount factor = exp(-discount rate * T) are calculated + discountr_factors : A real 'Tensor' of same shape and dtype as 'spots' The + discounting factor; discount_rates = -log(discount factor) * expiries + is_call_options : A boolean `Tensor` of a shape compatible with + `volatilities`. Indicates whether the option is a call (if True) or a put + (if False). If not supplied, call options are assumed. + dtype: supplied dtype but converted to tf.float64 + name: name of the function + + Returns + + option_prices: A `Tensor` of the same shape as `spots`. The Bachelier + price of the options. + + + + """ + with tf.compat.v1.name_scope( + name, + default_name='bachelier_option_price', + values=[ + spots, strikes, volatilities, expiries, discount_rates, + discount_factors, is_call_options + ]): + + spots = tf.convert_to_tensor(spots, dtype=tf.float64, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + volatilities = tf.convert_to_tensor(volatilities, tf.float64, name='volatilities') + expiries = tf.convert_to_tensor(expiries, tf.float64, name='expiries') + if (discount_rates != None and discount_factors != None): + raise ValueError('Either discount rates or discount factors have to be used.') + + if (discount_rates != None and discount_factors == None): + rates = tf.convert_to_tensor(discount_rates, tf.float64, name='rates') + df = tf.math.exp(-rates * expiries) + elif (discount_factors != None and discount_rates == None): + rates = -tf.math.log(tf.convert_to_tensor(discount_rates, tf.float64, name='rates')) / expiries + df = discount_factors + else: + rates = 0.0 + df = tf.convert_to_tensor(rates, dtype=tf.float64, name='discount_rates') + + + z = tf.zeros_like(strikes) + + #normal = tfp.distributions.Normal( + # loc=tf.zeros([], dtype=spots.dtype), scale=1) + + vt = volatilities * tf.math.sqrt(expiries) + + z = tf.where(rates == 0., (spots - strikes)/vt, + (spots - strikes * df) / (volatilities + * tf.math.sqrt(0.5 * (1.-tf.math.exp(-2. * rates*expiries)) / rates))) + + n1 = _ncdf(z) + n2 = _npdf(z) + calls = tf.where(rates==0., (spots - strikes) * n1 + vt * n2, + (spots - strikes * df) * n1 + + volatilities * tf.math.sqrt(0.5 * (1 - tf.math.exp(-2 * rates * expiries)) / rates)) + + + if is_call_options is None: + return calls + + puts = calls - spots + strikes * tf.math.exp(-rates * expiries) + + return tf.where(is_call_options, calls, puts) + + + +def dawson_option_price(forwards, + strikes, + volatilities, + expiries, + discount_rates = None, + discount_factors = None, + is_call_options=None, + dtype = None, + name = None): + + """Computes the Bachelier price for a batch of European options. + We assume a standard Brownian motion of the form + dS = r dt + sigma dW + for the underlying + + ## References: + [1] Dawson, P., Blake, D., Cairns, A. J. G. and Dowd, K.: Options on normal under- + lyings, CRIS Discussion Paper Series – 2007.VII, 2007. + + #### Examples + spots = np.array([0.03, 0.02]) + strikes = np.array([.02, .02]) + volatilities = np.array([.004, .005]) + expiries = 2.0 + expiries = 1.0 + computed_prices = dawson_option_price( + forwards, + strikes, + volatilities, + expiries, + dtype=tf.float64) + # Expected print output of computed prices: + # + + Args: + + forwards: A real `Tensor` of any shape. The current forward prices to + expiry. + strikes: A real `Tensor` of the same shape and dtype as `forwards`. The + strikes of the options to be priced. + volatilities: A real `Tensor` of same shape and dtype as `forwards`. The + volatility to expiry. + expiries : A real `Tensor` of same shape and dtype as `spots`. + discount_rates : rates from which discount factor via + discount factor = exp(-discount rate * T) are calculated + discount_factors : A real 'Tensor' of same shape and dtype as 'spots' The + discounting factor; discount_rates = -log(discount factor) * expiries + is_call_options : A boolean `Tensor` of a shape compatible with + `volatilities`. Indicates whether the option is a call (if True) or a put + (if False). If not supplied, call options are assumed. + dtype: supplied dtype but converted to tf.float64 + name: name of the function + + Returns: + option_prices: A `Tensor` of the same shape as `forwards`. The Bachelier + price of the options. + + + + """ + with tf.compat.v1.name_scope( + name, + default_name='dawson_option_price', + values=[ + forwards, strikes, volatilities, expiries, discount_factors, + discount_rates, is_call_options + ]): + + forwards = tf.convert_to_tensor(forwards, dtype=tf.float64, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + volatilities = tf.convert_to_tensor(volatilities, dtype=tf.float64, name='volatilities') + expiries = tf.convert_to_tensor(expiries, dtype=tf.float64, name='expiries') + + # check if discount rates or discount factor version is used + if (discount_rates != None and discount_factors != None): + raise ValueError('Either discount rates or discount factors have to be used.') + + if (discount_rates != None and discount_factors == None): + discount_factors = tf.math.exp(-tf.convert_to_tensor(discount_rates, tf.float64, name='discount factors')*expiries) + else: + if (discount_factors == None and discount_rates == None): + discount_factors = 1.0 + discount_factors = tf.convert_to_tensor(discount_factors, dtype=tf.float64, name='discount_factors') + + vt = volatilities * tf.math.sqrt(expiries) + + z = (forwards - strikes) / vt + + # calculate constituents of Bachelier formula + n1 = _ncdf(z) + n2 = _npdf(z) + undiscounted_calls = (forwards - strikes) * n1 + vt * n2 # Bachelier option price + + # check if calls or puts need to be considered + if is_call_options is None: + return discount_factors * undiscounted_calls + undiscounted_forward = forwards - strikes + undiscounted_puts = undiscounted_calls - undiscounted_forward + + # return call, resp. put prices + return discount_factors * tf.where(is_call_options, undiscounted_calls, + undiscounted_puts) + + diff --git a/tf_quant_finance/volatility/sabr_approx_hagan_tf_test.py b/tf_quant_finance/volatility/sabr_approx_hagan_tf_test.py new file mode 100644 index 000000000..58a1df8de --- /dev/null +++ b/tf_quant_finance/volatility/sabr_approx_hagan_tf_test.py @@ -0,0 +1,137 @@ +""" +Created on Fri Nov 22 15:22:13 2019 + +# Copyright 2020 Joerg Kienitz + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +@author: Joerg Kienitz +""" + +import sabr_approx_tf as sabr +import bachelier_tf as vptf +import volbachelier_tf as bvtf + +import numpy as np +import matplotlib.pyplot as plt + +# SABR parameters +# SABR parameters +f = 1.0 #0.00434015 +alpha_org = 0.16575423 +beta_org = .6#0.7#0.16415365 +nu_org = 0.2632859 +rho_org = -0.32978014 +T = 5 + +displacement_org = 0. #0.005 +kmin = -displacement_org +kmax = 10 +kval = np.arange(kmin, kmax, 0.01) +kval[0] = (kval[0] + kval[1])/2 +vol = np.zeros(len(kval)) + +alpha_vec = [0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.175, 0.2, 0.25, 0.3, 0.5, 0.75, 1., 1.5] +beta_vec = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +rho_vec = [-1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +nu_vec = [0.001, 0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.5, 0.75, 1.0, 1.5] +displacement_vec = [0.0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.05] + +print('alpha varies') +for alpha in alpha_vec: + print('alpha: ', alpha) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha, beta_org, rho_org,nu_org) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(alpha) + label2 = 'iv ' + str(alpha) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('alpha varies') +plt.legend() +plt.show() + + +print('beta varies') +for beta in beta_vec: + print('parameters: ', beta) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta, rho_org,nu_org) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(beta) + label2 = 'iv ' + str(beta) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('beta varies') +plt.legend() +plt.show() + +print('rho varies') +for rho in rho_vec: + print('parameters: ', rho) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho,nu_org) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(rho) + label2 = 'iv ' + str(rho) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('rho varies') +plt.legend() +plt.show() + +print('nu varies') +for nu in nu_vec: + print('parameters: ', nu) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho_org,nu) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(nu) + label2 = 'iv ' + str(nu) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('nu varies') +plt.legend() +plt.show() + +print('displacement varies') +for displacement in displacement_vec: + print('parameters: ', displacement) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho_org,nu_org) + label1 = 'approx ' + str(displacement) + plt.plot(kval,yval,label= label1) +plt.title('displacement varies') +plt.legend() +plt.show() + +# different approximation techniques for SABR and Mean Reverting SABR +kappa = 0.5 +cap = 3. + +yval1 = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho_org,nu_org) +yval2 = sabr.volsabr_mr_n_tf(f,kval,T,displacement_org, alpha_org, beta_org, rho_org, nu_org, kappa) +yval3 = sabr.volsabr_h_n_cap_tf(f,kval,T,displacement_org, alpha_org, beta_org, rho_org, nu_org, cap) + +label1 = 'Hagan approx ' +label2 = 'MR SABR approx ' +label3 = 'Capped SABR approx ' + +plt.plot(kval,yval1,label= label1) +plt.plot(kval,yval2,label= label2) +plt.plot(kval,yval3,label= label3) +plt.title('different SABR approximation') +plt.legend() +plt.show() + + + diff --git a/tf_quant_finance/volatility/sabr_approx_tf.py b/tf_quant_finance/volatility/sabr_approx_tf.py new file mode 100644 index 000000000..38fe13855 --- /dev/null +++ b/tf_quant_finance/volatility/sabr_approx_tf.py @@ -0,0 +1,375 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import tensorflow as tf + +# straight fwd implementation of the Bachelier pricing +# there is a version with just one call to exp !! + +def volsabr_h_n_tf( + forwards, + strikes, + expiries, + displacements, + alpha, + beta, + rho, + nu, + name = None): + """ + Standard Hagan SABR approximation for the Bachelier/Normal volatility + + We assume a SABR model dynamic given by + dS = r dt + alpha * (S+d)^beta sigma dW + dalpha = nu * alpha * dZ + = rho dt + S(0) = forward + alpha(0) = alpha + + ## References: + [1] Hagan, S. Patrick, Deep Kumar, Andrew Lesniewski and Diana E. Woodward, 2002. - “Managing Smile Risk”. Wilmott Magazine, September/October. + and on Researchgate + Link: https://www.researchgate.net/publication/235622441_Managing_Smile_Risk + [2] Patrick S. Hagan Deep Kumar Andrew S. Lesniewski Diana E. Woodward, " Universal Smiles" + Link: https://onlinelibrary.wiley.com/doi/abs/10.1002/wilm.10523 + [3] Kienitz, J. "Interest Rate Derivatives Explained I", Palgrave McMillan (2014) p.119 + Link: https://www.palgrave.com/gp/book/9781137360069 + [4] Kienitz, J. "Interest Rate Derivatives Explained II", Palgrave McMillan (2017), p. + Link: https://www.palgrave.com/gp/book/9781137360182 + + parameters: + forwards - forward rates + strikes - strikes as array + expiries - expiry/maturity + displacements - displacement + alpha - SABR parameter (initial vol) + beta - SABR parameter (CEV coefficient) + rho - SABR parameter (correlation) + nu - SABR parameter (vol of vol) + + For SABR we can always use f=1 and apply the scaling: + if f = f0 + we knew -> k/f + alphanew -> f**(beta-1) * alpha + ivol1 = sabrapprox.volsabr_h_n(f, kval, T, displacement, alpha, beta, rho, nu) + price1 = vp.vanilla_n(f,kval,implVolApprox1,0,0,T,1) + + ivol2 = f * sabrapprox.volsabr_h_n(1, knew, T, displacement, alphanew, beta, rho, nu) + ivol3 = sabrapprox.volsabr_h_n(1, knew, T, displacement, alphanew, beta, rho, nu) + price2 = vp.vanilla_n(f,kval,ivol2,0,0,T,1) + price3 = f * vp.vanilla_n(1,kvalnew,ivol3,0,0,T,1) + + price1 = price2 = price3 + ivol1 = ivol2 = f * ivol3 + + Returns + + vol_sabr_approx: A `Tensor` of the same shape as `forwards`. The implied + Bachelier volatility approximation for the SABR model. + + """ + with tf.compat.v1.name_scope( + name, + default_name='sabr_implied_vol_hagan', + values=[ + forwards, strikes, expiries, displacements, alpha, beta, rho, nu + ]): + # conversion maybe outside function!!!! + forwards = tf.convert_to_tensor(forwards, dtype=tf.float64, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + expiries = tf.convert_to_tensor(expiries, dtype=tf.float64, name='expiries') + displacements = tf.convert_to_tensor(displacements, dtype=tf.float64, name='displacement') + alpha = tf.convert_to_tensor(alpha, dtype=tf.float64, name='alpha') + beta = tf.convert_to_tensor(beta, dtype=tf.float64, name='beta') + rho = tf.convert_to_tensor(rho, dtype=tf.float64, name='rho') + nu = tf.convert_to_tensor(nu, dtype=tf.float64, name='nu') + + + # identify ATM and non ATM strikes + eps = tf.constant(1.0e-06, dtype = tf.float64) # small number + index_natm = (tf.math.abs(forwards-strikes) > eps) # itm/otm strikes + + # case rho = 1 may cause divide by zero problems + # case rho = 1 may cause divide by zero problems + # for rho == 1 + # xik = log(2) * log(x+1) -> x > -1 + # for rho == -1 + # xik = log(2) * log(x-1) -> x > 1 + # this means only certain strikes are feasible for calculating implied vols! + # to this end we only consider the SABR model for correlation values + # between rho = -0.99 and rho = 0.99 + rho = tf.where(rho == 1., 0.999, tf.where(rho == -1., -0.999, rho)) + + betam = tf.constant(1.0 - beta) # often used + + fa = forwards+displacements # account for displacement for forward + ka = strikes+displacements # account for displacement for strikes + + # different cases due to normal, cev, log-normal + if 0. < beta and beta < 1.: # case of true CEV SABR + gk = tf.zeros_like(strikes) + gk = tf.where(index_natm,(beta**2 - 2. * beta) / 24. + * fa**(-betam) * ka**(-betam)*alpha**2, gk) + xik = nu / alpha * (fa**betam - ka**betam) / betam + xxik = tf.math.log((tf.math.sqrt(1.0 - 2.0 * rho * xik + xik**2) - rho + xik) + / (1 - rho)) / nu + vol_sabr_approx = tf.where(index_natm, + (fa - ka) / xxik * (1. + (gk + 0.25 * rho * nu * alpha * beta * fa**(0.5 * (beta-1.)) * ka**(0.5 * (beta - 1.)) + + (2. - 3. * rho**2) / 24. * nu**2) * expiries), + alpha * fa**beta * (1 + (beta * (beta-2.) * alpha**2 / 24. / fa**(2. * betam) + + 0.25 * rho * nu * alpha * beta / fa**(betam) + + (2. - 3. * rho**2) / 24. * nu**2) * expiries) ) + + elif beta == 0.: # case of a Gaussian SV model (normal SABR) + xik = nu / alpha * (fa - ka) + xxik = tf.math.log((tf.math.sqrt(1.0 - 2.0 * rho * xik + xik**2) - rho + xik) / (1 - rho)) + vol_sabr_approx = tf.where(index_natm, alpha * xik / xxik * (1. + (1. / 24. * (2. - 3. * rho**2) * nu**2) * expiries), + alpha * (1. + (2. - 3. * rho**2) / 24. * nu**2 * expiries)) + + else: # case of log-normal SV model (log-normal SABR) + gk = - 1. / 24. * alpha**2 + xik = nu / alpha * tf.math.log(fa / ka) + sum2 = 0.25 * rho * nu * alpha + xxik = tf.math.log((tf.math.sqrt(1. - 2. * rho * xik + xik**2) - rho + xik) + / (1. - rho)) / nu + vol_sabr_approx = tf.where(index_natm, + (fa - ka) / xxik * (1 + (gk + sum2 + 1. / 24. * (2. - 3.*rho**2) * nu**2) * expiries), + alpha * fa * (1.+ (gk + sum2 + (2. - 3. * rho**2) / 24. * nu**2) * expiries)) + return vol_sabr_approx + +def volsabr_mr_n_tf(forwards, + strikes, + expiries, + displacements, + alpha, + beta, + rho, + nu, + kappa, + name = None): + """ computes the Bachelier implied volatility batch of European options. + We assume a SABR model dynamic given by + dS = alpha * (S+d)^beta sigma dW + dalpha = kappa (1-alpha) dt + nu * alpha * dZ + = rho dt + S(0) = forward + alpha(0) = alpha + + ## References: + [1] Kienitz, J. "Interest Rate Derivatives Explained I", Plagrave McMillan (2014) p.119 + Link: https://www.palgrave.com/gp/book/9781137360069 + [2] Terakado, Satoshi: On the Option Pricing Formula Based on the Bachelier Model + Link: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3428994 + + + Args: + + forwards : A real `Tensor` of any shape. The current spot prices to + expiry. + strikes : A real `Tensor` of the same shape and dtype as `forwards`. The + strikes of the options to be priced. + expiries : A real `Tensor` of the same shape and dtype as `forwards`. The + expiries of the options to be priced. + expiries : A real `Tensor` of same shape and dtype as `spots`. + displacements:A real `Tensor` of the same shape and dtype as `forwards`. The + displacements for the forwards. + alpha : A real `Tensor` of the same shape and dtype as `forwards`. Initial + volatility of the stochastic volatility component. + beta : A real `Tensor` of the same shape and dtype as `forwards`. CEV + parameter of the SABR model. + rho : A real `Tensor` of the same shape and dtype as `forwards`. correlation + of the Brownian motions driving the forward and the volatility. + nu : A real `Tensor` of the same shape and dtype as `forwards`. volatiltiy + of volatility. + kappa: A real `Tensor` of the same shape and dtype as `forwards`. mean + reversion that is applied. + name: name of the function + + Returns + + vol_sabr_mr: A `Tensor` of the same shape as `forwards`. The Bachelier + implied volatiltiy for the approximation in the mean reverting SABR model. + + """ + with tf.compat.v1.name_scope( + name, + default_name='sabr_implied_vol_hagan', + values=[ + forwards, strikes, expiries, displacements, alpha, beta, rho, nu + ]): + + forwards = tf.convert_to_tensor(forwards, dtype=tf.float64, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + expiries = tf.convert_to_tensor(expiries, dtype=tf.float64, name='expiries') + displacements = tf.convert_to_tensor(displacements, dtype=tf.float64, name='displacement') + alpha = tf.convert_to_tensor(alpha, dtype=tf.float64, name='alpha') + beta = tf.convert_to_tensor(beta, dtype=tf.float64, name='beta') + rho = tf.convert_to_tensor(rho, dtype=tf.float64, name='rho') + nu = tf.convert_to_tensor(nu, dtype=tf.float64, name='nu') + kappa = tf.convert_to_tensor(kappa, dtype=tf.float64, name='kappa') + + + if kappa > 0.: + bbar = 2. * rho * nu / alpha * (kappa * expiries - 1. + + tf.math.exp(-kappa * expiries)) / (kappa**2 * expiries**2) + cbar = 1.5 * nu ** 2 / alpha**2. * (1. + rho**2) * (1. + 2. * kappa * expiries - (2. - tf.math.exp(-kappa * expiries))**2) \ + / (kappa**3 * expiries**3) + 12. * rho**2 * nu**2 / alpha**2 * (kappa**2 * expiries**2 * tf.math.exp(-kappa * expiries) + - (1. - tf.math.exp(-kappa * expiries))**2) / (kappa**4 * expiries**4) + Gstar = (-0.5 * cbar + nu**2 / alpha**2 * (2. * kappa * expiries - 1. + tf.math.exp(-2. * kappa * expiries)) \ + / (4. * kappa**2 * expiries**2)) * alpha**2 * expiries + + # std case + astd = alpha * tf.math.exp(0.5 * Gstar) + rhostd = bbar / tf.math.sqrt(cbar) + nustd = alpha * tf.math.sqrt(cbar) + vol_sabr_mr = volsabr_h_n_tf(forwards, strikes, expiries, displacements, astd, beta, rhostd, nustd) + else: + vol_sabr_mr = volsabr_h_n_tf(forwards, strikes, expiries, displacements, alpha, beta, rho, nu) + + return vol_sabr_mr + + +def volsabr_h_n_cap_tf( + forwards, + strikes, + expiries, + displacements, + alpha, + beta, + rho, + nu, + cap, + name = None): + """ + Standard Hagan SABR approximation for the Bachelier/Normal volatility with + capping the volatility to reduce the right wing of the smile + + We assume a SABR model dynamic given by + dS = alpha * (S+d)^beta * sigma * dW + dalpha = nu * alpha * dZ + = rho dt + S(0) = forward + alpha(0) = alpha + + ## References: + [1] Patrick S. Hagan Deep Kumar Andrew S. Lesniewski Diana E. Woodward, " Universal Smiles" + Link: https://onlinelibrary.wiley.com/doi/abs/10.1002/wilm.10523 + [2] Kienitz, J. "Interest Rate Derivatives Explained II", Palgrave McMillan (2017) + https://www.palgrave.com/gp/book/9781137360182 + + parameters: + forwards - forward rates + strikes - strikes as array + expiries - expiry/maturity + displacements - displacement + alpha - SABR parameter (initial vol) + beta - SABR parameter (CEV coefficient) + rho - SABR parameter (correlation) + nu - SABR parameter (vol of vol) + cap - parameter for capping volatility + + Returns + + sabr_vol_capped : A `Tensor` of the same shape as `forwards'. The Bachelier + implied volatility for the capped SABR model. + """ + with tf.compat.v1.name_scope( + name, + default_name='sabr_implied_vol_hagan', + values=[ + forwards, strikes, expiries, displacements, alpha, beta, rho, nu + ]): + # conversion maybe outside function!!!! + forwards = tf.convert_to_tensor(forwards, dtype=tf.float64, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + expiries = tf.convert_to_tensor(expiries, dtype=tf.float64, name='expiries') + displacements = tf.convert_to_tensor(displacements, dtype=tf.float64, name='displacement') + alpha = tf.convert_to_tensor(alpha, dtype=tf.float64, name='alpha') + beta = tf.convert_to_tensor(beta, dtype=tf.float64, name='beta') + rho = tf.convert_to_tensor(rho, dtype=tf.float64, name='rho') + nu = tf.convert_to_tensor(nu, dtype=tf.float64, name='nu') + cap = tf.convert_to_tensor(cap, dtype=tf.float64, name='cap') + + eps = 1.0e-06 + index_atm = tf.math.abs(forwards - strikes) < eps + + sabr_vol_capped = tf.zeros_like(strikes) + fa = forwards + displacements + ka = strikes + displacements + + betam = 1. - beta + # atm but different cases due to normal, cev, log-normal + if 0. < beta and beta < 1.: # case of true CEV SABR + volatm = alpha * fa**beta * (1 + (beta * (beta - 2.) * alpha**2 / 24. / fa**(2. * betam) + + 0.25 * rho * nu * alpha * beta / fa**(betam) + + (2. - 3. * rho**2) / 24. * nu**2) * expiries) + + elif beta == 0.: # case of a Gaussian SV model (normal SABR) + volatm = alpha * (1. + (2. - 3. * rho**2) / 24. * nu**2 * expiries) + else: # case of log-normal SV model (log-normal SABR) + volatm = alpha * fa * (1.+ (-1. / 24. * alpha**2 + 0.25 * rho * nu * alpha + (2. - 3. * rho**2) / 24. * nu**2) * expiries) + + + if beta == 1: + beta = 0.999 + + rho = tf.where(rho == 1, 0.999, tf.where(rho == -1, -0.999,rho)) + + # the cap can only be applied if the term under sqrt is positive + term = tf.math.sqrt(tf.math.maximum(cap**2 - 1. + rho**2,0.)) + + xip = -rho + term + xim = -rho - term + + Yp = -tf.math.log((cap - term) / (1. - rho)) + Ym = -tf.math.log((cap + term) / (1. - rho)) + + # here we need ATM consideration + ic = ((ka)**(1. - beta) - (fa)**(1. - beta)) / (1. - beta) + + f0 = 0.5 * ((forwards + strikes) + displacements) # 2* displace? + + gamma = beta * f0**(beta-1.) + + Delta0 = (beta**2 - 2. * beta) * f0 ** (2. * beta-2.) + + xi = nu / alpha * ic + + Yxi = tf.where(xi > xip, + Yp + (xi - xip) / cap, + tf.where(xi < xim, + Ym + (xi -xim) / cap, + -tf.math.log((tf.math.sqrt(1. + 2. * rho * xi + xi ** 2) - rho - xi) / (1 - rho)) + ) + ) + + sK0 = tf.where(xi > xip, + nu**2 / (8. * alpha**2 * Yxi) * (-Yp + 3. * (-rho * cap + term) / cap) \ + + Delta0 / (16. * Yxi) * (2. * cap**2 * (Yxi - Yp) + (1 - rho**2) * Yp + cap * term - rho), + tf.where(xi < xim, + nu**2 / (8 * alpha**2 * Yxi) * (-Ym - 3. * (rho * cap + term) / cap) + Delta0 / (16. * Yxi) * (2. * cap**2 * (Yxi - Ym) + (1. - rho**2) * Ym - cap * term - rho), + nu**2 / (8. * alpha**2 * Yxi) * (-Yxi + 3. * (xi + rho - rho * tf.math.sqrt(1. + 2. * rho * xi + xi**2)) / tf.math.sqrt(1. + 2. * rho * xi + xi**2)) \ + + Delta0 / (16. * Yxi) * ((1. - rho**2) * Yxi + (xi + rho) * tf.math.sqrt(1. + 2. * rho * xi + xi**2) - rho) + ) + ) + + theta = -0.5 * rho * nu * alpha * gamma - 2. / 3. * alpha ** 2 * sK0 + highorder = tf.where(theta <0., + tf.math.sqrt(1. - theta * expiries), + 1. / tf.math.sqrt(1. + theta * expiries) + ) + sabr_vol_capped = tf.where(index_atm, volatm, nu * (strikes - forwards) / Yxi * highorder) + return sabr_vol_capped + \ No newline at end of file diff --git a/tf_quant_finance/volatility/test_sabr_hagan_tf.py b/tf_quant_finance/volatility/test_sabr_hagan_tf.py new file mode 100644 index 000000000..2ce3d8737 --- /dev/null +++ b/tf_quant_finance/volatility/test_sabr_hagan_tf.py @@ -0,0 +1,87 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# f = 1.0 is sufficient due to scaling property + +import volsabr_h_n_tf as sabr + +import matplotlib.pyplot as plt +#import numpy as np +import autograd.numpy as np +# SABR parameters +f = 1.0 #0.00434015 +alpha_org = 0.16575423 +beta_org = .6#0.7#0.16415365 +nu_org = 0.0632859 +rho_org = -0.32978014 +T = 5 + +displacement = 0 #0.005 +kmin = -displacement +kmax = 10 +kval = np.arange(kmin, kmax, 0.01) + +alpha_vec = [0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.175, 0.2, 0.25, 0.3, 0.5, 0.75, 1., 1.5] +beta_vec = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +rho_vec = [-1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +nu_vec = [0.001, 0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.5, 0.75, 1.0, 1.5] +displacement_vec = [0.0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.05] + +print('alpha varies') +for alpha in alpha_vec: + print('alpha:', alpha) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha, beta_org, rho_org, nu_org) + plt.plot(kval,yval,label=alpha) +plt.title('alpha varies') +plt.legend() +plt.show() + + +print('beta varies') +for beta in beta_vec: + print('beta:', beta) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta, rho_org, nu_org) + plt.plot(kval,yval,label=beta) +plt.title('beta varies') +plt.legend() +plt.show() + +print('rho varies') +for rho in rho_vec: + print('rho:', rho) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta_org, rho, nu_org) + plt.plot(kval,yval,label=rho) +plt.title('rho varies') +plt.legend() +plt.show() + +print('nu varies') +for nu in nu_vec: + print('nu:' nu) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha, beta, rho_org,nu) + plt.plot(kval,yval,label=nu) +plt.title('nu varies') +plt.legend() +plt.show() + +print('displacement varies') +for displacement in displacement_vec: + print('displacement:', displacement) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta_org, rho_org,nu_org) + plt.plot(kval,yval,label=displacement) +plt.title('displacement varies') +plt.legend() +plt.show() + diff --git a/tf_quant_finance/volatility/test_sabr_hagan_tf_1.py b/tf_quant_finance/volatility/test_sabr_hagan_tf_1.py new file mode 100644 index 000000000..84da792b1 --- /dev/null +++ b/tf_quant_finance/volatility/test_sabr_hagan_tf_1.py @@ -0,0 +1,90 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# f = 1.0 is sufficient due to scaling property + +import sabr_approx_tf as sabr + +import matplotlib.pyplot as plt +#import numpy as np +import autograd.numpy as np +# SABR parameters +f = 1.0 #0.00434015 +alpha_org = 0.16575423 +beta_org = .6#0.7#0.16415365 +nu_org = 0.0632859 +rho_org = -0.32978014 +T = 5 + +displacement = 0 #0.005 +kmin = -displacement +kmax = 10 +kval = np.arange(kmin, kmax, 0.01) + +alpha_vec = [0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.175, 0.2, 0.25, 0.3, 0.5, 0.75, 1., 1.5] +beta_vec = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +rho_vec = [-1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +nu_vec = [0.001, 0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.5, 0.75, 1.0, 1.5] +displacement_vec = [0.0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.05] + +print('alpha varies') +for alpha in alpha_vec: + print('alpha:', alpha) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha, beta_org, rho_org, nu_org) + plt.plot(kval,yval,label=alpha) +plt.title('alpha varies') +plt.legend() +plt.show() + + +print('beta varies') +for beta in beta_vec: + print('beta:', beta) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta, rho_org, nu_org) + plt.plot(kval,yval,label=beta) +plt.title('beta varies') +plt.legend() +plt.show() + +print('rho varies') +for rho in rho_vec: + print('rho:', rho) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta_org, rho, nu_org) + plt.plot(kval,yval,label=rho) +plt.title('rho varies') +plt.legend() +plt.show() + +print('nu varies') +for nu in nu_vec: + print('nu:', nu) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha, beta, rho_org,nu) + plt.plot(kval,yval,label=nu) +plt.title('nu varies') +plt.legend() +plt.show() + +print('displacement varies') +for d in displacement_vec: + print('displacement:', d) + kmin = -d + kmax = 10 + kval = np.arange(kmin, kmax, 0.01) + yval = sabr.volsabr_h_n_tf(f, kval, T, d, alpha_org, beta_org, rho_org,nu_org) + plt.plot(kval,yval,label=d) +plt.title('displacement varies') +plt.legend() +plt.show() + diff --git a/tf_quant_finance/volatility/test_sabr_hagan_tf_2.py b/tf_quant_finance/volatility/test_sabr_hagan_tf_2.py new file mode 100644 index 000000000..343ba8ec5 --- /dev/null +++ b/tf_quant_finance/volatility/test_sabr_hagan_tf_2.py @@ -0,0 +1,132 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import sabr_approx_tf as sabr +import bachelier_tf as vptf +import volbachelier_tf as bvtf + +import numpy as np +import matplotlib.pyplot as plt + +# SABR parameters +# SABR parameters +f = 1.0 #0.00434015 +alpha_org = 0.16575423 +beta_org = .6#0.7#0.16415365 +nu_org = 0.2632859 +rho_org = -0.32978014 +T = 5 + +displacement_org = 0. #0.005 +kmin = -displacement_org +kmax = 10 +kval = np.arange(kmin, kmax, 0.01) +kval[0] = (kval[0] + kval[1])/2 +vol = np.zeros(len(kval)) + +alpha_vec = [0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.175, 0.2, 0.25, 0.3, 0.5, 0.75, 1., 1.5] +beta_vec = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +rho_vec = [-1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +nu_vec = [0.001, 0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.5, 0.75, 1.0, 1.5] +displacement_vec = [0.0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.05] + +print('alpha varies') +for alpha in alpha_vec: + print('alpha: ', alpha) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha, beta_org, rho_org,nu_org) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(alpha) + label2 = 'iv ' + str(alpha) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('alpha varies') +plt.legend() +plt.show() + + +print('beta varies') +for beta in beta_vec: + print('parameters: ', beta) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta, rho_org,nu_org) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(beta) + label2 = 'iv ' + str(beta) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('beta varies') +plt.legend() +plt.show() + +print('rho varies') +for rho in rho_vec: + print('parameters: ', rho) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho,nu_org) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(rho) + label2 = 'iv ' + str(rho) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('rho varies') +plt.legend() +plt.show() + +print('nu varies') +for nu in nu_vec: + print('parameters: ', nu) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho_org,nu) + cval = vptf.bachelier_option_price(f,kval,yval,T,0.) + yval1 = bvtf.volbachelier_tf(1, kval, f, T, cval) + label1 = 'approx ' + str(nu) + label2 = 'iv ' + str(nu) + plt.plot(kval,yval,label= label1) + plt.plot(kval, yval1, label=label2) +plt.title('nu varies') +plt.legend() +plt.show() + +print('displacement varies') +for displacement in displacement_vec: + print('parameters: ', displacement) + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho_org,nu_org) + label1 = 'approx ' + str(displacement) + plt.plot(kval,yval,label= label1) +plt.title('displacement varies') +plt.legend() +plt.show() + +# different approximation techniques for SABR and Mean Reverting SABR +kappa = 0.5 +cap = 3. + +yval1 = sabr.volsabr_h_n_tf(f, kval, T, displacement_org, alpha_org, beta_org, rho_org,nu_org) +yval2 = sabr.volsabr_mr_n_tf(f,kval,T,displacement_org, alpha_org, beta_org, rho_org, nu_org, kappa) +yval3 = sabr.volsabr_h_n_cap_tf(f,kval,T,displacement_org, alpha_org, beta_org, rho_org, nu_org, cap) + +label1 = 'Hagan approx ' +label2 = 'MR SABR approx ' +label3 = 'Capped SABR approx ' + +plt.plot(kval,yval1,label= label1) +plt.plot(kval,yval2,label= label2) +plt.plot(kval,yval3,label= label3) +plt.title('different SABR approximation') +plt.legend() +plt.show() + + + diff --git a/tf_quant_finance/volatility/test_volbachelier_tf.py b/tf_quant_finance/volatility/test_volbachelier_tf.py new file mode 100644 index 000000000..be5359df2 --- /dev/null +++ b/tf_quant_finance/volatility/test_volbachelier_tf.py @@ -0,0 +1,80 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# f = 1.0 is sufficient due to scaling property + +import sabr_approx_tf as sabr + +import matplotlib.pyplot as plt +import vanilla_n_tf as vp +import volbachelier_tf as vb +import autograd.numpy as np + +# SABR parameters +f = 1.0 #0.00434015 +alpha_org = 0.16575423 +beta_org = .6#0.7#0.16415365 +nu_org = 0.0632859 +rho_org = -0.32978014 +T = 5 + +displacement = 0 #0.005 +kmin = .25 +kmax = 5 +kval = np.arange(kmin, kmax, 0.01) + +alpha_vec = [0.1, 0.15, 0.175, 0.2, 0.25, 0.3, 0.5, 0.75] +beta_vec = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] +rho_vec = [-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] +nu_vec = [0.001, 0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.5, 0.75, 1.0, 1.5] + + +for alpha in alpha_vec: + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha, beta_org, rho_org, nu_org) + oval = vp.option_price(f, kval, yval, T, 0) + ival = vb.volbachelier_tf(1., kval, f, T, oval) + plt.plot(kval, yval - ival) +plt.title('Bachelier implied vol vs SABR 1') +plt.legend() +plt.show() + +for beta in beta_vec: + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta, rho_org, nu_org) + oval = vp.option_price(f, kval, yval, T, 0) + ival = vb.volbachelier_tf(1., kval, f, T, oval) + plt.plot(kval, yval - ival) +plt.title('Bachelier implied vol vs SABR 2') +plt.legend() +plt.show() + + +for rho in rho_vec: + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta_org, rho, nu_org) + oval = vp.option_price(f, kval, yval, T, 0) + ival = vb.volbachelier_tf(1., kval, f, T, oval) + plt.plot(kval, yval - ival) +plt.title('Bachelier implied vol vs SABR 3') +plt.legend() +plt.show() + + +for nu in nu_vec: + yval = sabr.volsabr_h_n_tf(f, kval, T, displacement, alpha_org, beta_org, rho_org, nu) + oval = vp.option_price(f, kval, yval, T, 0) + ival = vb.volbachelier_tf(1., kval, f, T, oval) + plt.plot(kval, yval - ival) +plt.title('Bachelier implied vol vs SABR 4') +plt.legend() +plt.show() \ No newline at end of file diff --git a/tf_quant_finance/volatility/vanilla_n_tf.py b/tf_quant_finance/volatility/vanilla_n_tf.py new file mode 100644 index 000000000..8868bf4f0 --- /dev/null +++ b/tf_quant_finance/volatility/vanilla_n_tf.py @@ -0,0 +1,228 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +import tensorflow as tf +import tensorflow_probability as tfp + +# straight fwd implementation of the Bachelier pricing +# there is a version with just one call to exp !! + + +def option_price(spots, + strikes, + volatilities, + expiries, + rates, + is_call_options=None, + dtype = None, + name = None): + """ Compute the Bachelier price for a batch of European options. + + ## References: + [1] Kienitz, J. "Interest Rate Derivatives Explained I", Plagrave McMillan (2014) p.119 + [2] https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3428994 + + Parameters + ---------- + spots : A real `Tensor` of any shape. The current spot prices to + expiry. + strikes : A real `Tensor` of the same shape and dtype as `spots`. The + strikes of the options to be priced. + volatilities : A real `Tensor` of same shape and dtype as `spots`. The + volatility to expiry. + expiries : A real `Tensor` of same shape and dtype as `spots`. The expiry + for each option. The units should be such that `expiry * volatility**2` is + dimensionless. + rates : A real `Tensor` of the same shape and dtype as `spots`. The + rates of the options to be priced. + is_call_options: A boolean `Tensor` of a shape compatible with `forwards`. + Indicates whether to compute the price of a call (if True) or a put (if + False). If not supplied, it is assumed that every element is a call. + dtype: Optional `tf.DType`. If supplied, the dtype to be used for conversion + of any supplied non-`Tensor` arguments to `Tensor`. + Default value: None which maps to the default dtype inferred by TensorFlow + (float32). + name: str. The name for the ops created by this function. + Default value: None which is mapped to the default name `option_price`. + + Returns + ------- + option_prices: A `Tensor` of the same shape as `spots`. The Bachelier + price of the options. + + #### Examples + ```python + spots = np.array([0.03, 0.02]) + strikes = np.array([.02, .02]) + volatilities = np.array([.004, .005]) + expiries = 2.0 + rates = [0.02, 0.01] + computed_prices = option_price( + spots, + strikes, + volatilities, + expiries, + rates, + dtype=tf.float64) + # Expected print output of computed prices: + # + ``` + + """ + with tf.compat.v1.name_scope( + name, + default_name='option_price', + values=[ + spots, strikes, volatilities, expiries, rates, + is_call_options + ]): + + spots = tf.convert_to_tensor(spots, dtype=tf.float64, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + volatilities = tf.convert_to_tensor(volatilities, tf.float64, name='volatilities') + expiries = tf.convert_to_tensor(expiries, tf.float64, name='expiries') + rates = tf.convert_to_tensor(rates, tf.float64, name='rates') + + z = tf.zeros_like(strikes) + + normal = tfp.distributions.Normal( + loc=tf.zeros([], dtype=spots.dtype), scale=1) + + df = tf.math.exp(-rates*expiries) + vt = volatilities * tf.math.sqrt(expiries) + + z = tf.where(rates == 0., (spots - strikes)/vt, + (spots-strikes*df)/(volatilities + * tf.math.sqrt(0.5*(1.-tf.math.exp(-2.*rates*expiries))/rates))) + + n1 = normal.cdf(z) + n2 = normal.prob(z) + calls = tf.where(rates==0., (spots - strikes) * n1 + vt * n2, + (spots - strikes*df)*n1 + + volatilities*tf.math.sqrt(0.5*(1-tf.math.exp(-2*rates*expiries))/rates)) + + + if is_call_options is None: + return calls + + puts = calls - spots + strikes * tf.math.exp(-rates*expiries) + + return tf.where(is_call_options, calls, puts) + + + +def vanilla_n_dawson_tf(forwards, + strikes, + volatilities, + expiries, + discount_factors = None, + is_call_options=None, + dtype = None, + name = None): + + """Computes the Black Scholes price for a batch of European options. + + ## References: + [1] Dawson, P., Blake, D., Cairns, A. J. G. and Dowd, K.: Options on normal under- + lyings, CRIS Discussion Paper Series – 2007.VII, 2007. + + Args: + forwards: A real `Tensor` of any shape. The current forward prices to + expiry. + strikes: A real `Tensor` of the same shape and dtype as `forwards`. The + strikes of the options to be priced. + volatilities: A real `Tensor` of same shape and dtype as `forwards`. The + volatility to expiry. + expiries: A real `Tensor` of same shape and dtype as `forwards`. The expiry + for each option. The units should be such that `expiry * volatility**2` is + dimensionless. + discount_factors: A real `Tensor` of same shape and dtype as the `forwards`. + The discount factors to expiry (i.e. e^(-rT)). If not specified, no + discounting is applied (i.e. the undiscounted option price is returned). + Default value: None, interpreted as discount factors = 1. + is_call_options: A boolean `Tensor` of a shape compatible with `forwards`. + Indicates whether to compute the price of a call (if True) or a put (if + False). If not supplied, it is assumed that every element is a call. + dtype: Optional `tf.DType`. If supplied, the dtype to be used for conversion + of any supplied non-`Tensor` arguments to `Tensor`. + Default value: None which maps to the default dtype inferred by TensorFlow + (float32). + name: str. The name for the ops created by this function. + Default value: None which is mapped to the default name `option_price`. + + Returns: + option_prices: A `Tensor` of the same shape as `forwards`. The Bachelier + price of the options. + + + #### Examples + ```python + spots = np.array([0.03, 0.02]) + strikes = np.array([.02, .02]) + volatilities = np.array([.004, .005]) + expiries = 2.0 + expiries = 1.0 + computed_prices = option_price( + forwards, + strikes, + volatilities, + expiries, + dtype=tf.float64) + # Expected print output of computed prices: + # + ``` + """ + with tf.compat.v1.name_scope( + name, + default_name='option_price_dawson', + values=[ + forwards, strikes, volatilities, expiries, discount_factors, + is_call_options + ]): + + forwards = tf.convert_to_tensor(forwards, dtype=None, name='forwards') + strikes = tf.convert_to_tensor(strikes, dtype=None, name='strikes') + volatilities = tf.convert_to_tensor(volatilities, dtype=None, name='volatilities') + expiries = tf.convert_to_tensor(expiries, dtype=None, name='expiries') + + if discount_factors is None: + discount_factors = 1. + discount_factors = tf.convert_to_tensor( + discount_factors, dtype=dtype, name='discount_factors') + + vt = volatilities * tf.math.sqrt(expiries) + normal = tfp.distributions.Normal( + loc=tf.zeros([], dtype=forwards.dtype), scale=1) + + z = (forwards - strikes) / vt + + n1 = normal.cdf(z) + n2 = normal.prob(z) + undiscounted_calls = (forwards-strikes) * n1 + vt * n2 + + if is_call_options is None: + return discount_factors * undiscounted_calls + undiscounted_forward = forwards - strikes + undiscounted_puts = undiscounted_calls - undiscounted_forward + + return discount_factors * tf.where(is_call_options, undiscounted_calls, + undiscounted_puts) + + + \ No newline at end of file diff --git a/tf_quant_finance/volatility/volbachelier_tf.py b/tf_quant_finance/volatility/volbachelier_tf.py new file mode 100644 index 000000000..589a1081e --- /dev/null +++ b/tf_quant_finance/volatility/volbachelier_tf.py @@ -0,0 +1,173 @@ +# Lint as: python3 +# Copyright 2020 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + + +import math as m +import tensorflow.compat.v2 as tf + +def func1(t, + x, + z, + betaStart, + betaEnd, + cut_off2, + cut_off3, + bLFK4, + cLFK4, + dLFK4): + """ helper function + + """ + + u = -(tf.math.log(z) + betaStart) / (betaEnd - betaStart) + + num = tf.where(u < cut_off2, + bLFK4[0] + u * (bLFK4[1] +u * (bLFK4[2] +u * (bLFK4[3] + u * (bLFK4[4] +u * (bLFK4[5] +u * (bLFK4[6] + u * (bLFK4[7] +u * (bLFK4[8] +u * bLFK4[9] )))) )))), + tf.where(u < cut_off3, + cLFK4[0] + u * (cLFK4[1] +u * (cLFK4[2] +u * (cLFK4[3] +u * (cLFK4[4] +u * (cLFK4[5] +u * (cLFK4[6] + u * (cLFK4[7] + u * (cLFK4[8] + u * cLFK4[9] )))))))), + dLFK4[0] + u * (dLFK4[1] +u * (dLFK4[2] +u * (dLFK4[3] +u * (dLFK4[4] +u * (dLFK4[5] +u * (dLFK4[6] + u * (dLFK4[7] + u * (dLFK4[8] + u * dLFK4[9] )))))))) + ) + ) + + den = tf.where(u < cut_off2, + 1.0 + u * (bLFK4[10] + u * (bLFK4[11] + u * (bLFK4[12] + u * (bLFK4[13] + u * (bLFK4[14] + u * (bLFK4[15] + u * (bLFK4[16] ) ))))) ), + tf.where(u < cut_off3, + 1.0 + u * (cLFK4[10] + u * (cLFK4[11] + u * (cLFK4[12] + u * (cLFK4[13] + u * (cLFK4[14] + u * (cLFK4[15] + u * (cLFK4[16]))))))), + 1.0 + u * (dLFK4[10] + u * (dLFK4[11] + u * (dLFK4[12] + u * (dLFK4[13] + u * (dLFK4[14] + u * (dLFK4[15] + u * (dLFK4[16]))))))) + ) + ) + hz = num / den + return tf.math.abs(x) / (tf.math.sqrt(hz * t ) ) + +def func2(price, + t, + x, + aLFK4): + """ helper function + + """ + p = tf.where(x < 0., price - x, price) + u = eta(tf.math.abs(x) / p ) + num = aLFK4[0] + u * (aLFK4[1] + u * (aLFK4[2] + u * (aLFK4[3] + u * (aLFK4[4] + u * (aLFK4[5] + u * (aLFK4[6] + u * (aLFK4[7]))))))) + den = 1.0 + u * (aLFK4[8] + u * (aLFK4[9] + u * (aLFK4[10] + u * (aLFK4[11] + u * (aLFK4[12]))))) + return p * num / den / tf.math.sqrt(t) + +def eta(z): +# case for avoiding incidents of 0/0, z close to zero + return tf.where(z < 1e-2, + 1 -z *(0.5+ z * (1.0 / 12.0 + z * (1.0 / 24.0 + z * (19.0 / 720.0 + z * (3.0 / 160.0 + z * (863.0/60480.0 + z * (275.0/24192.0))))))), + -z / tf.math.log1p(-z)) + +def vol_atm(price, t): + # atm case + return price * tf.math.sqrt(2. * m.pi / t) + +def vol_iotm(x,betaStart,betaEnd,price, t, cut_off1,cut_off2,cut_off3, aLFK4,bLFK4,cLFK4,dLFK4): + # other cases (ITM/OTM) + z = tf.where(x>=0.,(price - x) / x, -price / x) + + return tf.where(z <= cut_off1, + func1(t,x,z,betaStart,betaEnd,cut_off2,cut_off3,bLFK4,cLFK4,dLFK4), + func2(price,t,x,aLFK4) + ) + + + + +def volbachelier_tf(signs, + strikes, + forwards, + expiries, + prices): + """ computes the Bachelier implied volatility for a batch of prices of + European Call or Put options. + We assume a standard Brownian motion of the form + dS = sigma dW + for the underlying. sigma is the implied volatility. + + sign : A `Tensor` of any shape. The current sign that specifies if the + prices are Call (sign=1) or Put (sign=-1). + strikes : A real `Tensor`. The strikes of the options to be priced. + forwards : A real `Tensor` of same shape and dtype as `strikes`. + expiries : A real `Tensor` of same shape and dtype as `strikes`. + prices : A real `Tensor` of same shape and dtype as `strikes`. The prices + of the European Call, resp. Put prices + + Returns + + implied_bachelier_vols: A `Tensor` of the same shape as `strieks`. + The Bachelier implied volatilties of the Call, resp. Put options. + + +References: + [1] Fabien Le Floc'h, "Fast and Accurate Analytic Basis Point Volatility" + Link: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2420757 + + +""" + sign = tf.convert_to_tensor(signs, dtype=tf.float64, name='sign') + strike = tf.convert_to_tensor(strikes, dtype=tf.float64, name='strikes') + forward = tf.convert_to_tensor(forwards, dtype=tf.float64, name='forward') + t = tf.convert_to_tensor(expiries, dtype=tf.float64, name='expiries') + price = tf.convert_to_tensor(prices, dtype=tf.float64, name='price') + + # for rational expansion + aLFK4 = tf.constant([0.06155371425063157,2.723711658728403,10.83806891491789, + 301.0827907126612,1082.864564205999, 790.7079667603721, + 109.330638190985, 0.1515726686825187, 1.436062756519326, + 118.6674859663193, 441.1914221318738, 313.4771127147156, + 40.90187645954703], dtype=tf.float64) + bLFK4 = tf.constant([0.6409168551974357, 788.5769356915809, 445231.8217873989, + 149904950.4316367, 32696572166.83277, 4679633190389.852, + 420159669603232.9, 2.053009222143781e+16, 3.434507977627372e+17, + 2.012931197707014e+16, 644.3895239520736, 211503.4461395385, + 42017301.42101825, 5311468782.258145, 411727826816.0715, + 17013504968737.03, 247411313213747.3], dtype=tf.float64) + cLFK4 = tf.constant([0.6421106629595358, 654.5620600001645, 291531.4455893533, + 69009535.38571493, 9248876215.120627, + 479057753706.175, 9209341680288.471, 61502442378981.76, + 107544991866857.5, 63146430757.94501, + 437.9924136164148, 90735.89146171122, 9217405.224889684, + 400973228.1961834, 7020390994.356452, + 44654661587.93606, 76248508709.85633], dtype=tf.float64) + dLFK4 = tf.constant([0.936024443848096, 328.5399326371301, 177612.3643595535, + 8192571.038267588, 110475347.0617102, + 545792367.0681282, 1033254933.287134, 695066365.5403566, + 123629089.1036043, 756.3653755877336, + 173.9755977685531, 6591.71234898389, 82796.56941455391, + 396398.9698566103, 739196.7396982114, + 493626.035952601, 87510.31231623856], dtype=tf.float64) + + #eps = tf.cast(tf.constant(np.finfo(float).tiny),tf.float32) + # constants + cut_off_atm = tf.constant(1.0e-10, dtype=tf.float64) # ATM Cutoff level #np.finfo(float).eps + cut_off1 = tf.constant(0.15, dtype=tf.float64) # cut-off for -C(x)/x + cut_off2 = tf.constant(0.0091, dtype=tf.float64) # 1st cut-off for tilde(eta) + cut_off3 = tf.constant(0.088, dtype=tf.float64) # 2nd cut-off for tilde(eta) + + betaStart = - tf.math.log(cut_off1) ; betaEnd = 708.3964185322641#- tf.math.log(machine eps) + + x = ( forward - strike) * sign # intrinsic value of the Call (sign=1) + # or Put (sign = -1) + implied_bachelier_vols = tf.where(tf.math.abs(x) < cut_off_atm, + vol_atm(price,t), + vol_iotm(x,betaStart,betaEnd,price, t, cut_off1,cut_off2,cut_off3, aLFK4,bLFK4,cLFK4,dLFK4) + ) + + # return the full tensor of implied Bachelier volatilities due to atm and itm/otm + return implied_bachelier_vols + +