I have only a little question.My problem is to compute the inverse of the cumulative normal distribution. You say that there are some methods, Boris Moro's and Peter Acklam's.Could you tell me please where I can find something about the implementation of the 2 algorithms?thank youandrea

Here is a inverse gaussian in 1-d which might not be the best, but I find it is pretty much OK. /*.....................................................**** function ppnd() **** Applied Statistics algorithm AS111, **** Beasley and Springer **** compute the inverse distribution function of a **** standard normal distribution that is set to the **** value of z satisfying P(Z<z) = p. therefore, **** 0 < p < 1. ifault = 0 means that all is well. ****.....................................................*/double ppnd( double p, int *ifault ){ double ppnd1, q, zero = 0., split = 0.42, half = 0.5, one = 1.; double r, a0 = 2.50662823884, a1 = -18.61500062529, a2 = 41.39119773534, a3 = -25.44106049637, b1 = -8.47351093090, b2 = 23.08336743743, b3 = -21.06224101826, b4 = 3.13082909833, c0 = -2.78718931138, c1 = -2.29796479134, c2 = 4.85014127135, c3 = 2.32121276858, d1 = 3.54388924762, d2 = 1.63706781897; *ifault = 0; q = p - half; if( fabs(q) > split ) { r = p; if ( q > zero ) r = one - p; if ( r <= zero ) { *ifault = 1; ppnd1 = zero; return(ppnd1); } else { r = sqrt(-log(r)); ppnd1 = ((( c3 * r + c2 ) * r + c1 ) * r + c0 ); ppnd1 = ppnd1 / (( d2 * r + d1 ) * r + one ); if ( q < 0 ) ppnd1 = -ppnd1; return(ppnd1); } } else { r = q * q; ppnd1 = q * ((( a3 * r + a2 ) * r + a1 ) * r + a0 ); ppnd1 = ppnd1 / (((( b4 * r + b3 ) * r + b2 ) * r + b1 ) * r + one ); return( ppnd1 ); }}

thank you!bye

I have some simulation results for you.Ordinary Monte Carlo squares off against the Sobol sequence.I would particurlarly like to hear from others what they think of the accuracy.A European call on a constant volatility asset was priced (S(0)=50, sigma=0.4, T=1)for 10 different strikes K=40,....,90.The accuracy declines markedly as the call moves more and more out of the money.N=1000 asset paths were generated.The problem can be set up as a problem in any dimension (the number of time steps to expiry).Graphs are porvided for dimensions d=1,2,3,4.....,17, 60,300,1600,2000;This simulation was not tweaked at all (no variance reduction or control variate techniques)How much better are you doing?

trc, why don't you test your Sobol' generator before using it for modelling ( L2 discrepancy in n dimensions for the first 2^16 points is a good one, see also P. Jaeckel's book ). It looks as if generator you are using is not a good one. You can get a noncommercial version of Sobol' generator from www.broda.co.uk

Broda only distributes a DLL for Windows. I am on Linux.Does Broda publish discrepancies for its Sobol sequence?What sort of accuracy would you expect?

trc, Yes, there are such results ( obtained by Peter Jaeckel. Probably, you can aslo find them in his book )

Audetto,in situations like these, please first try your favourite search engine. As for Boris Moro's Chebyshev expansion, you find it at http://www.mathfinance.de/FF/cpplib.html(found via typing "Boris Moro" into Google). As for Peter Acklam, whose web site has moved since it was mentioned in my book,Google found http://home.online.no/~pjacklam/notes/i ... index.html for me.Regards,pj

trc,I don't have the time to look at your data (or to download them, for that matter),but it looks like you would actually benefit from buying or borrowing my book. Thereare key points about the use of Sobol' numbers:- 1.) Construction and initialisation of Sobol' numbers. Even number theoreticians still struggle with Sobol' 197x paper(s) on the additional uniformity properties provided by careful initialisation. 2.) Effective dimensionality reduction. Use the Brownian bridge as the best compromise betwen effective dimensionality reduction and construction speed. You get a Brownian bridge class with my book.You mention up to 2000 dimensions for Sobol' numbers. Where did you get thecode to produce up to 2000 dimensions? To the best of my knowledge, the onlypublished source of that many (or more) primitive polynomials modulo two, whichare of paramount importance for Sobol' numbers, is my book. Note that particularlyfor high dimensions the difference between primitive and irreducible (and only thelatter is the equivalent to 'prime' for ordinary numbers) becomes important. You geta few million of those with my book.I am terribly sorry about the continuous advertising of my book, but it just so happensthat you ask all the questions that I had to struggle with in the past and whose solution,or at least the best I could find, I put down in that book. Honestly.Regards,pj

