I'm trying to use the Stehfest algorithm to numerically invert a Laplace transform (as in Sepp and Skachkov). I was wondering if anybody has already implemented this algo as I have a couple of questions. Namely: how is it possible that I have totally different results for different number of nodes/weights I chose.

QuoteOriginally posted by: meteor..algorithm to numerically invert a Laplace transform . I was wondering if anybody has already implemented this algo as I have a couple of questions. Namely: how is it possible that I have totally different results for different number of nodes/weights I chose.I would expect this erratic behaviour, yes.The inverse Laplace transform is an example of a Fredholm equation of the 1st kind. This is an ill-posed problem (not well-posed in the Hadamard sense). Thus, small perturbations (0.001) in the input can lead to wild changes in the solutin (10^6). Is this he kind of behaviour you are getting.This kind of problem is well-known and Tikonov regularisation techniques are needed.I am not familiar with the specific algorithm you mention but it needs to take ill-posedness into account. Of course, I might be completely wrong. P.S. do you use even and odd nodes? e.g. N = 50 and N = 51. Any difference in output?

Thank you for your replyI'm using even number of nodes (since it is not possible to use odd number) with this algo which is briefly is presented here (pg 57): http://www.wilmott.com/pdfs/060103_sepp.pdfWhat do you mean by "it needs to take ill-posedness into account"?

The Stehfest method is a pretty stable algorithm for numerical invertion of a Laplace transform. Number of ponts, N, should be somewhere in the range 10-30. I personally use 16 or 14. Since Q_{j} represent large positive and negative numbers, one needs to take care implementing the code.Note there was a typo in formula (B.2) in Sepp-Skachkov paper: the last term in denominator should be (2k-j)! instead of (2k-1)!. But the VBA code given in the paper is ok, and using it one can invert the call option value under the Black-Scholes with accuracy up to 6 digits.The nice property of the Stehfest algorithm is that since constants Q_{j} depend neither on t nor on x they can be tabulated and subsequently used in inversion formula. Below I report values of Q_{j} for N=14 using 8-digit accuracy:j Q_{j}1 0.002777782 -6.402777783 924.050000004 -34597.927777785 540321.111111116 -4398346.366666667 21087591.777777708 -63944913.044444409 127597579.5500000010 -170137188.0833330011 150327467.0333330012 -84592161.5000000013 27478884.7666666014 -3925554.96666666Let me know if you need more details

Thank you very much for the info about the typo, this was the cause of my problems.I'm actually using this algorithm in the factor copula setting (ie to compute the density of the factor).

How does Stehfest algorithm compare to Gaver-Stehfest algorithm? The latter is claimed to be more efficient, but I haven't seen any performance analysis of these two algorithms.

Does someone have a copy or link to Stehfest's original paper - or indeed any paper that explains how and _why_ it works?

QuoteOriginally posted by: ZmeiGorynychDoes someone have a copy or link to Stehfest's original paper - or indeed any paper that explains how and _why_ it works?Algorithms is proposed in - D. P. Gaver Jr., Observing Stochastic Processes, and Approximate Transform Inversion, Operations Research, Vol. 14, No. 3. (May - Jun., 1966), pp. 444-459.Stehfest's implementation is in - Stehfest, H. 1970. Remark on algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13, 10 (Oct. 1970)Both are available online. The explanation of algo and its reasoning are in Gaver's work. It's not super precise from what i can understandQuoteOriginally posted by: J2How does Stehfest algorithm compare to Gaver-Stehfest algorithm? The latter is claimed to be more efficient, but I haven't seen any performance analysis of these two algorithms.i thought it's the same algorithm 368. are there two different algos?

Gents,may I recommend "our" inversion algorithm, which is pretty fast and stable ? ( I have checked the Gammas (!!) of a call very close to maturity @ ATM - no oscillations at all )http://www.cardano.nl/pagina/Numerical_ ... ture.pdfIt is explained, why it works in the paper.With regards,

