cal 1971 2 12
all 2008:12
OPEN DATA "C: \ Desktop\ heterosc_option_Last_2.xls"
DATA(FORMAT=xls,ORG=COL) / f1 f2 x r
OPEN DATA "C: \ Desktop\ Data_T.xls"
DATA(FORMAT=xls,ORG=COL) / S T1 r2 K T2
com N = 3 , M = 61 , J = 600
dec vect[series]  ST(N) Call(N)
set temp1 = 0.
set temp2 = 0.
set temp3 = 0.
set temp4 = 0.
set temp5 = 0.
com mux = 0
nonlin nu  a20 a21 a22 beta10 beta11 beta12 beta20 beta21 beta22 lambda10 lambda11 lambda12 lambda20 lambda21 lambda22 $
alpha11 alpha12 alpha21 alpha22  alpha31 alpha32 alpha41 alpha42 gama11 gama12 gama21 gama22
frml d1 = beta10^2+(beta11^2)*(temp1{1}-gama11)^2+(beta12^2)*temp2{1}+(alpha11^2)*f1{1}^2+(alpha12^2)*f2{1}^2
frml d2 = beta20^2+(beta21^2)*(temp1{1}-gama12)^2+(beta22^2)*temp3{1}+(alpha21^2)*f1{1}^2+(alpha22^2)*f2{1}^2
frml k1 = lambda10^2+(lambda11^2)*(temp1{1}-gama21)^2+(lambda12^2)*temp4{1}+(alpha31^2)*f1{1}^2+(alpha32)^2*f2{1}^2
frml k2 = lambda20^2+(lambda21^2)*(temp1{1}-gama22)^2+(lambda22^2)*temp5{1}+(alpha41^2)*f1{1}^2+(alpha42^2)*f2{1}^2
frml sigmaTx1 = (beta10*lambda10)+(beta11*lambda11)*(temp1{1}-gama11)*(temp1{1}-gama21)+(beta12*lambda12)*sqrt(temp2{1})*sqrt(temp4{1})+(alpha11*alpha31)*f1{1}^2+(alpha12*alpha32)*f2{1}^2
frml sigmaTx2 = (beta20*lambda20)+(beta21*lambda21)*(temp1{1}-gama12)*(temp1{1}-gama22)+(beta22*lambda22)*sqrt(temp3{1})*sqrt(temp5{1})+(alpha21*alpha41)*f1{1}^2+(alpha22*alpha42)*f2{1}^2
frml muT = a20+a21*f1{1}+a22*f2{1}
frml m11 = exp(-nu*mux+nu^2*d1/2)
frml m22 = exp(-nu*mux+nu^2*d2/2)
frml q1 = sqrt(k1+d1-2*sigmaTx1)
frml q2 = sqrt(k2+d2-2*sigmaTx2)
frml m33 =  %cdf((muT-mux-nu*(sigmaTx1-d1))/q1)
frml m44 =  %cdf(-(muT-mux-nu*(sigmaTx2-d2))/q2)
frml m31 =  %cdf((muT-mux)/q1)
frml m41 =  %cdf(-(muT-mux)/q2)
frml c = 1.0/(m31+m41)
frml m1 = c*(m11*m33+m22*m44)
frml m12 = exp((1-nu)*mux+(1-nu)^2*d1/2)
frml m21 = exp((1-nu)*mux+(1-nu)^2*d2/2)
frml m13 =  %cdf((muT-mux+(1-nu)*(sigmaTx1-d1))/q1)
frml m14 =  %cdf(-(muT-mux+(1-nu)*(sigmaTx2-d2))/q2)
frml m2 = c*(m12*m13+m21*m14)
frml psinu0 = log(m2)
frml psinu1 = log(m1)
frml  Rt = r-psinu0+psinu1
frml e = x-Rt
frml g1 = (1.0/sqrt(d1))*%density(e/sqrt(d1))
frml q3 = sqrt(k1-d1^(-1)*sigmaTx1^2)
frml q4 = sqrt(k2+d2-2*sigmaTx2)
frml g21 = -(1+d1^(-1)*sigmaTx1)*e/q3
frml g22 = -(muT-mux)/q1
frml g2 = %cdf(g21+g22)
frml g3 =  %cdf((muT-mux)/q1)
frml g4 =  %cdf(-(muT-mux)/q2)
frml c1 = g3+g4
frml h13 =(1.0/sqrt(d2))*%density(e/sqrt(d2))
frml h21 =  (1+d2^(-1)*sigmaTx2)*e/q4
frml h22 = (muT-mux)/q2
frml h23 = %cdf(h21+h22)
frml L = (temp1=e), (temp2=d1), (temp3 = d2), (temp4 = k1), (temp5 = k2),  log((g1*g2)+(h13*h23))-log(c1)
com nu=0.5, a20=-1.0, a21=-1.0, a22=-0.2, beta10=0.01, beta11= 0.2, beta12=-0.5, beta20=0.01, beta21= -0.1 , beta22= -0.8, $
lambda10=-1.0, lambda11=0.1  , lambda12=-0.3, lambda20=-1.0 , lambda21=-1.0  , lambda22=-0.4 ,$
alpha11=0.1 ,alpha12=0.1, alpha21=-0.5, alpha22=0.01, alpha31=-0.5, alpha32= 0.6, alpha41=1.0, alpha42= 1.0, $
gama11=-0.1, gama12=1.0,  gama21=-1.0 , gama22=-1.0
maximize(pmethod=simplex, piters=13,  iters = 5000, method=bfgs, cvcrit=0.000001 ) L 8 *
set muxstr1 = mux -nu*d1
set muxstr2 = mux -nu*d2
set muTstr1 = muT -nu*sigmaTx1
set muTstr2 = muT -nu*sigmaTx2
set lambda01 = (muTstr1-muxstr1)/q3
set lambda02 = (muTstr2-muxstr2)/q4
set lambda1 = sqrt(d1)*(1+d1^(-1)*sigmaTx1)/q3
set lambda2 = sqrt(d2)*(1+d2^(-1)*sigmaTx2)/q4
set as = -lambda01/sqrt(1+lambda1^2)
set bs = -lambda02/sqrt(1+lambda2^2)
set cs = 1.0/sqrt(1+lambda1^2)
set ds = 1.0/sqrt(1+lambda2^2)
set g3s = exp(-nu*mux+nu^2*d1/2)*%cdf((muTstr1-muxstr1)/q1)
set g4s = exp(-nu*mux+nu^2*d2/2)*%cdf(-(muTstr2-muxstr2)/q2)
set pistr = g3s/(g3s+g4s)
do i = 1,N
set v  =  %ran(1.0)
set simul1 = %uniform(0,1)
set u1 = %invnormal(%cdf(as)+(1-%cdf(as))*simul1)
set u2 = %invnormal(%cdf(bs)*simul1)
set simul2 = %uniform(0,1)
set j1 = %sign(simul2-pistr)
set j2 = j1+1
set j3 = %sign(j2)
set z = cs*(lambda1*u1+v)*(1-j3) + ds*(lambda2*u2+v)*j3
set estr = (sqrt(d1)*z+muxstr1)*(1-j3) + (sqrt(d2)*z+muxstr2)*j3
set  temp22  = 0.
set temp33  = 0.
set temp44 = 0.
set temp55 = 0.
set d11 = beta10^2+(beta11^2)*(estr(t-1)-gama11)^2+(beta12^2)*temp22{1}+(alpha11^2)*f1{1}^2+(alpha12^2)*f2{1}^2
set d22 = beta20^2+(beta21^2)*(estr(t-1)-gama12)^2+(beta22^2)* temp33{1}+(alpha21^2)*f1{1}^2+(alpha22^2)*f2{1}^2
set k11 = lambda10^2+(lambda11^2)*(estr(t-1)-gama21)^2+(lambda12^2)* temp44{1}+(alpha31^2)*f1{1}^2+(alpha32)^2*f2{1}^2
set k22 = lambda20^2+(lambda21^2)*(estr(t-1)-gama22)^2+(lambda22^2)*temp55{1}+(alpha41^2)*f1{1}^2+(alpha42^2)*f2{1}^2
set sigmaTx11 = (beta10*lambda10)+(beta11*lambda11)*(estr(t-1)-gama11)*(estr(t-1)-gama21)+(beta12*lambda12)*sqrt(temp22{1})*sqrt(temp44{1})+(alpha11*alpha31)*f1{1}^2+(alpha12*alpha32)*f2{1}^2
set sigmaTx22 = (beta20*lambda20)+(beta21*lambda21)*(estr(t-1)-gama12)*(estr(t-1)-gama22)+(beta22*lambda22)*sqrt(temp33{1})*sqrt(temp55{1})+(alpha21*alpha41)*f1{1}^2+(alpha22*alpha42)*f2{1}^2
set m_11 = exp(-nu*mux+nu^2*d11/2)
set m_22 = exp(-nu*mux+nu^2*d22/2)
set q11 = sqrt(k11+d11-2*sigmaTx11)
set q22 = sqrt(k22+d22-2*sigmaTx22)
set m_33 = %cdf((muT-mux-nu*(sigmaTx11-d11))/ q11)
set m_44 = %cdf(-(muT-mux-nu*(sigmaTx22-d22))/q22)
set m_31 = %cdf((muT-mux)/q11)
set m_41 = %cdf(-(muT-mux)/q22)
set c2 = 1/(m_31+m_41)
set m_1 = c2*(m_11*m_33+m_22*m_44)
set m_12 = exp((1-nu)*mux+(1-nu)^2*d11/2)
set m_21 = exp((1-nu)*mux+(1-nu)^2*d22/2)
set m_13 =  %cdf((muT-mux+(1-nu)*(sigmaTx11-d11))/q11)
set m_14 =  %cdf(-(muT-mux+(1-nu)*(sigmaTx22-d22))/q22)
set m_2 = c2*(m_12*m_13+m_21*m_14)
set psinustr0 = log(m_2)
set psinustr1 = log(m_1)
do m = 1, M
     do j = 1,J
set Rtstr = r2(m)-psinustr0+psinustr1+estr
set dummy = %if(T2(j) > T1(m), Rtstr, %na)
set ST(i)     = (temp22=d11), (temp33 = d22), (temp44 = k11), (temp55 = k22), S(m)*exp(Rtstr)
set Call(i)   = (temp22=d11), (temp33 = d22), (temp44 = k11), (temp55 = k22), exp(-r2(m)*(T2(j)-T1(m))/365)*%max(0,ST(i)-K(j))
    end do j
end do m
end do i
copy(Format=xls,ORG=col) / ST Call
tab