dec real a10 a11 a12 a20 a21 a22 beta10 beta11 beta 12 beta20 beta21 beta 22 lambda10 lambda11 lambda12 lambda20 lambda21 lambda22
dec real alpha11 alpha12 alpha21 alpha22  alpha31 alpha32 alpha41 alpha42 gama11 gama12 gama21 gama22  gama31 gama32
OPEN DATA "C: \ Desktop\ heterosc.xls"
CALENDAR 1978 1 12
DATA(FORMAT=xls,ORG=COL) / f1 f2 x
set temp1 = 0.
set temp2 = 0.
set temp3 = 0.
set temp4 = 0.
set temp5 = 0.
set temp6 = 0.
set temp7 = 0.

nonlin a10 a11 a12 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  e = x-a10-a11*x{1}-a12*temp1{1}
frml mux = a10+a11*x{1}+a12*e{1}
frml muT = a20+a21*f1{1}+a22*f2{1}
frml d1 = beta10^2+(beta11^2)*(e{1}-gama11)^2+(beta12^2)*temp2{1}+(alpha11^2)*f1{1}^2+(alpha12^2)*f2{1}^2
frml d2 = beta20^2+(beta21^2)*(e{1}-gama12)^2+(beta22^2)*temp3{1}+(alpha21^2)*f1{1}^2+(alpha22^2)*f2{1}^2
frml k1 = lambda10^2+(lambda11^2)*(e{1}-gama21)^2+(lambda12^2)*temp4{1}+(alpha31^2)*f1{1}^2+(alpha32)^2*f2{1}^2
frml k2 = lambda20^2+(lambda21^2)*(e{1}-gama22)^2+(lambda22^2)*temp5{1}+(alpha41^2)*f1{1}^2+(alpha42^2)*f2{1}^2
frml sigmaTx1 = (beta10*lambda10)+(beta11*lambda11)*(e{1}-gama11)*(e{1}-gama21)+(beta12*lambda12)*temp6{1}+(alpha11*alpha31)*f1{1}^2+(alpha12*alpha32)*f2{1}^2
frml sigmaTx2 = (beta20*lambda20)+(beta21*lambda21)*(e{1}-gama12)*(e{1}-gama22)+(beta22*lambda22)*temp7{1}+(alpha21*alpha41)*f1{1}^2+(alpha22*alpha42)*f2{1}^2
frml L = $
             g1 = (1.0/sqrt(d1))*%density(e/sqrt(d1)),$
             g21 = -((1+((d1^(-1))*sigmaTx1))*e)/sqrt(k1-((d1^(-1))*(sigmaTx1^2))),$
             g22 = -(muT-mux)/sqrt(k1-((d1^(-1))*(sigmaTx1^2))),$
             g2 = %cdf(g21+g22),$
             g3 = %cdf((muT-mux)/sqrt(k1+d1-2*sigmaTx1)),$
             g4 = %cdf(-(muT-mux)/sqrt(k2+d2-2*sigmaTx2)),$
             h10= (1.0/sqrt(d2))*%density(e/sqrt(d2)),$
             h21 = ((1+((d2^(-1))*sigmaTx2))*e)/sqrt(k2-((d2^(-1))*(sigmaTx2^2))),$
             h22 = (muT-mux)/sqrt(k2-((d2^(-1))*(sigmaTx2^2))),$
             h20 = %cdf(h21+h22),$
             (temp1=e), (temp2=d1), (temp3=d2), (temp4=k1), (temp5=k2),(temp6=sigmaTx1), (temp7=sigmaTx2), log((g1*g2)+(h10*h20))-log(g3+g4)
             box(constant, ar=1,ma=1,noprint) x
             com a10=%beta(1), a11=%beta(2), a12=%beta(3), a20=%beta(1), a21=%beta(2), a22=%beta(3), $
             beta10=-5.0, beta11=-0.5, beta12 = 0.0, beta20=3.0, beta21= 0.1, beta22=-0.6, $
             lambda10=0.1, lambda11=-0.0, lambda12= 0.0, lambda20=0.0, lambda21=0.0, lambda22=0.0, $
             alpha11=-0.2, alpha12=0.3 , alpha21=-0.2, alpha22=0.2, $
            alpha31=0.0, alpha32=0.1, alpha41=0.1, alpha42=0.2, $
            gama11 =-0.6, gama12 =-0.6, gama21=-0.2, gama22 =0.2
            maximize(pmethod=simplex, piters=6, iters = 5000, method=bfgs, cvcrit=0.000001 ) L 12 *