I think the question where towards a (fast, simple) algorithmwhile my 'feeling' after scanning the paper is more towards amethod (but pls do not ask for a definition ...)What I do not see is, how singularities or piecewise stuff isactually handeled ... would you mind to sketch in brief again?The reason for the question:If integration is allowed I like the following: given F find areal constant a, such that F has no singularities in the complexhalf plane on the right side of a. Then inverse Laplace of F isf(t) = 2*exp(a*t)/Pi*Int(Re(F(a+u*I))*cos(u*t),u = 0 ... infinity).This follows from the proof of Abate & Whitt for 'Euler' and isa version of Bromwich's contour integral making 'everything' real.Hence one has to compute integrals with *constant* periodics tover the positive Reals and a standard solution is in Quadpack,for example contained in the GSL library as C code.The problem is to find the singularities automatically (and forthe choice of 'a' there is some trade off), especially if they arenot on the Reals.If possible some ready-to-run code for testing would be fine(but of course all would accept that's too much publicity).

What do you mean, too much publicity?"We got this terrific new algo, here's the code" is reproducible research and therefore a Good Thing."We got this terrific new algo, but we won't tell you the details of how it works" on the other hand would be a pretty pointless thing to say on this forum - if so, why put the paper here?Not so, GogolaAnita?

To some extent, yes - but there's always the issue that if you hand out the code to too many people, you're giving a starter's advantage to other academics/practitioners to come up with improvements. I do believe however that a paper should be self-explanatory and reproducible, so that anyone can code it up. But not all clever tricks/details have to be shared - at the end of the day, that's our job.Same dilemma in university - if students are doing a project, as the university you'd want them to code up a model themselves, and understand what they're doing, not for them to surf the web, download the best implementation of a model (or bug an author of a paper for their code), run it and report the results.

Two possibilities: either you're promoting your reseach as a commercial enterprise, looking for customers or at least to raise your profile as service provider. A fair thing to do, and in that case, you certainly won't give full disclosure, but these forums are not really the place for that kind of thing.Another possibility is you advertise your research as an academic - likewise a fair thing to do, but in that case I would expect full disclosure upon request, as reproducibility is an essential part of the scientific method. This is quite a fundamental principle of science and therefore academe I would say.Can't say I agree with your point about students either - if an implementation is freely available online, that a) gives the students who try to replicate it a benchmark to compare against (and checking whether they wrote their own, if you even care, is easy enough); b) allows them to play with several models without having to code them all up. Anyway, what matters is understanding (and that is testable), not how easy it is to grade homework (well that's a separate discussion right there). In any case, standards of scientific discourse do not in my opinion have anything to do with convenience of educating students.Not so?

> Two possibilities: either you're promoting your reseach as a commercial enterprise, looking for customers or at least to raise your> profile as service provider. A fair thing to do, and in that case, you certainly won't give full disclosure, but these forums are not really> the place for that kind of thing.Maybe these forums aren't, but I'm sure there are people writing papers with this kind of motivation - surely they would not disclose all their code?> Another possibility is you advertise your research as an academic - likewise a fair thing to do, but in that case I would expect full> disclosure upon request, as reproducibility is an essential part of the scientific method. This is quite a fundamental principle of science> and therefore academe I would say.Absolutely - reproducibility is essential, I fully agree on that. It's just my point of view that reproducibility does not necessarily require you to dish out all your code on request. I believe it's sufficient that you can adequately describe the algorithm such that anyone is able to implement it. But there are pro's - if you do disclose your code, more people might end up using it, just because it is readily available.> In any case, standards of scientific discourse do not in my opinion have anything to do with convenience of educating students.Point taken.

Well, I think the authors have a commercial employer and thus it is noteven obvious he agrees they publish. However if they only have done acommercial implementation it would be astonishing to find that code inpublic.So my feeling for the discussion is: it is honest to 'demand' academicpublicity (like the given link suggests), but it might not be aware ofindustrial needs.For students (and their profs) I would not care: either they have to doexercises (then a very long or sophisticated solution is not plausible)or they do some final project/thesis/... and then it would be too riskyfor them to be dishonest, as they have to aware that the prof in chargeknows the field quite well