# Inverse Laplace transform

Let $r=\mu = 0.15; \sigma = 0.05; T = 1; S_0 = 100; K = 95;$

Let $\nu:=\frac{2\mu}{\sigma^2}-1$

and $\eta \equiv\eta(\alpha):=-\frac{\nu}{2}+\frac{1}{2}\sqrt{\nu^2+\frac{8\alpha}{\sigma^2}}$.

Define $u$ as follows: $$F(\alpha)=\int_K^{\infty}\left(\frac{2S_0}{y\sigma^2}\right)^{\eta} \frac{\Gamma(\eta+2\mu/\sigma^2)}{\Gamma(2\eta +2\mu/\sigma^2)} M(\eta, 2\eta+\frac{2\mu}{\sigma^2}, -\frac{2S_0}{y\sigma^2}) \, dy$$ where $\Gamma$ is the gamma function and $M$ is the hypergeometric function ($M$ is the Hypergeometric1F1 function in Mathematica).

I am trying to find the inverse Laplace transform of $F(\alpha)$ $$V(t)=\lim_{z\to\infty}\frac{1}{2\pi i}\int_{\gamma- iz}^{\gamma +iz} e^{t\alpha}F(\alpha)d \alpha$$ at $t=1$.

Here is what I tried

Clear[S0, r, T, K, mu, sigma]
r = 0.15;
sigma = 0.05;
T = 1;
S0 = 100;
K = 95;
mu = r;
v0 = S0;
nu = 2*mu/sigma^2 - 1;
k1 = 8/sigma^2;
sigma2mu = 2*mu/sigma^2;

F[alpha_] := NIntegrate[
(2*v0/(a*sigma^2))^(-0.5*nu + 0.5*Sqrt[nu^2 + k1*alpha])*
Gamma[(-0.5*nu + 0.5*Sqrt[nu^2 + k1*alpha]) + sigma2mu]*
Hypergeometric1F1[(-0.5*nu + 0.5*Sqrt[nu^2 + k1*alpha]),
2*(-0.5*nu + 0.5*Sqrt[nu^2 + k1*alpha]) + sigma2mu,
-2*v0/(a*sigma^2)]/
Gamma[2*(-0.5*nu + 0.5*Sqrt[nu^2 + k1*alpha]) + sigma2mu],
{a, K, 10000}];


Edit: Here What I added after the Hint from Dr. Hintze. I use the Fourier series method to find the inverse numerically.

A = 18.4;
n = 38;
m = 11;
S = ConstantArray[0, m + 1];
B = 1;
For[k = 0, k <= m, k++,
S[[k + 1]] = Exp[A/2]/(2*B)*Re[F[A/(2*B)]];
For[j = 1, j <= n + k, j++,
S[[k + 1]] =
S[[k + 1]] + Exp[A/2]/B*(-1)^j*Re[F[(A + 2*j*Pi*I)/(2*B)]];
] ]
f = 0
For[k = 0, k <= m, k++,
f = f + Binomial[m, k]*S[[k + 1]]*2^(-m);
]
f


Still:It seems take forever to find V[1]. Using a different method, I know that exact answer is $11.094$. Could some one tell me where I am wrong ? Or a better way to compute V[1]? Thank you so much for your time. I truly appreciate it.

• (at)Dave Nguyen (sorry, at the moment the at-symbol is not accepted here): you can do the y-integral exactly in terms of HypergeometricPFQ, but I think you have to do the inverse Laplace transform numerically. – Dr. Wolfgang Hintze Dec 14 '15 at 10:06
• Dr. Hintze : Thank you for your time. I updated the question by addiing your comment. You meant I can find the closed-form of $u(\alpha)$? I used mathematica and cannot get the function HypergeometricPFQ. Could you please give me more details ? Thank you for your time ! – D. Nguyen Dec 14 '15 at 12:43
• (at)Dave Nguyen: I have explained it in an answer. Where did this problem arise from? – Dr. Wolfgang Hintze Dec 14 '15 at 12:50
• Dr. Hintze: Thank you for your answer. I truly appreciate it. The question comes from finding the Asian options using Laplace transform. Can you give more detail about a, b,c,p,q, and k ? did you do it by hand ? – D. Nguyen Dec 14 '15 at 12:51
• The parameters are defined through your own expression for u. For instance b = 2 S0/sigma^2 and in this manner you can identify them easily by yourself. The point for me here was that they do not depend on y, so that the y-integral could be done as shown. – Dr. Wolfgang Hintze Dec 14 '15 at 12:57

ClearAll["Global*"]; Remove["Global*"];

eq = (2*S0/(y*sigma^2))^nu*(Gamma[nu + (2*mu)/sigma^2]/
Gamma[2*nu + (2*mu)/sigma^2])*Hypergeometric1F1[nu, 2*nu + (2*mu)/sigma^2, -2*S0/(y*sigma^2)];
int = Integrate[eq, {y, K, Infinity}, GenerateConditions -> False]
int2 = int /. nu -> (-vu/2 + Sqrt[vu^2 + 8*alpha/sigma^2]/2) /.vu -> (2*mu/sigma^2 - 1)

