dec real b a10 a11 a12 a20 a21 a22 beta10 beta11 beta12 beta20 beta21 beta22 lambda10 lambda11 lambda12 lambda20 lambda21 lambda22
dec real alpha10 alpha11 alpha12 alpha20 alpha21 alpha22
OPEN DATA "C: \ Desktop\ heterosc.xls"
CALENDAR 1978 1 12
Allocate 2013:8
DATA(FORMAT=xls,ORG=COL) / f1 f2 x
set temp1 = 0.
set temp2 = 0.
set temp3 = 0.
com b=1.0/sqrt(2*%pi)
com a10= 2.125569 , a11=-0.320445, a12=0.266288, a20=143.615219, a21=-2.889713, a22=0.982924, $
beta10=-0.268960, beta11=0.681791, beta12=-0.000001, beta20=0.961134, beta21=-0.368144, beta22=0.000000, $
lambda10=107.886349, lambda11=-53.421775, lambda12=2.577898, lambda20=-8.229714, lambda21=-1.121924, lambda22=0.106505, $
alpha10 = 0.204003, alpha11= 0.483518, alpha12= 0.322963, alpha20=-9.731820, alpha21=1.050807, alpha22=0.045248

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 = exp(beta10+beta11*e{1}+beta12*(temp2{1}-1-(temp2{1}-1)^2/2))
set d2 = exp(beta20+beta21*e{1}+beta22*(temp3{1}-1-(temp3{1}-1)^2/2))
set k1 = exp(lambda10+lambda11*f1{1}+lambda12*f2{1})
set k2 = exp(lambda20+lambda21*f1{1}+lambda22*f2{1})
set rho1 = (2/(1+exp(alpha10+alpha11*f1{1}+alpha12*f2{1})))-1
set rho2 = (2/(1+exp(alpha20+alpha21*f1{1}+alpha22*f2{1})))-1
set s1 = sqrt(d1)*sqrt(k1)*rho1
set s2 = sqrt(d2)*sqrt(k2)*rho2
set k11 = %cdf((muT-mux)/sqrt(k1+d1-2*s1))
set k22 = %cdf(-(muT-mux)/sqrt(k2+d2-2*s2))
set c = 1.0/(k11+k22)
set sigmaT1str = sqrt(k1+d1-2*s1)
set sigmaT2str = sqrt(k2+d2-2*s2)
set sigmaTx1str = s1-d1
set sigmaTx2str = s2-d2
set muTstr = muT-mux
set h1 = (muTstr)/sigmaT1str
set h2 = (muTstr)/sigmaT2str
set t1 = %if(h1>=0,1,0)
set t2 = %if(h2>=0,1,0)

set i0str = 1/2*(1+%sign(h1)*%gammainc(h1^2/2,1/2))
set i0strstr = 1/2*(1+(-1)^(t2)*%gammainc(h2^2/2,1/2))
set i1str = b*(1+%sign(h1)*(-1)^(t1)*%gammainc(h1^2/2,1))
set i1strstr = b*(-1+%gammainc(h2^2/2,1))
set i2str = 1/2*(1+%sign(h1)*%gammainc(h1^2/2,3/2))
set i2strstr =1/2*(1+(-1)^(3*t2)*%gammainc(h2^2/2,3/2))
set i3str = 2*b*(1+%sign(h1)*(-1)^3*(t1)*%gammainc(h1^2/2,2))
set i3strstr = 2*b*(-1+%gammainc(h2^2/2,2))
set i4str = 3/2*(1+%sign(h1)*%gammainc(h1^2/2,5/2))
set i4strstr = 3/2*(1+(-1)^(5*t2)*%gammainc(h2^2/2,5/2))

set m1 = mux-(sigmaT1str^(-1))*sigmaTx1str*i1str+(sigmaT2str^(-1))*sigmaTx2str*i1strstr

