Hi, I need to calculate the Heston option price with any desired accuracy as a true price. For me, it seems the famous paper by Carr and Madan (1999) is only able to calculate the Heston European option price with very high accuracy, but not within any accuracy. So, how to find a true price? Thanks. I would be grateful if someone can show me the Matlab code. (BTW, I have the code of Carr and Madan (1999) on Heston model)

No, I think you don't appreciate the role of numeric precision in computations.Just as computing sin(x) requires specifying a precision, so does the result from an integration.In Mathematica, for example, both can be specified arbitrarily high, subject only to memory constraints.A related issue is that, in Mathematica, you can set a numerical integration cutoff to [$]\infty[$], so thatdoes not become a separate source of error. Here are some examples

Thanks, Alan. I am not quite familiar with the FFT, there are several values to choose, alpha, N, eta. I choose according to Carr an Madam, with the Matlab code below. It calcuates with very high accuracy. However, I have no idea how to get a price within any desired accuracy, say 4th decimal point, by changing these parameters? I find that just change N would help little. function CallValue = HestonFft(k,theta,sigma,rho,r ,v0,s0,strike,T)%Heston European Call Price%source code from: http://www2.math.uni-wuppertal.de/~rued ... Lars.pdf%k = rate of reversion%theta = long run variance%sigma = Volatility of volatility%v0 = initial Variance%rho = correlation%T = maturity%r = interest rate%s0 = initial asset pricex0 = log(s0);alpha = 1.25;N= 4096;;eta = 0.25;b =pi/eta;u = [0:N-1]*eta;lamda = 2*b/N;position = (log(strike) + b)/lamda + 1; %position of callv = u - (alpha+1)*i;zeta = -.5*(v.^2 +i*v);gamma = k - rho*sigma*v*i;PHI = sqrt(gamma.^2 - 2*sigma^2*zeta);A = i*v*(x0 + r*T);B = v0*((2*zeta.*(1-exp(-PHI.*T)))./(2*PHI - ...(PHI-gamma).*(1-exp(-PHI*T))));C = -k*theta/sigma^2*(2*log((2*PHI - ...(PHI-gamma).*(1-exp(-PHI*T)))./ ...(2*PHI)) + (PHI-gamma)*T);charFunc = exp(A + B + C);ModifiedCharFunc = charFunc*exp(-r*T)./(alpha^2 ...+ alpha - u.^2 + i*(2*alpha +1)*u);SimpsonW = 1/3*(3 + (-i).^[1:N] - [1, zeros(1,N-1)]);FftFunc = exp(i*b*u).*ModifiedCharFunc*eta.*SimpsonW;payoff = real(fft(FftFunc));CallValueM = exp(-log(strike)*alpha)*payoff/pi;CallValue = CallValueM(round(position));end

neither use fft or matlab.You probably need to both increase the cutoff and decrease the effective spacing of the points. If you want help from me, just code a simple, non-fft, simpson's rule integration and we candiscuss convergence in that. It looks like your code can be easily modified to that as I see simpson weights in there.Once you've got that working, you can go back to your fft stuff.

The Carr Madan paper is rightly famous but it was published 14 years ago - since then, there have been improvements to the implementation such as Lord & Kahl or, as Alan says, the Heston option value is simply a numerical integral - so worth looking at his paper on contour integration

Alan, sorry for my late response. I tried the formula from Carr and Madan (1999) for ATM option, but only replace FFT with simpson' rule. The result I calculate is quite strange, with complex number. I checked my code many times, but still cannot find out what goes wrong. I would thanks a lot if you are willing to help me.function call = HestonCallQuad(kappa,theta,sigma,rho,v0,r,T,s0,K)% Exact solution of European call option pricealpha = 1.25;k = log(K);N = 100;h = 0.0001;ts = zeros(1,N/h);call = 0;for i = 0:N/h-1 ts(i+1)=i*h; if i == 0; kroneckerDelta=1; else kroneckerDelta=0; end call = call+exp(-alpha*k)/pi*HestonIntegrand(ts(i+1),kappa,theta,sigma,rho,v0,r,T,s0,alpha,k)... *(h/3*(3+(-1)^(i+1)-kroneckerDelta));endendfunction Integ = HestonIntegrand(v,kappa,theta,sigma,rho,v0,r,T,s0,alpha,k)Integ = exp(-i*v*k).*exp(-r*T).*Hestonf(v-(alpha+1)*i,kappa,theta,sigma,rho,v0,r,T,s0)./(alpha^2+alpha-v.^2+i*(2*alpha+1).*v);endfunction f = Hestonf(phi,kappa,theta,sigma,rho,v0,r,T,s0)% characteristic function of logSt in Heston model% This code from: http://math.nyu.edu/~atm262/fall06/comp ... odley.pdfu = -0.5;b = kappa;a = kappa*theta; x = log(s0);d = sqrt((rho*sigma*phi.*i-b).^2-sigma^2*(2*u*phi.*i-phi.^2));g = (b-rho*sigma*phi*i + d)./(b-rho*sigma*phi*i - d);C = r*phi.*i*T + a/sigma^2.*((b- rho*sigma*phi*i + d)*T - 2*log((1-g.*exp(d*T))./(1-g)));D = (b-rho*sigma*phi*i + d)./sigma^2.*((1-exp(d*T))./(1-g.*exp(d*T)));f = exp(C + D*v0 + i*phi*x);endQuoteOriginally posted by: Alanneither use fft or matlab.You probably need to both increase the cutoff and decrease the effective spacing of the points. If you want help from me, just code a simple, non-fft, simpson's rule integration and we candiscuss convergence in that. It looks like your code can be easily modified to that as I see simpson weights in there.Once you've got that working, you can go back to your fft stuff.

What is (i) the cutoff value (vmax), (ii) the real part of your answer and (iii) the known answer?

QuoteOriginally posted by: AlanWhat is (i) the cutoff value (vmax), (ii) the real part of your answer and (iii) the known answer?In my code, the cutoff value vmax is N=100, and the step size h is 0.001. I test an european call with true price 34.9998, the second case in Broadie & Kaya2006, and however, the result I calcuate is 30.0562 - 7.2189i. The formula in my code is from Carr and Madan (1999), but I did not use FFT. By the way, I tried to use FFT as suggest by Carr and Madan (1999), with the code is also posted here, and find that by changing N =2000000, eta =0.01, the result is exactly the same as the true price, 34.9998. I also test the first case in Broadie & Kaya, it would have error around 0.0002. I guess it would be sufficient to get the 'true price' if I calculate european call with other parameters, but just not sure. How do you choose N, eta, and alpha, all parameters with the same meaning as in Carr and Madan (1999), for the true price?

I don't read matlab, but just looking at it,1. it look like an index i and the imaginary i are being used together?? 2. It looks like your integrand is the "bad" version with branch-cut crossings As spursfan suggested, your source material is obsolete. Try Rouah's book

Never mind, Alan. Thanks for the book you recommend.

By the way, you may want to try the Fractional FFT, which is a little more flexible and perhaps better suited to your needs. The FRFT is also coded in my book that Alan is referring to.

Thank you, Fabrice, I like your very useful book.