可不可以幫我看看我的code  ，當s1=40;s2=40;k=30;sig1=0.2;sig2=0.3;p12=0.5;r=0.05;t=0.5833;q1=0.02;q2=0.03; Cmax應該是要13.4800
可是怎麼寫都不對 ，可能是二元常態裡變數問題，我順便附上題目和原本PAPER的PDF檔給大家看看
謝謝各位^ ^






new;
cls;
s1=40;
s2=40;
k=30;
sig1=0.2;
sig2=0.3;
p12=0.5;
r=0.05;
t=0.5833;
q1=0.02;
q2=0.03;

{cal1}=bs(s1,k,r,t,sig1,q1);
print "cal1" cal1;
{cal2}=bs(s2,k,r,t,sig2,q2);
print "cal2" cal2;
{Cmin}=bsm(s1,s2,k,r,q1,q2,sig1,sig2,p12,t);
print "Cmin" Cmin;
Cmax=cal1+cal2-Cmin;
print "Cmax" Cmax;
end;
proc(1)=bs(s,k,r,t,sig,q);
local d1,d2,cal;
  d1=(ln(s./k)+(r+(sig.^2)./2).*t)./(sig.*sqrt(t));
  d2=d1-sig.*sqrt(t);
  cal=s.*cdfn(d1)-k.*(exp(-(r-q).*t)).*cdfn(d2);
retp(cal);
endp;
proc(1)=bsm(s1,s2,k,r,q1,q2,sig1,sig2,p12,t);
local d11,d22,vari,A1,A2,B1,B2,C1,C2,Cmin;
  d11=(ln(s1/k)+(r-q1-(sig1^2)/2)*t)/(sig1*sqrt(t));
  d22=(ln(s2/k)+(r-q2-(sig2^2)/2)*t)/(sig2*sqrt(t));
  vari=(sig1^2)+(sig2^2)-2*p12*sig1*sig2;
  A1=d11+sig1*sqrt(t);
  A2=d22+sig2*sqrt(t);
  B1=(ln(s1/s2)+(q1-q2-(vari/2))*t)/(sqrt(vari)*sqrt(t));
  B2=(ln(s2/s1)+(q2-q1-(vari/2))*t)/(sqrt(vari)*sqrt(t));
  C1=(p12*sig2-sig1)/sqrt(vari);
  C2=(p12*sig1-sig2)/sqrt(vari);
Cmin=s1*exp(-q1*t)*cdfbvn(A1,B1,C1)+s2*exp(-q2*t)*cdfbvn(A2,B2,C2)-k*exp(-r*t)*cdfbvn(d11,d22,C2);
retp(Cmin);
endp;