set a1 = (mux-m1)^2
set a2 = 2*(mux-m1)*((sigmaT1str)^(-1)*sigmaTx1str)*i1str
set a3 = 2*(mux-m1)*((sigmaT2str)^(-1)*sigmaTx2str)*i1strstr
set a4 = -((sigmaT1str)^(-2)*(sigmaTx1str^2)*i2str+(sigmaT2str)^(-2)*(sigmaTx2str^2)*i2strstr)
set a5 = -(d1-sigmaT1str^(-2)*(sigmaTx1str^2)*i0str)
set a6 = -(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i0strstr
set v1 = c*(a1+a2+a3+a4+a5+a6)

set aa1 = (mux-m1)^3
set aa2 = -3*(mux-m1)^2*((sigmaT1str)^(-1)*sigmaTx1str*i1str+(sigmaT2str)^(-1)*sigmaTx2str*i1strstr)
set aa3 = -3*(mux-m1)*(sigmaT1str)^(-2)*(sigmaTx1str^2)*i2str+(sigmaT2str^(-2))*(sigmaTx2str^2)*i2strstr
set aa4 = 3*(mux -m1)*(d1-(sigmaT1str^(-2))*(sigmaTx1str^2))*i0str
set aa5 = 3*(mux -m1)*(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i0strstr
set aa6 = (sigmaT1str)^(-3)*(sigmaTx1str)^3*i3str+(sigmaT2str)^(-3)*(sigmaTx2str)^3*i3strstr
set aa7 = (sigmaT1str)^(-1)*sigmaTx1str*(d1-(sigmaT1str^(-2))*(sigmaTx1str^2))*i1str
set aa8 = (sigmaT2str)^(-1)*sigmaTx2str*(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i1strstr
set aa9 = c*(aa1+aa2+aa3+aa4+aa5+aa6+aa7+aa8)
set beta1 = aa9/(v1)^(1.5)

set aaa1 = (mux-m1)^4
set aaa2 = -4*(mux -m1)^3*((sigmaT1str)^(-1)*sigmaTx1str*i1str+(sigmaT2str)^(-1)*sigmaTx2str*i1strstr)
set aaa3 = 6*(mux -m1)^2*((sigmaT1str)^(-2)*(sigmaTx1str^2)*i2str+(sigmaT2str^(-2))*(sigmaTx2str^2)*i2strstr)
set aaa4 = 6*(mux -m1)^2*(d1-(sigmaT1str^(-2))*(sigmaTx1str^2))*i0str
set aaa5 = 6*(mux -m1)^2*(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i0strstr
set aaa6 = -4*(mux -m1)*(sigmaT1str)^(-3)*(sigmaTx1str)^3*i3str+(sigmaT2str)^(-3)*(sigmaTx2str)^3*i3strstr
set aaa7 = -12*(mux -m1)*(sigmaT1str)^(-1)*sigmaTx1str*(d1-(sigmaT1str^(-2))*(sigmaTx1str^2))*i1str
set aaa8 = -12*(mux -m1)*(sigmaT2str)^(-1)*sigmaTx2str*(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i1strstr
set aaa9 = -(sigmaT1str)^(-4)*(sigmaTx1str)^4*i4str-(sigmaT2str)^(-4)*(sigmaTx2str)^4*i4strstr
set aaa10 = -6*(sigmaT1str)^(-2)*(sigmaTx1str^2)*(d1-(sigmaT1str^(-2))*(sigmaTx1str^2))*i2str
set aaa10 = -6*(sigmaT2str)^(-2)*(sigmaTx2str^2)*(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i2strstr
set aaa11 = -3*((d1-(sigmaT1str^(-2))*(sigmaTx1str^2))*i0str+(d2-(sigmaT2str^(-2))*(sigmaTx2str^2))*i0strstr)
set aaa12 = c*(aaa1+aaa2+aaa3+aaa4+aaa5+aaa6+aaa7+aaa8+aaa9+aaa10+aaa11+aaa12)
set theta1 = aaa12/v1^2

spgraph(hea='Graphs of the Four Conditional Moments', hfi=2, vfi=2)
gra(hea='Panel 1: Time path of mean') 1 ;  # m1
gra(hea='Panel 2: Time path of variance') 1 ; # v1
gra(hea='Panel 3: Time path of skewness') 1 ; # beta1
gra(hea='Panel 4: Time path of kurtosis') 1 ; # theta1
spgraph(done)