This sheet shows, how one can extract a "risk neutral density" from prices of traded options: theoretical one can price options through a probability function (describing the stock's behavior) and if prices would be given 'in a continuous way' one should recover that functions to get information on the stock (differentiating prices twice w.r.t. strike). The first step is a quite brute way: interpolate extracted volatility by a polynomial of low degree to get the observable part of such an density. Now extend to the unobservable left and right tails by selecting a Black-Scholes situation with reasonable parameters.

If one only wants numerical values in double precision for real arguments then Maple is quite slow, even for the quite typical case of index 1 or 0. A way is to use the GNU Scientific Library GSL as external source for this special cases or to compile the needed stuff into a DLL and use that through external calls. However I (once) wanted to have that as Maple code and translated Ooura's solution into Maple (which is not too complicated). The motivation was to have something fast for evaluating integrals over such Bessel functions. Putting the coefficients needed for the numericial approximation into global arrays this makes Maple as fast as GSL (almost, M9.5 did, M11 is somewhat slower for this approach) for this functions, if one uses hardware float evaluations 'evalhf' (=double precision).

/ ranting on I use Win XP home on an Acer, which has an Athlon (do not care for specifica). It installs with a blank space in the directory name as default (which may cause troubles if using Open Maple with MS compilers). Installing Watcom not even gives a warning to prior installed ones (so I hope it does not matter ... say having some modifications in it [well, a backup is never a bad idea before any install ...]). It does not adapt settings from M10 (I think 9.5 -> 10 did so), i.e. global settings and style sheet had to be done again. And: it does not associate M11 to *.mws (classical sheets are my default).

In a posting posting at http://www.math.utexas.edu/pipermail/maxima/2006/000126.html Fateman cites Gosper with an interesting approach to compute the hypergeometric function 2F1. I used that to produce a compiled version (double precision only), which can be achieved from Excel for example. The approach however can be used within Maple as well. The idea is, that with his way the linear transformations given in A&S need to be applied only once to cover the whole complex plane: up to 1 transformation one can use either Gosper's recursion or even lives in the unit circle with radius 1/2 (where the hypergeometric series already converges quite fast).

For those interested in financial Math: A classical in mathematical finance is evaluating option prices by binomial trees. This has many advantages (like easy, but coarse results for American options), but it is well known to quite inaccurate for various reasons - even for European options. The 'best' known improvement is due to Leisen-Reimer. They all suffer from low order convergence towards a continuous model. The standard reference model is CCR, the Cox-Ross-Rubinstein tree.