com M = 412,  N = 100
dec vect[series] ST(N) Call(N) v(M)  simul1(M)  simul2(M)  u1(M)  u2(M)  v1(M)  v2(M)  j1(M)  j2(M)  j3(M) s(M)  estr(M)
dec vect[series] psinustr1(M) psinustr2(M) psinustr(M) d11(M) d22(M) k11(M) k22(M) sigmaTx11(M) sigmaTx22(M)
dec vect[series] mstr1(M) mstr2(M) mstr3(M) mstr4(M) mstr31(M) mstr41(M) Rtstr(M)
all M
com K = 40.0, S0 = 40.0
OPEN DATA "C: \ Desktop\ heterosc_option.xls"
CALENDAR 1978 1 12
DATA(FORMAT=xls,ORG=COL) / f1 f2 y r
set temp1 = 0.
set temp2 = 0.
set temp3 = 0.
set temp4 = 0.
set temp5 = 0.
com mux = 0
set x = 100*log(y/y{1})
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 m33 =  %cdf((muT-mux-nu*(sigmaTx1-d1))/sqrt(k1+d1-2*sigmaTx1))
frml m44 =  %cdf(-(muT-mux-nu*(sigmaTx2-d2))/sqrt(k2+d2-2*sigmaTx2))
frml m31 =  %cdf((muT-mux)/sqrt(k1+d1-2*sigmaTx1))
frml m41 =  %cdf(-(muT-mux)/sqrt(k2+d2-2*sigmaTx2))
frml c = 1/(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))/sqrt(k1+d1-2*sigmaTx1))
frml m14 =  %cdf(-(muT-mux+(1-nu)*(sigmaTx2-d2))/sqrt(k2+d2-2*sigmaTx2))
frml m31 =  %cdf((muT-mux)/sqrt(k1+d1-2*sigmaTx1))
frml m41 =  %cdf(-(muT-mux)/sqrt(k2+d2-2*sigmaTx2))
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 g21 =  -(1+d1^(-1)*sigmaTx1)*e/sqrt(k1-d1^(-1)*sigmaTx1^2)
frml g22 = -(muT-mux)/sqrt(k1-d1^(-1)*sigmaTx1^2)
frml g2 = %cdf(g21+g22)
frml g3 =  %cdf((muT-mux)/sqrt(k1+d1-2*sigmaTx1))
frml g4 =  %cdf(-(muT-mux)/sqrt(k2+d2-2*sigmaTx2))
frml c1 = g3+g4
frml h13 =(1.0/sqrt(d2))*%density(e/sqrt(d2))
frml h21 =  (1+d2^(-1)*sigmaTx2)*e/sqrt(k2-d2^(-1)*sigmaTx2^2)
frml h22 = (muT-mux)/sqrt(k2-d2^(-1)*sigmaTx2^2)
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=-0.1, a22=-0.2, beta10=-0.01, beta11= 0.2, beta12=-0.3, 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.01 ,alpha12=-0.01, alpha21=-0.05, 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=50,  iters = 5000, method=bfgs, cvcrit=0.000001 ) L 12 *
set muxstr1 = mux -nu*d1
set muxstr2 = mux -nu*d2
set lambda01 = (muT-muxstr1)/sqrt(k1-d1^(-1)*sigmaTx1^2)
set lambda02 = (muT-muxstr2)/sqrt(k2-d2^(-1)*sigmaTx2^2)
set lambda1 = sqrt(d1)*(1+d1^(-1)*sigmaTx1)/sqrt(k1-d1^(-1)*sigmaTx1^2)
set lambda2 = sqrt(d2)*(1+d2^(-1)*sigmaTx2)/sqrt(k2-d2^(-1)*sigmaTx2^2)
set as = -lambda01/sqrt(1+lambda1^2)
set bs = -lambda02/sqrt(1+lambda2^2)
set cs = 1/sqrt(1+lambda1^2)
set ds = 1/sqrt(1+lambda2^2)
set g3s = %cdf((muT-muxstr1)/sqrt(k1+d1-2*sigmaTx1))
set g4s =%cdf(-(muT-muxstr2)/sqrt(k2+d2-2*sigmaTx2))
set pis = g3s/(g3s+g4s)
set pistr11 = pis*exp(-nu*mux+nu^2*d1/2)
set pistr22 = (1-pis)*exp(-nu*mux+nu^2*d2/2)
set g33s =  %cdf((muT-mux-nu*(sigmaTx1-d1))/sqrt(k1+d1-2*sigmaTx1))
set g44s =  %cdf(-(muT-mux-nu*(sigmaTx2-d2))/sqrt(k2+d2-2*sigmaTx2))
set pistr1 = pistr11/(pistr11*g33s/g3s+pistr22*g44s/g4s)
set pistr2 = pistr22/(pistr11*g33s/g3s+pistr22*g44s/g4s)


