Option Pricing Models Driven by the Space-Time Fractional Diffusion: Series Representation and Applications
Jean-Philippe Aguilar Jan Korbel










§ INTRODUCTION


The pricing of derivatives, and notably of options, is a central subject in mathematical finance. It allows the market practitioner to estimate the value of its portfolio, and to construct appropriate hedging strategies. The most popular option pricing model is the one introduced by Black and Scholes <cit.>, because of its simplicity (e.g., the option price can be expressed in terms of simple mathematical functions), and can be used to imply the market parameters, such as volatility surfaces, from the observation of traded market prices.

On the other hand, the simplicity of the Black-Scholes (BS) model is also its main limitation. The dynamics of the underlying asset is described by a geometric Brownian motion, so the resulting option price is given by the Gaussian distribution. This assumption does not describe well extreme market events such as sudden jumps of asset prices, which happen far more frequently than expected in the Gaussian world (see for instance the influential book by Taleb <cit.>). This makes the Black-Scholes formula less reliable in abnormal conditions or in illiquid markets.

During the past years, several models have been introduced to describe the market dynamics more realistically; one can mention regime switching multifractal models <cit.>, stochastic volatility models <cit.> or jump (Lévy-stable) processes <cit.>. More recently, a model based on space-time fractional diffusion has been introduced <cit.>, and can be regarded as a generalization of the Lévy-stable model; analytic resolution of this model has been provided in <cit.> under the form of a series representation for its pricing formula. In the present article, we briefly recall these analytic results (and notably how they recover previously known models) and test their efficiency in real market applications. We also discuss the related topics as at-the-money approximation or implied volatility.   


The paper is organized as follows. In the following section, we introduce some fundamental concepts in option pricing, and the main models that, under risk-neutral approach, can be reduced to a space-time fractional diffusion problem. This includes the BS model (which reduces to the classical heat equation), the Lévy-stable model (which reduces to the space-fractional diffusion equation) and the generic space-time (or double) fractional model. In section 3, we briefly recall the analytic solution to the space-time fractional model. In section 4 we present various applications of the theoretical results: we calculate the call prices and compare it with the real data, we introduce at-the-money volatility and discuss the construction of volatility smile.  The last section is dedicated to conclusions.




§ OPTION PRICING


The price of an option of strike K and maturity T, is a function of market parameters such as an underlying asset price S, a risk-free interest rate r and a market volatility σ. We will denote this price by V(S,K,r,σ,t). In the case of an European option, it is characterized by its payoff, that its, its value at the exercise time T; for an European call, this value is equal to 

    V(S,K,r,σ,T)   =  max{ S-K, 0 }  :=   [S-K]^+

For a  put option, the corresponding payoff is [K-S]^+. Now we recall the principles of option pricing, that is, the way of determining V(S,K,r,σ,t).



 §.§ The risk-neutral approach


The risk-neutral, or risk-free approach is based on the idea that one can construct a portfolio where the (market) risk can be totally eliminated <cit.>. Schematically, it consists in buying an option and selling a certain quantity Δ (to be determined) of the underlying price, so that the total value of the portfolio reads Π  =   V - Δ S and therefore:

    Π  =   V - Δ S

On the other hand, the markets are assumed to offer no arbitrage opportunity, that is, any risk-less portfolio will have the same yield as if it were capitalized at the risk-free interest rate:

    Π  =   rΠ t

Equalizing (<ref>) and (<ref>) and making an appropriate choice for Δ transforms the option pricing problem into the resolution of a partial differential equation with terminal condition.

From a more theoretical point of view, within risk-neutral approach, the price can be formulated as the discounted expectations of the terminal payoff <cit.>:

    V(S,K,r,μ,τ)   =   e^-rτ 𝔼^ℚ[ [S-K]^+ ]

where we have introduced the time-to-maturity τ = T-t. The expectations are to be taken under the risk-neutral measure ℚ, which is associated to the original probability measure ℙ via the Radon-Nikodym derivative:

    ℚ_t/ℙ_t  =   e^S_t - μ t

The risk-neutral parameter μ can be expressed as

    μ  =   - log𝔼^ℙ[ e^S_t=1]

More details can be found in <cit.>.




 §.§ Black-Scholes model

