oops, inverse of a cumulative bivariate doesn't make sense.

The Goursat FDM generates the matrix of values in one sweep. Then use interpolation for in-between values.
On a 50X50 grid we can easily get 7 digits accuracy.

I admit I don't know what this is about, I briefly looked at this thread just now and I saw 50x50 grid and 7 digits accuracy in one sentence and I don't think I've ever seen that. I'm assuming this is a discretized 2D PDE solved with FDM you're talking about. Do you use a sixth or eighth order discretization?
EDIT: OK, I just saw extrapolation, if the base scheme is 4th order then OK. If it's 2nd then again this sounds strange, one extrapolation would get you to 4th (in theory) and a double one to 6th, but that would need a dense grid to actually deliver that theoretical order even if we're talking about a smooth function.

Looking at the trials with 70x80-ish grid the accuracy for 50x50 would be approx 2 digits? 

However I think these reported values are flawed. If you look at the trial 1,2 and 3 results you see they have different parameters and grid sizes, but the maximum error is reported to be *exactly* the same for all three cases:

I think Cuch is just keeping track of the max error up to that trial. So until a larger one appears in a subsequent trial the max error isn't updated.

Is the error estimated based on comparing against Genz West? They claim near double precision accuracy.

What's interesting to note is that the central motive for fast and high precision algorithms in Genz West paper is that conventional lower precision methods are to be avoided in option pricing (see section 4 "Option pricing disasters" in https://lyle.smu.edu/~aleskovs/emis/sqc ... umnorm.pdf )

I think Cuch is just keeping track of the max error up to that trial. So until a larger one appears in a subsequent trial the max error isn't updated.

Indeed. The list of all errors is printed at the end.
I had to get used to it as well. max stays max until a bigger one comes along.

The Goursat FDM generates the matrix of values in one sweep. Then use interpolation for in-between values.
On a 50X50 grid we can easily get 7 digits accuracy.

I admit I don't know what this is about, I briefly looked at this thread just now and I saw 50x50 grid and 7 digits accuracy in one sentence and I don't think I've ever seen that. I'm assuming this is a discretized 2D PDE solved with FDM you're talking about. Do you use a sixth or eighth order discretization?
EDIT: OK, I just saw extrapolation, if the base scheme is 4th order then OK. If it's 2nd then again this sounds strange, one extrapolation would get you to 4th (in theory) and a double one to 6th, but that would need a dense grid to actually deliver that theoretical order even if we're talking about a smooth function.

I am using the 2nd order one-step box scheme for Goursat (it's old as Metusalah) on grids N, 2N, 4N until catastrophic cancellation creeps in.
Some output is (I compare against the posse)
Haug 2007 value =  0.273 122 222 6 DREWES2
Steps NX, NY: 200, 200
Goursat FDM Classico: 0.2731193155592519
Goursat FDM extrapolate h/2: 0.273 122 2267525446

Goursat FDM extrapolate h/4: 0.273 122 0325492174 (cancellation issues??)

What does extrapolate h/2 mean? I'm assuming h is the uniform grid spacing, so you use a 400x400 grid (since your base is 200x200) and then get Extrapolated solution= (4/3)*V(h/2) - (1/3)*V(h)?
And h/4 ?

What does extrapolate h/2 mean? I'm assuming h is the uniform grid spacing, so you use a 400x400 grid (since your base is 200x200) and then get Extrapolated solution= (4/3)*V(h/2) - (1/3)*V(h)?
And h/4 ?

Yes to 1st
h/4 under assumption of even powers
(16 v2 - v1)/15
(I have not proven this mathematically)

1. Is the error estimated based on comparing against Genz West? They claim near double precision accuracy.

2. What's interesting to note is that the central motive for fast and high precision algorithms in Genz West paper is that conventional lower precision methods are to be avoided in option pricing (see section 4 "Option pricing disasters" in https://lyle.smu.edu/~aleskovs/emis/sqc ... umnorm.pdf )

1. Yes; the error is between my FDM and Genz/West as benchmark.
2. I have tested Goursat against the two option prices in Haug 2007 and I get the same answer.
A penetrating question is how accurate must M(a,b,rho) be to give a given accuracy in a two-asset option pricing formula?

oops, inverse of a cumulative bivariate doesn't make sense.

No unique solution; it is a level set.

OK, so that's using 400x400(plus 200x200), but you cannot get such accuracy from a 50x50 plus 25x25 grid. Extrapolation can deliver very high accuracy but only when the grid is rather fine and you're well into the asymptotic range. 50x50 is not there I'm pretty sure. (plus of course a well behaved scheme with constant convergence order and a smooth function, which you have)

oops, inverse of a cumulative bivariate doesn't make sense.

No unique solution; it is a level set.

Indeed.

So the sweetspot of your method would (IMO) be low precision grids?

If you do h/4 extrapolation the speed will drop with a factor 16 and at that point the near double precision Genz West method would be both faster and more accurate -even on a grid-.

The Genz West paper mentions option pricing issues, perhaps (just guessing!) negative probabilities that might not arise in your method?

I have tested some option, Haug page 224 (Haug gives 2.7939).. Some nice accuracy and efficiciency results on an awfully rough grid.
Steps NX, NY: 50, 50

Call best cash-or-nothing Goursat: 2.79451280353896
Elapsed time: 0.007

Call best cash-or-nothing Goursat extrapolation: 2.793935288536937
Elapsed time: 0.008

Call best cash-or-nothing Genz/West: 2.793935852582142
Elapsed time: 0.006

Call best cash-or-nothing Tanh second opinion: 2.793935850289368
Elapsed time: 0.007

// And on a even rougher grid
Steps NX, NY: 15, 15
 
Call best cash-or-nothing Goursat: 2.800610790093465
Call best cash-or-nothing Goursat extrapolation: 2.793862604891467
Call best cash-or-nothing Genz/West: 2.793935852582142
Call best cash-or-nothing Tanh second opinion: 2.793935850289368

The Genz West paper mentions option pricing issues, perhaps (just guessing!) negative probabilities that might not arise in your method?

Good spot. 
The point about Goursat FDM is that it satisfies the maximum principle (Holy Grail). See theorem 2.1
http://www.ams.org/journals/mcom/1964-1 ... 0337-7.pdf
So, IMO there is no way negative values can occur.

//
For that case (section 4.3) I get 0.017979, GW get 0.0180005.
For rho = -.95 I get 5.36249e-8 , GW gets 5.45857e-8.

// Problems in section 4.4. can be solved with multiprecision. Just 'mechanical' problems.