set e = x-a10-a11*x{1}-a12*temp1{1}
set mux = a10+a11*x{1}+a12*e{1}
set muT = a20+a21*f1{1}+a22*f2{1}
set d1 = beta10^2+(beta11^2)*(e{1}-gama11)^2+(beta12^2)*temp2{1}+(alpha11^2)*f1{1}^2+(alpha12^2)*f2{1}^2
set et sigmax1 = (temp1=e), (temp2=d1),  sqrt(d1)
set d2 = beta20^2+(beta21^2)*(e{1}-gama12)^2+(beta22^2)*temp3{1}+(alpha21^2)*f1{1}^2+(alpha22^2)*f2{1}^2
set sigmax2 = (temp1 = e), (temp3 = d2), sqrt(d2)
set k1 = lambda10^2+(lambda11^2)*(e{1}-gama21)^2+(lambda12^2)*temp4{1}+(alpha31^2)*f1{1}^2+(alpha32)^2*f2{1}^2
set sigmaT1 = (temp1 = e), (temp4 = k1), sqrt(k1)
set k2 = lambda20^2+(lambda21^2)*(e{1}-gama22)^2+(lambda22^2)*temp5{1}+(alpha41^2)*f1{1}^2+(alpha42^2)*f2{1}^2
set sigmaT2 = (temp1 = e), (temp5 = k2), sqrt(k2)
set sigmaTx1 =  (beta10*lambda10)+(beta11*lambda11)*(e{1}-gama11)*(e{1}-gama21)+(beta12*lambda12)*temp6{1}+(alpha11*alpha31)*f1{1}^2+(alpha12*alpha32)*f2{1}^2
set sigmaTx2 = (beta20*lambda20)+(beta21*lambda21)*(e{1}-gama12)*(e{1}-gama22)+(beta22*lambda22)*temp7{1}+(alpha21*alpha41)*f1{1}^2+(alpha22*alpha42)*f2{1}^2
set rho1 = (temp1=e), (temp6 = sigmaTx1), sigmaTx1/(sigmax1*sigmaT1)
set rho2 = (temp1=e), (temp7 = sigmaTx2), sigmaTx2/(sigmax2*sigmaT2)
set k11 = %cdf((muT-mux)/sqrt(k1+sigmax1^2-2*sigmaTx1))
set k22 = %cdf(-(muT-mux)/sqrt(k2+sigmax2^2-2*sigmaTx2))
set c = 1.0/(k11+k22)
set sigmaT1str = sqrt(k1+sigmax1^2-2*sigmaTx1)
set sigmaT2str = sqrt(k2+sigmax2^2-2*sigmaTx2)
set sigmaTx1str = sigmaTx1-sigmax1^2
set sigmaTx2str = sigmaTx2-sigmax2^2
set muTstr = muT-mux
set h1 = muTstr/sigmaT1str
set h2 = muTstr/sigmaT2str
set t1 = %if(h1>=0.0,1.0,0.0)
set t2 = %if(h2>=0.0,1.0,0.0)
com b=1.0/sqrt(2.0*%pi)

set i0str = 0.5*(1.0+%sign(h1)*%gammainc(h1^2/2,0.5))
set i0strstr = 0.5*(1.0+(-1)^(-t2)*%gammainc(h2^2/2,0.5))
set i1str = b*(1.0+%sign(h1)*(-1)^(t1)*%gammainc(h1^2/2,1.0))
set i1strstr = b*(-1.0+%gammainc(h2^2/2,1.0))
set i2str = b*sqrt(2.0)*%gamma(3.0/2)*(1+%sign(h1)*%gammainc(h1^2/2,3.0/2))
set i2strstr = b*sqrt(2.0)*%gamma(3.0/2)*(1+(-1)^(t2)*%gammainc(h2^2/2,3.0/2))
set i3str = b*2.0*%gamma(2.0)*(1+%sign(h1)*(-1)^(t1)*%gammainc(h1^2/2,2))
set i3strstr = b*2.0*%gamma(2.0)*(-1+%gammainc(h2^2/2,2))
set i4str = b*2*sqrt(2.0)*%gamma(5.0/2)*(1+%sign(h1)*%gammainc(h1^2/2,5.0/2))
set i4strstr = b*2*sqrt(2.0)*%gamma(5.0/2)*(1+(-1)^(t2)*%gammainc(h2^2/2,5.0/2))

set m11 = ((sigmaT1str^(-2))*sigmaTx1str)*(sigmaT1str*i1str+muTstr*i0str)
set m12 = ((sigmaT2str^(-2))*sigmaTx2str)*(sigmaT2str*i1strstr+muTstr*i0strstr)
set m1 = mux+(m11+m12)/(i0str+i0strstr)