PJ,I have your book and I like it.<<I am terribly sorry about the continuous advertising of my book>>I will gladly advertise the book for you and can wholeheartedly endorse it.Anyone who has any interest at all in Monte Carlo methods should buy it.The code on the CD is very valuable and will save everyone a huge amount of work.<< Where did you get the code to produce up to 2000 dimensions?>>In fact I did code up the Sobol generator using your primitive polynomials and the table ofinitialization numbers (these only up to dimension 13) and generally following the algorithm in the book.I don't know why Dostojewski believes the Sobol generator is not good.I checked it in the following admittedly unscientific manner: I plotted projections on 2 dimensional planesparallel to coordinate axes for a number of dimensions (i,j) up to dimension 3000 and never saw any structureor regularity. I was very happy with that.The results of my simulations are that ordinary Monte Carlo has anyhwere from 5-90% probabilityof beating the Sobol sequence depending on dimension and also the nature of the function being integrated.Interestingly the nature of the function plays a big role. The problem was to value a vanilla European callon the usual constant volatlity asset. The results were highly variable with the call strike.I have a feeling that in more than one dimension the transfomation uniform x -> multinormal z plays a significant role in the accuracy. It came out that coordinatewise inverse CDF is much worse than Box-Muller on pairsof coordinates. Dostojewski,I assume your belief that the Sobol generator is not good results from your dissatisfation with the acuracy.There could be other issues such as the transformation from uniform to multinormal.I would like to hear what sort of accuracy you and others want to get (or are in fact getting) out of a path computation that simulates only 1000 paths Are you looking for relative errors1. less than 0.1%2. less than 1/2%3. less than 5%Since I coded up the Sobol generator following Peter Jaeckel's book (with optimized initialization numbers up to dimension 13 not 32 as in the book)I would assume that my L2-discrepancy should be similar to his. That's why it would be nice to have some other standards to compare against.