mu = 15/100;
sigma = 5/100;
S0 = 100;
K = 95;
T = 1;

fun[alpha_] := int2;


Let's try a Stehfest.H method to calculate inverse Laplace transform:

Vi[n_, i_] :=
Vi[n, i] = (-1)^(i + n/2) Sum[
k^(n/2) (2 k)! /( (n/2 - k)! k! (k - 1)! (i - k)! (2 k -
i)! ), { k, Floor[ (i + 1)/2 ], Min[ i, n/2] } ] // N
NStehfest[F_, s_, t_, n_] := NStehfest[F, s, t, n] =
Log[2]/t Sum[ Vi[n, i]* F /. s -> i Log[2]/t , {i, 1, n} ] // N
n=16;(*n=16 is a optimal value and must be even*)
V1[t_] := NStehfest[fun[alpha], alpha, t, n]; V1[T]
(*Out : 93.4984*)


Let's try a second Durbin, F method:

NDurbin[F_, s_, t_, n_, a_, e_] :=
NDurbin[F, s, t, n, a, e] = Module[{k}, ((1/4) Exp[a t + e]/
t) ( ((1/2) Re[ F] /. s -> a + e/t) +
Sum[Re[ Exp[(1/4) I k Pi] (F /.
s -> a + e/t + (1/4) I k Pi/t )], {k, 1, n}] )] // N
n=200;(*increase n  to improve accuracy,default is n=2000*)
V2[t_] := NDurbin[fun[alpha], alpha, t, n, 0, 2]; V2[T]
(*Out: 116.31*)


Lets try Week W.T method:

 Tn[n_, T_] := Tn[n, T] = T / n // N
sig[T_, a_] := sig[T, a] = If[ (a + 1 / T) >= 0, a + 1 / T, 0 ] // N
theta[n_, j_] := theta[n, j] = 0.5 Pi (2 j + 1) / (n + 1) // N
Hj[F_, s_, n_, T_, a_, j_] :=
Hj[F, s, n, T, a, j] = ( 0.5 / Tn[n, T] ) (
Re[ F ] - Cot[ 0.5 theta[n, j]] *Im[ F ] ) /.
s -> sig[T, a] + ( 0.5 I /Tn[n, T] ) Cot[0.5 theta[n, j]] // N
Ak[F_, s_, n_, T_, a_, k_] := Ak[F, s, n, T, a, k] =
2 / (n + 1) Sum[
Hj[F, s, n, T, a, j] Cos[ k theta[n, j] ], {j, 0, n} ] // N
Ak[F_, s_, n_, T_, a_, 0] := Ak[F, s, n, T, a, 0] =
1 / (n + 1) Sum[Hj[F, s, n, T, a, j], {j, 0, n} ] // N
NWeeks[F_, s_, t_, n_, a_] := NWeeks[F, s, t, n, a] =
Exp[ ( sig[ t, a ] - ( 0.5 / Tn[n, t] ) ) t]*
Sum[ Ak[F, s, n, t, a, k] LaguerreL[ k, t / Tn[n, t] ], {k, 0,
n} ]  // N

n=200;(*increase n  to improve accuracy,default is n=100*)
V3[t_] := NWeeks[fun2[alpha], alpha, t, n, 0]; V3[1]
(*Out: 115.394*)


Lets try GWR package thanks to xzczd user:

 << NumericalLaplaceInversion
<< FixedTalbotNumericalLaplaceInversion
funG[alpha_] = int2;(* or: funG[alpha_] = int2/alpha. See below !!! *)
N[GWR[funG, T, 31], 20]

(*Out: 116.6460181143383387 *)


OP answers is 11.094 , but mine is about 116.The first Stehfest.H method is not very accurate. I do not know why but the result is different.I checked a couple of times but did not find any personally mistakes.

EDITED 3:

When changing to fun[alpha_] := int2/alpha; because OP wanted to,we have:

1. Stehfest.H method at t=1and n=16 is: 13.3339
2. Durbin, F method at t=1and n=100 is: 12.8816
3. Week W.T method at t=1and n=100 is: 12.3234
4. GWR package at t=1and n=31 is: 13.0424

APPENDIX:

Let's try simple case at t=1: $$f(s)=\frac{1}{(s-1)^2+1}$$

f[s_] := 1/((s - 1)^2 + 1)
Re[InverseLaplaceTransform[f[s], s, t = 1]] // N
(*2.28736*)


A compiled Durbin method (is 25 time faster than no-compiled method).

But here is the problem.Can't get speed up with this "Hypergeometric1F1Regularized" function,because it can not compilable.It very slow compared with no-compiled version in some cases

CompileCompilerFunctions[] // Sort


