QuoteOriginally posted by: AlanThanks, Daniel. Posting the code would be very helpful, too.If I have questions about it, I'll bug you off-line.Sure, any time

hi everybody, I am currently implementing on matlab 3 models of stochastic volatility, and I have trouble with CEV implementation. I don't know why but the summation of the product of the density gamma distribution with the complementary cumulative does not work. I tried to fix it but I don't know what is wrong. Here is the code of the matlab function I created, if anyone could help me out on this, it will be very appreciatedfunction callCEV= CEV(N,alpha,sigma,r,q,T,s,E)% E=strike/N=upperbound of the summation/q=dividend% proportional to S%%CEV model: ds = s*r*dt+ s^(alpha/2)*sigma * dwt%%Parameter used in the gamma functionsk = (2*(r-q))/(sigma^2*(2-alpha)*(exp((r-q)*(2-alpha)*T)-1));x = k * s^(2-alpha)*exp((r-q)*(2-alpha)*T);v = k * E^(2-alpha);Som1=0;Som2=0;m1=0;m2=0;for n=0:1:N m1=n+1+1/(2-alpha); m2=n+1; Som1 = Som1 + gampdf(x,m2,1) * (1-gamcdf(v,m1,1)); Som2 = Som2 + gampdf(x,m1,1) * (1-gamcdf(v,m2,1));endcallCEV= s*exp(-q*T)*Som1 - E*exp(-r*T)*Som2;

I am not a matlab user and don't know the algorithm.But, since no one has replied, I will make one suggestion basedon my Mathematica experience with this problem. In Mathematica,I first relied upon a built-in cdf (non-central chisquared), but thisled to problems, I corrected the problems by coding the cdfdirectly from the density (involving a Bessel function). Similarly, youmight investigate if the source of your problems is the gamcdf(a built-in?), in which case you can fix it the same way I did.

thank you Alan. I will code my own cumulative and compare it with the matlab built-in function. However, I found another formula for the square root CEV which involves polar coordinates and does not require to code gamma function. Even with this, my code does not work (returning O). I was wondering if the term defined by Cox k = 2*r / (sigma^2 * ( exp(r*T) - 1 )) takes very high values and therefore impact my results. I will try to code my own gamma function and I will see. Thank you anyway

Quote1) 0.3412572) 0.357819I use a direct computation of the noncentral chi function as a series of incomplete gamma function.How do you determine the effects of round-off errors in all of these computations? Are the results really correct or do solutions give the same round off errors? Strictly speaking we should used interval arithmetic. Is the extra computational effort in getting accuracy fro 10^-3 to 10^-5 worth it? This is meant as a serious question.

QuoteOriginally posted by: CuchulainnQuote1) 0.3412572) 0.357819I use a direct computation of the noncentral chi function as a series of incomplete gamma function.How do you determine the effects of round-off errors in all of these computations? Are the results really correct or do solutions give the same round off errors? Strictly speaking we should used interval arithmetic. Is the extra computational effort in getting accuracy fro 10^-3 to 10^-5 worth it? This is meant as a serious question.I use a floating bound to truncate the series plus some reshuffling of the terms. As far as I can see, my results agree with those posted by Alan, who most likely used a Mathematica software to get them, so no bad for my own implementation of the noncentral chi function. Since it is an analytical formula assuming no discretization error, so I would say that you must get accuracy of 5-6 digits.

QuoteOriginally posted by: CuchulainnYou're welcome, Mike. I use a slightly different notation from Schoder, but it worked out fine. BTW, I also tried to get an answer on MC and Euler, but the answer is not accurate (S can become negative). There was a thread on this issue when r = 0.//It might be interesting to test the boost's C++ NCHSquared distribution with this algorithm. Is it documented somewhere?I used a bit the boost non_central_chi_squared cdf.it is fairly easy to get negatives call prices for deep OTM calls..

QuoteOriginally posted by: ApprenticeQuoteOriginally posted by: CuchulainnYou're welcome, Mike. I use a slightly different notation from Schoder, but it worked out fine. BTW, I also tried to get an answer on MC and Euler, but the answer is not accurate (S can become negative). There was a thread on this issue when r = 0.//It might be interesting to test the boost's C++ NCHSquared distribution with this algorithm. Is it documented somewhere?I used a bit the boost non_central_chi_squared cdf.it is fairly easy to get negatives call prices for deep OTM calls.. Ouch! Is the problem Schroder or the boost algo?Which values did you use? I would like to try on ADE. Thanks.