Ok, trc, I did look at your file after all.First, let me reiterate: NEVER USE INCREMENTAL PATH GENERATION, PLEASE. It does make a massive difference if you assign the least discrepancy to the most important dimensions. For Brownian motion, the best compromise method is the Brownian bridge.Of course, the discrepancy advantages of low-discrepancy numberstail off in very high effective dimensions. We practically never have very higheffective dimensions in finance, only very high apparent dimensions.If you try to compute the volume of a 1200-dimensional unit sphere incartesian coordinates, low-discrepancy numbers will be of no use becauseall dimensions have equal importance. Mind you, this problem is so difficult thatyou'll hardly like the results that you get from conventional Monte Carlo either,so perhaps this statement is a bit harsh. However, if you transform this problemto polar coordinates, low-discrepancy numbers will work like a charm.Low-discrepancy numbers are not a panacea, though.As for your statement that the Box-Muller transform is better than theinverse cumulative normal, I am completely at loss what you base this on.I cannot see any evidence of it in your diagrams, just a random set ofsome of the black curves being higher for Box-Muller, and some of thembeing higher for N^-1. As we know, the Box-Muller method is a variatetransformation, and thus guaranteed to generate the exact distributionwithin the usual numerical error propagation limits. The inverse cdf methodof Peter Acklam also gives the correct distribution within the numerical limitsgiven by the machine accuracy. Both methods produce an accuracy wellbeyond the numerical resolution of your experiments. They do, of course,result in different sets of numbers, and thus you could only really make astatement about their relative benefits by having an incredibly large numberof individual simulations such that you can make statistically relevant statementsabout some assumed null-hypothesis. Luckily, this particular experiment (primarilydue to the use of the Mersenne-Twister) is not susceptible to the Neave effect,but with number generators like Ran0 or IBM's original RANDU, you might havesome serious mispricing in the wings (for far out of the money options).As for the worse performance for out-of-the money options, this is also a wellknown effect: In this case, the yield of your simulation, i.e. the ratio of the numberof paths that result in a non-zero contribution to the average is very low. If you knowbeforehand that your running this kind of risk, use the importance sampling methodin addition to other variance reduction techniques. Horses for courses. Still, I have neverseen MC outperform low-dimensional Sobol'-QMC even for extremely far out of the moneyoptions precisely because of their regularity features. Sobol' numbers tend to sample theextreme events more reliably, albeit, of course, with the same average frequency.I hope this helped.Regards,pj

Wow, that's what I call a quick turnaround.First of all: thank you very much for your kind endorsement.As for the Brownian bridge: If I am not mistaken, it will actuallyrender all of the intermediate steps as irrelevant since yourexperiment is to price plain vanilla options (is this right?). Wether you have2000 intermediate steps or none, any one path will end in the samepoint. The entire experiment is reduced to its 1-dimensional nature.You may say this is cheating, but the identification of the most importantdimensions does make a tremendous difference. For almost all financialsimulations involving Wiener processes I found it beneficial to first constructthe underlying Wiener path(s) using the Brownian bridge, and then possiblychopping it apart again if I actually need increments (because I might beEuler-integrating some nasty SDE driven by an underlying Wiener process).And yes, perhaps there is a little more roundoff error in your multivariatetransformation to Gaussians. Did you use the method posted earlier in this streamor Boris Moro's or Peter Acklam's ?As for your relative error question, the short answer is: this depends on the problem.I'd never use the result of a simulation with 1000 paths for anything that has anythingto do with my company's money. For the pricing of exotic deals with computation of Greekswe hardly ever trust anything with less than 2^20 - 1 paths. You do always use a Mersennenumber of paths with Sobol' sequences, don't you?Regards,pj(signing out for tonight, this is London and the City is expecting me for work in a few hours...)

PJ,If you are still there could you please take a look at my message directed to you on the thread about your bookand respond to me at spyqqqdia@yahoo.com at your leisure.I believe the graphs really do show Box-Muller beats coordinate wise inverse normal CDF.Look at the blue curves (relative accuracy of the Sobol price).Also if the black curve was higher Monte carlo had a better chance of beating the accuracy of the Sobol priceand that can only happen if the accuracy of the Sobol price was worse (algorithm details, you have to believe me on this).Yes, the Brownian bridge will reduce this to a one dimensional problem.I am mostly interested in Libor market models where the dimension might be anythingbetween 10 and 120.I'll play around with this some more and will post how much better importance sampling fares than unrefined QMC.Thanks,

trc,I'll see what I can do. Of course, one would naturally wish to look at all possible causes for such an unexpected difference, andthe definition of 'higher probability of default' and the means to assess this probability would be part of the complete picture,but, as I say, for incremental path generation with many steps I can believe that Sobol' numbers can appear to be equal or worsethan pseudo-random numbers. Use the Brownian bridge. The class given in the book might just work for you, give it a go!Regards,pj

Congratulations Dostojewski,I suspect a problem with my Sobol generator too.Quite possibly I made a mistake implementing the generator.I will check everything thoroughly.Are you an expert on the subject?