gives you a list of compilable functions, Hypergeometric1F1Regularized is not one of them.

In Module function is there a e parameter,typically is 1 < e < 4. I use e=2.Parameter a is typicaly 0, but not always.a is optimal number greater than the real parts of all singularities of F(s) if any.We note in particular that the function f[s] has a pole at s = 1. The proper inversion requires the use of a = 1 Parametr a can be calculated with:

Re@TransferFunctionPoles[
TransferFunctionModel[Unevaluated[{{f[s]}}], s,
SamplingPeriod ->None, SystemsModelLabels -> None]]
(*{{1,1}}*)


Durbin's method complied version:

NInverseLaplaceT =
Compile[{{t, _Real}, {n, _Integer}},
Module[{e = 2, a = 1,
k}, ((1/4) Exp[a t + e]/ t) ( ((1/2) Re[f[s]] /. s -> a + e/t) +
Sum[Re[ Exp[(1/4) I k Pi] (f[s] /.
s -> a + e/t + (1/4) I k Pi/t )], {k, 1, n}] )],
RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"];

t = 1;
n = 50000;
NInverseLaplaceT[t, n]
(*2.28736*)


Outputs is the same.

Sources:

(1) Durbin, F., Numerical inversion of Laplace transforms: an
efficient improvement to Dubner and Abate's method'', Computer J., v. 17, 371-376, 1974.

(2) Stehfest, H., Numerical inversion of Laplace transforms'', Comm.
ACM, v. 13, 47-49 and 624, 1970.

(3) Weeks, W.T., Numerical inversion of Laplace transforms using
Laguerre functions,'' J. ACM., v. 13, 419-426, 1966.

A pdf file Approximate Inversion of the Laplace Transform in this book provided five approximate inversion algorithms (Stehfest, Papoulis, Durbin-Crump, Weeks, Piessens).You can find a Mathematica package here.

• Dear I_Marius: Thank you SO MUCH for your work. I really appreciate it. How long did it take you to get the answer ? the answer 11.094 is correct, i checked it using PDE method and Monter carlo method. – D. Nguyen Dec 14 '15 at 16:55
• few hours :) strange results..... – Mariusz Iwaniuk Dec 14 '15 at 17:14
• : )) Thanks for your time! – D. Nguyen Dec 14 '15 at 17:22
• It works. You have to change from e=1 to e=2,(e is a special parameter-typically 1<e<4),like this :NDurbin[fun[alpha], alpha, t, n, 10^-5, 2]. – Mariusz Iwaniuk Dec 16 '15 at 12:02
• Here is another package for numerical Laplace inversion. – xzczd Dec 18 '15 at 3:25

I was asked to explain how I arrived at the statement of my comment, regarding the first step of the complete solution. Here it is:

The y-integral is of the form

u0 := c Integrate[y^-a Hypergeometric1F1[p, q, -b/y], {y, k, \[Infinity]}]


where a, b, c, p, q, and k are parameters independent of y.

Letting

y = b/z, dy = - b/z^2 dz


u0 becomes

u1 := c1 Integrate[z ^ a1 Hypergeometric1F1[p, q, -z], {z, 0, k1}]


with some other constants c1 and k1, derived from the previous ones.

Explicitly we have

p = eta
q = 2 eta + 2 mu/sigma^2
a = eta
b = 2 S0/sigma^2
c = (2 S0/sigma^2)^eta Gamma(eta + 2 mu/sigma^2)/Gamma(2 eta + 2 mu/sigma^2)


and, in u2 these are

a1 = a - 2
c1 = c*b^(eta-1)
k1 = b/K


The dependence on alpha is via

eta = - nu/2 + nu/2 Sqrt(nu^2 - 8 alpha/sigma^2)


It appears in all parameters except a1 and k1.

Now Mathematica can do the z-integral

u1

(*
Out[14]= ConditionalExpression[(
c1 k1^(1 + a1) HypergeometricPFQ[{1 + a1, p}, {2 + a1, q}, -k1])/(1 + a1),
Re[k1] > 0 && Re[a1] > -1]
*)


For the values given in the OP the conditions are fulfilled. Hence we have

u2 = u1[[1]]

(*
Out[15]= (c1 k1^(1 + a1) HypergeometricPFQ[{1 + a1, p}, {2 + a1, q}, -k1])/(1 + a1)
*)


The dependence of u2 on $\alpha$ is given via the dependence of the parameters on $\alpha$ which is rather complex.

I don't believe that a simple InverseLaplaceTransform[]` ($\mathcal{L}_{\alpha }^{-1}[u](t)$) will work.

Hence the next step is to analyze the singularities of u2 in the complex $\alpha$-plane. This then permits to determine $\gamma$, i.e. to locate the complex path integral to the right of all those singularities, and to end up with an integral between $-\infty$ and $\infty$.

This then should be attempted numerically, the more as it is only required at t = 1.