set temp22 = 0.
set temp33 = 0.
set temp44 = 0.
set temp55 = 0.
report(action = define, hlabels= || 'ST', 'Call'||)
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-pistr1)
set j2 = j1+1
set j3 = %sign(j2)
set s = cs*(lambda1*u1+v)*(1-j3) + ds*(lambda2*u2+v)*j3
set estr =  (sqrt(d1)*s+muxstr1)*(1-j3) + (sqrt(d2)*s+muxstr2)*j3
set d11 = beta10^2+(beta11^2)*(estr-gama11)^2+(beta12^2)*temp22{1}+(alpha11^2)*f1{1}^2+(alpha12^2)*f2{1}^2
set d22 = beta20^2+(beta21^2)*(estr-gama12)^2+(beta22^2)* temp33{1}+(alpha21^2)*f1{1}^2+(alpha22^2)*f2{1}^2
set k11 = lambda10^2+(lambda11^2)*(estr-gama21)^2+(lambda12^2)* temp44{1}+(alpha31^2)*f1{1}^2+(alpha32)^2*f2{1}^2
set k22 = lambda20^2+(lambda21^2)*(estr-gama22)^2+(lambda22^2)*temp55{1}+(alpha41^2)*f1{1}^2+(alpha42^2)*f2{1}^2
set sigmaTx11 = (beta10*lambda10)+(beta11*lambda11)*(estr-gama11)*(estr-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-gama12)*(estr-gama22)+(beta22*lambda22)*sqrt(temp33{1})*sqrt(temp55{1})+(alpha21*alpha41)*f1{1}^2+(alpha22*alpha42)*f2{1}^2
set mstr1 = exp(-nu*muxstr1+nu^2*d11/2)
set mstr2 = exp(-nu*muxstr2+nu^2*d22/2)
set mstr3 =  %cdf((muT-mux-(1-nu)*(sigmaTx11-d11))/sqrt(k11+d11-2*sigmaTx11))
set mstr4 =  %cdf(-(muT-mux-(1-nu)*(sigmaTx22-d22))/sqrt(k22+d22-2*sigmaTx22))
set mstr31 =  %cdf((muT-mux)/sqrt(k11+d11-2*sigmaTx11))
set mstr41 =  %cdf(-(muT-mux)/sqrt(k22+d22-2*sigmaTx22))
set psinustr1 = mstr1*mstr3*exp(-muxstr1+d11/2)
set psinustr2 = mstr2*mstr4*exp(-muxstr2+d22/2)
set psinustr = log(pistr1*psinustr1/mstr31+pistr2*psinustr2/mstr41)
set Rtstr = r-psinustr+estr
set ST(i) = (temp22=d11), (temp33 = d22), (temp44 = k11), (temp55 = k22), S0*exp(Rtstr)
set Call(i) = (temp22=d11), (temp33 = d22), (temp44 = k11), (temp55 = k22), exp(-r)*%max(0,ST(i)-K)
              report(row=current, atcol=1, align=decimal) ST(i)
              report(row=current, atcol=2, align = decimal) Call(i)
             report(action = format, align = center, width = 9)
             report(row = new, atcol = 1) ST Call
             end do
copy(Format=xls,ORG=Columns) / ST Call