[有偿编程] 老师们，为什么我尝试复现的这个Asian Option公式程序算出来是NaN [推广有奖]

0关注
1粉丝

高中生

25%

还不是VIP/贵宾

-

0%

威望: 0 级
论坛币: 10 个
通用积分: 0.6000
学术水平: 0 点
热心指数: 0 点
信用等级: 0 点
经验: 69 点
帖子: 5
精华: 0
在线时间: 44 小时
注册时间: 2018-7-22
最后登录: 2023-1-25

楼主

kuluomu 发表于 2022-3-3 00:14:14 |AI写论文

300论坛币

各位老师们，为什么我这个代码算出来的结果和中途积分的结果都是NaN，可以正常运行函数代码复现的是这篇文章的Asian option定价公式（定理4.1）

Kim, B. & Wee, I. 2014, "Pricing of geometric Asian options under Heston's stochastic volatility model", Quantitative finance, vol. 14, no. 10, pp. 1795-1809.

函数代码如下

function call=KimCallQuad(nu,theta,s0,r,T,t,sigma,rho,kappa,w,n,K)
call=exp(-r*(T-t))*KimB(nu,theta,s0,r,T,t,sigma,rho,kappa,w,n,K);
function ret=KimB(nu,theta,s0,r,T,t,sigma,rho,kappa,w,n,K)
global psi10
ret=((psi10-K)/2)+(1/pi)*quadl(@KimIntegral,0,10^5,[],[],...
nu,theta,s0,r,T,t,sigma,rho,kappa,w,n,K);
function ret=KimIntegral(nu,theta,s0,r,T,t,sigma,rho,kappa,xi,w,n,K)
ret=real((Kimpsi(nu,theta,s0,r,T,t,sigma,rho,kappa,xi,w,n,3)-...
Kimpsi(nu,theta,s0,r,T,t,sigma,rho,kappa,xi,w,n,2)*K)*...
(exp(-i*xi*log(K)))/(i*xi));

复制代码

function psi=Kimpsi(nu,theta,s0,r,T,t,sigma,rho,kappa,xi,w,n,type)
if type == 1
xi=1;
else
if type == 2
xi=1i*xi;
else
type == 3
xi=1+1i*xi;
end
end
a1=2*nu/(sigma^2);
a2=2*kappa*theta/(sigma^2);
a3=((T-t)/T)*log(s0)+((r*sigma-kappa*theta*rho)*(T-t)^2)/(2*sigma*T)-(rho*(T-t)/(sigma*T))*nu;
a4=log(s0)-(rho/sigma)*nu+(r-(rho*kappa*theta/sigma))*(T-t);
a5=(kappa*nu+(kappa^2)*theta*(T-t))/(sigma^2);
psi=exp(-a1*(KimHtilde(T,t,sigma,rho,kappa,xi,w,n,type)/KimH(T,t,sigma,rho,kappa,xi,w,n,type))-...
a2*log(KimH(T,t,sigma,rho,kappa,xi,w,n,type))+...
a3*xi+a4*w+a5);
function H=KimH(T,t,sigma,rho,kappa,xi,w,n,type)
if type == 1
xi=1;
else
if type == 2
xi=1i*xi;
else
type == 3
xi=1+1i*xi;
end
end
h=[];
h(1)=0; % h-2
h(2)=0;% h-1
h(3)=1;% h0
h(4)=((T-t)*(kappa-w*rho*sigma))/2; % h1
for j=5:n+4 % h for n bigger than or equal 2
h(j)=(((T-t)^2)/ (4*n*(n-1)*T^2))*(-1*(xi^2)*(sigma^2)*(1-(rho^2))*((T-t)^2)*h(j-4)+...
(xi*sigma*T*(sigma-2*rho*kappa)-2*xi*w*sigma^2*T*(1-(rho^2)))*(T-t)*h(j-3)+...
T*((kappa^2)*T-2*xi*rho*sigma-w*(2*rho*kappa-sigma)*sigma*T-(w^2)*(1-(rho^2)*(sigma^2)*T)*h(j-2)));
end
H=sum(h(3:n+4));
function Htilde=KimHtilde(T,t,sigma,rho,kappa,xi,w,n,type)
if type == 1
xi=1;
else
if type == 2
xi=1i*xi;
else
type == 3
xi=1+1i*xi;
end
end
h=[];
h(1)=0; % h-2
h(2)=0; % h-1
h(3)=1; % h0
h(4)=((T-t)*(kappa-w*rho*sigma))/2; % h1
for j=5:n+4 % h for n bigger than or equal 2
h(j)=(((T-t)^2)/ (4*n*(n-1)*T^2))*(-1*(xi^2)*(sigma^2)*(1-(rho^2))*((T-t)^2)*h(j-4)+...
(xi*sigma*T*(sigma-2*rho*kappa)-2*xi*w*sigma^2*T*(1-(rho^2)))*(T-t)*h(j-3)+...
T*((kappa^2)*T-2*xi*rho*sigma-w*(2*rho*kappa-sigma)*sigma*T-(w^2)*(1-(rho^2)*(sigma^2)*T)*h(j-2)));
hh(j)=(j-4)/(T-t)*h(j);
end
Htilde=sum(hh(5:n+4))+1/(T-t)*h(4);