QuoteOriginally posted by: CuchulainnQuoteOriginally posted by: ApprenticeQuoteOriginally posted by: CuchulainnYou're welcome, Mike. I use a slightly different notation from Schoder, but it worked out fine. BTW, I also tried to get an answer on MC and Euler, but the answer is not accurate (S can become negative). There was a thread on this issue when r = 0.//It might be interesting to test the boost's C++ NCHSquared distribution with this algorithm. Is it documented somewhere?I used a bit the boost non_central_chi_squared cdf.it is fairly easy to get negatives call prices for deep OTM calls.. Ouch! Is the problem Schroder or the boost algo?Which values did you use? I would like to try on ADE. Thanks.I think the problem comes from the boost implementation (to confirm yet)I used the formula from Schroeder paper for the call price , but with no drift (so just doing CEV on the forward directly, I am more looking at rates actually)for the following set of params, I get some issues :type = callforward = 2%maturity = 1delta/sigma = 2%beta = 1.6strike = 1.4% => price is lower than the intrinsic valuestrikes = from 2.8% to 4% => negative pricethe pb I think comes from the fact that boost actually uses 3 methods to compute the nc chisquare cdf, depending on the respective values of x , k (d.f.) and l (n.c param)so for a given application of the call formula, the two uses of the nc chi square cdf end up actually using two different approximations, hence breaking basic cdf properties it seems.

QuoteOriginally posted by: ApprenticeQuoteOriginally posted by: CuchulainnQuoteOriginally posted by: ApprenticeQuoteOriginally posted by: CuchulainnYou're welcome, Mike. I use a slightly different notation from Schoder, but it worked out fine. BTW, I also tried to get an answer on MC and Euler, but the answer is not accurate (S can become negative). There was a thread on this issue when r = 0.//It might be interesting to test the boost's C++ NCHSquared distribution with this algorithm. Is it documented somewhere?I used a bit the boost non_central_chi_squared cdf.it is fairly easy to get negatives call prices for deep OTM calls.. Ouch! Is the problem Schroder or the boost algo?Which values did you use? I would like to try on ADE. Thanks.I think the problem comes from the boost implementation (to confirm yet)I used the formula from Schroeder paper for the call price , but with no drift (so just doing CEV on the forward directly, I am more looking at rates actually)for the following set of params, I get some issues :type = callforward = 2%maturity = 1delta/sigma = 2%beta = 1.6strike = 1.4% => price is lower than the intrinsic valuestrikes = from 2.8% to 4% => negative pricethe pb I think comes from the fact that boost actually uses 3 methods to compute the nc chisquare cdf, depending on the respective values of x , k (d.f.) and l (n.c param)so for a given application of the call formula, the two uses of the nc chi square cdf end up actually using two different approximations, hence breaking basic cdf properties it seems.Just want to get the notation right:dF = delta*F^(beta/2)dW or is delta detemined bydelta * F^(beta/2 -1) = 0.02?And forward = 2% is the strike?

dF = delta*F^(beta/2)dW is spot on for the notation.forward F_0 = 2%then I price options with different absolute strikes { 1.4% , 2.8% , etc.. } where I observe negative time value for calls.

See my comment in this thread, Sat Nov 22, 08 08:20 PM -- you may need to do the same.

QuoteOriginally posted by: ApprenticedF = delta*F^(beta/2)dW is spot on for the notation.forward F_0 = 2%then I price options with different absolute strikes { 1.4% , 2.8% , etc.. } where I observe negative time value for calls.I tested the model using both FDM (ADE) and Monte Carlo (good old Euler)F_0 = .02T = 1delta = .02beta = 1.6 / 2 = 0.8K = 0.014C (FDM) = 0.006, C (MC) = 0.00599807K = 0.028For both cases I get C = 0Is this the expected result?

Here's what I get for these two cases:K=0.014 Call = 0.006104906313K=0.028 Call = 0.0000665295511 (using a new experimental representation for the formula I have been playing around with)Likely all the digits are good.

QuoteOriginally posted by: AlanHere's what I get for these two cases:K=0.014 Call = 0.006104906313(using a new experimental representation for the formula I have been playing around with)Likely all the digits are good.Now I tried it with NY = NT = 4000 with ADE and implicit Euler; the value remains 0.006. The values on either side are 0.00598004, 0.00602004. Same with NY = NT = 6000 and the interval on each side is [0.00598667, 0.00601334].