set a1 = (mux-m1)^2*(i0str+i0strstr)
set a2 = 2*(mux-m1)*((sigmaT1str)^(-2)*sigmaTx1str)*(sigmaT1str*i1str+muTstr*i0str)
set a3 = 2*(mux-m1)*((sigmaT2str)^(-2)*sigmaTx2str)*(sigmaT2str*i1strstr+muTstr*i0strstr)
set a4 = (sigmaT1str^(-2)*sigmaTx1str)^2
set a5 = (sigmaT1str^2)*i2str+2*muTstr*sigmaT1str*i1str+(muTstr^2)*i0str
set a6 = (sigmaT2str^(-2)*sigmaTx2str)^2
set a7 = (sigmaT2str^2)*i2strstr+2*muTstr*sigmaT2str*i1strstr+(muTstr^2)*i0strstr
set a8 = (sigmax1^2-sigmaT1str^(-2)*sigmaTx1str^2)*i0str
set a9 = (sigmax2^2-sigmaT2str^(-2)*sigmaTx2str^2)*i0strstr
set v1 = c*(a1+a2+a3+(a4*a5)+(a6*a7)+a8+a9)
set d = sqrt(v1)

set aa1 = (mux-m1)^3*(i0str+i0strstr)
set aa2 = 3*(mux-m1)^2*(sigmaT1str^(-2))*sigmaTx1str*(sigmaT1str*i1str+muTstr*i0str)
set aa3 = 3*(mux-m1)^2*(sigmaT2str^(-2))*sigmaTx2str*(sigmaT2str*i1strstr+muTstr*i0strstr)
set aa4 = 3*(mux-m1)*(sigmaT1str^(-2)*sigmaTx1str)^2
set aa5 = sigmaT1str^2*i2str+2*muTstr*sigmaT1str*i1str+muTstr^2*i0str
set aa6 = 3*(mux -m1)*(sigmax1^2-sigmaT1str^(-2)*sigmaTx1str^2)*i0str
set aa7 = 3*(mux-m1)*(sigmaT2str^(-2)*sigmaTx2str)^2
set aa8 = (sigmaT2str^2)*i2strstr+2*muTstr*sigmaT2str*i1strstr+(muTstr^2)*i0strstr
set aa9 = 3*(mux -m1)*(sigmax2^2-sigmaT2str^(-2)*sigmaTx2str^2)*i0strstr
set aa10 = 3*(sigmaT1str^(-2))*sigmaTx1str*(sigmax1^2-sigmaT1str^(-2)*sigmaTx1str^2)
set aa11 = sigmaT1str*i1str+muTstr*i0str
set aa12 = 3*(sigmaT2str^(-2))*sigmaTx2str*(sigmax2^2-sigmaT2str^(-2)*sigmaTx2str^2)
set aa13 = sigmaT2str*i1strstr+muTstr*i0strstr
set aa14 = (sigmaT1str^(-2)*sigmaTx1str)^3
set aa15 = sigmaT1str^3*i3str+3*muTstr*sigmaT1str^2*i2str+3*muTstr^2*sigmaT1str*i1str+muTstr^3*i0str
set aa16 = ((sigmaT2str^(-2))*sigmaTx2str)^3
set aa17 = sigmaT2str^3*i3strstr+3*muTstr*sigmaT2str^2*i2strstr+3*muTstr^2*sigmaT2str*i1strstr+muTstr^3*i0strstr
set aa18 = c*(aa1+aa2+aa3+(aa4*aa5)+aa6+(aa7*aa8)+aa9+(aa10*aa11)+(aa12*aa13)+(aa14*aa15)+(aa16*aa17))
set beta1 = aa18/(v1^1.5)

