• C++ Programming for Financial Engineering
    Highly recommended by thousands of MFE students. Covers essential C++ topics with applications to financial engineering. Learn more Join!
    Python for Finance with Intro to Data Science
    Gain practical understanding of Python to read, understand, and write professional Python code for your first day on the job. Learn more Join!
    An Intuition-Based Options Primer for FE
    Ideal for entry level positions interviews and graduate studies, specializing in options trading arbitrage and options valuation models. Learn more Join!

Problem with R code, with option pricing

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 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+alpha-v**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,((N-1)*h),by=h)
eta=h
b=pi/eta
lambda=(2*pi)/(N*eta)
u=seq(1,N,by=1)
k_u=-b+lambda*(u-1)
f=function(i){
sum=0
for(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^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!
 
Top