

I don't know if this is the right place to post this, or if anyone can help with this, but I am having problems with my R code not producing accurate results. I am trying to implement the CarrMadan approach to option pricing, using the BlackScholes model. The formula can be found in equation (6)on page 64 here:http://engineering.nyu.edu/files/jcfpub.pdf. I am using the FFT with simpson weights (eq(22)) to try and approximate the integral, however the results aren't as near as they should be. My question is can you tell me where I have gone wrong with my code, or if I have made an error?
R code:
#The characteristic function of BlackScholes model cf=function(S_0,mu,sigma,T,u){ im=complex(real=0,imaginary=1) x=exp(im*u*(log(S_0)+(mu((sigma**2)/2))*T)(((sigma**2)*T*(u**2))/2)) x } #The psi function (eq(4) in link) psim=function(r,S_0,mu,sigma,T,alpha,v){ im=complex(real=0,imaginary=1) u=v(alpha+1)*im x=(exp(r*T)*cf(S_0,mu,sigma,T,u))/(alpha**2+alphav**2+im*(2*alpha+1)*v) x } #function used in the simpson weights kdf=function(x){ if(x==0){ result=1 } else{ result=0 } result } #function used in simpson weights kdfn=function(x,N){ if(x==N){ result=3 } else{ result=0 } result } #FFT of the eq(6) in the link callcm=function(r,S_0,mu,sigma,T,alpha,N,h){ im=complex(real=0,imaginary=1) v=seq(0,((N1)*h),by=h) eta=h b=pi/eta lambda=(2*pi)/(N*eta) u=seq(1,N,by=1) k_u=b+lambda*(u1) f=function(i){ sum=0 for(j in 1:N){ sum=sum+((exp((im*lambda*eta*(j1)*(i1))+im*v[j]*b)* psim(r,S_0,mu,sigma,T,alpha,v[j]))*(eta/3)*(3+(1)**(j)kdf(j1)kdfn(j,N))) } sum } result=sapply(u,f) result=((exp(alpha*k_u))/(pi))*Re(result) result }
I've been using parameter values: N=2^11 h=0.01 S_0=210.59 T=4/365 r=0.002175 alpha=1 mu=r sigma=0.1404
But have at best been getting results that are correct to 2 dp, whereas I am led to believe I should be getting results as accurate as say 10^14. So any help on where you can see errors in my formulation or the code would be much appreciated thank you! 




goldorak


Total Posts: 986 
Joined: Nov 2004 


May be the first thing to do is to quit "programming" in R.... or at least quit the way you "program" in R...

If you are not living on the edge you are taking up too much space. 


 

goldorak


Total Posts: 986 
Joined: Nov 2004 


2nd post already. Now you can at least pretend that you are contributing.
Seriously young man... did you think you were connecting to backoffice_helpdesk_111_please_help_me_i_am_a_lazy_student_and_cannot_get_answer_on_snapshat.com?

If you are not living on the edge you are taking up too much space. 



haHA, good one. You got me, I've done nothing and have asked people to do it for me, teehee. Thanks for all the help, and thanks for reading the question. cya soon ;) 




goldorak


Total Posts: 986 
Joined: Nov 2004 

 


Nice speaking to you, try to keep that blood pressure down though. :) 




rod


Total Posts: 370 
Joined: Nov 2006 


http://quant.stackexchange.com 




Thank you, I've got a post on there as well, I've improved on the function as well: callcm=function(r,S_0,mu,sigma,T,alpha,N,h){ im=complex(real=0,imaginary=1) v=seq(0,((N1)*h),by=h) eta=h lambda=(2*pi)/(N*eta) b=N*lambda/2 u=seq(1,N,by=1) k_u=b+lambda*(u1) w=vector() for(i in 1:N){ w[i]=(3+(1)^(i)) } w=wc(1,rep(0,N1)) w=1/3*w f=function(v){ psim(r,S_0,mu,sigma,T,alpha,v) } result=exp(im*b*v)*(sapply(v,f)*eta)*w result=Re(fft(result)) result=((exp(alpha*k_u))/(pi))*result result }
Using R's fft function. 





I have managed to obtain accurate results(error of order 10^9) but with single strike prices, using the pracma package in R, specifically using the adapted Simpsons numerical approximation method.
callp=function(r,S_0,mu,sigma,T,alpha,k){ im=complex(real=0,imaginary=1) fr=function(v){ Re(exp(im*v*k)*psim(r,S_0,mu,sigma,T,alpha,v)) } resultr=simpadpt(fr,0,2^11,tol=0.00001) result=resultr*(exp(alpha*k))/(pi) result } 