set aaa1 = (mux-m1)^4*(i0str+i0strstr)
set aaa2 = 4*(mux -m1)^3*((sigmaT1str)^(-2)*sigmaTx1str)*(sigmaT1str*i1str+muTstr*i0str)
set aaa3 = 4*(mux -m1)^3*((sigmaT2str)^(-2)*sigmaTx2str)*(sigmaT2str*i1strstr+muTstr*i0strstr)
set aaa4 = 6*(mux -m1)^2*((sigmaT1str)^(-2)*sigmaTx1str)^2
set aaa5 = (sigmaT1str^2)*i2str+2*muTstr*sigmaT1str*i1str+(muTstr^2)*i0str
set aaa6 = 6*(mux -m1)^2*(sigmax1^2-(sigmaT1str^(-2))*(sigmaTx1str^2))*i0str
set aaa7 = 6*(mux -m1)^2*((sigmaT2str)^(-2)*sigmaTx2str)^2
set aaa8 = (sigmaT2str^2)*i2strstr+2*muTstr*sigmaT2str*i1strstr+(muTstr^2)*i0strstr
set aaa9 = 6*(mux -m1)^2*(sigmax2^2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i0strstr
set aaa10 = 4*(mux -m1)*(((sigmaT1str^(-2))*sigmaTx1str))^3*(sigmaT1str^3*i3str+3*muTstr*sigmaT1str^2*i2str+3*muTstr^2*sigmaT1str*i1str+muTstr^3*i0str)
set aaa11 = 12*(mux -m1)*sigmaT1str^(-2)*sigmaTx1str
set aaa12 = (sigmax1^2-sigmaT1str^(-2)*sigmaTx1str^2)*(sigmaT1str*i1str+muTstr*i0str)
set aaa13 = 4*(mux -m1)*((sigmaT2str^(-2))*sigmaTx2str)^3*(sigmaT2str^3*i3strstr+3*muTstr*sigmaT2str^2*i2strstr+3*muTstr^2*sigmaT2str*i1strstr+muTstr^3*i0strstr)
set aaa14 = 12*(mux -m1)*sigmaT2str^(-2)*sigmaTx2str
set aaa15 = (sigmax2^2-sigmaT2str^(-2)*sigmaTx2str^2)*(sigmaT2str*i1strstr+muTstr*i0strstr)
set aaa16 = (sigmaT1str^(-2)*sigmaTx1str)^4
set aaa17 = sigmaT1str^4*i4str+4*muTstr*sigmaT1str^3*i3str+6*muTstr^2*sigmaT1str^2*i2str+4*muTstr^3*sigmaT1str*i1str+muTstr^4*i0str
set aaa18 = (sigmaT2str^(-2)*sigmaTx2str)^4
set aaa19 = sigmaT2str^4*i4strstr+4*muTstr*sigmaT2str^3*i3strstr+6*muTstr^2*sigmaT2str^2*i2strstr+4*muTstr^3*sigmaT2str*i1strstr+muTstr^4*i0strstr
set aaa20 = 6*(sigmaT1str^(-2)*sigmaTx1str)^2
set aaa21 = (sigmax1^2-sigmaT1str^(-2)*sigmaTx1str^2)*(sigmaT1str^2*i2str+2*muTstr*sigmaT1str*i1str+muTstr^2*i0str)
set aaa22 = 6*(sigmaT2str^(-2)*sigmaTx2str)^2
set aaa23 = (sigmax2^2-sigmaT2str^(-2)*sigmaTx2str^2)*(sigmaT2str^2*i2strstr+2*muTstr*sigmaT2str*i1strstr+muTstr^2*i0strstr)
set aaa24 = 3*(sigmax1^2-sigmaT1str^(-2)*sigmaTx1str^2)^2*i0str
set aaa25 = 3*(sigmax2^2-sigmaT2str^(-2)*sigmaTx2str^2)^2*i0strstr
set aaa26 = c*(aaa1+aaa2+aaa3+(aaa4*aaa5)+aaa6+(aaa7*aaa8)+aaa9+aaa10+(aaa11*aaa12)+aaa13+(aaa14*aaa15)+(aaa16*aaa17)+(aaa18*aaa19)+(aaa20*aaa21)+(aaa22*aaa23)+aaa24+aaa25)
set theta1 = aaa26/(v1^2)

spgraph(hea='Graphs of the Four Conditional Moments', hfi=5, vfi=2)
gra(hea='Panel 1: Time path of return') 1 ; # x
gra(hea='Panel 2: Time path of Conditional Volatility') 1 ; # d
gra(hea='Panel 3: Time path of standard Deviation_x1') 1 ; # sigmax1
gra(hea='Panel 4: Time path of standard Deviation_x2') 1;  # sigmax2
gra(hea='Panel 5: Time path of standard Deviation_tau1') 1 ; # sigmaT1
gra(hea='Panel 6: Time path of standard Deviation_tau2') 1;  # sigmaT2
gra(hea='Panel 7: Time path of Correlation_rho1') 1;  # rho1
gra(hea='Panel 8: Time path of Correlation_rho2') 1;  # rho2
gra(hea='Panel 9: Time path of Conditional Skewness') 1;  # beta1
gra(hea='Panel 10: Time path of Conditional Kurtosis') 1;  # theta1
spgraph(done)