Bivariate Normal Integral paper

jambodev · March 30th, 2011, 11:08 pm

HiAnyone has Drezner, and Wesolowsky's 1989 paper: "On the Computation of the Bivariate Normal Integral" that are willing to share?appreciate it.

drone · July 18th, 2011, 11:44 am

You could try using Alan Genz's - 'Numerical Computation of Multivariate Normal Probabilities' (if that is what you are looking for!)?

Cuchulainn · May 18th, 2015, 1:19 pm

Another way is to write the integral as a Goursat PDE. Initial tests look very promising, e.g. 2-asset options.I do not see a reason why this FDM approach cannot also be applied to trivariate cumulative normal function. http://www.wilmott.com/messageview.cfm? ... adid=23637

AVt · May 18th, 2015, 7:01 pm

Yes, Genz. Or Vasicek. I have my own. Or http://arxiv.org/abs/1004.3616 - which I would prefer.

cm27874 · May 19th, 2015, 4:27 am

QuoteOriginally posted by: AVtYes, Genz. Or Vasicek. I have my own. Or http://arxiv.org/abs/1004.3616 - which I would prefer.Many thanks... :)An extended version has been published in JSS: http://www.jstatsoft.org/v52/i10

Cuchulainn · May 20th, 2015, 5:32 pm

Quote http://arxiv.org/abs/1004.3616Some preliminary remarks.I have studied this and related papers. In general I prefer mathematics that works for all parameter ranges without having to tweak the parameters ..Maybe I am missing something (it's not my area as such) but I have some questions:1. The article is implicitly suggesting that not all what-if cases have been taken care of; is there a nasty parameter value lurking in the wings?2. The algorithm is iterative (== recursive?) so convergence is not deterministic nor always assured? see the remarks on page 11.3. Assumptions are implicit/unflagged, e.g. in the C code if (rho > 0.99).4. How robust is the code? (IEEE 754).What I have done is to realize that the integral can be posed as the solution of a characteristic boundary value problem for a Goursate PDE and then solve by FDM. In this way I avoid asymptotics, cancellation errors etc.I downloaded the C code from the late Dr. West's site. His code is ~20% faster than mine but I create a matrix of values for later lookup. The core code is about 6 lines and it took me 30 minutes to write. I can easily change the parameters (2nd gear or fifth gear speed/accuracy).The stress tests were done and the results were in agreement (1.0e-16) with those of Graeme. Finally, I used M(a,b,rho) in the 2-aseet option formulae in the Collector's bok.

cm27874 · May 21st, 2015, 5:22 am

Quote1. The article is implicitly suggesting that not all what-if cases have been taken care of; is there a nasty parameter value lurking in the wings?No nasty parameters that I am aware of.Quote2. The algorithm is iterative (== recursive?) so convergence is not deterministic nor always assured? see the remarks on page 11.Yes, "iterative" might be the better word. I used "recursive" because Marsaglia did so in his paper on Phi. Convergence is discussed in the JSS version (the referee rightly demanded it)Quote3. Assumptions are implicit/unflagged, e.g. in the C code if (rho > 0.99).Yes, there are some cutoff/branching points, and in the paper I mentioned that they might be optimized.Quote4. How robust is the code? (IEEE 754).For the paper, I tried it in quad-double precision (using the QD library). Mike Staunton also sent me a VBA/Excel implementation. In general, if you re-compute the constants and adapt the check against the upper bound, it should work in other settings as well.Of course, your algorithm might be more flexible, and faster. At the time I wrote the paper, my main motivation was to gain deeper understanding of the bivariate normal distribution (it occurs a lot in CreditMetrics-type models and their parameterization, often in disguise). The algorithm popped up as a bonus.

Cuchulainn · May 21st, 2015, 6:00 am

Thanks for the feedback. I also stumbled on the bivariate <-> Goursat association which just seemed to work in one go. Normally numerical methods need a lot of tweaking. I can avoid many possible numerical issues because the FDM uses only +,- and a called to bivariate pdf. The code works right up to +- 0.99999.. without any cut off. I am at a lost to know what to do when abs(pdf) =1. In general, I would like to avoid trigonometric functions.I have not seen this approach before, so my interest in knowing if others have used it. Thinking out loud: solve trivariate case as u_xyz = pdf(x,y,z).Here was my V1 get it working code:

Cuchulainn · May 21st, 2015, 6:10 am

QuoteFor the paper, I tried it in quad-double precision (using the QD library). Just wondering how it works with Boost multiprecision.

logos01 · May 23rd, 2015, 1:54 pm

Seemed to me that the Genz (popularized by West) algorithm was quite good in terms of accuracy vs speed the last time I checked. There was a minor accuracy issue a while ago, I believe it's now fixed in Genz & Quantlib implementations:Trapped by the Tails of the Bivariate Normal Distribution

Cuchulainn · May 24th, 2015, 7:40 am

QuoteOriginally posted by: logos01Seemed to me that the Genz (popularized by West) algorithm was quite good in terms of accuracy vs speed the last time I checked. There was a minor accuracy issue a while ago, I believe it's now fixed in Genz & Quantlib implementations:Trapped by the Tails of the Bivariate Normal DistributionAt face value, I would expect problems with Genz when attempting to compute greeks in this way. Essentially it is numerical differentiation applied to M(a,b,r), which is inherently unstable as is well-known. And one never knows when new bugs will crop up as we tackle other pricing problems?For the PDE approach, the bespoke Aziz and Hubbard article give 2nd-order accuracy for deltas(x=S1,y=S2) and C_xy. I have not tried it get but my hunch is that it will be OK. Divided differences for PDE tend to be more stable.BTW the PDE code only needs to be written once (each time it is just a different forcing function). So, for example it is very easy to compute bivariate (NC)Chi^2, etc.

Cuchulainn · May 29th, 2015, 1:44 pm

The SSRN article by Harvey Stein sums up the idiosyncracies of numerical differentiation well. Equations (11) and (12) are asking for problems IMO.Alternatively:The original PDE is u_xy = f(x,y,rho)We can transform by t = x +y, z = x - y to get the hyperbolic PDEu_tt - u_zz = f (1)We can reduce (1) to a 1st order system to compute option price and the 2 deltas: define v = u_z, w = u_t. Thenw_t - v_z = fu_t = wv_t = w_zu = option price, v = delta_1, w = delta_2.which is a standard PDE/FDM problem.//For the Appendix B Sample Data I get 2.51517...E-13 with many mesh points,

Cuchulainn · May 31st, 2015, 6:53 pm

I did a test on the cumulative bivariate centred t distribution using the Goursat PDE approach. M(a, b; dof, rho)Results look promising:1. The integral in whole domain is 1.2. The bivariate normal is a special case.3. Checked against Dunnett/Sobel asymptotic expansions.4. A classical solution entails a horrendous series of gamma and incomplete gamma functions. In the PDE case it's just a matter of define a new pdf in the PDEu_xy = pdf(x,y,...)The PDE scheme uses no scary functions and just +, - and *.

Cuchulainn · May 31st, 2015, 6:56 pm

The trivariate case:I solve u_xyz = pdf()by PDE and I use boost multi-array to hold the data. I took the correlation matrix == I for the moment. Even now, NX = NY = NZ= 50 steps gives very accurate results. So, the data is stored and can be reused to compute stuff, like greeks and percentage points, for example.// For 4 dimensions it is just a boost multi-array again.

Cuchulainn · June 2nd, 2015, 5:28 pm

How would Euler have solved this problem?