Inverse Laplace transform

Question

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.

Answer 1

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

Compile`CompilerFunctions[] // 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.

Answer 2

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.