7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ (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. $\endgroup$ – Dr. Wolfgang Hintze Dec 14 '15 at 10:06
  • $\begingroup$ 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 ! $\endgroup$ – Duy Nguyen Dec 14 '15 at 12:43
  • $\begingroup$ (at)Dave Nguyen: I have explained it in an answer. Where did this problem arise from? $\endgroup$ – Dr. Wolfgang Hintze Dec 14 '15 at 12:50
  • $\begingroup$ 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 ? $\endgroup$ – Duy Nguyen Dec 14 '15 at 12:51
  • $\begingroup$ 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. $\endgroup$ – Dr. Wolfgang Hintze Dec 14 '15 at 12:57
10
$\begingroup$
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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Duy Nguyen Dec 14 '15 at 16:55
  • $\begingroup$ few hours :) strange results..... $\endgroup$ – Mariusz Iwaniuk Dec 14 '15 at 17:14
  • $\begingroup$ : )) Thanks for your time! $\endgroup$ – Duy Nguyen Dec 14 '15 at 17:22
  • 1
    $\begingroup$ 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]. $\endgroup$ – Mariusz Iwaniuk Dec 16 '15 at 12:02
  • 1
    $\begingroup$ Here is another package for numerical Laplace inversion. $\endgroup$ – xzczd Dec 18 '15 at 3:25
8
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.