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 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!

goldorak


Total Posts: 1000
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: 1000
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: 1000
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: 375
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
}
Previous Thread :: Next Thread 
Page 1 of 1