Help with calculating a complicated integral

Edit: The question is solved : )) Thank you for your attention !

I am sorry for mass questions, but the question bothers me a lot. So I would like to ask this question here.

This problem arises when I tried to compute the valua of Asian call otions using the inverse Laplace transform method. The method I use here is very new since it computes the Asian call options in very new way ! george2079 asked me whether the theoretical part of the"new method" has a problem. I would like to confirm that: two of us have checked the theoretial part more than 10 times and we confirm that it is correct.

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

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

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

Define $F$ as follows: $$F(\alpha)=\frac{1}{\alpha} \frac{\Gamma(\eta+2\mu/\sigma^2)}{\Gamma(2\eta +2\mu/\sigma^2)} \int_K^{\infty}\left(\frac{2S_0}{y\sigma^2}\right)^{\eta}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 a confluent hypergeometric function

($M$ is the Hypergeometric1F1 function in Mathematica).

I am trying to find the inverse Laplace transform of $F(\alpha)$ numerically.

$$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$.

Note that if we let $b=\left(\frac{2S_0}{\sigma^2}\right)$, and $z=b/y$ then we have $$F(\alpha)=\frac{1}{\alpha} \left(\frac{2S_0}{\sigma^2}\right)\frac{\Gamma(\eta+2\mu/\sigma^2)}{\Gamma(2\eta +2\mu/\sigma^2)} \int_0^{\frac{b}{K}}z^{\eta-2}M(\eta, 2\eta+\frac{2\mu}{\sigma^2}, -z) \, dz$$

So using the ideas suggested by Dr. Hintze and I_Mariusz and the package written by Prof Valko, I write a Mathematica program to compute it. Please see the Mathematica code here

Clear["Global*"]

(*The program used to find the Inverse Laplace transform numerically*)

GWR[F_, t_, M_: 32, precin_: 0] :=
Module[
{M1, G0, Gm, Gp, best, expr, τ = Log[2]/t, Fi, broken, prec},
If[precin <= 0, prec = 21 M/10, prec = precin];
If[prec <= $MachinePrecision, prec =$MachinePrecision];
broken = False;
If[Precision[τ] < prec, τ = SetPrecision[τ, prec]];
Do[Fi[i] = N[F[i τ], prec], {i, 1, 2 M}];
M1 = M;
Do[
G0[n - 1] = τ (2 n)!/(n! (n - 1)!) Sum[
Binomial[n, i] (-1)^i  Fi[n + i], {i, 0, n}];
If[Not[NumberQ[G0[n - 1]]], M1 = n - 1; G0[n - 1] =.;
Break[]];
, {n, 1, M}];
Do[Gm[n] = 0, {n, 0, M1}];
best = G0[M1 - 1];

Do[
Do[
expr = G0[n + 1] - G0[n];

If[Or[Not[NumberQ[expr]], expr == 0], broken = True; Break[]];
expr = Gm[n + 1] + (k + 1)/expr;
Gp[n] = expr;
If[OddQ[k],
If[n == M1 - 2 - k, best = expr]
];
, {n, M1 - 2 - k, 0, -1}];
If[broken, Break[]];
Do[Gm[n] = G0[n]; G0[n] = Gp[n], {n, 0, M1 - k}];
, {k, 0, M1 - 2}];
best
]

(*Model parameters*)
mu = 0.15;
sigma = 0.20;
S0 = 100;
K = 95;
T = 1;

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

F[alpha_] := F[alpha] = int2[alpha]/alpha;
GWR[F, 1]
TimeUsed[] (*Time used to compute the value*)


13.4836
2.184

Before I ask my question, I would like to make a remark that the routine: GWR[] has been tested with different functions and it gives correct answers up to 15 decimal places for those tested functions. It gives me the answers in less than a minute.

However, when I tested with F[], it gives me a wrong answer:

• when K=95, it gives me 13.4836. The correct anwer is 15.643.

• when K=105, it gives me 7.94016. The correct answer is : 3.560

It gives me the answers in less than 1 minute.

So here is my question:

Is there any thing which is bad for the function F[alpha_]` that I don't know? If there is such a bad singulariy, how to remove it.

• I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. – m_goldberg Dec 17 '15 at 13:51
• @m_goldberg -- I think the question has merit. Wouldn't the OP and this site better benefit from improving the question rather then closing it. The OP has done a lot of work, trying to solve his problem. That alone merits praise. As it stands the question seems focused on the problem of translating the computation into code. I recommend the poster analyze the code and step through the calculations to identify the specific part of his code where the calculation goes wrong. This may point him to a Mma problem or a translation problem. – Jagra Dec 17 '15 at 14:07
• @DuyNguyen -- Note how I formatted the functions and output you describe in the last section of your post so that it becomes clearer. It now better shows what comes from theory and what you've done in Mathematica. Even this little bit of formatting can help the participants in the forum, focus better on the parts of the problem where in they can best help. – Jagra Dec 17 '15 at 14:31
• Maybe This be better: arxiv.org/pdf/1304.2505.pdf – Mariusz Iwaniuk Dec 17 '15 at 18:34
• "the singularities of Hypergeometric1F1" - the Kummer function is holomorphic (no poles or branch cuts), so your bet is way off in this case. – J. M. will be back soon Dec 18 '15 at 1:48