Hi, I am calculating exact prices for volswaps under the Heston model in Matlab using the Laplace transform of square root function (see formula 11.6 in Gatheral's "Volatility Surface"). In particular I am trying to reproduce Carr & Lee's value of 19.02 (in their paper "Robust replication of volatility derivatives"). Although I manage to get quite close (18.98) I am not sure I should be satisfied with the slight difference, especially since sometimes I get NaN, when for instance I take b too large in the integration domain [0,b]. Does anyone know why the NaN?The codeI am using is below. Thanks.-------------------------------------------------function ret = vs()ret = quadgk(@vsint,0,6500000)-------------------------------------------------function result = vsint(x)%Heston parameterslmd = 1.15;tht = 0.04;eta = 0.39;T = .5;v0 = 0.04;%%%%%%%%%%%%%%%%%%phi = sqrt(lmd^2 + 2*eta^2*x);A1 = 2*phi.*exp(0.5*(phi + lmd)*T);A2 = (phi + lmd).*(exp(phi*T)-1) + 2*phi;B1 = 2*(exp(phi*T)-1);B2 = (phi + lmd).*(exp(phi*T)-1) + 2*phi;A = (A1./A2).^(2*lmd*tht/eta^2);B = B1./B2;E = A.*exp(-x.*v0.*B);result = (1 - E)./(x.^(3/2).*2.*sqrt(pi.*T));

Don't know Matlab, but in Mathematica, it's poor practice to write (0,b) as an integration range for large b. Instead, oneshould always write (0,1,b), forcing the integrator to sample a point where the integrand is O(1). This will cure some errors in that system anyway. In Mathematica, with xmax = 10^nTable[{n,NIntegrate[result[x],{x,0,1,10^n}]},{n,1,14}]1 0.0482362 0.1151233 0.1649364 0.1821845 0.1876396 0.1893657 0.189918 0.1900839 0.19013710 0.19015411 0.1901612 0.19016213 0.19016214 0.190162so apparently you need [$]x_{max} > 10^{12}[$] for 6 good digits.This is not surprising, as the same would be true of [$]\int_1^{x_{max}} x^{-3/2} \, dx[$]. A simpler approach in Mathematica is to use the keyword Infinity as the integration limit, which prompts the use of a suitable transformation method:NIntegrate[result[x],{x,0,1,Infinity}]0.190162 I will guess Matlab offers something similar.

You need to debug your code.NaN can occur whenQuoteOperations generating NaN[edit]There are three kinds of operations that can return NaN:[4]Operations with a NaN as at least one operand.Indeterminate forms The divisions 0/0 and ±∞/±∞The multiplications 0×±∞ and ±∞×0The additions ∞ + (−∞), (−∞) + ∞ and equivalent subtractionsThe standard has alternative functions for powers: The standard pow function and the integer exponent pown function define 00, 1∞, and ∞0 as 1.The powr function defines all three indeterminate forms as invalid operations and so returns NaN.Real operations with complex results, for example: The square root of a negative number.The logarithm of a negative numberThe inverse sine or cosine of a number that is less than −1 or greater than +1.I don't know if Matlab has function to test NaN and to catch NaNs?Is the NaN quiet or signallig NaN? QNaN don't crash but can give the wrong answer; sNaN crashes.

Thanks Alan & Cuch! I still need / would like to figure out what caused the NaN, but indeed if I let Matlab integrate over (0,inf) instead of (0,x_max) I get the exact value of 0.1902. Strange, but interesting to find out. If / when I do will post it here.

Frolloos,I see Matlab does support isnanMaybe try to put in some statements to pinpoint. In any case sqrt(x) ==> NaN for x < 0.

QuoteOriginally posted by: CuchulainnFrolloos,I see Matlab does support isnanMaybe try to put in some statements to pinpoint. In any case sqrt(x) ==> NaN for x < 0.Yes I'll have to debug. Thought I was out of the woods with (0,inf), but when I ran the code with T > 0.5 I get more NaNs. Matlab says infinity encountered.

QuoteOriginally posted by: frolloosQuoteOriginally posted by: CuchulainnFrolloos,I see Matlab does support isnanMaybe try to put in some statements to pinpoint. In any case sqrt(x) ==> NaN for x < 0.Yes I'll have to debug. Thought I was out of the woods with (0,inf), but when I ran the code with T > 0.5 I get more NaNs. Matlab says infinity encountered.It's a start the first NaN is a symptom of something more serious? I once solved a problem that only crashed for vol = 62.32... But it was tip of ye olde ijsberg. Quoteresult = (1 - E)./(x.^(3/2).*2.*sqrt(pi.*T)); 11.6

QuoteOriginally posted by: CuchulainnQuoteOriginally posted by: frolloosQuoteOriginally posted by: CuchulainnFrolloos,I see Matlab does support isnanMaybe try to put in some statements to pinpoint. In any case sqrt(x) ==> NaN for x < 0.Yes I'll have to debug. Thought I was out of the woods with (0,inf), but when I ran the code with T > 0.5 I get more NaNs. Matlab says infinity encountered.It's a start the first NaN is a symptom of something more serious? I once solved a problem that only crashed for vol = 62.32... But it was tip of ye olde ijsberg.Looked at the Matlab documentation but didn't find anything particularly illuminating. Not sure what is meant with 'decays sufficiently rapidly': (will contact them tomorrow and ask them 'watskebeurd' )"If the interval is infinite, [a,Inf), then for the integral of fun(x) to exist, fun(x) must decay as x approaches infinity, and quadgk requires it to decay rapidly. Special methods should be used for oscillatory functions on infinite intervals, but quadgk can be used if fun(x) decays fast enough.The quadgk function will integrate functions that are singular at finite endpoints if the singularities are not too strong. For example, it will integrate functions that behave at an endpoint c like log|x-c| or |x-c|p for p >= -1/2. If the function is singular at points inside (a,b), write the integral as a sum of integrals over subintervals with the singular points as endpoints, compute them with quadgk, and add the results"

Do you get NaN's if the lower integration limit is 1?

QuoteOriginally posted by: AlanDo you get NaN's if the lower integration limit is 1?a = 0b = 6 500 000ret = quadgk(@vsint,0,ret = quadgk(@vsint,0,6500000)b looks kind of big. If the integrand decays fast then indeed a smaller b is OK? And move a away from 0? Have you tried other methods, e.g. Clenshaw-Curtis, even midpoint?

QuoteOriginally posted by: AlanDo you get NaN's if the lower integration limit is 1?Yes, still getting "Not a Number" if 1 is lower limit. If integration range is (0,1e4) I don't have any problems at all with various maturities and parameter choices. I think it's because the integrand does not decay fast enough that gives rise to the issues for (0,inf) or (0,b) with b very large. I'll try other methods than gauss quadrature, maybe that works

In some cases if the integrand falls off exponentially, one could try a change of variables t = exp(-x). We then have a finite range [0,1] of integrand - f(-log t) / t Hopefully a removable singularity at t = 0.

A colleague of mine came up with a nice 'renormalization' of the integral - I will post it later tonight - which so far makes the thing work. The integral blows up for large z because z occurs in phi, and phi in exp(phi). Not sure if only Matlab has problems with this integrand / integral or other languages/packages as well. So the workaround I'll post later may be of interest to others as well. EDIT: Although as mentioned above re-scaling of the integral is possible, there is actually a more straighforward solution I found out. Matlab calculates the denominator and numerator separately in for instance A1. The exp(phi*T) can become so large that Matlab treats it as infinite and returns NaN. The solution is to divide everything by this exponential term. Posted the modified code below for those interested. No issues anymore. Some values for the volswap:T = 0.5, 0.1902T = 1, 0.1874T = 3, 0.1888T = 5, 0.1912T = 10, 0.1945T = 50 (out of curiosity), 0.1986-------------------------------------------------------------------function result = vsint(x)%Heston parameterslmd = 1.15;tht = 0.04;eta = 0.39;T = 5;v0 = 0.04;%Integrandphi = sqrt(lmd^2 + 2*eta^2*x);%A1 = 2*phi.*exp(0.5*(phi + lmd)*T);%A2 = (phi + lmd).*(exp(phi*T)-1) + 2*phi;%B1 = 2*(exp(phi*T)-1);%B2 = (phi + lmd).*(exp(phi*T)-1) + 2*phi;A1 = 2*phi.*exp(0.5*lmd*T);A2 = (phi + lmd).*(exp(0.5*phi*T)-exp(-0.5*phi*T)) + 2*phi.*exp(-0.5*phi*T);B1 = 2;B2 = (phi + lmd) + 2*phi.*(exp(phi*T)-1).^-1;A = (A1./A2).^(2*lmd*tht/eta^2);B = B1./B2;E = A.*exp(-x.*v0.*B);result = (1 - E)./(x.^(3/2).*2.*sqrt(pi.*T));

Nice.I have tested exp(Big) in C++ and it also gives NaN.

Hi,using log() might do the job... You can do all calculations using logs and finally taking the exp, then, at least in intermediate calculation you reduce the appearance of very large numbers that you have to pass to exp()...Best, Lapsi