Package CharFuncPricing

Status	Coverage

Julia package to provide reference European option prices for stochastic volatility models with a known characteristic function, such as the Heston stochastic volatility model. For the Heston and Schobel-Zhu models, a function provides the first, second and fourth cumulants through analytical formulas.

The code is not meant for production purpose and does not always cater for corner cases. It however supports arbitrary precision via the Nemo package.

Installation

In a Julia REPL, enter pkg mode (by pressing ]) and run:

(v1.0) pkg> add CharFuncPricing

Julia version 1.2 or higher is required.

Cos method

Float64 Usage

Start by creating a HestonParams structure, which represents the parameters of the Heston model: v0, κ, θ, ρ, σ.

params = HestonParams(v0, κ, θ, ρ, σ)

Then make a CosCharFuncPricer structure via makeCosCharFuncPricer. This function will store the relevant m values of the characteristic function for the range [a,b] defined by l deviations using the cumulants rule a = c1 - l * sqrt(c2+sqrt(c4)), b = c1 + l * sqrt(c2+sqrt(c4)).

cf = DefaultCharFunc(params)
pricer = makeCosCharFuncPricer(cf, τ, m, l)  

And price vanilla call and puts of a given strike as following

priceEuropean(pricer, false, strike, forward, τ, df)

The second parameter specifies whether we want to price a call (true) or a put (false). The last parameter specifies the discount factor to maturity.

The first, second and fourth cumulants are given by

c1,c2,c4 = computeCumulants(params, τ)

Nemo Usage

The only difference is to make sure the parameters are of ArbField type (real arbitrary precision), the function makeCosCharFuncPricer should also be called on the AcbField type.

using Nemo
R = ArbField(256)
CC = AcbField(256)
cf = NemoCharFunc{HestonParams{arb}}(params, CC)
pricer = makeCosCharFuncPricer(cf, τ, m, l)

Float64 Example

Here is how to price a put option with parameters

r=0.01; q=0.02
κ=4.0; θ=0.25; σ=1.0; ρ=-0.5; v0=0.04
τ = 1.0
spot = 100.0; strike = 80.0
spot *= exp((r - q) * τ)
df = exp(-r * τ)
params = HestonParams(v0, κ, θ, ρ, σ)
cf = DefaultCharFunc(params)
l = 32; m = 1024
pricer = makeCosCharFuncPricer(cf, τ, m, l)
priceEuropean(pricer, false, strike, spot, τ, df)

The result is 7.95887811325676.

### Nemo Example
```julia
using Nemo
using CharFuncPricing

R = ArbField(256)
CC = AcbField(256)
r = R(QQ(1, 100)); q = R(QQ(2, 100))
κ = R(4.0); θ = R(QQ(1, 4)); σ = R(1.0); v0 = R(QQ(4, 100)); ρ = R(-0.5);
τ = R(1.0)
spot = R(100); strike = R(80)
isCall = false
spot *= exp((r - q) * τ)
df = exp(-r * τ)
params = HestonParams(v0, κ, θ, ρ, σ)
cf = NemoCharFunc{HestonParams{arb}}(params, CC)
l = 48; m = 1024*4
pricer = makeCosCharFuncPricer(cf, τ, m, l)
priceEuropean(pricer, isCall, strike, spot, τ, df)

The result is:

7.95887811325676828521326060761429303089865693725960319205094095681790030 +/- 4.83e-72.

With l=64; m=1024*8, we obtain:

7.95887811325676828521326060761429303089865693725960319205094095681878397 +/- 3.71e-72.

Andersen-Lake

Float64 Example

It starts similarly as for the Cos method, the only change is how to build the pricer variable.

pricer = ALCharFuncPricer(cf)
priceEuropean(pricer, false, strike, spot, τ, df)

The result is 7.9588781132567705.

BigFloat Example

We use BigFloat for the Heston parameters and option characteristics, as well as in the quadrature tolerance.

r=BigFloat("0.01"); q=BigFloat("0.02")
κ=BigFloat(4.0); θ=BigFloat("0.25"); σ=BigFloat(1.0); ρ=BigFloat("-0.5"); v0=BigFloat("0.04")
τ = BigFloat(1.0)
spot = BigFloat(100.0); strike = BigFloat(80.0)
spot *= exp((r - q) * τ)
df = exp(-r * τ)
params = HestonParams(v0, κ, θ, ρ, σ)
cf = DefaultCharFunc{HestonParams{BigFloat},Complex{BigFloat}}(params)
quad = TanhSinhQuadrature(800, BigFloat(1e-200))
pricer = ALCharFuncPricer(cf,quad)
priceEuropean(pricer, false, strike, spot, τ, df)

The result is

7.95887811325676828521326060761429303089865693725960319205094095681918541918632

Adaptive Flinn

This is the adaptive Flinn quadrature using the transformation to (-1, 1) interval. No truncation is involved.

Float64 Example

With a quadrature tolerance of 1e-8:

pricer = AdaptiveFlinnCharFuncPricer(cf, τ, qTol = 1e-8)
priceEuropean(pricer, false, strike, spot, τ, df)

The result is 7.958878112874899

BigFloat Example

The adaptive Flinn pricer works with high accuracy, but does not perform very well then. It is more intended for the calculation of prices with a absolute error tolerance of around 1e-8 or 1e-10.

pricer = AdaptiveFlinnCharFuncPricer(cf, τ, qTol = BigFloat(1e-24))
priceEuropean(pricer, false, strike, spot, τ, df)

The result is

7.958878113256768285213257572750089190600415520655780637746847607110529890012863

and the effective accuracy is 3e-24. In total, length(pricer.kcos[1,:])=31511 points are used, compared to 295 for an tolerance of 1e-10. For a tolerance of 1e-32, 450799 points are used. This means that algorithm is asymptotically linear on this example.

Testing

In a Julia REPL session, enter pkg mode and run test CharFuncPricing.

Unit tests verify the call and put option prices against the reference prices of Alan Lewis in double and arbitrary precision. In fact, the implementation here gives more precise results (minimum accuracy of 1e-60 while Alan Lewis numbers are only accurate up to 1e-25).

Cumulants are checked against a Taylor series algorithmic differentiation.

References

Andersen, L.B.G. and Lake, M. (2018) [Robust high-precision option pricing by Fourier transforms: Contour deformations and double-exponential quadrature](SSRN 3231626)

Fang, F. and Oosterlee, C. W. (2008) A novel pricing method for European options based on Fourier-cosine series expansions

Healy, J. (2021) Applied Quantitative Finance for Equity Derivatives

Le Floc'h, F. (2018) More Robust Pricing of European Options Based on Fourier Cosine Series Expansions

Le Floc'h, F. (2020) An adaptive Filon quadrature for stochastic volatility models