Forums > Pricing & Modelling > Problem with R code, with option pricing

 Page 1 of 1
Display using:
 quithating Total Posts: 6 Joined: Feb 2017
 Posted: 2017-02-18 21:07 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 Carr-Madan approach to option pricing, using the Black-Scholes 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 Black-Scholes modelcf=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)*imx=(exp(-r*T)*cf(S_0,mu,sigma,T,u))/(alpha**2+alpha-v**2+im*(2*alpha+1)*v)x}#function used in the simpson weightskdf=function(x){if(x==0){result=1}else{result=0}result}#function used in simpson weightskdfn=function(x,N){if(x==N){result=3}else{result=0}result}#FFT of the eq(6) in the linkcallcm=function(r,S_0,mu,sigma,T,alpha,N,h){im=complex(real=0,imaginary=1)v=seq(0,((N-1)*h),by=h)eta=hb=pi/etalambda=(2*pi)/(N*eta)u=seq(1,N,by=1)k_u=-b+lambda*(u-1)f=function(i){sum=0for(j in 1:N){sum=sum+((exp((-im*lambda*eta*(j-1)*(i-1))+im*v[j]*b)*psim(r,S_0,mu,sigma,T,alpha,v[j]))*(eta/3)*(3+(-1)**(j)-kdf(j-1)-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^11h=0.01S_0=210.59T=4/365r=0.002175alpha=1mu=rsigma=0.1404But 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: 1041 Joined: Nov 2004
 Posted: 2017-02-19 16:44 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.
 quithating Total Posts: 6 Joined: Feb 2017
 Posted: 2017-02-19 17:15 Thanks for the top tip!
 goldorak Total Posts: 1041 Joined: Nov 2004
 Posted: 2017-02-19 17:22 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.
 quithating Total Posts: 6 Joined: Feb 2017
 Posted: 2017-02-19 17:25 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: 1041 Joined: Nov 2004
 Posted: 2017-02-19 17:26 21st century... these guys will never find jobs If you are not living on the edge you are taking up too much space.
 quithating Total Posts: 6 Joined: Feb 2017
 Posted: 2017-02-19 17:28 Nice speaking to you, try to keep that blood pressure down though. :)
 rod Total Posts: 376 Joined: Nov 2006
 Posted: 2017-02-19 20:58 http://quant.stackexchange.com
 quithating Total Posts: 6 Joined: Feb 2017
 Posted: 2017-02-19 21:44 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,((N-1)*h),by=h) eta=h lambda=(2*pi)/(N*eta) b=N*lambda/2 u=seq(1,N,by=1) k_u=-b+lambda*(u-1) w=vector() for(i in 1:N){ w[i]=(3+(-1)^(i)) } w=w-c(1,rep(0,N-1)) 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.
 quithating Total Posts: 6 Joined: Feb 2017
 Posted: 2017-02-20 18:54 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}
 Page 1 of 1