In the BS model, the underlying asset price S is assumed to be described by a geometric Brownian motion:

    S_t  =   r   S  t  +  σ  S   W_t

It follows from Itô's lemma <cit.> that the total differential of the option price is:

    V   =  ∂ V/∂ t  t   +  ∂ V/∂ S  S   +  1/2σ^2 S^2 ∂^2 V/∂ S^2  t

Choosing Δ=∂ V/∂ S, using (<ref>) and equalizing (<ref>) and (<ref>), we have shown that the call price satisfies the famous Black-Scholes equation, which is a partial differential equation (PDE) with terminal condition:

    {   ∂ V/∂ t  +  1/2σ^2 S^2∂^2 V/∂ S^2  +   r S ∂ V/∂ S  -   rV   =   0   t∈[0,T]  
        V(S,K,r,σ,t=T)   =   [S-K]^+
    .

It is known that, with the change of variables

    {    x   :=  log S   +   (r-σ^2/2)  τ
       τ  :=   T-t 
        V(S,K,r,σ,t)    :=   e^-rτW(x,K,r,σ,τ)
    .

then the Black-Scholes PDE (<ref>) resumes to the diffusion (or heat) equation

    ∂ W/∂τ  -  σ^2/2∂^2 W/∂ x^2  =   0

which is a particular case of the double fractional diffusion (<ref>) with time fractionality γ=1 and space fractionality α = 2. It is well known 
that the Green function for (<ref>) is the heat kernel

    g(x,K,r,σ,τ)   =  1/σ√(2πτ)  e^-x^2/2σ^2τ

and therefore, by the method of Green functions and turning back to the initial variables, we obtain the solution for the Black-Scholes PDE (in the call case):

    V(S,K,r,σ,τ)   =   e^-rτ ∫_-∞^+∞ [Se^(r-σ^2/2)τ+y-K]^+  g(y,K,r,σ,τ)  𝕀 y

Basic manipulations on the integral (<ref>) yield

    V(S,K,r,σ,τ)   =   SN(d_+) - Ke^-rτ N(d_-)  d_± = 1/σ√(τ)(logS/K+rτ) ±1/2σ√(τ)

where N(.) is the normal distribution function; formula (<ref>) is the celebrated Black-Scholes formula for the European call. The corresponding risk-neutral parameter is therefore 

    μ_BS = - σ^2/2




 §.§ Finite-Moment Lévy-stable model


An interesting generalization of the BS model is the so-called Finite Moment Lévy (or Log) Stable (FMLS) model; it was introduced in <cit.> and assumes that the underlying asset price S_t is described by:

    S_t   =   r S_t  t   +  σ S_t  L_α,β(t)

where L_α,β(t) is the Lévy process <cit.>. α∈[0,2] and β∈[-1,1] are the so-called stability and asymmetry parameters and determine the decay of the tails and the asymmetry of the probability distributions g_α,β(x,t). Under the (strong) hypothesis that β=-1 (maximal negative asymmetry hypothesis) then the distribution g_α,-1:=g_α possesses one heavy-tail in the negative axis, and another tail in the positive axis with exponential decay as soon as α>1, and finite exponential moments 

    𝔼^ℙ [ e^- λ S_t]   =   e^-λ^α( σ/√(2))^α/cosπα/2


These particular Lévy distributions are sometimes called Lévy-Pareto distributions; their relevance in financial modelling has been known since the works of Mandelbrot and Fama in the 1960s <cit.>. They are known to satisfy the space-fractional equation:

    ∂ g_α(x,t)/∂τ  +  μ  D^α g_α(x,t) = 0

where D^α:=   ^α-2 D ^α is a particular case of the Riesz-Feller operator (<ref>) for θ = α-2. The condition θ = α -2 turns out to be the fractional analogue to the probabilistic condition β = -1. Note that equation (<ref>) degenerates into the the reduced BS equation (<ref>) when α=2; the corresponding call option price is then

    V_α(S,K,r,μ,τ)   =   e^-rτ ∫_-∞^+∞ [Se^(r+μ)τ+y-K]^+  g_α(y,τ)  𝕀 y

where the risk-neutral parameter μ follows from (<ref>):

    μ  =  ( σ/√(2))^α/cosπα/2

and reduces to μ=-σ^2/2 in the Gaussian case (α=2). An analytic resolution of the FMLS model has been provided in <cit.>, under the form of a quickly convergent series representation for the call price (<ref>). The proof is based on the Mellin-Barnes representation for the solutions of the space fractional equation (<ref>) (see <cit.>): if x>0 then

    g_α(x,τ)   =  1/α x  ∫_c_1-i∞^c_1+ i ∞ Γ(1-t_1)/Γ(1-t_1/α) ( x/(-μτ)^1/α)^t_1 𝕀 t_1/2iπ  0 < c_1 < 1




 §.§ Space-time option pricing model




Let us discuss option pricing models based on space-time (double)-fractional diffusion equation, which can be expressed as

    (^∗_0 𝒟^γ_t +  μ [^θ𝒟^α_x] ) g(x,t) = 0

where  α∈ (0,2], γ∈ (0,α]. Asymmetry parameter θ is defined in the so-called Feller-Takayasu diamond |θ| ≤min{α, 2-α}.  ^∗_0 𝒟^γ_t denotes the Caputo fractional derivative, which is defined as

    ^∗_t_0𝒟^ν_t f(t) =
    1/Γ(⌈ν⌉ - ν)∫_t_0^t f^⌈ν⌉(τ)/(t - τ)^ν +1-⌈ν⌉τ

and ^θ𝒟^α_x denotes the Riesz-Feller fractional derivative, which is usually defined via its Fourier image as

    ℱ[^θ𝒟^ν_x
    f(x)](k) = - ^θψ^ν (k)F[f(x)](k) = - μ |k|^ν e^i(sign k) θπ /2ℱ[f(x)](k)


Let us describe the various financial models that are included in (<ref>).
The FMLS model, although far more generic than the BS one, can still be regarded as too restrictive; this is because the maximal negative asymmetry hypothesis β = -1, or equivalently θ = α - 2, does not describe well all capital markets (in particular illiquid ones, where financial assets often exhibit an almost symmetric heavy-tail). Nevertheless, it is not a priori possible to relax the maximal negative asymmetry hypothesis, because when β≠-1 the expectations (<ref>) are known to diverge <cit.>. The fact that the risk-neutral parameter is infinite in this case traduces the fact that the risk cannot completely be eliminated from this class of Lévy processes. Risk-minimal (instead of risk-neutral) approach has been introduced (see <cit.>) in this case; an interesting possibility, to generalize the the FMLS model and to remain within the risk-neutral framework, is to allow the time derivative to be also fractional (in the Caputo sense) in eq. (<ref>):

    ^∗_0 𝒟^γ_t g_α,γ(x,t)   +  μ_γ [^2-α  D^α] g_α,γ(x,t) = 0

The corresponding call option price now reads

    V_α,γ(S,K,r,μ_γ,τ)   =   e^-rτ ∫_-∞^+∞ [Se^(r+μ_γ)τ+y-K]^+  g_α,γ(y,τ)  𝕀 y

The Green functions are also known under the form of a Mellin-Barnes line integral <cit.>:

    g_α,γ(x,τ)   =  1/α x  ∫_c_1-i∞^c_1+ i ∞ Γ(1-t_1)/Γ(1-γ/αt_1) ( x/(-μ_γτ)^1/α)^t_1 𝕀 t_1/2iπ  0 < c_1 < 1

for any x>0. The main difference with the Lévy-stable price (<ref>) is that the risk-neutral parameter μ_γ now depends on the time-fractionality γ, and is not known analytically like in the Lévy stable case (<ref>). In <cit.>, an efficient and simple series expansion is derived for the risk neutral parameter, as well as a fast converging series expansion for the call price (<ref>). In the next section, we discuss these results in detail, and test them in the real market conditions.



§ SERIES REPRESENTATION OF THE PRICING FORMULAS UNDER THE SPACE-TIME FRACTIONAL DIFFUSION


Now, let us provide an analytic pricing formulas for the call options driven by the fractional diffusion (<ref>) (details of the proofs can be found in <cit.>). We assume that 1<α≤ 2 and 0< γ≤α. 



 §.§ Risk-neutral parameter


The expectations in definition (<ref>) over the probability measure ℙ can be expressed in terms of its probability densities g_α,γ(y,τ), that is:

    μ_γ  =  -log∫_-∞^∞ e^y g_α,γ(y,τ=1)   y

It is possible to bring the calculation back to the non time-fractional case, by writing (see details in <cit.>):

    g_α,γ(x,τ) = ∫_0^∞ g_γ(τ,l) g_α(l,x)  l

where g_γ and g_α^θ are solutions of single-fractional diffusion equations

    ∂ g_γ(t,l)/∂ l   =   ^∗_0 𝒟^γ_t   g_γ(t,l) 
    ∂ g_α(l,x)/∂ l   =   𝒟^α_x   g_α(l,x)






  §.§.§ Mellin-Barnes representation of the risk-neutral parameter

The solution to the Caputo equation (<ref>) is known to be <cit.>:

    g_γ(τ,l)   =  1/τ^γ  M_γ( l/τ^γ)

where M_ν(z) is a function of Wright type, admitting the following Mellin-Barnes representation <cit.>:

    M_ν(z)    =   ∫_c-i∞^c+i∞ Γ(s)/Γ(ν s + 1 - ν)   z^-s  s/2iπ c>0  .

Inserting (<ref>) and (<ref>) in (<ref>), interverting the integrals and using (<ref>) we obtain a Mellin-Barnes representation for the risk-neutral parameter:

    μ_γ  =   - log[ 1/α ∫_c-i∞^c+i∞ Γ(s)Γ(1-s/α)/Γ(γ s + 1 - γ)  μ_1^s-1/α  s/2iπ]  0 < c < 1  .

where μ_1=( σ/√(2)) ^απα/2 is the risk-neutral parameter in the Lévy-stable case (γ=1), cf. eq. (<ref>).



  §.§.§ Series representation of the risk-neutral parameter

It is possible to express the integral over a vertical line (<ref>) as a sum of residues of its analytic continuation, on the condition that the integrand decreases sufficiently fast at infinity. For Gamma function, this condition is determined by the well-known Stirling approximation <cit.>

    Γ(x   +   i y)  |x| →∞
     |y| →∞∼ √(2π)  |y|^x-1/2 e^-π|y|/2  x,y  ∈ ℝ

It follows from (<ref>) that, for a Gamma function of linear arguments of the type Γ(a s +b),   a,b∈ℝ, its behavior at infinity depends on the sign of a, namely:

    {    a > 0  |Γ(as+b)||s|→∞⟶ 0   s  ∈  (3π/2 , -π/2 )  
        a < 0  |Γ(as+b)||s|→∞⟶ 0   s  ∈  (-π/2 , π/2 )
    .

(<ref>) easily generalizes to a ratio of products of Gamma functions of linear arguments. Let us assume that a function f admits a Mellin transform f^* of the form:

    f^* (s)   =  Γ(a_1 s + b_1)  … Γ(a_n s + b_n)/Γ(c_1 s + d_1)  … Γ(c_m s + d_m) a_j, b_j, c_k, d_k ∈ℝ

and that this Mellin transform converges on some non-empty strip { Re(s)  ∈  (c_1,c_2)}, so that the Mellin inversion formula holds:

    f(x)   =  ∫_c-i∞^c+i ∞  f^*(s)   x^-s  s/2iπ c ∈ (c_1,c_2)

Introduce the characteristic quantity Δ:

    Δ  =  ∑_j=0^n   a_j   -  ∑_k=0^m   c_k

It follows from (<ref>) that Δ governs the asymptotic behavior of f^*(s):

    {   Δ > 0  |f^*(s)| |s|→∞⟶ 0   s  ∈  (3π/2 , -π/2 ) 
       Δ < 0  |f^*(s)| |s|→∞⟶ 0   s  ∈  (-π/2 , π/2 )
    .

Therefore, applying the residue theorem to the inversion formula (<ref>) yields:

    {   Δ > 0  f(x)   =  ∑_Re(s) < c[ f^*(s)x^-s] 
       Δ < 0  f(x)   =   -∑_Re(s) > c[ f^*(s)x^-s]
    .

where the choice { Re(s) > c } or { Re(s) < c } is determined by the fact that |x|^-s goes to 0 at infinity in the chosen half-plane. In the case of the  Mellin-Barnes representation (<ref>), the characteristic quantity is 

    Δ  =    1 - 1/α  -  γ

and is negative as soon as:

    γ  >   1 - 1/α

It follows from rule (<ref>) that one can express (<ref>) as the sum of the residues in the right half-plane (we choose it because μ_1^s-1/α goes to 0 at infinity in this half plane as soon as |μ_1|<1, which is the case in all financial applications). These poles are induced by the singularities of the Γ(1-s/α) term, which arise at every negative integer value -n of its argument, that is at every point of the type s=1+α n,   n∈ℕ; it is well-known the residue of the Gamma function at a negative integer -n is (-1)^n/n! <cit.>, and therefore we obtain:

    μ_γ  =   - log ∑_n=0^∞ (-1)^n Γ(1+α n)/n!Γ(1+γα n)μ_1^n

as soon as the condition (<ref>) is fulfilled. An interesting approximation of (<ref>) can be easily derived from the Taylor approximation log(1+u) ≃ u:

    μ_γ    =     -log[  1   -  Γ(1+α)/Γ(1+γα )μ_1   +  …] 
       ≃   Γ(1+α)/Γ(1+γα )μ_1

which, as expected, coincides with the Lévy-stable risk-neutral parameter μ_1 when γ=1. As a particular case, we obtain a nice approximation for the risk-neutral parameter in the fractional Black-Scholes model (α = 2):

    μ_γ ≃  -σ^2/Γ(1+2γ)

which resumes to the well-known gaussian parameter -σ^2/2 when γ=1.




In Fig. <ref>, we plot the evolution of μ in function of the parameters γ, σ and α; thanks to the exponential convergence of the series (<ref>), it suffices to consider only the very first few terms of the series to obtain an excellent level of precision.








 §.§ Option price

In all the following we will use the notation [log]:=logS/K+rτ, so that the payoff in (<ref>) can be written:

    [Se^(r+μ_γ)τ+y-K]^+   =   K[e^[log]+μ_γτ+y-1]^+




  §.§.§ Mellin-Barnes representation of the option price

The call price (<ref>) can be expressed as a double Mellin-Barnes integral. First, one introduces in (<ref>) the Mellin-Barnes representation (<ref>) for the Green function; second, one writes a Mellin-Barnes representation for the exponential term in (<ref>):

    e^[log]+μ_γτ + y  =  ∫_c_2-i∞^c_2+i∞  (-1)^-t_2Γ(t_2) ([log] + μ_γτ + y )^-t_2 dt_2/2iπ c_2 > 0

The integral over the Green variable y becomes a particular case of a Bêta integral, which is straightforward to calculate, and one obtains the representation for the option price:

    V_α,γ (S,K,r,μ_γ,τ)   =  K e^-r τ/α∫_c_1-i∞^c_1+i∞∫_c_2-i∞^c_2+i∞  (-1)^-t_2Γ(t_2)Γ(1-t_2)Γ(-1-t_1+t_2)/Γ(1-γ/αt_1) 
    ×
    (-[log]-μ_γτ)^1+t_1-t_2(-μ_γτ^γ)^-t_1/α t_1/2iπ∧ t_2/2iπ

The vector c:= [ c_1 , c_2 ] is an element of the ℂ^2-polyhedra P:={ (t_1,t_2)∈ℂ^2,   0<Re(t_2)<1 ,   Re(t_2-t_1)>0  }, which generalizes the notion of convergence strip for one-dimensional Mellin transform.



  §.§.§ Series representation of the option price

The double Mellin-Barnes integral (<ref>) can also be expressed as a sum of residues. In one dimension, it is usual to sum the residues right or left to the convergence strip of the Mellin transform, like we have done for the risk-neutral parameter, where the integral (<ref>) has been computed by right-summing the residues to obtain the series (<ref>). In two dimensions, this procedure generalizes to a summation to a subregion of ℂ^2, determined by a characteristic vector associated to the integrand. The incoming residues are computed by the two-dimensional analogue to the Cauchy formula:

    [ f(t_1,t_2)   t_1/2iπ t_1∧ t_2/2iπ t_2]   =   f(0,0)

This procedure has been introduced in <cit.>. Namely, the characteristic quantity (<ref>) generalizes to a characteristic vector, which, in the case of the double Mellin-Barnes integral (<ref>) reads

    Δ  =  [ -1 + γ/α;        1 ]

The rule (<ref>) generalizes to

    f(x_1,x_2)   =  ∑_[t_1,t_2]∈Π_Δ [ f^*(t_1,t_2)   x_1^-t_1  x_2^-t_2 ]

where Π_Δ is the subset of ℂ^2 defined by

    Π_Δ  =  {t:=[t_1,t_2]∈ℂ^2   ,  Re( Δ . s )   <   Re ( Δ . c )  }

in the sense of the euclidean scalar product. In the plane Re  ( ℂ^2 ), Π_Δ is therefore the region located under the line

    Re(t_2)   =   (1 - γ/α ) (Re(t_1) - c_1)   +   c_2

whose slope is positive because by hypothesis γ≤α. In this region, poles come from functions Γ(-1-t_1+t_2) and Γ(t_2) which are singular at every negative integer value of their argument (see Fig. <ref>).



From the singular behavior of the Gamma function around a singularity <cit.> and the Cauchy formula (<ref>), we obtain the series for the call price under double-fractional model:

    V_α,γ (S,K,r,μ_γ,τ)   =  Ke^-rτ/α ∑_n = 0 
     m = 1^∞ (-1)^n/n!Γ(1-γn-m/α) (-[log]-μ_γτ)^n(-μ_γτ^γ)^m-n/α


Full details of this calculation can be found in <cit.>.




§ APPLICATIONS


Let us discuss several applications of the series formula for the space-time fractional option prices. We show that it can be used for estimating the market parameters of the option prices. The calculation of the option price is very quick compared to the other methods (Mellin-Barnes representation, numerical estimation, …). We also briefly discuss the applications to implied volatility. 



 §.§ Call price


When fixing an upper bound for the n (resp. m) summation in the double series (<ref>), we are left with a simple series whose m-(resp. n) partial sums converges very quickly to the option price (see fig.  <ref>, where the parameters are S=3800,K=4000,r=1%,σ=20%, α=1.7, γ = 0.9). We may observe that the convergence of the m-sums are monotone, while the n-sums oscillate around the final price.



In the graphs in fig. <ref>, we study the evolution of the option price (<ref>) in function of different parameters. In the first graph we fix S=3800, K=4000, r=1% σ= 20% and we plot the evolution of the price in function of γ, for different stability parameters α∈[1.5,2]; we choose to consider only γ>0.33 so that the condition γ>1-1/α is satisfied for all stabilities, and observe that the prices are a decreasing function of the time fractionality. In graph 2 we let α vary between 1 and 2 and note that when γ≤1 then the prices are always a decreasing function of the stability, while for γ>1 they possess a maximum. In graph 3 (resp. 4) we plot the evolution of the option price in function of the spot price S (resp. market volatility σ) for various time fractionality γ, and with fixed stability α=1.7; note that the prices are, as expected, always a monotonous (growing) function of the spot and of the volatility, which is coherent with the non-arbitrage hypothesis of financial markets. Finally, we show the estimated parameters of the three option pricing models for the real options of S&P 500 options. The results are presented in Tab. <ref>, which has been taken from <cit.>.









 §.§ Implied volatility



The process of implying the market volatility consists in finding for which volatility σ_I a model-driven option price coincides with the observable price C, that is when

    V(S,K,r,σ_I,τ)   =   C

A typical procedure is to imply a Black-Scholes volatility (by using the Black-Scholes formula for the price and solving (<ref>) by means of numerical methods, such as a Newton-Raphson algorithm, see for instance <cit.>) and use it as an input parameter in a more sophisticated model. Let us show how the analytic series (<ref>) allows to imply a market volatility, and compare with the Gaussian one.



  §.§.§ At-the-money volatility



When the asset is "at-the-money forward", that is when

    S   =   K e^-rτ

then there exists an approximation for the Black-Scholes formula <cit.>

    V(S,K,r,σ,τ)  ≃ S/√(2π)σ√(τ)

and therefore the solution to the implied volatility equation (<ref>) reads

    σ_I  ≃ C/S √(2π/τ)

Such an approximation can also be derived in the double-fractional Black-Scholes model (α=2): note that, with our notations, the ATM-forward hypothesis (<ref>) reads [log]=0 and therefore in this case the pricing formula (<ref>) becomes a power series (i.e., with only positive powers of μ_γ and τ):

    V_2,γ(S,K,r,μ,τ)   =  S/2 [  1/Γ(1+γ/2)√(-μ_γτ^γ)  +   O(-μτ^γ)  ]

Using approximation (<ref>) for the risk-neutral parameter

    μ_γ  =   - σ^2/Γ(1+2γ)

in the first order term of the power series (<ref>), we obtain the implied fractional Black-Scholes volatility (in the ATM forward case):

    σ_I  ≃  2  C/S Γ(1+γ/2)  √(Γ(1+2γ)/τ^γ)

Let us remark that the formula (<ref>) resumes to the Black-Scholes implied volatility formula (<ref>) when γ=1 (recall that Γ(3/2) = √(π)/2). In graph <ref> we plot the evolution of formula (<ref>) in function of γ for a time to maturity τ = 1.027 and various exercise and call prices (see market datas in table <ref>).





  §.§.§ Volatility smile




In table <ref>, we provide observable market bid (offered) prices for S & P 500  call options with several exercise (strike) prices traded end 2008, quotation date = 03 nov. 2008, expiry = 17 jan. 2009 (source: eurexchange.com). We compute the implied Black-Scholes volatility as well as the implied fractional Black-Scholes volatility for various time fractionalities. They are obtained via a truncation of the series (<ref>) to n,m=4 and the approximation (<ref>) for the parameter μ_γ.




In fig <ref> we plot the implied volatilities obtained in table <ref> for and for 0.8 ≤γ≤ 1.1. We observe that the usual volatility smile (that is, the existence of a minimum around the spot price) is preserved, although less smooth when γ>1. Interestingly, when γ≤ 1, the minimal implied volatility is attained for the same strike price (independently of γ).
































§ CONCLUSIONS

In this paper, we have discussed the application of space-time fractional diffusion in option pricing and its relation to Black-Scholes model and Finite moments Lévy stable model. Models based on fractional diffusion enable to model the risk redistribution in order to incorporate large drops, memory effects and abnormal periods. We have briefly introduced all aforementioned models and described their main properties. Additionally, we have presented the series representation for all models, which is based on Mellin-Barnes integral representation of the option price and residue summation in ℂ^2. This mathematical techniques can overcome the technical difficulties of the fractional models, which is caused by the fact that the resulting prices are normally expressed in terms of integral transforms (Fourier, Laplace, Mellin) and the practical calculation is time consuming and understandable only to people trained in fractional calculus. The resulting series representation can be easily grasped by any financial practitioner. We have also applied the formulas to real financial data in order to demonstrate fast convergence and stability of the method. We have particularly shown numerical estimations of model parameters from the real data, applications to implied volatility and presence of volatility smile. 

Fractional models provide a fruitful field for further investigations of  financial systems, including portfolio management, derivative pricing, commodity pricing and many other possible applications. Naturally, in these applications it is necessary to carefully define the proper fractional derivatives and boundary conditions. In some cases, as e.g. in the case of fractional geometric Brownian motion, it is also necessary to overcome the mathematical issues, as non-existence of moments, etc.  Some of these topics will be addressed in the future research.  



§ ACKNOWLEDGEMENTS

J. K. acknowledges support from the Austrian Science Fund, Grant No. I 3073-N32., and from the Czech Science Foundation, Grant No. 17–33812L.

99

Abramowitz72
M. Abramowitz and I. Stegun, Handbook of mathematical functions. Dover Publications (1972).

ACK17-Levy
J.-Ph. Aguilar, C. Coste, J. Korbel, Non-Gaussian analytic option pricing: a closed formula for the Lévy-stable model, 2017, 	arXiv:1609.00987, submitted to SIAM Journal of Finance

ACK17
J.-Ph. Aguilar, C. Coste, J. Korbel, Series representation of the pricing formula for the European option driven by space-time fractional diffusion, 2017, arXiv:1712.04990, submitted to FCAA

BS73 
F. Black and M. Scholes, The pricing of options and corporate liabilities, 1973, Journal of Political Economy, 81, 637

BS94
M. Brenner and M. G. Subrahmanyam, A simple approach to option valuation and hedging in the Black-Scholes model, 1994, Financial Analysts Journal, 25-28

Bouchaud_Sornette94
J.-Ph. Bouchaud and D. Sornette, The Black-Scholes option pricing problem in mathematical finance: generalization and extensions for a large class of stochastic processes, 1994, J. Phys. I France, 4, 863–881

Calvet08
L. Calvet and A. Fisher, Multifractal Volatility: Theory, Forecasting, and Pricing. Academic Press Advanced Finance, Elsevier (2008).

Carr03
P. Carr and L. Wu, The Finite Moment Log Stable Process and Option Pricing, J. Finance  58 (2003), 753–778; doi:10.1111/1540-6261.00544.

Fama65
E. F. Fama,  The behavior of stock market prices, Journal of Business, 1965, 38, 34–105

Gerber93
H. Gerber, U. Hans and E. Shiu, Option Pricing by Esscher Transforms, HEC Ecole des hautes études commerciales (1993)

Gorenflo99
R. Gorenflo, Yu. Luchko. and F. Mainardi, Analytical properties and applications of the Wright function, Fract. Calc. Appl. Anal. 2, No 4 (1999). 383–414

Zatloukal14
H. Kleinert and V. Zatloukal, Green function of the double-fractional Fokker-Planck equation: Path integral and stochastic differential equations. Phys. Rev. E 88 (2013), Paper ID 052106; doi:10.1103/PhysRevE.88.052106.

KK16
H. Kleinert and J. Korbel, Option pricing beyond Black-Scholes based on double-fractional diffusion. Physica A 449 (2016), 200–214; doi:10.1016/j.physa.2015.12.125.

Korbel16
J. Korbel and Yu. Luchko. Modeling of financial processes with a space-time fractional diffusion equation of varying order. Fract. Calc. Appl. Anal. 19, No 6 (2016), 1414–1433; doi:10.1515/fca-2016-0073; 

Heston93
S. L. Heston, A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Rev. Financial Stud. 6, No 2 (1993), 327–343; doi:10.1093/rfs/6.2.327.

Mainardi01
F. Mainardi, Yu. Luchko and G. Pagnini,
The fundamental solution of the space-time fractional diffusion equation. Fract. Calc. Appl. Anal. 4, No 2 (2001), 153–192.

Mainardi10
F. Mainardi, A. Mura and G. Pagnini, The M-Wright function in time-fractional diffusion precesses: a tutorial survey, Int. J. Diff. Eq. 2010 (2010), Paper ID: 104505; doi:10.1155/2010/104505. 

Mandelbrot63
B. Mandelbrot, The variation of certain speculative prices, Journal of Business, 1963, 36, 394–419

Oksendal00
B. Øksendal, Stochastic differential equations: Introduction with applications, 2000, Springer (5th version).

Privault14
N. Privault, Stochastic Finance, Chapman & Hall, 2014

Taleb10
N. N. Taleb, The Black Swan: The Impact of the Highly Improbable, 2010, Random House Publishing Group

Tankov03
P. Tankov and R. Cont, Financial Modelling with Jump Processes, Chapman & Hall/CRC Financial Mathematics Series, Taylor &
Francis (2003).

Tsikh94
M. Passare, A. Tsikh and O. Zhdanov, A multidimensional Jordan residue lemma with an application to Mellin-Barnes integrals. In: 
Contributions to Complex Analysis and Analytic Geometry. Aspects of Mathematics E26. Vieweg+Teubner Verlag, Wiesbaden (1994), 233–241;
doi:10.1007/978-3-663-14196-9_8.

Tsikh97
M. Passare, A. Tsikh and A. A. Cheshel, Multiple Mellin-Barnes integrals as periods of Calabi-Yau manifolds with several moduli, Theor. Math. Phys.  109, No 3 (1997), 1544–1555; doi:10.1007/BF02073871.

Wilmott06 P. Wilmott, Paul Wilmott on Quantitative Finance, Wiley & Sons, 2006.

Zolotarev86 V. M. Zolotarev, One-dimensional stable distributions, 1986, American Mathematical Society.



