/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ OUTLINE OF THE MAIN PROCEDURES: * 1. estimate model using first t obs 2. simulate x_sim(t+tao) N times for tao-step ahead prediction, based on parameters estimated step 1. 3. calculate simulated conditional distribution prob(u1=v[i+rows(tao),1]; V_sup=V_sup|v[i,1]; v1=v1|v[i,2]; v2=v2|v[i,3]; else; V_sup=V_sup|v[i+rows(tao),1]; v1=v1|v[i+rows(tao),2]; v2=v2|v[i+rows(tao),3]; endif; endfor; print "The CS V_sup~v1~v2 for tao=1,2,3,4,5, 6, 12 are: " V_sup~v1~v2; // doing b times to get the Bootstrap critical value for CS option=1; //using the same parameters as from X(t) //option=2; // using different parameters in Bootstrap mod_ind1=1;//CIR as the benchmark lval=5; //block size V_sup={}; SimInd=3;//dont simulate latent for boot(1,100,1); //bootstrap replications xxt_boot= boot1(xxt,lval,see); v_boot={}; if SimInd==1; // skip this becoz it is tedious to simulate a long series for SV each time each step elseif SimInd==2; v_boot= v_boot~(CS_statBoot2(xxt_boot,R,N,tao,S,mod_ind1,1,hh,u_bar,svSim1,mvSim1,svSim2,mvSim2,option)); //CIR, recenter,P=(1-a)*T, so sqrt(P/T)=sqrt(1-a) v_boot= v_boot~(CS_statBoot2(xxt_boot,R,N,tao,S,mod_ind1,3,hh,u_bar,svSim1,mvSim1,svSim2,mvSim2,option));//SVJ elseif SimInd==3; v_boot= v_boot~(CS_statBoot3(xxt_boot,R,N,tao,S,mod_ind1,2,hh,u_bar,option)); //SV, recenter,P=(1-a)*T, so sqrt(P/T)=sqrt(1-a) v_boot= v_boot~(CS_statBoot3(xxt_boot,R,N,tao,S,mod_ind1,3,hh,u_bar,option));//SVJ endif; V_sup=V_sup|(maxc(v_boot'))'; print "V_sup for tao=1,2,4,12 are: " V_sup; print "boot=" boot; see=see+100;// change see to change the rndu in boot endfor; print " the bootsrap CS stat are: " V_sup; print " the 5%, 10%,15%,20% critical values for tao=1,2 4,12 are : " ; for loop(1,7,1); V_b=sortc(V_sup,loop); v_b[95,loop]~v_b[90,loop]~v_b[85,loop]~v_b[80,loop]; endfor; //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ " for u_bar=(meanc(xxt)-1*stdc(xxt))|(meanc(xxt)+1*stdc(xxt))"; //the interval u_bar see=12343; u_bar=(meanc(xxt)-1*stdc(xxt))|(meanc(xxt)+1*stdc(xxt));//the interval u_bar // CS model selection test statistic // CS model selection test statistic mod_ind1=1;//CIR as the benchmark v={};//hold vp~v1~v2; vp={};// D(k,P,N), hold test statistic for all 2 alternative models v1={};// D1(k,P,N), v2={};// D2(k,P,N), V_sup={};// sup(D(k,P,N)) /* */ if SimInd==1; v= v|CS_stat(xxt,R,N,tao,S,mod_ind1,1,hh,u_bar); //Against CIR v= v|CS_stat(xxt,R,N,tao,S,mod_ind1,3,hh,u_bar); //Against SVJ elseif SimInd==2; v= v|CS_stat2(xxt,R,N,tao,S,mod_ind1,1,hh,u_bar,svSim1,mvSim1,svSim2,mvSim2); //against CIR v= v|CS_stat2(xxt,R,N,tao,S,mod_ind1,3,hh,u_bar,svSim1,mvSim1,svSim2,mvSim2); //Against SVJ elseif SimInd==3; v= v|CS_stat3(xxt,R,N,tao,S,mod_ind1,2,hh,u_bar); //against SV v= v|CS_stat3(xxt,R,N,tao,S,mod_ind1,3,hh,u_bar); //Against SVJ endif; print " The CS stat D=D1-D2 for tao=1,2,3,4,5, 6, 12 are:" v; for i(1,rows(tao),1); if v[i,1]>=v[i+rows(tao),1]; V_sup=V_sup|v[i,1]; v1=v1|v[i,2]; v2=v2|v[i,3]; else; V_sup=V_sup|v[i+rows(tao),1]; v1=v1|v[i+rows(tao),2]; v2=v2|v[i+rows(tao),3]; endif; endfor; print "The CS V_sup~v1~v2 for tao=1,2,3,4,5, 6, 12 are: " V_sup~v1~v2; // doing b times to get the Bootstrap critical value for CS option=1; //using the same parameters as from X(t) //option=2; // using different parameters in Bootstrap mod_ind1=1;//CIR as the benchmark lval=5; //block size V_sup={}; SimInd=3;//dont simulate latent for boot(1,100,1); //bootstrap replications xxt_boot= boot1(xxt,lval,see); v_boot={}; if SimInd==1; // skip this becoz it is tedious to simulate a long series for SV each time each step elseif SimInd==2; v_boot= v_boot~(CS_statBoot2(xxt_boot,R,N,tao,S,mod_ind1,1,hh,u_bar,svSim1,mvSim1,svSim2,mvSim2,option)); //CIR, recenter,P=(1-a)*T, so sqrt(P/T)=sqrt(1-a) v_boot= v_boot~(CS_statBoot2(xxt_boot,R,N,tao,S,mod_ind1,3,hh,u_bar,svSim1,mvSim1,svSim2,mvSim2,option));//SVJ elseif SimInd==3; v_boot= v_boot~(CS_statBoot3(xxt_boot,R,N,tao,S,mod_ind1,2,hh,u_bar,option)); //SV, recenter,P=(1-a)*T, so sqrt(P/T)=sqrt(1-a) v_boot= v_boot~(CS_statBoot3(xxt_boot,R,N,tao,S,mod_ind1,3,hh,u_bar,option));//SVJ endif; V_sup=V_sup|(maxc(v_boot'))'; print "V_sup for tao=1,2,4,12 are: " V_sup; print "boot=" boot; see=see+100;// change see to change the rndu in boot endfor; print " the bootsrap CS stat are: " V_sup; print " the 5%, 10%,15%,20% critical values for tao=1,2 4,12 are : " ; for loop(1,7,1); V_b=sortc(V_sup,loop); v_b[95,loop]~v_b[90,loop]~v_b[85,loop]~v_b[80,loop]; endfor; y2=time; print "time is :" y2-y1; end; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Lden: returns the Kernal density * ************************************************************************************************* * Inputs: arg - the values at which the density is to evaluated * * v - The realised values of the series for which density is to be estimates * * Output: p - density * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc Lden(arg,v); local Narg,r,N,h,P,j,d,t1; Narg = rows(arg); r=1/(4+cols(v)); N=rows(v); h=(1/(N^r))*stdc(v)'; P=zeros(Narg,1); j=0; do while j < Narg;j=j+1; d=(arg[j,.]-V)./h; t1=pdfn(d)./h; t1=prodc(t1'); P[j]=meanc(t1); endo; retp(P); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_CJ: Returns the obj func for CHEN-JUMP Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_CJ(b); local moms,g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,gsubT,gsubTall,NWlags,NW,invNW,flagnw,g11; moms=CJ_moms(b); g1=meanc(xt[2:rows(xt)-3]-moms[1:rows(moms)-4,1]); g2=meanc(xt[2:rows(xt)-3]^2-moms[1:rows(moms)-4,2]); g3=meanc(xt[2:rows(xt)-3]^3-moms[1:rows(moms)-4,3]); g4=meanc(xt[2:rows(xt)-3]^4-moms[1:rows(moms)-4,4]); g5=meanc(xt[2:rows(xt)-3].*xt[3:rows(xt)-2]-moms[1:rows(moms)-4,5]); g6=meanc(xt[2:rows(xt)-3].*xt[4:rows(xt)-1]-moms[1:rows(moms)-4,6]); g7=meanc(xt[2:rows(xt)-3].*xt[5:rows(xt)]-moms[1:rows(moms)-4,7]); g8=meanc((xt[2:rows(xt)-3]^2).*xt[3:rows(xt)-2]-moms[1:rows(moms)-4,8]); g9=meanc((xt[2:rows(xt)-3]^2).*xt[4:rows(xt)-1]-moms[1:rows(moms)-4,9]); g10=meanc(xt[2:rows(xt)-3].*(xt[3:rows(xt)-2]^2)-moms[1:rows(moms)-4,10]); g11=meanc(xt[2:rows(xt)-3].*(xt[4:rows(xt)-1]^2)-moms[1:rows(moms)-4,11]); gsubT=g1|g2|g3|g4|g5|g6|g7|g8|g9|g10|g11; gsubTall=((xt[2:rows(xt)-3]-moms[1:rows(moms)-4,1]) ~ (xt[2:rows(xt)-3]^2-moms[1:rows(moms)-4,2]) ~(xt[2:rows(xt)-3]^3-moms[1:rows(moms)-4,3]) ~ (xt[2:rows(xt)-3]^4-moms[1:rows(moms)-4,4]) ~(xt[2:rows(xt)-3].*xt[3:rows(xt)-2]-moms[1:rows(moms)-4,5]) ~ (xt[2:rows(xt)-3].*xt[4:rows(xt)-1]-moms[1:rows(moms)-4,6]) ~(xt[2:rows(xt)-3].*xt[5:rows(xt)]-moms[1:rows(moms)-4,7]) ~ ((xt[2:rows(xt)-3]^2).*xt[3:rows(xt)-2]-moms[1:rows(moms)-4,8]) ~((xt[2:rows(xt)-3]^2).*xt[4:rows(xt)-1]-moms[1:rows(moms)-4,9]) ~ (xt[2:rows(xt)-3].*(xt[3:rows(xt)-2]^2)-moms[1:rows(moms)-4,10]) ~(xt[2:rows(xt)-3].*(xt[4:rows(xt)-1]^2)-moms[1:rows(moms)-4,11]) ); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*inv(NW)*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc CJ_moms(b); local moms,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11; m1=(exp((-((b[1]+b[10]+b[8]))).*(1/52)).*((exp(b[1]*(1/52))))^((b[10]+b[8])/b[1]).*((b[1]*xt-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*(1/52))))/b[1]; m2=(1/(b[1]^2*b[5]^3)).*((exp((-((2*b[1]+b[5]+b[10]+b[8]))).*(1/52)).*((exp(b[1]*(1/52))))^((b[10]+b[8])/b[1]).*((exp(((2*b[1]+b[5])).*(1/52))*b[5]^3*((((b[11]*b[10]-b[9]*b[8]))^2+b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))))+b[1]^4*b[6]*b[7]^2-2*exp(((b[1]+b[5])).*(1/52))*b[5]^3*((b[11]*b[10]-b[9]*b[8])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(b[5]*(1/52)).*((b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^3*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^3*xt*b[6]*(1/52)+b[1]^2*b[5]^3*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52))))))))))); m3=(1/(2*b[1]^3*b[5]^6)).*((exp((-((3*b[1]+3*b[5]+b[10]+b[8]))).*(1/52)).*((exp(b[1]*(1/52))))^((b[10]+b[8])/b[1]).*(((-2)*exp(3*b[5]*(1/52)).*(((-1)+exp(b[1]*(1/52))))^3*b[5]^6*((b[11]*b[10]-b[9]*b[8]))^3-6*exp(3*b[5]*(1/52)).*(((-1)+exp(b[1]*(1/52))))^2*b[1]*b[5]^6*((b[11]*b[10]-b[9]*b[8])).*((xt*(((-b[11])*b[10]+b[9]*b[8]))+((1+exp(b[1]*(1/52)))).*((b[11]^2*b[10]+b[9]^2*b[8]))))+2*exp(3*b[5]*(1/52))*b[1]^3*b[5]^6*((xt^3+3*(((-1)+exp(2*b[1]*(1/52)))).*((b[11]^2*b[10]+b[9]^2*b[8]))*b[6]*(1/52)+3*xt*((b[3]-2*(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))*b[6])).*(1/52)))+6*exp(2*b[5]*(1/52))*b[1]^5*b[5]^3*xt*b[6]*((exp(b[5]*(1/52))*b[5]^3*b[6]*(1/52)^2+b[7]^2*((1+exp(b[5]*(1/52)).*(((-1)+b[5]*(1/52)))))))+2*exp(2*b[5]*(1/52))*b[1]^6*b[5]*b[6]*((exp(b[5]*(1/52))*b[5]^5*b[6]^2*(1/52)^3+3*b[7]^4*((2+b[5]*(1/52)+exp(b[5]*(1/52)).*(((-2)+b[5]*(1/52)))))+3*b[5]^2*b[6]*b[7]^2*(1/52).*((1+exp(b[5]*(1/52)).*(((-1)+b[5]*(1/52)))))))-3*b[1]^4*b[6]*(((-(((-1)+exp(b[1]*(1/52))))).*(((-1)+exp(b[5]*(1/52))))^2*((1+exp(b[5]*(1/52)))).*((b[5]^2))^(3/2).*((b[11]*b[10]-b[9]*b[8]))*b[7]^2+(((-1)+exp(b[1]*(1/52)))).*(((-1)+exp(3*b[5]*(1/52))))*b[5]^4*((b[11]*b[10]-b[9]*b[8]))*b[7]^2*(1/52)-(((-1)+exp(3*b[5]*(1/52))))*b[5]^6*(1/52).*((xt^2+((b[3]-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))*b[6])).*(1/52)))-((1+exp(3*b[5]*(1/52))))*b[5]^6*(1/52).*((xt^2+((b[3]-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))*b[6])).*(1/52)))+(((-1)+exp(b[1]*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8]))*b[7]^2*((1-exp(b[5]*(1/52))+exp(2*b[5]*(1/52))+b[5]*(1/52)+exp(3*b[5]*(1/52)).*(((-1)+b[5]*(1/52)))))))+2*exp(3*b[5]*(1/52))*b[1]^2*b[5]^6*(((-3).*(((-1)+exp(b[1]*(1/52))))*xt^2*((b[11]*b[10]-b[9]*b[8]))+3*(((-1)+exp(2*b[1]*(1/52))))*xt*((b[11]^2*b[10]+b[9]^2*b[8]))-(((-1)+exp(b[1]*(1/52)))).*((2*((1+exp(b[1]*(1/52))+exp(2*b[1]*(1/52))))*b[11]^3*b[10]-3*(((-1)+exp(b[1]*(1/52))))*b[11]^2*b[10]^2*b[6]*(1/52)+3*b[11]*b[10]*((b[3]+2*(((-1)+exp(b[1]*(1/52))))*b[9]*b[8]*b[6])).*(1/52)-b[9]*b[8]*((2*((1+exp(b[1]*(1/52))+exp(2*b[1]*(1/52))))*b[9]^2+3*b[3]*(1/52)+3*(((-1)+exp(b[1]*(1/52))))*b[9]*b[8]*b[6]*(1/52))))))))))); m4=(1/b[1]^4).*((exp((-((4*b[1]+b[10]+b[8]))).*(1/52)).*((exp(b[1]*(1/52))))^((b[10]+b[8])/b[1]).*((((b[1]*xt-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*(1/52)))^4+(1/b[5]^5).*((4*exp((-b[5]).*(1/52))*b[1]^2*((b[1]*xt-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*(1/52))).*(((-2)*exp(((3*b[1]+b[5])).*(1/52))*b[5]^5*((b[11]^3*b[10]-b[9]^3*b[8]))+3*b[1]^4*b[6]*b[7]^4*((2+b[5]*(1/52)))+exp(b[5]*(1/52)).*((2*b[5]^5*((b[11]^3*b[10]-b[9]^3*b[8]))-6*b[1]^4*b[6]*b[7]^4+3*b[1]^4*b[5]*b[6]*b[7]^4*(1/52)))))))+(1/b[5]^3).*((6*exp((-b[5]).*(1/52))*b[1]*((b[1]*xt-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*(1/52)))^2*((exp(((2*b[1]+b[5])).*(1/52))*b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]))+b[1]^3*b[6]*b[7]^2-exp(b[5]*(1/52)).*((b[1]^3*b[6]*b[7]^2-b[1]^3*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]-b[1]*b[3]*(1/52)))))))))+(1/b[5]^6).*((3*exp((-2)*b[5]*(1/52))*b[1]^2*((exp(((2*b[1]+b[5])).*(1/52))*b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]))+b[1]^3*b[6]*b[7]^2-exp(b[5]*(1/52)).*((b[1]^3*b[6]*b[7]^2-b[1]^3*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]-b[1]*b[3]*(1/52)))))))^2))+(1/(2*b[2]^3*b[5]^7)).*((3*exp((-((b[2]+2*b[5]))).*(1/52))*b[1]^3*((4*exp(((4*b[1]+b[2]+2*b[5])).*(1/52))*b[2]^3*b[5]^7*((b[11]^4*b[10]+b[9]^4*b[8]))+2*exp(2*b[5]*(1/52))*b[1]*b[5]^7*b[3]*b[4]^2+exp(b[2]*(1/52))*b[1]^5*b[2]^3*b[6]*b[7]^6+4*exp(((b[2]+b[5])).*(1/52))*b[1]^5*b[2]^3*b[6]*b[7]^6*((7+5*b[5]*(1/52)+b[5]^2*(1/52)^2))-exp(((b[2]+2*b[5])).*(1/52)).*((2*b[1]*b[5]^7*b[3]*b[4]^2-2*b[1]*b[2]*b[5]^7*b[3]*b[4]^2*(1/52)+b[2]^3*((4*b[5]^7*((b[11]^4*b[10]+b[9]^4*b[8]))+29*b[1]^5*b[6]*b[7]^6-10*b[1]^5*b[5]*b[6]*b[7]^6*(1/52))))))))))))); m5=((-((1/(b[1]^2*b[5]^3)).*((exp((-2)*b[1]*(1/52)-b[5]*(1/52)-3*b[1]*(1/52)-2*b[5]*(1/52)).*(((-exp(2*b[1]*(1/52)+b[5]*(1/52)+3*b[1]*(1/52)+2*b[5]*(1/52)))*b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+exp(b[5]*(((1/52)+2*(1/52)))+b[1]*(((1/52)+3*(1/52))))*b[1]*b[5]^3*(((-b[11]^2)*b[10]+b[1]*(1/52)*b[11]*b[10]*b[6]-b[9]*b[8]*((b[9]+b[1]*(1/52)*b[6]))))-exp(((b[1]+b[5])).*(((1/52)+(1/52))))*b[1]^4*b[6]*b[7]^2+exp(2*b[1]*(((1/52)+(1/52)))+b[5]*(((1/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(((b[1]+b[5])).*(((1/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8]-b[1]^2*(1/52)*b[6])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))-exp(b[1]*(((1/52)+(1/52)))+b[5]*(((1/52)+2*(1/52)))).*((b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^3*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^3*xt*b[6]*(1/52)+b[1]^2*b[5]^3*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))); m6=((-((1/(b[1]^2*b[5]^3)).*((exp((-2)*b[1]*(2/52)-b[5]*(2/52)-3*b[1]*(1/52)-2*b[5]*(1/52)).*(((-exp(2*b[1]*(2/52)+b[5]*(2/52)+3*b[1]*(1/52)+2*b[5]*(1/52)))*b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+exp(b[5]*(((2/52)+2*(1/52)))+b[1]*(((2/52)+3*(1/52))))*b[1]*b[5]^3*(((-b[11]^2)*b[10]+b[1]*(2/52)*b[11]*b[10]*b[6]-b[9]*b[8]*((b[9]+b[1]*(2/52)*b[6]))))-exp(((b[1]+b[5])).*(((2/52)+(1/52))))*b[1]^4*b[6]*b[7]^2+exp(2*b[1]*(((2/52)+(1/52)))+b[5]*(((2/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(((b[1]+b[5])).*(((2/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8]-b[1]^2*(2/52)*b[6])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))-exp(b[1]*(((2/52)+(1/52)))+b[5]*(((2/52)+2*(1/52)))).*((b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^3*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^3*xt*b[6]*(1/52)+b[1]^2*b[5]^3*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))); m7=((-((1/(b[1]^2*b[5]^3)).*((exp((-2)*b[1]*(3/52)-b[5]*(3/52)-3*b[1]*(1/52)-2*b[5]*(1/52)).*(((-exp(2*b[1]*(3/52)+b[5]*(3/52)+3*b[1]*(1/52)+2*b[5]*(1/52)))*b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+exp(b[5]*(((3/52)+2*(1/52)))+b[1]*(((3/52)+3*(1/52))))*b[1]*b[5]^3*(((-b[11]^2)*b[10]+b[1]*(3/52)*b[11]*b[10]*b[6]-b[9]*b[8]*((b[9]+b[1]*(3/52)*b[6]))))-exp(((b[1]+b[5])).*(((3/52)+(1/52))))*b[1]^4*b[6]*b[7]^2+exp(2*b[1]*(((3/52)+(1/52)))+b[5]*(((3/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(((b[1]+b[5])).*(((3/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8]-b[1]^2*(3/52)*b[6])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))-exp(b[1]*(((3/52)+(1/52)))+b[5]*(((3/52)+2*(1/52)))).*((b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^3*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^3*xt*b[6]*(1/52)+b[1]^2*b[5]^3*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))); m8=((-((1/(b[1]^3*b[5]^5)).*((exp((-3)*b[1]*(1/52)-b[5]*(1/52)-5*b[1]*(1/52)-3*b[5]*(1/52)).*((exp(3*b[1]*(1/52)+b[5]*(1/52)+5*b[1]*(1/52)+3*b[5]*(1/52))*b[5]^5*((b[11]*b[10]-b[9]*b[8])).*((((b[11]*b[10]-b[9]*b[8]))^2+b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))))-exp(2*b[1]*(1/52)+b[5]*(1/52)+5*b[1]*(1/52)+3*b[5]*(1/52))*b[1]*b[5]^5*(((-2)*b[11]^3*b[10]^2+2*b[11]^2*b[9]*b[10]*b[8]-2*b[11]*b[9]^2*b[10]*b[8]+2*b[9]^3*b[8]^2+b[1]^2*(1/52).*((b[11]^2*b[10]+b[9]^2*b[8]))*b[6]+b[1]*(((-2)*b[11]^3*b[10]+(1/52)*b[11]^2*b[10]^2*b[6]-2*(1/52)*b[11]*b[9]*b[10]*b[8]*b[6]+b[9]^2*b[8]*((2*b[9]+(1/52)*b[8]*b[6]))))))+exp(3*b[1]*(((1/52)+(1/52)))+b[5]*(((1/52)+2*(1/52))))*b[1]^4*b[5]^2*((b[11]*b[10]-b[9]*b[8]))*b[6]*b[7]^2-exp(2*b[1]*(1/52)+b[5]*(1/52)+3*b[1]*(1/52)+2*b[5]*(1/52))*b[1]^4*b[5]^2*b[6]*(((-2)*b[11]*b[10]+2*b[9]*b[8]+b[1]^2*(1/52)*b[6]))*b[7]^2-2*exp(3*b[1]*(1/52)+b[5]*(1/52)+4*b[1]*(1/52)+3*b[5]*(1/52))*b[5]^5*((b[11]*b[10]-b[9]*b[8]))^2*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(2*b[1]*(((1/52)+2*(1/52)))+b[5]*(((1/52)+3*(1/52))))*b[5]^5*(((-((b[11]*b[10]-b[9]*b[8]))^2)-3*b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))+2*b[1]^2*(1/52).*((b[11]*b[10]-b[9]*b[8]))*b[6])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(3*b[1]*(((1/52)+(1/52)))+b[5]*(((1/52)+3*(1/52))))*b[5]*((b[11]*b[10]-b[9]*b[8])).*((b[5]^4*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^4*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^4*xt*b[6]*(1/52)+b[1]^2*b[5]^4*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[5]*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))+exp(2*b[1]*(1/52)+b[5]*(1/52)+3*b[1]*(1/52)+3*b[5]*(1/52))*b[5]*((2*b[11]*b[10]-2*b[9]*b[8]-b[1]^2*(1/52)*b[6])).*((b[5]^4*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^4*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^4*xt*b[6]*(1/52)+b[1]^2*b[5]^4*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[5]*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))-3*exp(2*b[1]*(((1/52)+(1/52)))+b[5]*(((1/52)+2*(1/52))))*b[1]^4*b[6]*b[7]^2*((b[1]*b[5]^2*xt+b[5]^2*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*((b[5]^2*b[6]*(1/52)+b[7]^2*((2+b[5]*(1/52)))))))-exp(2*b[1]*(((1/52)+(1/52)))+b[5]*(((1/52)+3*(1/52)))).*((b[5]^5*((b[11]*b[10]-b[9]*b[8]))^3*-3*b[1]*b[5]^5*((b[11]*b[10]-b[9]*b[8])).*(((-xt)*b[11]*b[10]+b[11]^2*b[10]+xt*b[9]*b[8]+b[9]^2*b[8]))+b[1]^3*b[5]^5*((xt^3-3*((b[11]^2*b[10]+b[9]^2*b[8]))*b[6]*(1/52)+3*xt*((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+3*b[1]^5*b[5]^2*xt*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))+b[1]^6*b[6]*((b[5]^5*b[6]^2*(1/52)^3+3*b[7]^4*(((-2)+b[5]*(1/52)))+3*b[5]^2*b[6]*b[7]^2*(1/52).*(((-1)+b[5]*(1/52)))))+b[1]^2*b[5]^5*((2*b[11]^3*b[10]+3*xt^2*((b[11]*b[10]-b[9]*b[8]))-3*xt*((b[11]^2*b[10]+b[9]^2*b[8]))+3*b[11]^2*b[10]^2*b[6]*(1/52)+3*b[11]*b[10]*((b[3]-2*b[9]*b[8]*b[6])).*(1/52)+b[9]*b[8]*(((-2)*b[9]^2-3*b[3]*(1/52)+3*b[9]*b[8]*b[6]*(1/52)))))+3*b[1]^4*b[5]^2*b[6]*((b[5]^3*(1/52).*((xt^2+b[3]*(1/52)))+b[11]*b[10]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))-b[9]*b[8]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))))); m9=((-((1/(b[1]^3*b[5]^5)).*((exp((-3)*b[1]*(2/52)-b[5]*(2/52)-5*b[1]*(1/52)-3*b[5]*(1/52)).*((exp(3*b[1]*(2/52)+b[5]*(2/52)+5*b[1]*(1/52)+3*b[5]*(1/52))*b[5]^5*((b[11]*b[10]-b[9]*b[8])).*((((b[11]*b[10]-b[9]*b[8]))^2+b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))))-exp(2*b[1]*(2/52)+b[5]*(2/52)+5*b[1]*(1/52)+3*b[5]*(1/52))*b[1]*b[5]^5*(((-2)*b[11]^3*b[10]^2+2*b[11]^2*b[9]*b[10]*b[8]-2*b[11]*b[9]^2*b[10]*b[8]+2*b[9]^3*b[8]^2+b[1]^2*(2/52).*((b[11]^2*b[10]+b[9]^2*b[8]))*b[6]+b[1]*(((-2)*b[11]^3*b[10]+(2/52)*b[11]^2*b[10]^2*b[6]-2*(2/52)*b[11]*b[9]*b[10]*b[8]*b[6]+b[9]^2*b[8]*((2*b[9]+(2/52)*b[8]*b[6]))))))+exp(3*b[1]*(((2/52)+(1/52)))+b[5]*(((2/52)+2*(1/52))))*b[1]^4*b[5]^2*((b[11]*b[10]-b[9]*b[8]))*b[6]*b[7]^2-exp(2*b[1]*(2/52)+b[5]*(2/52)+3*b[1]*(1/52)+2*b[5]*(1/52))*b[1]^4*b[5]^2*b[6]*(((-2)*b[11]*b[10]+2*b[9]*b[8]+b[1]^2*(2/52)*b[6]))*b[7]^2-2*exp(3*b[1]*(2/52)+b[5]*(2/52)+4*b[1]*(1/52)+3*b[5]*(1/52))*b[5]^5*((b[11]*b[10]-b[9]*b[8]))^2*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(2*b[1]*(((2/52)+2*(1/52)))+b[5]*(((2/52)+3*(1/52))))*b[5]^5*(((-((b[11]*b[10]-b[9]*b[8]))^2)-3*b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))+2*b[1]^2*(2/52).*((b[11]*b[10]-b[9]*b[8]))*b[6])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(3*b[1]*(((2/52)+(1/52)))+b[5]*(((2/52)+3*(1/52))))*b[5]*((b[11]*b[10]-b[9]*b[8])).*((b[5]^4*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^4*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^4*xt*b[6]*(1/52)+b[1]^2*b[5]^4*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[5]*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))+exp(2*b[1]*(2/52)+b[5]*(2/52)+3*b[1]*(1/52)+3*b[5]*(1/52))*b[5]*((2*b[11]*b[10]-2*b[9]*b[8]-b[1]^2*(2/52)*b[6])).*((b[5]^4*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^4*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^4*xt*b[6]*(1/52)+b[1]^2*b[5]^4*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[5]*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))-3*exp(2*b[1]*(((2/52)+(1/52)))+b[5]*(((2/52)+2*(1/52))))*b[1]^4*b[6]*b[7]^2*((b[1]*b[5]^2*xt+b[5]^2*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*((b[5]^2*b[6]*(1/52)+b[7]^2*((2+b[5]*(1/52)))))))-exp(2*b[1]*(((2/52)+(1/52)))+b[5]*(((2/52)+3*(1/52)))).*((b[5]^5*((b[11]*b[10]-b[9]*b[8]))^3*-3*b[1]*b[5]^5*((b[11]*b[10]-b[9]*b[8])).*(((-xt)*b[11]*b[10]+b[11]^2*b[10]+xt*b[9]*b[8]+b[9]^2*b[8]))+b[1]^3*b[5]^5*((xt^3-3*((b[11]^2*b[10]+b[9]^2*b[8]))*b[6]*(1/52)+3*xt*((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+3*b[1]^5*b[5]^2*xt*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))+b[1]^6*b[6]*((b[5]^5*b[6]^2*(1/52)^3+3*b[7]^4*(((-2)+b[5]*(1/52)))+3*b[5]^2*b[6]*b[7]^2*(1/52).*(((-1)+b[5]*(1/52)))))+b[1]^2*b[5]^5*((2*b[11]^3*b[10]+3*xt^2*((b[11]*b[10]-b[9]*b[8]))-3*xt*((b[11]^2*b[10]+b[9]^2*b[8]))+3*b[11]^2*b[10]^2*b[6]*(1/52)+3*b[11]*b[10]*((b[3]-2*b[9]*b[8]*b[6])).*(1/52)+b[9]*b[8]*(((-2)*b[9]^2-3*b[3]*(1/52)+3*b[9]*b[8]*b[6]*(1/52)))))+3*b[1]^4*b[5]^2*b[6]*((b[5]^3*(1/52).*((xt^2+b[3]*(1/52)))+b[11]*b[10]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))-b[9]*b[8]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))))); m10=((1/(2*b[1]^2*b[5]^5)).*((2*exp((-2)*b[1]*(1/52)-3*b[1]*(1/52)-b[5]*(1/52))*b[1]*(((-2)*exp(((3*b[1]+b[5])).*(1/52))*b[5]^5*((b[11]^3*b[10]-b[9]^3*b[8]))+3*b[1]^4*b[6]*b[7]^4*((2+b[5]*(1/52)))+exp(b[5]*(1/52)).*((2*b[5]^5*((b[11]^3*b[10]-b[9]^3*b[8]))-6*b[1]^4*b[6]*b[7]^4+3*b[1]^4*b[5]*b[6]*b[7]^4*(1/52)))))+4*exp((-2)*b[1]*(1/52)-3*b[1]*(1/52)-b[5]*(1/52))*b[5]^2*((b[1]*xt-(((-1)+exp(b[1]*(((1/52)+(1/52)))))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(1/52)+(1/52))))).*((exp(((2*b[1]+b[5])).*(1/52))*b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]))+b[1]^3*b[6]*b[7]^2-exp(b[5]*(1/52)).*((b[1]^3*b[6]*b[7]^2-b[1]^3*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]-b[1]*b[3]*(1/52)))))))+(1/b[1]).*((2*exp((-2).*((b[5]*(((1/52)+(1/52)))+b[1]*((3*(1/52)+2*(1/52))))))*b[5]^2*((b[1]*xt-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*(1/52))).*((exp(2*b[5]*(((1/52)+(1/52)))+3*b[1]*((2*(1/52)+(1/52))))*b[5]^3*((((b[11]*b[10]-b[9]*b[8]))^2+b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))))-2*exp(2*b[5]*(((1/52)+(1/52)))+b[1]*((5*(1/52)+3*(1/52))))*b[1]^2*b[5]^3*(1/52).*((b[11]*b[10]-b[9]*b[8]))*b[6]+exp(4*b[1]*(1/52)+2*b[5]*(1/52)+b[1]*(1/52)+b[5]*(1/52))*b[1]^4*b[6]*b[7]^2+exp(4*b[1]*(1/52)+b[5]*(1/52)+3*b[1]*(1/52)+2*b[5]*(1/52))*b[1]^4*b[6]*b[7]^2+exp(2*b[5]*(((1/52)+(1/52)))+b[1]*((4*(1/52)+3*(1/52))))*b[1]^2*((b[5]^3*(1/52).*((b[3]+b[1]^2*(1/52)*b[6]^2))-b[1]^2*b[6]*b[7]^2+b[1]^2*b[5]*(1/52)*b[6]*b[7]^2))-2*exp(2*b[5]*(((1/52)+(1/52)))+b[1]*((5*(1/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+2*exp(2*b[5]*(((1/52)+(1/52)))+2*b[1]*((2*(1/52)+(1/52))))*b[1]^2*b[5]^3*(1/52)*b[6]*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(2*b[5]*(((1/52)+(1/52)))+b[1]*((4*(1/52)+(1/52)))).*((b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^3*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^3*xt*b[6]*(1/52)+b[1]^2*b[5]^3*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))); m11=((1/(2*b[1]^2*b[5]^5)).*((2*exp((-2)*b[1]*(2/52)-3*b[1]*(1/52)-b[5]*(1/52))*b[1]*(((-2)*exp(((3*b[1]+b[5])).*(1/52))*b[5]^5*((b[11]^3*b[10]-b[9]^3*b[8]))+3*b[1]^4*b[6]*b[7]^4*((2+b[5]*(1/52)))+exp(b[5]*(1/52)).*((2*b[5]^5*((b[11]^3*b[10]-b[9]^3*b[8]))-6*b[1]^4*b[6]*b[7]^4+3*b[1]^4*b[5]*b[6]*b[7]^4*(1/52)))))+4*exp((-2)*b[1]*(2/52)-3*b[1]*(1/52)-b[5]*(1/52))*b[5]^2*((b[1]*xt-(((-1)+exp(b[1]*(((2/52)+(1/52)))))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(2/52)+(1/52))))).*((exp(((2*b[1]+b[5])).*(1/52))*b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]))+b[1]^3*b[6]*b[7]^2-exp(b[5]*(1/52)).*((b[1]^3*b[6]*b[7]^2-b[1]^3*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((b[11]^2*b[10]+b[9]^2*b[8]-b[1]*b[3]*(1/52)))))))+(1/b[1]).*((2*exp((-2).*((b[5]*(((2/52)+(1/52)))+b[1]*((3*(2/52)+2*(1/52))))))*b[5]^2*((b[1]*xt-(((-1)+exp(b[1]*(1/52)))).*((b[11]*b[10]-b[9]*b[8]))+b[1]^2*b[6]*(1/52))).*((exp(2*b[5]*(((2/52)+(1/52)))+3*b[1]*((2*(2/52)+(1/52))))*b[5]^3*((((b[11]*b[10]-b[9]*b[8]))^2+b[1]*((b[11]^2*b[10]+b[9]^2*b[8]))))-2*exp(2*b[5]*(((2/52)+(1/52)))+b[1]*((5*(2/52)+3*(1/52))))*b[1]^2*b[5]^3*(2/52).*((b[11]*b[10]-b[9]*b[8]))*b[6]+exp(4*b[1]*(2/52)+2*b[5]*(2/52)+b[1]*(1/52)+b[5]*(1/52))*b[1]^4*b[6]*b[7]^2+exp(4*b[1]*(2/52)+b[5]*(2/52)+3*b[1]*(1/52)+2*b[5]*(1/52))*b[1]^4*b[6]*b[7]^2+exp(2*b[5]*(((2/52)+(1/52)))+b[1]*((4*(2/52)+3*(1/52))))*b[1]^2*((b[5]^3*(2/52).*((b[3]+b[1]^2*(2/52)*b[6]^2))-b[1]^2*b[6]*b[7]^2+b[1]^2*b[5]*(2/52)*b[6]*b[7]^2))-2*exp(2*b[5]*(((2/52)+(1/52)))+b[1]*((5*(2/52)+2*(1/52))))*b[5]^3*((b[11]*b[10]-b[9]*b[8])).*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+2*exp(2*b[5]*(((2/52)+(1/52)))+2*b[1]*((2*(2/52)+(1/52))))*b[1]^2*b[5]^3*(2/52)*b[6]*((b[1]*xt+b[11]*b[10]-b[9]*b[8]+b[1]^2*b[6]*(1/52)))+exp(2*b[5]*(((2/52)+(1/52)))+b[1]*((4*(2/52)+(1/52)))).*((b[5]^3*((b[11]*b[10]-b[9]*b[8]))^2+b[1]*b[5]^3*((2*xt*b[11]*b[10]-b[11]^2*b[10]-2*xt*b[9]*b[8]-b[9]^2*b[8]))+2*b[1]^3*b[5]^3*xt*b[6]*(1/52)+b[1]^2*b[5]^3*((xt^2+((b[3]+2*b[11]*b[10]*b[6]-2*b[9]*b[8]*b[6])).*(1/52)))+b[1]^4*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))))))))))); moms=m1~m2~m3~m4~m5~m6~m7~m8~m9~m10~m11; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_C: Returns the obj func for CHEN Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_C(b); local moms,g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,gsubT,gsubTall,NWlags,NW,invNW,flagnw,g11,g12,g13,g14; moms=C_moms(b); g1=meanc(xt[2:rows(xt)-1]-moms[1:rows(moms)-2,1]); g2=meanc(xt[2:rows(xt)-1]^2-moms[1:rows(moms)-2,2]); g3=meanc(xt[2:rows(xt)-1]^3-moms[1:rows(moms)-2,3]); g4=meanc(xt[2:rows(xt)-1]^4-moms[1:rows(moms)-2,4]); g5=meanc(xt[2:rows(xt)-1].*xt[3:rows(xt)] -moms[1:rows(moms)-2,5]); g6=meanc(xt[2:rows(xt)-1].*(xt[3:rows(xt)]^2) -moms[1:rows(moms)-2,6]); g7=meanc((xt[2:rows(xt)-1]^2).*(xt[3:rows(xt)]) -moms[1:rows(moms)-2,7]); gsubT=g1|g2|g3|g4|g5|g6|g7; gsubTall=((xt[2:rows(xt)-1]-moms[1:rows(moms)-2,1]) ~ (xt[2:rows(xt)-1]^2-moms[1:rows(moms)-2,2]) ~(xt[2:rows(xt)-1]^3-moms[1:rows(moms)-2,3]) ~ (xt[2:rows(xt)-1]^4-moms[1:rows(moms)-2,4]) ~(xt[2:rows(xt)-1].*xt[3:rows(xt)] -moms[1:rows(moms)-2,5]) ~(xt[2:rows(xt)-1].*(xt[3:rows(xt)]^2) -moms[1:rows(moms)-2,6]) ~((xt[2:rows(xt)-1]^2).*(xt[3:rows(xt)]) -moms[1:rows(moms)-2,7]) ); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*inv(NW)*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc C_moms(b); local moms,m1,m2,m3,m4,m5,m6,m7; m1=exp((-b[1]).*(1/52)).*((xt+b[1]*b[6]*(1/52))); m2=(exp((-((2*b[1]+b[5]))).*(1/52)).*((b[1]^2*b[6]*b[7]^2+exp(b[5]*(1/52)).*(((-b[1]^2)*b[6]*b[7]^2+b[1]^2*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52))))))))))/b[5]^3; m3=(1/b[5]^5).*((exp((-((3*b[1]+b[5]))).*(1/52)).*((3*b[1]^2*b[6]*b[7]^2*((2*b[1]*b[7]^2+b[1]*b[5]*b[7]^2*(1/52)+b[5]^2*((xt+b[1]*b[6]*(1/52)))))+exp(b[5]*(1/52)).*(((-6)*b[1]^3*b[6]*b[7]^4+3*b[1]^3*b[5]*b[6]*b[7]^4*(1/52)-3*b[1]^2*b[5]^2*b[6]*b[7]^2*((xt+b[1]*b[6]*(1/52)))+3*b[1]^2*b[5]^3*b[6]*b[7]^2*(1/52).*((xt+b[1]*b[6]*(1/52)))+b[5]^5*((xt+b[1]*b[6]*(1/52))).*((xt^2+3*b[3]*(1/52)+2*b[1]*xt*b[6]*(1/52)+b[1]^2*b[6]^2*(1/52)^2)))))))); m4=(1/(2*b[2]^3*b[5]^7)).*((exp((-((4*b[1]+b[2]+2*b[5]))).*(1/52)).*((6*exp(2*b[5]*(1/52))*b[5]^7*b[3]*b[4]^2+3*exp(b[2]*(1/52))*b[1]^4*b[2]^3*b[6]*b[7]^4*((2*b[5]*b[6]+b[7]^2))+12*exp(((b[2]+b[5])).*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2*((7*b[1]^2*b[7]^4+2*b[1]*b[5]^3*b[7]^2*(1/52).*((xt+b[1]*b[6]*(1/52)))+b[1]^2*b[5]*b[7]^2*(((-b[6])+5*b[7]^2*(1/52)))+b[5]^4*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52)))))+b[1]*b[5]^2*b[7]^2*((4*xt+b[1]*(1/52).*((5*b[6]+b[7]^2*(1/52)))))))+exp(((b[2]+2*b[5])).*(1/52)).*(((-6)*b[5]^7*b[3]*b[4]^2+6*b[2]*b[5]^7*b[3]*b[4]^2*(1/52)+b[2]^3*(((-87)*b[1]^4*b[6]*b[7]^6-12*b[1]^3*b[5]^2*b[6]*b[7]^4*((4*xt+5*b[1]*b[6]*(1/52)))+6*b[1]^3*b[5]^3*b[6]*b[7]^4*(1/52).*((4*xt+5*b[1]*b[6]*(1/52)))+6*b[1]^4*b[5]*b[6]*b[7]^4*((b[6]+5*b[7]^2*(1/52)))-12*b[1]^2*b[5]^4*b[6]*b[7]^2*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52)))))+12*b[1]^2*b[5]^5*b[6]*b[7]^2*(1/52).*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52)))))+2*b[5]^7*((xt^4+4*b[1]*xt^3*b[6]*(1/52)+6*xt^2*(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52)))+4*b[1]*xt*b[6]*(1/52)^2*((3*b[3]+b[1]^2*b[6]^2*(1/52)))+(1/52)^2*((3*b[3]^2+6*b[1]^2*b[3]*b[6]^2*(1/52)+b[1]^4*b[6]^4*(1/52)^2)))))))))))); m5=(((1/b[5]^3).*((exp((-b[1]).*(1/52)-2*b[1]*(1/52)-b[5]*(1/52)).*((b[1]^2*b[6]*b[7]^2+exp(((b[1]+b[5])).*(1/52))*b[1]*b[5]^3*(1/52)*b[6]*((xt+b[1]*b[6]*(1/52)))+exp(b[5]*(1/52)).*(((-b[1]^2)*b[6]*b[7]^2+b[1]^2*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52))))))))))))); m6=(((1/b[5]^5).*((exp((-4)*b[1]*(1/52)-b[5]*(1/52)-6*b[1]*(1/52)-3*b[5]*(1/52)).*((exp(3*b[1]*(1/52)+b[5]*(1/52)+4*b[1]*(1/52)+2*b[5]*(1/52))*b[1]^3*b[5]^2*(1/52)*b[6]^2*b[7]^2+3*exp(3*b[1]*(((2/52)))+b[5]*(((1/52)+2*(1/52))))*b[1]^2*b[6]*b[7]^2*((2*b[1]*b[7]^2+b[1]*b[5]*b[7]^2*(1/52)+b[5]^2*((xt+b[1]*b[6]*(1/52)))))+exp(3*b[1]*(1/52)+b[5]*(1/52)+4*b[1]*(1/52)+3*b[5]*(1/52))*b[1]*b[5]^2*(1/52)*b[6]*(((-b[1]^2)*b[6]*b[7]^2+b[1]^2*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52)))))))+exp(3*b[1]*(((2/52)))+b[5]*(((1/52)+3*(1/52)))).*(((-6)*b[1]^3*b[6]*b[7]^4+3*b[1]^3*b[5]*b[6]*b[7]^4*(1/52)-3*b[1]^2*b[5]^2*b[6]*b[7]^2*((xt+b[1]*b[6]*(1/52)))+3*b[1]^2*b[5]^3*b[6]*b[7]^2*(1/52).*((xt+b[1]*b[6]*(1/52)))+b[5]^5*((xt+b[1]*b[6]*(1/52))).*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((3*b[3]+b[1]^2*b[6]^2*(1/52))))))))))))); m7=(((1/b[5]^3).*((2*exp((-4)*b[1]*(1/52)-3*b[1]*(1/52)).*(((3*exp(2*b[1]*(1/52)-b[5]*(1/52))*b[1]^3*b[6]*b[7]^4*((2+b[5]*(1/52)+exp(b[5]*(1/52)).*(((-2)+b[5]*(1/52))))))/(2*b[5]^2)+exp(2*b[1]*(1/52)-b[5]*(1/52)).*((xt+b[1]*b[6]*((exp(b[1]*(1/52)).*(2/52))))).*((exp(b[5]*(1/52))*b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*((1+exp(b[5]*(1/52)).*(((-1)+b[5]*(1/52)))))))+1/2*exp((-2)*b[5]*(((2/52)))).*((xt+b[1]*b[6]*(1/52))).*((exp(2*b[1]*(1/52)+b[5]*((2*(2/52))))*b[1]^2*b[6]*b[7]^2+exp(2*b[1]*(((2/52)))+b[5]*(((1/52)+2*(1/52))))*b[1]^2*b[6]*b[7]^2+exp(2*((b[1]+b[5])).*(((2/52)))).*((b[5]^3*(1/52).*((b[3]+b[1]^2*(1/52)*b[6]^2))-b[1]^2*b[6]*b[7]^2+b[1]^2*b[5]*(1/52)*b[6]*b[7]^2))+2*exp(2*b[5]*(((2/52)))+b[1]*((2*(2/52))))*b[1]*b[5]^3*(1/52)*b[6]*((xt+b[1]*b[6]*(1/52)))+exp(2*b[1]*(1/52)+2*b[5]*(((2/52)))).*(((-b[1]^2)*b[6]*b[7]^2+b[1]^2*b[5]*b[6]*b[7]^2*(1/52)+b[5]^3*((xt^2+2*b[1]*xt*b[6]*(1/52)+(1/52).*((b[3]+b[1]^2*b[6]^2*(1/52))))))))))))))); moms=m1~m2~m3~m4~m5~m6~m7; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_Feed: Returns the obj func for Feed Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_feed(b); local moms,g1,g2,g3,g4,g5,g6,g7,g8,g9,gsubT,gsubTall,NWlags,NW,invNW,flagnw; moms=feed_moms(b); g1=meanc(xt[2:rows(xt)-2]-moms[1:rows(moms)-3,1]); g2=meanc(xt[2:rows(xt)-2]^2-moms[1:rows(moms)-3,2]); g3=meanc(xt[2:rows(xt)-2]^3-moms[1:rows(moms)-3,3]); g4=meanc(xt[2:rows(xt)-2]^4-moms[1:rows(moms)-3,4]); g5=meanc(xt[2:rows(xt)-2].*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,5]); g6=meanc(xt[2:rows(xt)-2].*(xt[3:rows(xt)-1]^2) -moms[1:rows(moms)-3,6]); g7=meanc((xt[2:rows(xt)-2]^2).*(xt[3:rows(xt)-1]) -moms[1:rows(moms)-3,7]); g8=meanc(xt[2:rows(xt)-2].*(xt[4:rows(xt)]^2) -moms[1:rows(moms)-3,8]); g9=meanc((xt[2:rows(xt)-2]^2).*(xt[4:rows(xt)]) -moms[1:rows(moms)-3,9]); gsubT=g1|g2|g3|g4|g5|g6|g7|g8|g9; gsubTall=( (xt[2:rows(xt)-2]-moms[1:rows(moms)-3,1]) ~ (xt[2:rows(xt)-2]^2-moms[1:rows(moms)-3,2]) ~(xt[2:rows(xt)-2]^3-moms[1:rows(moms)-3,3]) ~ (xt[2:rows(xt)-2]^4-moms[1:rows(moms)-3,4]) ~(xt[2:rows(xt)-2].*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,5]) ~ (xt[2:rows(xt)-2].*(xt[3:rows(xt)-1]^2) -moms[1:rows(moms)-3,6]) ~((xt[2:rows(xt)-2]^2).*(xt[3:rows(xt)-1]) -moms[1:rows(moms)-3,7]) ~ (xt[2:rows(xt)-2].*(xt[4:rows(xt)]^2) -moms[1:rows(moms)-3,8]) ~((xt[2:rows(xt)-2]^2).*(xt[4:rows(xt)]) -moms[1:rows(moms)-3,9]) ); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*inv(NW)*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc feed_moms(b); local moms,m1,m2,m3,m4,m5,m6,m7,m8,m9; m1=(exp((-b[1])*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))))/b[1]; m2=(1/(b[1]^2*b[2]^3*b[5]^3))*((exp((-((2*b[1]+b[2]+b[5])))*(1/52))*((exp(((2*b[1]+b[2]+b[5]))*(1/52))*b[8]^2*b[2]^3*b[5]^3*b[3]^2+exp(b[5]*(1/52))*b[1]^2*b[8]^2*b[5]^3*b[3]*b[4]^2+exp(b[2]*(1/52))*b[1]^4*b[2]^3*b[6]*b[7]^2+2*exp(((b[1]+b[2]+b[5]))*(1/52))*b[8]*b[2]^3*b[5]^3*b[3]*(((-b[8])*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))+exp(((b[2]+b[5]))*(1/52))*((b[8]^2*b[2]^3*b[5]^3*b[3]^2+2*b[1]^3*b[2]^3*b[5]^3*b[6]*(1/52)*((xt-b[8]*b[3]*(1/52)))+2*b[1]*b[8]*b[2]^3*b[5]^3*b[3]*(((-xt)+b[8]*b[3]*(1/52)))+b[1]^4*b[2]^3*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))+b[1]^2*b[5]^3*(((-b[8]^2)*b[3]*b[4]^2+b[8]^2*b[2]*b[3]*b[4]^2*(1/52)+b[2]^3*((xt^2-2*b[8]*xt*b[3]*(1/52)+b[3]*(1/52)*((1-2*b[8]*b[6]+b[8]^2*b[3]*(1/52)))))))))))))-(2*exp((-((2*b[1]+b[2])))*(1/52))*b[8]*b[3]*b[9]*b[4]*((1+exp(b[2]*(1/52))*(((-1)+b[2]*(1/52))))))/b[2]^2; m3=exp((-3)*b[1]*(1/52))*(((((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))^3/b[1]^3+(1/(2*b[2]^5*b[5]^5))*((3*exp((-((5*b[2]+b[5])))*(1/52))*((exp(b[5]*(1/52))*b[8]*b[2]*(((-b[2])+b[2]))*b[5]^5*b[3]*b[4]^2-exp(((b[2]+b[5]))*(1/52))*b[8]*b[2]*(((-b[2])+b[2]))*b[5]^5*b[3]*b[4]^2+2*exp(5*b[2]*(1/52))*b[1]^3*b[2]^5*b[6]*b[7]^4*((2+b[5]*(1/52)))+exp(((5*b[2]+b[5]))*(1/52))*(((-2)*b[8]^3*b[5]^5*b[3]*b[4]^4*(((-2)+b[2]*(1/52)))+b[8]*b[2]*b[5]^5*b[3]*b[4]^2*((b[2]+b[2]-2*b[2]^2*(1/52)))+2*b[1]^3*b[2]^5*b[6]*b[7]^4*(((-2)+b[5]*(1/52)))))-exp(((4*b[2]+b[5]))*(1/52))*b[8]*b[5]^5*b[3]*b[4]^2*((b[2]^2+4*b[8]^2*b[4]^2+b[2]*((b[2]+2*b[8]^2*b[4]^2*(1/52)))))))))+(1/(b[1]*b[2]^3*b[5]^3))*((3*exp((-((b[2]+b[5])))*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))*((exp(b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*b[4]^2+exp(b[2]*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(((b[2]+b[5]))*(1/52))*((b[8]^2*b[5]^3*b[3]*b[4]^2*(((-1)+b[2]*(1/52)))+b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52)))))))))))))+(1/(b[1]*b[2]^4))*((3*exp((-3)*b[1]*(1/52)-2*b[2]*(1/52))*b[3]*b[9]*b[4]*(((-2)*exp(((b[1]+b[2]))*(1/52))*b[8]^2*b[2]^2*b[3]-2*exp(((b[1]+2*b[2]))*(1/52))*b[8]^2*b[2]^2*b[3]*(((-1)+b[2]*(1/52)))+exp(2*b[2]*(1/52))*((2*b[8]^2*b[2]^2*b[3]*(((-1)+b[2]*(1/52)))-2*b[1]^2*b[8]*b[2]^2*b[6]*(1/52)*(((-1)+b[2]*(1/52)))+b[1]*(((-6)*b[8]^2*b[4]^2+b[2]^3*(1/52)*((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)))+b[8]*b[2]*b[4]*((4*b[9]+3*b[8]*b[4]*(1/52)))-b[2]^2*((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)+2*b[8]*b[9]*b[4]*(1/52)))))))+exp(b[2]*(1/52))*((2*b[8]^2*b[2]^2*b[3]-2*b[1]^2*b[8]*b[2]^2*b[6]*(1/52)+b[1]*((6*b[8]^2*b[4]^2+b[8]*b[2]*b[4]*(((-4)*b[9]+3*b[8]*b[4]*(1/52)))+b[2]^2*((1+2*b[8]^2*b[3]*(1/52)-2*b[8]*((xt+b[9]*b[4]*(1/52))))))))))))); m4=exp((-4)*b[1]*(1/52))*(((((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52) +b[1]*((xt-b[8]*b[3]*(1/52)))))^4/b[1]^4+(1/(b[1]*b[2]^5*b[5]^5))*((6*exp((-((5*b[2]+b[5])))*(1/52)) *(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))*((exp(b[5]*(1/52)) *b[8]*b[2]*(((-b[2])+b[2]))*b[5]^5*b[3]*b[4]^2-exp(((b[2]+b[5]))*(1/52))*b[8]*b[2]*(((-b[2])+b[2]))*b[5]^5*b[3] *b[4]^2+2*exp(5*b[2]*(1/52))*b[1]^3*b[2]^5*b[6]*b[7]^4*((2+b[5]*(1/52)))+exp(((5*b[2]+b[5]))*(1/52))*(((-2) *b[8]^3*b[5]^5*b[3]*b[4]^4*(((-2)+b[2]*(1/52)))+b[8]*b[2]*b[5]^5*b[3]*b[4]^2*((b[2]+b[2]-2*b[2]^2*(1/52))) +2*b[1]^3*b[2]^5*b[6]*b[7]^4*(((-2)+b[5]*(1/52)))))-exp(((4*b[2]+b[5]))*(1/52))*b[8]*b[5]^5*b[3]*b[4]^2 *((b[2]^2+4*b[8]^2*b[4]^2+b[2]*((b[2]+2*b[8]^2*b[4]^2*(1/52)))))))))+(1/(2*b[2]^7*b[5]^7))*((3*exp((-2) *((b[2]+b[5]))*(1/52))*((exp(2*b[5]*(1/52))*b[8]^4*b[5]^7*b[3]*b[4]^6+exp(2*b[2]*(1/52))*b[1]^4*b[2]^7*b[6] *b[7]^6+4*exp(((2*b[2]+b[5]))*(1/52))*b[1]^4*b[2]^7*b[6]*b[7]^6*((7+5*b[5]*(1/52)+b[5]^2*(1/52)^2)) +exp(2*((b[2]+b[5]))*(1/52))*(((-2)*b[2]^4*b[5]^7*b[3]*b[4]^2-24*b[8]^2*b[2]^2*b[5]^7*b[3]*b[4]^4 -29*b[8]^4*b[5]^7*b[3]*b[4]^6+2*b[2]^5*b[5]^7*b[3]*b[4]^2*(1/52)+12*b[8]^2*b[2]^3*b[5]^7*b[3]*b[4]^4*(1/52) +10*b[8]^4*b[2]*b[5]^7*b[3]*b[4]^6*(1/52)+b[1]^4*b[2]^7*b[6]*b[7]^6*(((-29)+10*b[5]*(1/52)))))+2*exp(((b[2] +2*b[5]))*(1/52))*b[5]^7*b[3]*b[4]^2*((b[2]^4+14*b[8]^4*b[4]^4+6*b[8]^2*b[2]^3*b[4]^2*(1/52) +10*b[8]^4*b[2]*b[4]^4*(1/52)+2*b[8]^2*b[2]^2*b[4]^2*((6+b[8]^2*b[4]^2*(1/52)^2)))))))) +(1/(b[1]^2*b[2]^3*b[5]^3))*((6*exp((-((b[2]+b[5])))*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3] +b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))^2*((exp(b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*b[4]^2 +exp(b[2]*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(((b[2]+b[5]))*(1/52))*((b[8]^2*b[5]^3*b[3]*b[4]^2*(((-1) +b[2]*(1/52)))+b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))) +(1/(b[2]^6*b[5]^6))*((3*exp((-2)*((b[2]+b[5]))*(1/52))*((exp(b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*b[4]^2 +exp(b[2]*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(((b[2]+b[5]))*(1/52))*((b[8]^2*b[5]^3*b[3]*b[4]^2*(((-1) +b[2]*(1/52)))+b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52)))))))))^2)))) -(1/(b[1]^2*b[2]^6*b[5]^3))*((6*exp((-((4*b[1]+2*b[2]+b[5])))*(1/52))*b[3]*b[9]*b[4]*((2*exp(((2*b[1] +b[2]+b[5]))*(1/52))*b[8]^3*b[2]^4*b[5]^3*b[3]^2+exp(b[5]*(1/52))*b[1]^2*b[8]^2*b[5]^3*b[4]*(((-b[2])*b[9] +b[8]*b[4]))*((2*b[2]*b[3]+b[4]^2))+2*exp(b[2]*(1/52))*b[1]^4*b[8]*b[2]^4*b[6]*b[7]^2+2*exp(((2*b[1]+2*b[2] +b[5]))*(1/52))*b[8]^3*b[2]^4*b[5]^3*b[3]^2*(((-1)+b[2]*(1/52)))+2*exp(2*b[2]*(1/52))*b[1]^4*b[8]*b[2]^4*b[6] *b[7]^2*(((-1)+b[2]*(1/52)))+2*exp(((b[1]+2*b[2]+b[5]))*(1/52))*b[8]*b[2]^2*b[5]^3*b[3]*(((-2)*b[8]^2*b[2]^2 *b[3]*(((-1)+b[2]*(1/52)))+2*b[1]^2*b[8]*b[2]^2*b[6]*(1/52)*(((-1)+b[2]*(1/52)))+b[1]*((6*b[8]^2*b[4]^2-b[2]^3 *(1/52)*((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)))-b[8]*b[2]*b[4]*((4*b[9]+3*b[8]*b[4]*(1/52)))+b[2]^2*((1-2*b[8]*xt +2*b[8]^2*b[3]*(1/52)+2*b[8]*b[9]*b[4]*(1/52)))))))-2*exp(((b[1]+b[2]+b[5]))*(1/52))*b[8]*b[2]^2*b[5]^3*b[3] *((2*b[8]^2*b[2]^2*b[3]-2*b[1]^2*b[8]*b[2]^2*b[6]*(1/52)+b[1]*((6*b[8]^2*b[4]^2+b[8]*b[2]*b[4]*(((-4)*b[9] +3*b[8]*b[4]*(1/52)))+b[2]^2*((1+2*b[8]^2*b[3]*(1/52)-2*b[8]*((xt+b[9]*b[4]*(1/52)))))))))+exp(((2*b[2]+b[5])) *(1/52))*((2*b[8]^3*b[2]^4*b[5]^3*b[3]^2*(((-1)+b[2]*(1/52)))+2*b[1]^4*b[8]*b[2]^4*b[6]*(((-1)+b[2]*(1/52))) *((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))+2*b[1]*b[8]*b[2]^2*b[5]^3*b[3]*(((-6)*b[8]^2*b[4]^2 +b[2]^3*(1/52)*((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)))+b[8]*b[2]*b[4]*((4*b[9]+3*b[8]*b[4]*(1/52)))-b[2]^2 *((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)+2*b[8]*b[9]*b[4]*(1/52)))))-2*b[1]^3*b[2]^2*b[5]^3*b[6]*(1/52)*(((-6) *b[8]^2*b[4]^2+b[2]^3*(1/52)*((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)))+b[8]*b[2]*b[4]*((4*b[9]+3*b[8]*b[4]*(1/52))) -b[2]^2*((1-2*b[8]*xt+2*b[8]^2*b[3]*(1/52)+2*b[8]*b[9]*b[4]*(1/52)))))+b[1]^2*b[5]^3*(((-29)*b[8]^3*b[4]^4 +2*b[2]^5*(1/52)*((b[8]*xt^2+b[8]*b[3]*(1/52)*((2-2*b[8]*b[6]+b[8]^2*b[3]*(1/52)))-xt*((1+2*b[8]^2*b[3] *(1/52)))))+2*b[2]^3*b[4]*((2*b[8]*b[9]^2*b[4]*(1/52)+b[8]*b[4]*(1/52)*((3-3*b[8]*xt+4*b[8]^2*b[3]*(1/52))) +b[9]*((2-4*b[8]*xt+6*b[8]^2*b[3]*(1/52)))))+b[8]^2*b[2]*b[4]^2*((35*b[9]*b[4]+2*b[8]*((b[3]+5*b[4]^2 *(1/52)))))-2*b[8]*b[2]^2*b[4]*((6*((1+b[9]^2))*b[4]+8*b[8]^2*b[3]*b[4]*(1/52)+b[8]*((b[3]*b[9]-6*xt*b[4] +6*b[9]*b[4]^2*(1/52)))))-2*b[2]^4*((b[8]*xt^2-xt*((1+2*b[8]^2*b[3]*(1/52)+2*b[8]*b[9]*b[4]*(1/52))) +(1/52)*((2*b[8]*b[3]+b[9]*b[4]+b[8]^3*b[3]^2*(1/52)+b[8]^2*b[3]*(((-2)*b[6]+3*b[9]*b[4]*(1/52))))))))))) +exp(((b[2]+b[5]))*(1/52))*((2*b[8]^3*b[2]^4*b[5]^3*b[3]^2+2*b[1]^4*b[8]*b[2]^4*b[6]*((b[5]^3*b[6]*(1/52)^2 +b[7]^2*(((-1)+b[5]*(1/52)))))+2*b[1]*b[8]*b[2]^2*b[5]^3*b[3]*((6*b[8]^2*b[4]^2+b[8]*b[2]*b[4]*(((-4)*b[9] +3*b[8]*b[4]*(1/52)))+b[2]^2*((1+2*b[8]^2*b[3]*(1/52)-2*b[8]*((xt+b[9]*b[4]*(1/52)))))))-2*b[1]^3*b[2]^2 *b[5]^3*b[6]*(1/52)*((6*b[8]^2*b[4]^2+b[8]*b[2]*b[4]*(((-4)*b[9]+3*b[8]*b[4]*(1/52)))+b[2]^2*((1+2*b[8]^2 *b[3]*(1/52)-2*b[8]*((xt+b[9]*b[4]*(1/52)))))))+b[1]^2*b[5]^3*((28*b[8]^3*b[4]^4+2*b[8]^2*b[2]*b[4]^2 *(((-17)*b[9]*b[4]-2*b[8]*((b[3]-5*b[4]^2*(1/52)))))+4*b[8]*b[2]^2*b[4]*((3*((1+b[9]^2))*b[4]+b[8]^2*b[4] *(1/52)*((4*b[3]+b[4]^2*(1/52)))+b[8]*((b[3]*b[9]-3*b[4]*((xt+2*b[9]*b[4]*(1/52)))))))+b[2]^3*b[4]*((8*b[8] *b[9]^2*b[4]*(1/52)+6*b[8]*b[4]*(1/52)*((1-b[8]*xt+b[8]^2*b[3]*(1/52)))-b[9]*((4-8*b[8]*xt+b[8]^2*(1/52) *((12*b[3]+5*b[4]^2*(1/52)))))))+2*b[2]^4*((b[8]*xt^2+xt*(((-1)-2*b[8]^2*b[3]*(1/52)+2*b[8]*b[9]*b[4]*(1/52))) +(1/52)*(((-b[9])*b[4]+b[8]^3*b[3]^2*(1/52)-2*b[8]^2*b[3]*((b[6]+b[9]*b[4]*(1/52)))+b[8]*((2*b[3]+b[9]^2*b[4]^2 *(1/52))))))))))))))); m5=(((1/(2*b[1]^2*b[2]^3*b[5]^3)).*((exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*(((1/52)+2*(1/52)))).*((2*exp(2*((b[2]+b[5])).*(1/52)+b[1]*(((1/52)+2*(1/52))))*b[8]^2*b[2]^3*b[5]^3*b[3]^2+2*exp(2*((b[1]+b[2]+b[5])).*(1/52))*b[1]*b[8]*b[2]^3*b[5]^3*(1/52)*b[3]*(((-b[8])*b[3]+b[1]*b[6]))+exp(2*b[5]*(1/52))*b[1]^2*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[1]^2*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^4*b[2]^3*b[6]*b[7]^2+2*exp(2*((b[2]+b[5])).*(1/52)+b[1]*(((1/52)+(1/52))))*b[8]*b[2]^3*b[5]^3*b[3]*(((-b[8])*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))+2*exp(((b[1]+2*((b[2]+b[5])))).*(1/52))*b[2]^3*b[5]^3*((b[8]*((b[3]-b[1]*(1/52)*b[3]))+b[1]^2*(1/52)*b[6])).*(((-b[8])*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))+exp(2*((b[2]+b[5])).*(1/52)).*((2*b[8]^2*b[2]^3*b[5]^3*b[3]^2+4*b[1]^3*b[2]^3*b[5]^3*b[6]*(1/52).*((xt-b[8]*b[3]*(1/52)))+4*b[1]*b[8]*b[2]^3*b[5]^3*b[3]*(((-xt)+b[8]*b[3]*(1/52)))+2*b[1]^4*b[2]^3*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))+b[1]^2*b[5]^3*((b[8]^2*b[3]*((b[4]^2-3*b[4]^2))-4*b[8]*b[2]^2*b[3]*b[9]*b[4]*(1/52)+2*b[8]*b[2]*b[3]*b[4]*((2*b[9]+b[8]*b[4]*(1/52)))+2*b[2]^3*((xt^2-2*b[8]*xt*b[3]*(1/52)+b[3]*(1/52).*((1-2*b[8]*b[6]+b[8]^2*b[3]*(1/52))))))))))))))); m6=(1/4*((((1/(b[1]*b[2]^3*b[5]^3)).*((4*exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*((2*(1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[8]*b[3]+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(1/52)+(1/52)))+b[1]*((xt-b[8]*b[3]*((exp(b[1]*(1/52)).*(1/52)+(1/52))))))).*((exp(2*b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(2*((b[2]+b[5])).*(1/52)).*(((-4)*b[8]*b[2]*b[5]^3*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))+b[8]^2*b[5]^3*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)))))+2*b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))))+((1/(b[2]^5*b[5]^5))*((2*exp((-3).*((b[2]+b[5])).*(1/52)-b[1]*((2*(1/52)+3*(1/52)))).*((exp(3*b[5]*(1/52))*b[8]^3*b[5]^5*b[3]*((b[4]^4-3*b[4]^2*b[4]^2+2*b[4]^4))+6*exp(3*b[2]*(1/52)+2*b[5]*(1/52))*b[1]^3*b[2]^5*b[6]*b[7]^4*((2+b[5]*(1/52)))-3*exp(((b[2]+3*b[5])).*(1/52))*b[8]*b[5]^5*b[3]*((b[4]^2-b[4]^2)).*((b[8]^2*b[4]^2+2*b[8]*b[2]*b[4]*(((-b[9])+b[8]*b[4]*(1/52)))+b[2]^2*((1-2*b[8]*b[9]*b[4]*(1/52)))))+exp(3*((b[2]+b[5])).*(1/52)).*((6*b[2]^3*(((-2)*b[1]^3*b[2]^2*b[6]*b[7]^4+b[1]^3*b[2]^2*b[5]*b[6]*b[7]^4*(1/52)+b[5]^5*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))))-b[8]^3*b[5]^5*b[3]*((b[4]^4+3*b[4]^2*b[4]^2+2*b[4]^4*(((-8)+3*b[2]*(1/52)))))+6*b[8]^2*b[2]*b[5]^5*b[3]*b[9]*b[4]*((b[4]^2+b[4]^2*(((-7)+3*b[2]*(1/52)))))-3*b[8]*b[2]^2*b[5]^5*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)+4*b[9]^2*(((-2)+b[2]*(1/52)))))))))+3*exp(2*b[2]*(1/52)+3*b[5]*(1/52))*b[5]^5*b[3]*((2*b[2]^3*b[9]*b[4]-2*b[8]^2*b[2]*b[9]*b[4]*((b[4]^2-4*b[4]^2)).*((2+b[2]*(1/52)))+b[8]^3*((b[4]^2-2*b[4]^2)).*((b[4]^2+b[4]^2*((3+2*b[2]*(1/52)))))-2*b[8]*b[2]^2*(((-b[4]^2)+2*b[4]^2*((1+b[9]^2*((2+b[2]*(1/52))))))))))))))-((1/b[1])*((2*exp((-b[1]).*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))).*(((-((2*b[8]^2*b[3]^2)/b[1]^2))+(4*exp((-b[1]).*(1/52))*b[8]*(1/52)*b[3]*((b[8]*b[3]-b[1]*b[6])))/b[1]-(exp((-2).*((b[1]+b[2])).*(1/52))*b[8]^2*b[3]*((b[4]^2-b[4]^2)))/b[2]^3-(exp((-2).*((b[2]*(1/52)+b[1]*(((1/52)+(1/52))))))*b[8]^2*b[3]*((b[4]^2-b[4]^2)))/b[2]^3+(2*exp((-((2*b[1]+b[2]))).*(1/52))*b[8]*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2)))))/b[2]^3+(2*exp((-b[2]).*(1/52)-2*b[1]*(((1/52)+(1/52))))*b[8]*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2)))))/b[2]^3-(2*exp((-((2*b[1]+b[5]))).*(1/52))*b[1]^2*b[6]*b[7]^2)/b[5]^3-(2*exp((-b[5]).*(1/52)-2*b[1]*(((1/52)+(1/52))))*b[1]^2*b[6]*b[7]^2)/b[5]^3+((1/(b[2]^3*b[5]^3))*((exp((-2)*b[1]*(1/52)).*((4*b[8]*b[2]^2*b[5]^3*(1/52)*b[3]*b[9]*b[4]-2*b[8]*b[2]*b[5]^3*b[3]*b[4]*((2*b[9]+b[8]*(1/52)*b[4]))-b[8]^2*b[5]^3*b[3]*((b[4]^2-3*b[4]^2))-2*b[2]^3*((b[5]^3*(1/52).*((b[3]+b[8]^2*(1/52)*b[3]^2-2*b[1]*b[8]*(1/52)*b[3]*b[6]+b[1]^2*(1/52)*b[6]^2))-b[1]^2*b[6]*b[7]^2+b[1]^2*b[5]*(1/52)*b[6]*b[7]^2)))))))+(4*exp((-b[1]).*(((1/52)+(1/52))))*b[8]*b[3]*(((-b[1])*xt+b[8]*b[3]+b[1]*b[8]*b[3]*(1/52)-b[1]^2*b[6]*(1/52))))/b[1]^2-(4*exp((-b[1]).*((2*(1/52)+(1/52)))).*(1/52).*(((-b[8])*b[3]+b[1]*b[6])).*(((-b[8])*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))))/b[1]+((1/(b[1]^2*b[2]^3*b[5]^3))*((exp((-2)*b[1]*(((1/52)+(1/52)))).*(((-2)*b[8]^2*b[2]^3*b[5]^3*b[3]^2+4*b[1]*b[8]*b[2]^3*b[5]^3*b[3]*((xt-b[8]*b[3]*(1/52)))+4*b[1]^3*b[2]^3*b[5]^3*b[6]*(1/52).*(((-xt)+b[8]*b[3]*(1/52)))-2*b[1]^4*b[2]^3*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))-b[1]^2*b[5]^3*((b[8]^2*b[3]*((b[4]^2-3*b[4]^2))-4*b[8]*b[2]^2*b[3]*b[9]*b[4]*(1/52)+2*b[8]*b[2]*b[3]*b[4]*((2*b[9]+b[8]*b[4]*(1/52)))+2*b[2]^3*((xt^2-2*b[8]*xt*b[3]*(1/52)+b[3]*(1/52).*((1-2*b[8]*b[6]+b[8]^2*b[3]*(1/52)))))))))))))))))))); m7=(1/2*(((2*exp((-b[1]).*(((1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))^2 .*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[8]*b[3]+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(1/52)+(1/52)))+b[1]*((xt-b[8]*b[3]*((exp(b[1]*(1/52)).*(1/52)+(1/52))))))))/b[1]^3+((1/(b[1]*b[2]^3*b[5]^3))*((exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*(((1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[8]*b[3]+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(1/52)+(1/52)))+b[1]*((xt-b[8]*b[3]*((exp(b[1]*(1/52)).*(1/52)+(1/52))))))).*((exp(2*b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(2*((b[2]+b[5])).*(1/52)).*(((-4)*b[8]*b[2]*b[5]^3*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))+b[8]^2*b[5]^3*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)))))+2*b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))))+((1/(b[1]*b[2]^3*b[5]^3))*((2*exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*(((1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))).*((exp(2*b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^2*((b[2]^2))^(3/2)*b[6]*b[7]^2+exp(2*((b[2]+b[5])).*(1/52)).*(((-4)*b[8]*b[2]*b[5]^3*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))+b[8]^2*b[5]^3*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)))))+2*b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))))+((1/(b[2]^5*b[5]^5))*((exp((-3).*((b[2]+b[5])).*(1/52)-b[1]*((2*(1/52)+7*(1/52)))).*((exp(3*b[5]*(1/52)+b[1]*(((1/52)+4*(1/52))))*b[8]^3*b[5]^5*b[3]*((b[4]^4-3*b[4]^2*b[4]^2+2*b[4]^4))+6*exp(b[1]*(1/52)+4*b[1]*(1/52)+3*b[2]*(1/52)+2*b[5]*(1/52))*b[1]^3*b[2]^5*b[6]*b[7]^4*((2+b[5]*(1/52)))-3*exp(b[1]*(1/52)+4*b[1]*(1/52)+b[2]*(1/52)+3*b[5]*(1/52))*b[8]*b[5]^5*b[3]*((b[4]^2-b[4]^2)).*((b[8]^2*b[4]^2+2*b[8]*b[2]*b[4]*(((-b[9])+b[8]*b[4]*(1/52)))+b[2]^2*((1-2*b[8]*b[9]*b[4]*(1/52)))))+exp(3*((b[2]+b[5])).*(1/52)+b[1]*(((1/52)+4*(1/52)))).*((6*b[2]^3*(((-2)*b[1]^3*b[2]^2*b[6]*b[7]^4+b[1]^3*b[2]^2*b[5]*b[6]*b[7]^4*(1/52)+b[5]^5*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))))-b[8]^3*b[5]^5*b[3]*((b[4]^4+3*b[4]^2*b[4]^2+2*b[4]^4*(((-8)+3*b[2]*(1/52)))))+6*b[8]^2*b[2]*b[5]^5*b[3]*b[9]*b[4]*((b[4]^2+b[4]^2*(((-7)+3*b[2]*(1/52)))))-3*b[8]*b[2]^2*b[5]^5*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)+4*b[9]^2*(((-2)+b[2]*(1/52)))))))))+3*exp(b[1]*(1/52)+4*b[1]*(1/52)+2*b[2]*(1/52)+3*b[5]*(1/52))*b[5]^5*b[3]*((2*b[2]^3*b[9]*b[4]-2*b[8]^2*b[2]*b[9]*b[4]*((b[4]^2-4*b[4]^2)).*((2+b[2]*(1/52)))+b[8]^3*((b[4]^2-2*b[4]^2)).*((b[4]^2+b[4]^2*((3+2*b[2]*(1/52)))))-2*b[8]*b[2]^2*(((-b[4]^2)+2*b[4]^2*((1+b[9]^2*((2+b[2]*(1/52))))))))))))))))); m8=(1/4*((((1/(b[1]*b[2]^3*b[5]^3)).*((4*exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*((2*(2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[8]*b[3]+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(2/52)+(1/52)))+b[1]*((xt-b[8]*b[3]*((exp(b[1]*(1/52)).*(2/52)+(1/52))))))).*((exp(2*b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(2*((b[2]+b[5])).*(1/52)).*(((-4)*b[8]*b[2]*b[5]^3*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))+b[8]^2*b[5]^3*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)))))+2*b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))))+((1/(b[2]^5*b[5]^5))*((2*exp((-3).*((b[2]+b[5])).*(1/52)-b[1]*((2*(2/52)+3*(1/52)))).*((exp(3*b[5]*(1/52))*b[8]^3*b[5]^5*b[3]*((b[4]^4-3*b[4]^2*b[4]^2+2*b[4]^4))+6*exp(3*b[2]*(1/52)+2*b[5]*(1/52))*b[1]^3*b[2]^5*b[6]*b[7]^4*((2+b[5]*(1/52)))-3*exp(((b[2]+3*b[5])).*(1/52))*b[8]*b[5]^5*b[3]*((b[4]^2-b[4]^2)).*((b[8]^2*b[4]^2+2*b[8]*b[2]*b[4]*(((-b[9])+b[8]*b[4]*(1/52)))+b[2]^2*((1-2*b[8]*b[9]*b[4]*(1/52)))))+exp(3*((b[2]+b[5])).*(1/52)).*((6*b[2]^3*(((-2)*b[1]^3*b[2]^2*b[6]*b[7]^4+b[1]^3*b[2]^2*b[5]*b[6]*b[7]^4*(1/52)+b[5]^5*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))))-b[8]^3*b[5]^5*b[3]*((b[4]^4+3*b[4]^2*b[4]^2+2*b[4]^4*(((-8)+3*b[2]*(1/52)))))+6*b[8]^2*b[2]*b[5]^5*b[3]*b[9]*b[4]*((b[4]^2+b[4]^2*(((-7)+3*b[2]*(1/52)))))-3*b[8]*b[2]^2*b[5]^5*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)+4*b[9]^2*(((-2)+b[2]*(1/52)))))))))+3*exp(2*b[2]*(1/52)+3*b[5]*(1/52))*b[5]^5*b[3]*((2*b[2]^3*b[9]*b[4]-2*b[8]^2*b[2]*b[9]*b[4]*((b[4]^2-4*b[4]^2)).*((2+b[2]*(1/52)))+b[8]^3*((b[4]^2-2*b[4]^2)).*((b[4]^2+b[4]^2*((3+2*b[2]*(1/52)))))-2*b[8]*b[2]^2*(((-b[4]^2)+2*b[4]^2*((1+b[9]^2*((2+b[2]*(1/52))))))))))))))-((1/b[1])*((2*exp((-b[1]).*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))).*(((-((2*b[8]^2*b[3]^2)/b[1]^2))+(4*exp((-b[1]).*(2/52))*b[8]*(2/52)*b[3]*((b[8]*b[3]-b[1]*b[6])))/b[1]-(exp((-2).*((b[1]+b[2])).*(2/52))*b[8]^2*b[3]*((b[4]^2-b[4]^2)))/b[2]^3-(exp((-2).*((b[2]*(1/52)+b[1]*(((2/52)+(1/52))))))*b[8]^2*b[3]*((b[4]^2-b[4]^2)))/b[2]^3+(2*exp((-((2*b[1]+b[2]))).*(2/52))*b[8]*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2)))))/b[2]^3+(2*exp((-b[2]).*(1/52)-2*b[1]*(((2/52)+(1/52))))*b[8]*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2)))))/b[2]^3-(2*exp((-((2*b[1]+b[5]))).*(2/52))*b[1]^2*b[6]*b[7]^2)/b[5]^3-(2*exp((-b[5]).*(1/52)-2*b[1]*(((2/52)+(1/52))))*b[1]^2*b[6]*b[7]^2)/b[5]^3+((1/(b[2]^3*b[5]^3))*((exp((-2)*b[1]*(2/52)).*((4*b[8]*b[2]^2*b[5]^3*(2/52)*b[3]*b[9]*b[4]-2*b[8]*b[2]*b[5]^3*b[3]*b[4]*((2*b[9]+b[8]*(2/52)*b[4]))-b[8]^2*b[5]^3*b[3]*((b[4]^2-3*b[4]^2))-2*b[2]^3*((b[5]^3*(2/52).*((b[3]+b[8]^2*(2/52)*b[3]^2-2*b[1]*b[8]*(2/52)*b[3]*b[6]+b[1]^2*(2/52)*b[6]^2))-b[1]^2*b[6]*b[7]^2+b[1]^2*b[5]*(2/52)*b[6]*b[7]^2)))))))+(4*exp((-b[1]).*(((2/52)+(1/52))))*b[8]*b[3]*(((-b[1])*xt+b[8]*b[3]+b[1]*b[8]*b[3]*(1/52)-b[1]^2*b[6]*(1/52))))/b[1]^2-(4*exp((-b[1]).*((2*(2/52)+(1/52)))).*(2/52).*(((-b[8])*b[3]+b[1]*b[6])).*(((-b[8])*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))))/b[1]+((1/(b[1]^2*b[2]^3*b[5]^3))*((exp((-2)*b[1]*(((2/52)+(1/52)))).*(((-2)*b[8]^2*b[2]^3*b[5]^3*b[3]^2+4*b[1]*b[8]*b[2]^3*b[5]^3*b[3]*((xt-b[8]*b[3]*(1/52)))+4*b[1]^3*b[2]^3*b[5]^3*b[6]*(1/52).*(((-xt)+b[8]*b[3]*(1/52)))-2*b[1]^4*b[2]^3*b[6]*((b[5]^3*b[6]*(1/52)^2+b[7]^2*(((-1)+b[5]*(1/52)))))-b[1]^2*b[5]^3*((b[8]^2*b[3]*((b[4]^2-3*b[4]^2))-4*b[8]*b[2]^2*b[3]*b[9]*b[4]*(1/52)+2*b[8]*b[2]*b[3]*b[4]*((2*b[9]+b[8]*b[4]*(1/52)))+2*b[2]^3*((xt^2-2*b[8]*xt*b[3]*(1/52)+b[3]*(1/52).*((1-2*b[8]*b[6]+b[8]^2*b[3]*(1/52)))))))))))))))))))); m9=(1/2*(((2*exp((-b[1]).*(((2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52)))))^2 .*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[8]*b[3]+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(2/52)+(1/52)))+b[1]*((xt-b[8]*b[3]*((exp(b[1]*(1/52)).*(2/52)+(1/52))))))))/b[1]^3+((1/(b[1]*b[2]^3*b[5]^3))*((exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*(((2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[8]*b[3]+b[1]^2*b[6]*((exp(b[1]*(1/52)).*(2/52)+(1/52)))+b[1]*((xt-b[8]*b[3]*((exp(b[1]*(1/52)).*(2/52)+(1/52))))))).*((exp(2*b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^2*b[2]^3*b[6]*b[7]^2+exp(2*((b[2]+b[5])).*(1/52)).*(((-4)*b[8]*b[2]*b[5]^3*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))+b[8]^2*b[5]^3*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)))))+2*b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))))+((1/(b[1]*b[2]^3*b[5]^3))*((2*exp((-2).*((b[2]+b[5])).*(1/52)-b[1]*(((2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[8]*b[3]+b[1]^2*b[6]*(1/52)+b[1]*((xt-b[8]*b[3]*(1/52))))).*((exp(2*b[5]*(1/52))*b[8]^2*b[5]^3*b[3]*((b[4]^2-b[4]^2))-2*exp(((b[2]+2*b[5])).*(1/52))*b[8]*b[5]^3*b[3]*((2*b[2]*b[9]*b[4]+b[8]*((b[4]^2-2*b[4]^2))))+2*exp(((2*b[2]+b[5])).*(1/52))*b[1]^2*((b[2]^2))^(3/2)*b[6]*b[7]^2+exp(2*((b[2]+b[5])).*(1/52)).*(((-4)*b[8]*b[2]*b[5]^3*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))+b[8]^2*b[5]^3*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)))))+2*b[2]^3*((b[5]^3*b[3]*(1/52)+b[1]^2*b[6]*b[7]^2*(((-1)+b[5]*(1/52))))))))))))+((1/(b[2]^5*b[5]^5))*((exp((-3).*((b[2]+b[5])).*(1/52)-b[1]*((2*(2/52)+7*(1/52)))).*((exp(3*b[5]*(1/52)+b[1]*(((2/52)+4*(1/52))))*b[8]^3*b[5]^5*b[3]*((b[4]^4-3*b[4]^2*b[4]^2+2*b[4]^4))+6*exp(b[1]*(2/52)+4*b[1]*(1/52)+3*b[2]*(1/52)+2*b[5]*(1/52))*b[1]^3*b[2]^5*b[6]*b[7]^4*((2+b[5]*(1/52)))-3*exp(b[1]*(2/52)+4*b[1]*(1/52)+b[2]*(1/52)+3*b[5]*(1/52))*b[8]*b[5]^5*b[3]*((b[4]^2-b[4]^2)).*((b[8]^2*b[4]^2+2*b[8]*b[2]*b[4]*(((-b[9])+b[8]*b[4]*(1/52)))+b[2]^2*((1-2*b[8]*b[9]*b[4]*(1/52)))))+exp(3*((b[2]+b[5])).*(1/52)+b[1]*(((2/52)+4*(1/52)))).*((6*b[2]^3*(((-2)*b[1]^3*b[2]^2*b[6]*b[7]^4+b[1]^3*b[2]^2*b[5]*b[6]*b[7]^4*(1/52)+b[5]^5*b[3]*b[9]*b[4]*(((-1)+b[2]*(1/52)))))-b[8]^3*b[5]^5*b[3]*((b[4]^4+3*b[4]^2*b[4]^2+2*b[4]^4*(((-8)+3*b[2]*(1/52)))))+6*b[8]^2*b[2]*b[5]^5*b[3]*b[9]*b[4]*((b[4]^2+b[4]^2*(((-7)+3*b[2]*(1/52)))))-3*b[8]*b[2]^2*b[5]^5*b[3]*((b[4]^2+b[4]^2*(((-3)+2*b[2]*(1/52)+4*b[9]^2*(((-2)+b[2]*(1/52)))))))))+3*exp(b[1]*(2/52)+4*b[1]*(1/52)+2*b[2]*(1/52)+3*b[5]*(1/52))*b[5]^5*b[3]*((2*b[2]^3*b[9]*b[4]-2*b[8]^2*b[2]*b[9]*b[4]*((b[4]^2-4*b[4]^2)).*((2+b[2]*(1/52)))+b[8]^3*((b[4]^2-2*b[4]^2)).*((b[4]^2+b[4]^2*((3+2*b[2]*(1/52)))))-2*b[8]*b[2]^2*(((-b[4]^2)+2*b[4]^2*((1+b[9]^2*((2+b[2]*(1/52))))))))))))))))); moms=m1~m2~m3~m4~m5~m6~m7~m8~m9; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_SVJ: Returns the obj func for SV Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_SVJ(b); local moms,g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,gsubT,gsubTall,NWlags,NW,invNW,flagnw,g11,g12,g13,g14; moms=SVJ_moms(b); g1=meanc(xt[2:rows(xt)-2]-moms[1:rows(moms)-3,1]); g2=meanc(xt[2:rows(xt)-2]^2-moms[1:rows(moms)-3,2]); g3=meanc(xt[2:rows(xt)-2]^3-moms[1:rows(moms)-3,3]); g4=meanc(xt[2:rows(xt)-2]^4-moms[1:rows(moms)-3,4]); g5=meanc(xt[2:rows(xt)-2].*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,5]); g6=meanc(xt[2:rows(xt)-2].*xt[4:rows(xt)] -moms[1:rows(moms)-3,6]); g7=meanc(xt[2:rows(xt)-2].*(xt[3:rows(xt)-1]^2) -moms[1:rows(moms)-3,7]); g8=meanc((xt[2:rows(xt)-2]^2).*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,8]); g9=meanc(xt[2:rows(xt)-2].*(xt[4:rows(xt)]^2) -moms[1:rows(moms)-3,9]); g10=meanc((xt[2:rows(xt)-2]^2).*xt[4:rows(xt)] -moms[1:rows(moms)-3,10]); gsubT=g1|g2|g3|g4|g5|g6|g7|g8|g9|g10; gsubTall=((xt[2:rows(xt)-2]-moms[1:rows(moms)-3,1]) ~ (xt[2:rows(xt)-2]^2-moms[1:rows(moms)-3,2]) ~(xt[2:rows(xt)-2]^3-moms[1:rows(moms)-3,3]) ~ (xt[2:rows(xt)-2]^4-moms[1:rows(moms)-3,4]) ~(xt[2:rows(xt)-2].*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,5]) ~ (xt[2:rows(xt)-2].*xt[4:rows(xt)] -moms[1:rows(moms)-3,6]) ~(xt[2:rows(xt)-2].*(xt[3:rows(xt)-1]^2) -moms[1:rows(moms)-3,7]) ~((xt[2:rows(xt)-2]^2).*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,8]) ~(xt[2:rows(xt)-2].*(xt[4:rows(xt)]^2) -moms[1:rows(moms)-3,9])~((xt[2:rows(xt)-2]^2).*xt[4:rows(xt)] -moms[1:rows(moms)-3,10]) ); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*inv(NW)*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc SVJ_moms(b); local moms,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,phi,c,d,Kapar,R_bar,Kapa_v,v_bar,Sigma_v,Rho,iota,xi,r1,r2,Lamdau,NUu,Lamdad,NUd,u; m1=exp((-b[1]).*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]+ xt))-(((((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))/(b[1]*exp(b[1]*(1/52))); m2=(1/b[1]^2).*((exp((-2)*b[1]*(1/52)).*(((((-1)+exp(b[1]*(1/52))))^2*b[2]^2*b[1]^2+xt^2*b[1]^2-b[10]^2*b[1]*b[9]+exp(2*b[1]*(1/52))*b[10]^2*b[1]*b[9]+b[10]^2*b[9]^2-2*exp(b[1]*(1/52))*b[10]^2*b[9]^2+exp(2*b[1]*(1/52))*b[10]^2*b[9]^2-b[8]^2*b[1]*b[7]+exp(2*b[1]*(1/52))*b[8]^2*b[1]*b[7]-2*b[10]*b[8]*b[9]*b[7]+4*exp(b[1]*(1/52))*b[10]*b[8]*b[9]*b[7]-2*exp(2*b[1]*(1/52))*b[10]*b[8]*b[9]*b[7]+b[8]^2*b[7]^2-2*exp(b[1]*(1/52))*b[8]^2*b[7]^2+exp(2*b[1]*(1/52))*b[8]^2*b[7]^2-2*(((-1)+exp(b[1]*(1/52))))*xt*b[1]*((b[10]*b[9]-b[8]*b[7]))+2*(((-1)+exp(b[1]*(1/52))))*b[2]*b[1]*((xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))+b[4]*b[1]^2*(1/52))))); m3=(1/(2*b[1]^3*b[3]^2)).*((exp((-((3*b[1]+3*b[3]))).*(1/52)).*(((-6)*exp(2*b[1]*(1/52)+3*b[3]*(1/52))*b[3]^2*((b[2]*b[1]-xt*b[1]-b[10]*b[9]+b[8]*b[7])).*((b[2]^2*b[1]^2+b[10]^2*b[9]*((b[1]+b[9]))-2*b[10]*b[8]*b[9]*b[7]+b[8]^2*b[7]*((b[1]+b[7]))+b[2]*(((-2)*b[10]*b[1]*b[9]+2*b[8]*b[1]*b[7]))))+2*exp(3*((b[1]+b[3])).*(1/52))*b[3]^2*((b[2]^3*b[1]^3-b[10]^3*b[9]*((2*b[1]^2+3*b[1]*b[9]+b[9]^2))+3*b[10]^2*b[8]*b[9]*((b[1]+b[9]))*b[7]-3*b[10]*b[8]^2*b[9]*b[7]*((b[1]+b[7]))+3*b[2]^2*b[1]^2*(((-b[10])*b[9]+b[8]*b[7]))+b[8]^3*b[7]*((2*b[1]^2+3*b[1]*b[7]+b[7]^2))+3*b[2]*b[1]*((b[10]^2*b[9]*((b[1]+b[9]))-2*b[10]*b[8]*b[9]*b[7]+b[8]^2*b[7]*((b[1]+b[7]))))))+6*exp(2*b[3]*(1/52))*b[4]*b[1]^3*b[6]*b[5]+6*exp(((b[1]+3*b[3])).*(1/52))*b[3]^2*((b[2]*b[1]-b[10]*b[9]+b[8]*b[7])).*((b[2]^2*b[1]^2+xt^2*b[1]^2-b[10]^2*b[1]*b[9]+b[10]^2*b[9]^2-b[8]^2*b[1]*b[7]-2*b[10]*b[8]*b[9]*b[7]+b[8]^2*b[7]^2+2*xt*b[1]*((b[10]*b[9]-b[8]*b[7]))-2*b[2]*b[1]*((xt*b[1]+b[10]*b[9]-b[8]*b[7]))+b[4]*b[1]^2*(1/52)))+2*exp(3*b[3]*(1/52)).*(((-b[2]^3)*b[1]^3*b[3]^2+xt^3*b[1]^3*b[3]^2+2*b[10]^3*b[1]^2*b[3]^2*b[9]-3*b[10]^3*b[1]*b[3]^2*b[9]^2+b[10]^3*b[3]^2*b[9]^3-2*b[8]^3*b[1]^2*b[3]^2*b[7]+3*b[10]^2*b[8]*b[1]*b[3]^2*b[9]*b[7]-3*b[10]*b[8]^2*b[1]*b[3]^2*b[9]*b[7]-3*b[10]^2*b[8]*b[3]^2*b[9]^2*b[7]+3*b[8]^3*b[1]*b[3]^2*b[7]^2+3*b[10]*b[8]^2*b[3]^2*b[9]*b[7]^2-b[8]^3*b[3]^2*b[7]^3+3*xt^2*b[1]^2*b[3]^2*((b[10]*b[9]-b[8]*b[7]))+3*b[2]^2*b[1]^2*b[3]^2*((xt*b[1]+b[10]*b[9]-b[8]*b[7]))-3*b[4]*b[1]^3*b[6]*b[5]+3*b[4]*b[10]*b[1]^2*b[3]^2*b[9]*(1/52)-3*b[4]*b[8]*b[1]^2*b[3]^2*b[7]*(1/52)+3*b[4]*b[1]^3*b[3]*b[6]*b[5]*(1/52)+3*xt*b[1]*b[3]^2*((b[10]^2*b[9]*(((-b[1])+b[9]))-2*b[10]*b[8]*b[9]*b[7]+b[8]^2*b[7]*(((-b[1])+b[7]))+b[4]*b[1]^2*(1/52)))-3*b[2]*b[1]*b[3]^2*((xt^2*b[1]^2+b[10]^2*b[9]*(((-b[1])+b[9]))-b[8]^2*b[1]*b[7]-2*b[10]*b[8]*b[9]*b[7]+b[8]^2*b[7]^2+2*xt*b[1]*((b[10]*b[9]-b[8]*b[7]))+b[4]*b[1]^2*(1/52))))))))); m4=((1/b[1]^4).*((4^(-((b[4]*b[3])/b[5]^2))*exp((-((4*b[1]))).*(1/52)).*((4^((b[4]*b[3])/b[5]^2).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))^4+3*2^(1+(2*b[4]*b[3])/b[5]^2)*b[1]*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))^2*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52)))+3*4^((b[4]*b[3])/b[5]^2)*b[1]^2*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52)))^2+(1/b[3]^2).*((4^(1+(b[4]*b[3])/b[5]^2)*exp((-b[3]).*(1/52))*b[1]^2*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*(((-2)*exp(((3*b[1]+b[3])).*(1/52))*b[3]^2*((b[10]^3*b[9]-b[8]^3*b[7]))+3*b[4]*b[1]*b[6]*b[5]+exp(b[3]*(1/52)).*((2*b[10]^3*b[3]^2*b[9]-2*b[8]^3*b[3]^2*b[7]+3*b[4]*b[1]*b[6]*b[5]*(((-1)+b[3]*(1/52)))))))))+(1/b[3]^3).*((3*4^((b[4]*b[3])/b[5]^2)*exp((-b[3]).*(1/52))*b[1]^3*((2*exp(((4*b[1]+b[3])).*(1/52))*b[3]^3*((b[10]^4*b[9]+b[8]^4*b[7]))+b[4]*b[1]*b[5]^2*((1+4*b[6]^2*((2+b[3]*(1/52)))))-exp(b[3]*(1/52)).*((2*b[10]^4*b[3]^3*b[9]+2*b[8]^4*b[3]^3*b[7]+b[4]*b[1]*b[5]^2*((1-b[3]*(1/52)+b[6]^2*((8-4*b[3]*(1/52)))))))))))))))); m5=(((exp((-b[1]).*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]+ xt))-(((((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))/(b[1]*exp(b[1]*(1/52))))).*((exp((-b[1]).*(((1/52)+(1/52)))).*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[2]+ xt))-(((((-1)+exp(b[1]*(((1/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7]))))/(b[1]*exp(b[1]*(((1/52)+(1/52)))))))-(exp((-b[1]).*(((1/52)+2*(1/52)))).*(((-(((-1)+exp(2*b[1]*(1/52)))))*b[10]^2*b[9]-(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]-b[4]*b[1]*(1/52))))/b[1]); m6=(((exp((-b[1]).*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]+ xt))-(((((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))/(b[1]*exp(b[1]*(1/52))))).*((exp((-b[1]).*(((2/52)+(1/52)))).*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[2]+ xt))-(((((-1)+exp(b[1]*(((2/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7]))))/(b[1]*exp(b[1]*(((2/52)+(1/52)))))))-(exp((-b[1]).*(((2/52)+2*(1/52)))).*(((-(((-1)+exp(2*b[1]*(1/52)))))*b[10]^2*b[9]-(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]-b[4]*b[1]*(1/52))))/b[1]); m7=((-2)*exp((-((b[9]+b[7]))).*(((1/52)+(1/52)))).*((exp(1/52)))^(b[9]+b[7]).*((exp(1/52)))^(b[9]+b[7]).*(((exp((-2)*b[1]*(((1/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-b[1]).*((2*(1/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-2).*(1/52)*b[1]-3*b[1]*(1/52)).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]-(exp((-2)*b[1]*(((1/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*((2*(1/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-2).*(1/52)*b[1]-3*b[1]*(1/52)).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-2).*(1/52)*b[1]-3*b[1]*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((1/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7]))))^2)./(2*b[1]^3)+(3*exp((-2).*(1/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((1-exp(b[3]*(1/52))+b[3]*(1/52))))./(2*b[3]^2)-(exp((-2).*(1/52)*b[1]-3*b[1]*(1/52)).*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((1/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52))))/b[1]^2-(exp((-2).*(1/52)*b[1]-3*b[1]*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*((exp(2*b[1]*(1/52)).*(1/52)*b[4]*b[1]-b[10]^2*b[9]-b[8]^2*b[7]+exp(2*b[1]*(((1/52)+(1/52)))).*((b[10]^2*b[9]+b[8]^2*b[7]))+b[4]*b[1]*(1/52))))./(2*b[1]^2)-(3*exp((-2).*(1/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((2+b[3]*(1/52)+exp(b[3]*(1/52)).*(((-2)+b[3]*(1/52))))))./(2*b[3]^2)-(exp((-2).*(1/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2-(exp((-2).*(1/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*(1/52).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]+(exp((-2).*(1/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((1+b[3]*(1/52))).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2))); m8=((-2)*exp((-((b[9]+b[7]))).*(((1/52)+(1/52)))).*((exp(1/52)))^(b[9]+b[7]).*((exp(1/52)))^(b[9]+b[7]).*(((exp((-b[1]).*(((1/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-b[1]).*(((1/52)+2*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-b[1]).*(((1/52)+3*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]-(exp((-b[1]).*(((1/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*(((1/52)+2*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*(((1/52)+3*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*(((1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))^2 .*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((1/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7])))))./(2*b[1]^3)+(3*exp((-b[3]).*(1/52)-b[1]*(((1/52)+3*(1/52))))*b[4]*b[6]*b[5]*((1-exp(b[3]*(1/52))+b[3]*(1/52))))./(2*b[3]^2)-(exp((-b[1]).*(((1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52))))/b[1]^2-(exp((-b[1]).*(((1/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((1/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52))))./(2*b[1]^2)-(3*exp((-(1/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*((2+b[3]*(1/52)+exp(b[3]*(1/52)).*(((-2)+b[3]*(1/52))))))./(2*b[3]^2)-(exp((-(1/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2-(exp((-(1/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*(1/52).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]+(exp((-(1/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*((1+b[3]*(1/52))).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2))); m9=((-2)*exp((-((b[9]+b[7]))).*(((2/52)+(1/52)))).*((exp(2/52)))^(b[9]+b[7]).*((exp(1/52)))^(b[9]+b[7]).*(((exp((-2)*b[1]*(((2/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-b[1]).*((2*(2/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-2).*(2/52)*b[1]-3*b[1]*(1/52)).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]-(exp((-2)*b[1]*(((2/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*((2*(2/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-2).*(2/52)*b[1]-3*b[1]*(1/52)).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-2).*(2/52)*b[1]-3*b[1]*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((2/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7]))))^2)./(2*b[1]^3)+(3*exp((-2).*(2/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((1-exp(b[3]*(1/52))+b[3]*(1/52))))./(2*b[3]^2)-(exp((-2).*(2/52)*b[1]-3*b[1]*(1/52)).*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((2/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52))))/b[1]^2-(exp((-2).*(2/52)*b[1]-3*b[1]*(1/52)).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*((exp(2*b[1]*(1/52)).*(2/52)*b[4]*b[1]-b[10]^2*b[9]-b[8]^2*b[7]+exp(2*b[1]*(((2/52)+(1/52)))).*((b[10]^2*b[9]+b[8]^2*b[7]))+b[4]*b[1]*(1/52))))./(2*b[1]^2)-(3*exp((-2).*(2/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((2+b[3]*(1/52)+exp(b[3]*(1/52)).*(((-2)+b[3]*(1/52))))))./(2*b[3]^2)-(exp((-2).*(2/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2-(exp((-2).*(2/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*(1/52).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]+(exp((-2).*(2/52)*b[1]-((b[3]+3*b[1])).*(1/52))*b[4]*b[6]*b[5]*((1+b[3]*(1/52))).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2))); m10=((-2)*exp((-((b[9]+b[7]))).*(((2/52)+(1/52)))).*((exp(2/52)))^(b[9]+b[7]).*((exp(1/52)))^(b[9]+b[7]).*(((exp((-b[1]).*(((2/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-b[1]).*(((2/52)+2*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]+(exp((-b[1]).*(((2/52)+3*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[10]^3*b[9])/b[1]-(exp((-b[1]).*(((2/52)+(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*(((2/52)+2*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*(((2/52)+3*(1/52)))).*(((-1)+exp(b[1]*(1/52))))*b[8]^3*b[7])/b[1]-(exp((-b[1]).*(((2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7]))))^2 .*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((2/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7])))))./(2*b[1]^3)+(3*exp((-b[3]).*(1/52)-b[1]*(((2/52)+3*(1/52))))*b[4]*b[6]*b[5]*((1-exp(b[3]*(1/52))+b[3]*(1/52))))./(2*b[3]^2)-(exp((-b[1]).*(((2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(1/52))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(1/52)))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52))))/b[1]^2-(exp((-b[1]).*(((2/52)+3*(1/52)))).*(((((-1)+exp(b[1]*(((2/52)+(1/52))))))*b[2]*b[1]+xt*b[1]-(((-1)+exp(b[1]*(((2/52)+(1/52)))))).*((b[10]*b[9]-b[8]*b[7])))).*(((((-1)+exp(2*b[1]*(1/52))))*b[10]^2*b[9]+(((-1)+exp(2*b[1]*(1/52))))*b[8]^2*b[7]+b[4]*b[1]*(1/52))))./(2*b[1]^2)-(3*exp((-(2/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*((2+b[3]*(1/52)+exp(b[3]*(1/52)).*(((-2)+b[3]*(1/52))))))./(2*b[3]^2)-(exp((-(2/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2-(exp((-(2/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*(1/52).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]+(exp((-(2/52))*b[1]-b[3]*(1/52)-3*b[1]*(1/52))*b[4]*b[6]*b[5]*((1+b[3]*(1/52))).*((1+exp(b[3]*(1/52)).*(((-1)+b[3]*(1/52))))))/b[3]^2))); moms=m1~m2~m3~m4~m5~m6~m7~m8~m9~m10; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_SV: Returns the obj func for SV Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_SV(b); local moms,g1,g2,g3,g4,g5,g6,gsubT,gsubTall,NWlags,NW,invNW,flagnw; moms=SV_moms(b); g1=meanc(xt[2:rows(xt)-2]-moms[1:rows(moms)-3,1]); g2=meanc(xt[2:rows(xt)-2]^2-moms[1:rows(moms)-3,2]); g3=meanc(xt[2:rows(xt)-2]^3-moms[1:rows(moms)-3,3]); g4=meanc(xt[2:rows(xt)-2]^4-moms[1:rows(moms)-3,4]); g5=meanc(xt[2:rows(xt)-2].*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,5]); g6=meanc(xt[2:rows(xt)-2].*xt[4:rows(xt)] -moms[1:rows(moms)-3,6]); gsubT=g1|g2|g3|g4|g5|g6; gsubTall=((xt[2:rows(xt)-2]-moms[1:rows(moms)-3,1]) ~ (xt[2:rows(xt)-2]^2-moms[1:rows(moms)-3,2]) ~ (xt[2:rows(xt)-2]^3-moms[1:rows(moms)-3,3]) ~ (xt[2:rows(xt)-2]^4-moms[1:rows(moms)-3,4]) ~ (xt[2:rows(xt)-2].*xt[3:rows(xt)-1] -moms[1:rows(moms)-3,5]) ~ (xt[2:rows(xt)-2].*xt[4:rows(xt)] -moms[1:rows(moms)-3,6]) ); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*inv(NW)*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc SV_moms(b); local moms,m1,m2,m3,m4,m5,m6,phi,c,d,Kapar,R_bar,Kapa_v,v_bar,Sigma_v,Rho,iota,xi,r1,r2,u; m1=exp((-b[1])*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[2]+xt)); m2=exp((-2)*b[1]*(1/52))*(((((-1)+exp(b[1]*(1/52))))^2*b[2]^2+2*(((-1)+exp(b[1]*(1/52))))*b[2]*xt+xt^2+b[4]*(1/52))); m3=(1/b[3]^2)*((exp((-3)*((2*b[1]+b[3]))*(1/52))*((exp(3*((2*b[1]+b[3]))*(1/52))*b[2]^3*b[3]^2-3*exp(5*b[1]*(1/52)+3*b[3]*(1/52))*b[2]^2*((b[2]-xt))*b[3]^2+3*exp(3*b[1]*(1/52)+2*b[3]*(1/52))*b[4]*b[6]*b[5]+3*exp(4*b[1]*(1/52)+3*b[3]*(1/52))*b[2]*b[3]^2*((b[2]^2-2*b[2]*xt+xt^2+b[4]*(1/52)))+exp(3*((b[1]+b[3]))*(1/52))*(((-b[2]^3)*b[3]^2+3*b[2]^2*xt*b[3]^2+xt^3*b[3]^2+3*xt*b[4]*b[3]^2*(1/52)-3*b[2]*b[3]^2*((xt^2+b[4]*(1/52)))+3*b[4]*b[6]*b[5]*(((-1)+b[3]*(1/52))))))))); m4=(1/b[3]^3)*((exp((-4)*((2*b[1]+b[3]))*(1/52))*((exp(4*((2*b[1]+b[3]))*(1/52))*b[2]^4*b[3]^3-4*exp(7*b[1]*(1/52)+4*b[3]*(1/52))*b[2]^3*((b[2]-xt))*b[3]^3+12*exp(5*b[1]*(1/52)+3*b[3]*(1/52))*b[2]*b[4]*b[3]*b[6]*b[5]+6*exp(6*b[1]*(1/52)+4*b[3]*(1/52))*b[2]^2*b[3]^3*((b[2]^2-2*b[2]*xt+xt^2+b[4]*(1/52)))+3*exp(4*b[1]*(1/52)+3*b[3]*(1/52))*b[4]*b[5]*(((-4)*b[2]*b[3]*b[6]+4*xt*b[3]*b[6]+b[5]+8*b[6]^2*b[5]+4*b[3]*b[6]^2*b[5]*(1/52)))-4*exp(5*b[1]*(1/52)+4*b[3]*(1/52))*b[2]*b[3]*((b[2]^3*b[3]^2-3*b[2]^2*xt*b[3]^2-xt^3*b[3]^2-3*xt*b[4]*b[3]^2*(1/52)+3*b[2]*b[3]^2*((xt^2+b[4]*(1/52)))-3*b[4]*b[6]*b[5]*(((-1)+b[3]*(1/52)))))+exp(4*((b[1]+b[3]))*(1/52))*((b[2]^4*b[3]^3-4*b[2]^3*xt*b[3]^3+xt^4*b[3]^3+6*xt^2*b[4]*b[3]^3*(1/52)+6*b[2]^2*b[3]^3*((xt^2+b[4]*(1/52)))+12*xt*b[4]*b[3]*b[6]*b[5]*(((-1)+b[3]*(1/52)))-4*b[2]*b[3]*((xt^3*b[3]^2+3*xt*b[4]*b[3]^2*(1/52)+3*b[4]*b[6]*b[5]*(((-1)+b[3]*(1/52)))))+3*b[4]*((b[4]*b[3]^3*(1/52)^2+b[5]^2*(((-1)+b[3]*(1/52)+4*b[6]^2*(((-2)+b[3]*(1/52))))))))))))); m5=(((((exp((-b[1])*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[2]+xt)))).*((exp((-b[1])*(((1/52)+(1/52))))*(((((-1)+exp(b[1]*(((1/52)+(1/52))))))*b[2]+xt))))+exp((-b[1])*((2*(1/52)+(1/52))))*b[4]*(1/52)))); m6=(((((exp((-b[1])*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[2]+xt)))).*((exp((-b[1])*(((1/52)+(2/52))))*(((((-1)+exp(b[1]*(((1/52)+(2/52))))))*b[2]+xt))))+exp((-b[1])*((2*(1/52)+(2/52))))*b[4]*(1/52)))); moms=m1~m2~m3~m4~m5~m6; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_SM: Returns the obj func for SM Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_SM(b); local moms,g1,g2,g3,g4,g5,gsubT,gsubTall,NWlags,NW,invNW,flagnw,gTall; moms=SM_moms(b); g1=meanc(xt[2:rows(xt)-1]-moms[1:rows(moms)-2,1]); g2=meanc(xt[2:rows(xt)-1]^2-moms[1:rows(moms)-2,2]); g3=meanc(xt[2:rows(xt)-1]^3-moms[1:rows(moms)-2,3]); g4=meanc(xt[2:rows(xt)-1]^4-moms[1:rows(moms)-2,4]); g5=meanc(xt[2:rows(xt)-1].*xt[3:rows(xt)] -moms[1:rows(moms)-2,5]); gsubT=g1|g2|g3|g4|g5; gsubTall=((xt[2:rows(xt)-1]-moms[1:rows(moms)-2,1]) ~ (xt[2:rows(xt)-1]^2-moms[1:rows(moms)-2,2]) ~ (xt[2:rows(xt)-1]^3-moms[1:rows(moms)-2,3]) ~ (xt[2:rows(xt)-1]^4-moms[1:rows(moms)-2,4]) ~ (xt[2:rows(xt)-1].*xt[3:rows(xt)] -moms[1:rows(moms)-2,5]) ); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*invNW*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc SM_moms(b); local moms,m1,m2,m3,m4,m5; m1=exp((-b[1]).*(1/52)).*((xt+b[1]*b[4].*(1/52))); m2=(1/(2*b[3]^3*b[1])).*((exp((-((b[3]+2*b[1]))).*(1/52)).*((exp(((b[3]+2*b[1])).*(1/52))*b[3]^3*b[2]^2+2*b[1]^3*b[4]*b[5]^2+exp(b[3].*(1/52)).*(((-2)*b[1]^3*b[4]*b[5]^2+2*b[3]*b[1]^3*b[4]*b[5]^2*(1/52)+b[3]^3*((2*xt^2*b[1]-b[2]^2+4*xt*b[1]^2*b[4].*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2)))))))); m3=(1/(2*b[3]^5*b[1])).*((exp((-((b[3]+3*b[1]))).*(1/52)).*((3*exp(((b[3]+2*b[1])).*(1/52))*b[3]^5*b[2]^2*((xt+b[1]*b[4].*(1/52)))+6*b[1]^3*b[4]*b[5]^2*((2*b[1]*b[5]^2+b[3]*b[1]*b[5]^2*(1/52)+b[3]^2*((xt+b[1]*b[4].*(1/52))))) +exp(b[3].*(1/52)).*(((-12)*b[1]^4*b[4]*b[5]^4+6*b[3]*b[1]^4*b[4]*b[5]^4*(1/52)-6*b[3]^2*b[1]^3*b[4]*b[5]^2*((xt +b[1]*b[4].*(1/52)))+6*b[3]^3*b[1]^3*b[4]*b[5]^2*(1/52).*((xt+b[1]*b[4].*(1/52)))+b[3]^5*((xt+b[1]*b[4].*(1/52))).*((2*xt^2*b[1]-3*b[2]^2+4*xt*b[1]^2*b[4].*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2)))))))); m4=(1/(4*b[3]^7*b[1]^2)).*((exp((-2).*((b[3]+2*b[1])).*(1/52)).*((3*exp(2*((b[3]+2*b[1])).*(1/52))*b[3]^7*b[2]^4+12*exp(((b[3]+2*b[1])).*(1/52))*b[3]^4*b[1]^3*b[4]*b[2]^2*b[5]^2+6*b[1]^6*b[4]*b[5]^4*((2*b[3]*b[4]+b[5]^2))+6*exp(2*((b[3]+b[1])).*(1/52))*b[3]^4*b[2]^2*(((-2)*b[1]^3*b[4]*b[5]^2+2*b[3]*b[1]^3*b[4]*b[5]^2*(1/52)+b[3]^3*((2*xt^2*b[1]-b[2]^2+4*xt*b[1]^2*b[4].*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2))))+12*exp(b[3].*(1/52))*b[1]^3*b[4]*b[5]^2*((14*b[1]^3*b[5]^4+4*b[3]^3*b[1]^2*b[5]^2*(1/52).*((xt+b[1]*b[4].*(1/52)))-2*b[3]*b[1]^3*b[5]^2*((b[4]-5*b[5]^2*(1/52)))+b[3]^4*((2*xt^2*b[1]-b[2]^2+4*xt*b[1]^2*b[4].*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2))+2*b[3]^2*b[1]^2*b[5]^2*((4*xt+b[1].*(1/52).*((5*b[4]+b[5]^2*(1/52))))))) +exp(2*b[3].*(1/52)).*(((-174)*b[1]^6*b[4]*b[5]^6-24*b[3]^2*b[1]^5*b[4]*b[5]^4*((4*xt+5*b[1]*b[4].*(1/52))) +12*b[3]^3*b[1]^5*b[4]*b[5]^4*(1/52).*((4*xt+5*b[1]*b[4].*(1/52)))+12*b[3]*b[1]^6*b[4]*b[5]^4*((b[4]+5*b[5]^2*(1/52)))-12*b[3]^4*b[1]^3*b[4]*b[5]^2*((2*xt^2*b[1]-b[2]^2+4*xt*b[1]^2*b[4].*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2))+12*b[3]^5*b[1]^3*b[4]*b[5]^2*(1/52).*((2*xt^2*b[1]-b[2]^2+4*xt*b[1]^2*b[4].*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2)) +b[3]^7*((4*xt^4*b[1]^2+3*b[2]^4+16*xt^3*b[1]^3*b[4].*(1/52)-12*b[1]^3*b[4]^2*b[2]^2*(1/52)^2+4*b[1]^6*b[4]^4*(1/52)^4+12*xt^2*(((-b[1])*b[2]^2+2*b[1]^4*b[4]^2*(1/52)^2))+8*xt*(((-3)*b[1]^2*b[4]*b[2]^2*(1/52)+2*b[1]^5*b[4]^3*(1/52)^3)))))))))); m5=((1/(2*((b[3]^2))^(3/2)*b[1]))*((exp((-((b[3]+b[1])))*(((1/52)+2*(1/52))))*((exp(2*b[1]*(1/52)+b[3]*(((1/52)+2*(1/52))))*b[3]^3*b[2]^2+2*exp(b[3]*(((1/52)+(1/52))))*b[1]^3*b[4]*b[5]^2+2*exp(b[1]*(1/52)+b[3]*(((1/52)+2*(1/52))))*b[3]^3*(1/52)*b[1]^2*b[4]*((xt+b[1]*b[4]*(1/52)))+exp(b[3]*(((1/52)+2*(1/52))))*(((-2)*b[1]^3*b[4]*b[5]^2+2*b[3]*b[1]^3*b[4]*b[5]^2*(1/52)+b[3]^3*((2*xt^2*b[1]-b[2]^2+4*xt*b[1]^2*b[4]*(1/52)+2*b[1]^3*b[4]^2*(1/52)^2))))))))); moms=m1~m2~m3~m4~m5; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * GMM_CIR: Returns the obj func for CIR Model * ************************************************************************************************* * Inputs: b - starting values * * Output: - the objective function to be minimised * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc gmm_CIR(b); local moms,g1,g2,g3,g4,g5,g6,gsubT,gsubTall,NWlags,NW,invNW,flagnw; moms=CIR_moms(b); g1=meanc(xt[2:rows(xt)-1]-moms[1:rows(moms)-2,1]); g2=meanc(xt[2:rows(xt)-1]^2-moms[1:rows(moms)-2,2]); g3=meanc(xt[2:rows(xt)-1].*xt[3:rows(xt)] -moms[1:rows(moms)-2,3]); gsubT=g1|g2|g3; gsubTall=((xt[2:rows(xt)-1]-moms[1:rows(moms)-2,1]) ~ (xt[2:rows(xt)-1]^2-moms[1:rows(moms)-2,2]) ~ (xt[2:rows(xt)-1].*xt[3:rows(xt)] -moms[1:rows(moms)-2,3])); NWlags=int(rows(xt)^(1/6)); NW=nwywest(gsubTall,NWlags); trap 1; invNW=inv(NW); flagnw = scalerr(invNW); if flagnw == 0; retp(gsubT'*inv(NW)*gsubT); else; retp(gsubT'*eye(rows(NW))*gsubT); endif; endp; proc CIR_moms(b); local moms,m1,m2,m12; m1=exp((-b[1])*(1/52))*(((((-1)+exp(b[1]*(1/52))))*b[2]+xt)); m2=(exp((-2)*b[1]*(1/52))*((2*b[1]*(((((-1)+exp(b[1]*(1/52))))*b[2]+xt))^2+(((-1)+exp(b[1]*(1/52))))*(((((-1)+exp(b[1]*(1/52))))*b[2]+2*xt))*b[3]^2)))/(2*b[1]); m12=(exp((-3)*b[1]*(1/52))*((2*b[1]*(((((-1)+exp(b[1]*(1/52))))^2*((1+exp(b[1]*(1/52))))*b[2]^2+(((-2)+exp(b[1]*(1/52))+exp(2*b[1]*(1/52))))*b[2]*xt+xt^2))+(((-1)+exp(b[1]*(1/52))))*(((((-1)+exp(b[1]*(1/52))))*b[2]+2*xt))*b[3]^2)))/(2*b[1]); moms=m1~m2~m12; retp(moms); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * ineq_One: One Factor Models Nonlinear Inequality constraint * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc ineq_one(b); retp(2*b[1]*b[2]-b[3]^2); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * ineq_two: Two Factor Models Nonlinear Inequality constraint * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc ineq_two(b); retp(2*b[3]*b[4]-b[5]^2); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * ineq_three: Three Factor Models Nonlinear Inequality constraint * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc ineq_three(b); retp((2*b[2]*b[3]-b[4]^2)|(2*b[5]*b[6]-b[7]^2)); endp; proc Nwywest(hhat,pp); local N,df,m,Omega0,jj,Shat,Omegaj,titl,flag1; hhat = hhat - meanc(hhat)'; N = rows(hhat); m = minc( pp|(N-2) ); Omega0 = hhat'hhat/N; jj = 1; Shat = Omega0; do until jj > pp; Omegaj = hhat[jj+1:N,.]'hhat[1:N-jj,.]/N; Shat = Shat + ( 1 - jj/(m+1) ) * (Omegaj + Omegaj'); jj = jj + 1; endo; retp(Shat); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * DGP_CIR: Data Generation as CIR process * * dp(t)=phi*(p_bar-p(t))*dt+sig*sqrt(p(t)*dW(t) * ************************************************************************************************* * Inputs: T - Length of time series * * h - discretization interval * * phi - mean reversion parameter Process * * p_bar - mean level * * sig1 - variance term * * see - seed of rndns * * Output: dat - time series * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc DGP_CIR(T,h,phi,p_bar,sig1,see); local TN,dat,Pt,i,X1,ee; TN=round(T/h); ee=sqrt(h)*rndns(TN,1,see); dat={}; Pt=p_bar;// here we use the estimated parameter as start value i=0; do while i < TN; i=i+1; X1=Pt+phi*(p_bar-Pt)*h+sig1*sqrt(pt)*ee[i] -0.5*(sig1*sqrt(pt))*(0.5*sig1/sqrt(pt))*h +0.5*(sig1*sqrt(pt))*(0.5*sig1/sqrt(pt))*(ee[i]^2); /*Theory implies that this condition will never be reached as when process approaches zero the volatility is switched off, that result depends on h going to zero here h is very small indeed I am including the condition as a final safe guard, here I am switching off the volatility manually something that should happens assymptotically*/ if X1 < 0; X1=Pt+phi*(p_bar-Pt)*h; print " same $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"; endif; if i%(h^-1)==0; dat=dat|X1; endif; Pt=X1; endo; retp(dat); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * DGP_CIRsim: Data Generation as CIR process for simulation, The only difference with above is that we are given the start value of X(0) * * dp(t)=phi*(p_bar-p(t))*dt+sig*sqrt(p(t)*dW(t) * *************************************************************************************************/ proc DGP_CIRsim(T,h,phi,p_bar,sig1,see,start); local TN,dat,Pt,i,X1,ee; TN=round(T/h); ee=sqrt(h)*rndns(TN,1,see); dat={}; Pt=start; // here we use the given X(0) as start value i=0; do while i < TN; i=i+1; X1=Pt+phi*(p_bar-Pt)*h+sig1*sqrt(pt)*ee[i] -0.5*(sig1*sqrt(pt))*(0.5*sig1/sqrt(pt))*h +0.5*(sig1*sqrt(pt))*(0.5*sig1/sqrt(pt))*(ee[i]^2); /*Theory implies that this condition will never be reached as when process approaches zero the volatility is switched off, that result depends on h going to zero here h is very small indeed I am including the condition as a final safe guard, here I am switching off the volatility manually something that should happens assymptotically*/ if X1 < 0; X1=Pt+phi*(p_bar-Pt)*h; endif; if i%((h^-1))==0; dat=dat|X1; endif; Pt=X1; endo; retp(dat); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data Generation * * DGP: SM * * dr(t) = kapa_r*(theta(t)-r(t))*dt+sig1*dW1(t) * * dTheta(t) = kapa_Theta*(theta_bar-theta(t))*dt+(sig3*sqrt(theta(t)))*dW3(t) * * W's are mulually independent Brownian motions * ************************************************************************************************* * Inputs: T: Length of time series * * h: discretization interval * * kapa_r: mean reversion parameter Interest Rate Process * * sig1: variance term interest rate Process * * kapa_theta: mean reversion parameter mean process * * theta_bar: mean Level, also used as the starting value of the mean process * * sig3: variance term mean Process * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc(2)= DGP_SM(T,h,kapa_r,sig1,kapa_theta,theta_bar,sig3,see); local TN,e1,e2,dat,Pt,Mt,i,zi1,zi2,r,r2,rhop,I21,M1,XX,test; test={}; dat={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); Pt=theta_bar; Mt=theta_bar; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+sig1*e1[i]; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e2[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; test=test|M1; endif; Pt=XX; Mt=M1; endo; retp(dat,test); endp; proc(2)= DGP_SMsim(T,h,kapa_r,sig1,kapa_theta,theta_bar,sig3,see,start1,start2); local TN,e1,e2,dat,Pt,Mt,i,zi1,zi2,r,r2,rhop,I21,M1,XX,test; test={}; dat={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); Pt=start1; // here we use the given X(0) as start value Mt=start2; // here we use the given seta(0) as start value i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+sig1*e1[i]; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e2[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; test=test|M1; endif; Pt=XX; Mt=M1; endo; retp(dat,test); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data Generation * * DGP:SV * * dr(t) = kapa_r*(r_bar-r(t))*dt+(sqrt(V(t)))*dW1(t) * * dV(t) = kapa_v*(v_bar-V(t))*dt+(sig2*sqrt(V(t)))*dW2(t) * ************************************************************************************************* * Inputs: T: Length of time series * * h: discretization interval * * kapa_r: mean reversion parameter Interest Rate Process * * r_bar: mean Level, also used as the starting value of the mean process * * kapa_v: mean reversion parameter volatility process * * v_bar: mean Volatility Level, also used as the starting value of the vol process * * sig2: variance term Volatility Process * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc(2)= DGP_SV(T,h,kapa_r,r_bar,kapa_v,v_bar,sig2,rho,see); local p,TN,e1,e2,dat,Pt,Vt,i,zi1,zi2,r,r2,rhop,I21,V1,XX,test; test={}; p=1/h; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); dat={}; Pt=r_bar; Vt=v_bar; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(r_bar-Pt)*h+sqrt(Vt)*(e1[i]*((1-rho^2)^0.5)+rho*e2[i]); V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; if i%(h^-1)==0; dat=dat|XX; test=test|v1; endif; Pt=XX; Vt=V1; endo; retp(dat,test); endp; proc(2)= DGP_SVsim(T,h,kapa_r,r_bar,kapa_v,v_bar,sig2,rho,see,start1,start2); local p,TN,e1,e2,dat,Pt,Vt,i,zi1,zi2,r,r2,rhop,I21,V1,XX,test; test={}; p=1/h; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); dat={}; Pt=start1; Vt=start2; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(r_bar-Pt)*h+sqrt(Vt)*(e1[i]*((1-rho^2)^0.5)+rho*e2[i]); V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; if i%(h^-1)==0; dat=dat|XX; test=test|v1; endif; Pt=XX; Vt=V1; endo; retp(dat,test); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data Generation * * DGP:SVj * ************************************************************************************************* * Inputs: T: Length of time series * * h: discretization interval * * kapa_r: mean reversion parameter Interest Rate Process * * r_bar: mean Level, also used as the starting value of the mean process * * kapa_v: mean reversion parameter volatility process * * v_bar: mean Volatility Level, also used as the starting value of the vol process * * sig2: variance term Volatility Process * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc(2)= DGP_SVJ(T,h,kapa_r,r_bar,kapa_v,v_bar,sig2,rho,lamdau,NUu,lamdad,NUd,see); local p,TN,e1,e2,dat,Pt,Vt,i,zi1,zi2,r,r2,rhop,I21,V1,XX,test,e3,jumpu,p_jumpu,jumpd,p_jumpd,jump; test={}; p=1/h; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3 = rndus(TN,4,see); jumpu=-NUu*ln(e3[.,2]); p_jumpu=exp(-lamdau*h); p_jumpu=e3[.,1].>p_jumpu; jumpu=jumpu.*p_jumpu; jumpd=-NUd*ln(e3[.,4]); p_jumpd=exp(-lamdad*h); p_jumpd=e3[.,3].>p_jumpd; jumpd=jumpd.*p_jumpd; jump=jumpu-jumpd; dat={}; Pt=r_bar; Vt=v_bar; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(r_bar-Pt)*h+sqrt(Vt)*(e1[i]*((1-rho^2)^0.5)+rho*e2[i])+jump[i]; V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; /*Theory implies that this condition will never be reached as when process approaches zero the volatility is switched off, that result depends on h going to zero. Here h is very small indeed, I am including the condition as a final safe-guard, here I am switching off the volatility manually as a precaution something that should happens assymptotically*/ if i%(h^-1)==0; dat=dat|XX; test=test|v1; endif; Pt=XX; Vt=V1; endo; retp(dat,test); endp; proc(2)= DGP_SVJsim(T,h,kapa_r,r_bar,kapa_v,v_bar,sig2,rho,lamdau,NUu,lamdad,NUd,see,start1,start2); local p,TN,e1,e2,dat,Pt,Vt,i,zi1,zi2,r,r2,rhop,I21,V1,XX,test,e3,jumpu,p_jumpu,jumpd,p_jumpd,jump; test={}; p=1/h; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3 = rndus(TN,4,see); jumpu=-NUu*ln(e3[.,2]); p_jumpu=exp(-lamdau*h); p_jumpu=e3[.,1].>p_jumpu; jumpu=jumpu.*p_jumpu; jumpd=-NUd*ln(e3[.,4]); p_jumpd=exp(-lamdad*h); p_jumpd=e3[.,3].>p_jumpd; jumpd=jumpd.*p_jumpd; jump=jumpu-jumpd; //calculate Jump; dat={}; Pt=start1; Vt=start2; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(r_bar-Pt)*h+sqrt(Vt)*(e1[i]*((1-rho^2)^0.5)+rho*e2[i])+jump[i]; V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; /*Theory implies that this condition will never be reached as when process approaches zero the volatility is switched off, that result depends on h going to zero. Here h is very small indeed, I am including the condition as a final safe-guard, here I am switching off the volatility manually as a precaution something that should happens assymptotically*/ if i%(h^-1)==0; dat=dat|XX; test=test|v1; endif; Pt=XX; Vt=V1; endo; retp(dat,test); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data Generation * * DGP: CHEN * * dr(t) = kapa_r*(theta(t)-r(t))*dt+(sqrt(V(t)))*dW1(t) * * dV(t) = kapa_v*(v_bar-V(t))*dt+(sig2*sqrt(V(t)))*dW2(t) * * dTheta(t) = kapa_Theta*(theta_bar-theta(t))*dt+(sig3*sqrt(theta(t)))*dW3(t) * * W's are mulually independent Brownian motions * ************************************************************************************************* * Inputs: T: Length of time series * * h: discretization interval * * kapa_r: mean reversion parameter Interest Rate Process * * kapa_v: mean reversion parameter volatility process * * v_bar: mean Volatility Level, also used as the starting value of the vol process * * sig2: variance term Volatility Process * * kapa_theta: mean reversion parameter mean process * * theta_bar: mean Level, also used as the starting value of the mean process * * sig3: variance term mean Process * * sim_ind: indicator for generating the errors * * sim_ind = 0 if errors are to be generated randomly * * sim_ind = 1 if errors are to be be taken from a global variable (SGMM) * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc(3)= DGP_chen(T,h,kapa_r,kapa_v,v_bar,sig2,kapa_theta,theta_bar,sig3,see); local theta_fact,v_fact,TN,e1,e2,e3,dat,pt,vt,mt,i,xx,v1,m1; v_fact={};theta_fact={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3= sqrt(h)*rndns(TN,1,see); dat={}; Pt=theta_bar; Vt=v_bar; Mt=theta_bar; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+sqrt(Vt)*e1[i]; V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e3[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; v_fact=v_fact|v1; theta_fact=theta_fact|M1; endif; Pt=XX; Vt=V1; Mt=m1; endo; retp(dat,v_fact,theta_fact); endp; proc(3)= DGP_chensim(T,h,kapa_r,kapa_v,v_bar,sig2,kapa_theta,theta_bar,sig3,see,start1,start2,start3); local theta_fact,v_fact,TN,e1,e2,e3,dat,pt,vt,mt,i,xx,v1,m1; v_fact={};theta_fact={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3= sqrt(h)*rndns(TN,1,see); dat={}; Pt=start1; Vt=start2; Mt=start3;// change start value for diffusion process i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+sqrt(Vt)*e1[i]; V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e3[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; v_fact=v_fact|v1; theta_fact=theta_fact|M1; endif; Pt=XX; Vt=V1; Mt=m1; endo; retp(dat,v_fact,theta_fact); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data Generation * * DGP: feed * ************************************************************************************************* * Inputs: T: Length of time series * * h: discretization interval * * kapa_r: mean reversion parameter Interest Rate Process * * kapa_r_v: parameter capturig the effect of volatility on conditional mean * * rho: conditional correlation between short rate and its stochastic volatility * * kapa_v: mean reversion parameter volatility process * * v_bar: mean Volatility Level, also used as the starting value of the vol process * * sig2: variance term Volatility Process * * kapa_theta: mean reversion parameter mean process * * kapa_theta_v: parameter capturig the effect of volatility on mean theta * * theta_bar: mean Level, also used as the starting value of the mean process * * sig3: variance term mean Process * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc(3)= DGP_Feed(T,h,kapa_r,kapa_v,v_bar,sig2,kapa_theta,theta_bar,sig3,kapa_r_v,rho,see); local theta_fact,v_fact,TN,e1,e2,e3,dat,pt,vt,mt,i,xx,v1,m1; v_fact={};theta_fact={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3= sqrt(h)*rndns(TN,1,see); dat={}; Pt=theta_bar; Vt=v_bar; Mt=theta_bar; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+kapa_r_v*(v_bar-Vt)*h+sqrt(Vt)*(e1[i]*((1-rho^2)^0.5)+rho*e2[i]); V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e3[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; v_fact=v_fact|v1; theta_fact=theta_fact|M1; endif; Pt=XX; Vt=V1; Mt=m1; endo; retp(dat,v_fact,theta_fact); endp; proc(3)= DGP_Feedsim(T,h,kapa_r,kapa_v,v_bar,sig2,kapa_theta,theta_bar,sig3,kapa_r_v,rho,see,start1,start2,start3); local theta_fact,v_fact,TN,e1,e2,e3,dat,pt,vt,mt,i,xx,v1,m1; v_fact={};theta_fact={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3= sqrt(h)*rndns(TN,1,see); dat={}; Pt=start1; Vt=start2; Mt=start3;// change start value for diffusion process i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+kapa_r_v*(v_bar-Vt)*h+sqrt(Vt)*(e1[i]*((1-rho^2)^0.5)+rho*e2[i]); V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e3[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; v_fact=v_fact|v1; theta_fact=theta_fact|M1; endif; Pt=XX; Vt=V1; Mt=m1; endo; retp(dat,v_fact,theta_fact); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data Generation * * DGP: CHEN_JUMP * ************************************************************************************************* * Inputs: T: Length of time series * * h: discretization interval * * kapa_r: mean reversion parameter Interest Rate Process * * kapa_v: mean reversion parameter volatility process * * v_bar: mean Volatility Level, also used as the starting value of the vol process * * sig2: variance term Volatility Process * * kapa_theta: mean reversion parameter mean process * * theta_bar: mean Level, also used as the starting value of the mean process * * sig3: variance term mean Process * * sim_ind: indicator for generating the errors * * sim_ind = 0 if errors are to be generated randomly * * sim_ind = 1 if errors are to be be taken from a global variable (SGMM) * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc(3)= DGP_chenJ(T,h,kapa_r,kapa_v,v_bar,sig2,kapa_theta,theta_bar,sig3,lamdau,NUu,lamdad,NUd,see); local theta_fact,v_fact,TN,e1,e2,e3,dat,pt,vt,mt,i,xx,v1,m1,e4,jumpu,p_jumpu,jumpd,p_jumpd,jump; v_fact={};theta_fact={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3= sqrt(h)*rndns(TN,1,see); dat={}; e4 = rndus(TN,4,see); jumpu=-NUu*ln(e4[.,2]); p_jumpu=exp(-lamdau*h); p_jumpu=e4[.,1].>p_jumpu; jumpu=jumpu.*p_jumpu; jumpd=-NUd*ln(e4[.,4]); p_jumpd=exp(-lamdad*h); p_jumpd=e4[.,3].>p_jumpd; jumpd=jumpd.*p_jumpd; jump=jumpu-jumpd; Pt=theta_bar; Vt=v_bar; Mt=theta_bar; i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+sqrt(Vt)*e1[i]+jump[i]; V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e3[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; v_fact=v_fact|v1; theta_fact=theta_fact|M1; endif; Pt=XX; Vt=V1; Mt=m1; endo; retp(dat,v_fact,theta_fact); endp; proc(3)= DGP_chenJsim(T,h,kapa_r,kapa_v,v_bar,sig2,kapa_theta,theta_bar,sig3,lamdau,NUu,lamdad,NUd,see,start1,start2,start3); local theta_fact,v_fact,TN,e1,e2,e3,dat,pt,vt,mt,i,xx,v1,m1,e4,jumpu,p_jumpu,jumpd,p_jumpd,jump; v_fact={};theta_fact={}; TN=round(T/h); e1= sqrt(h)*rndns(TN,1,see); e2= sqrt(h)*rndns(TN,1,see); e3= sqrt(h)*rndns(TN,1,see); dat={}; e4 = rndus(TN,4,see); jumpu=-NUu*ln(e4[.,2]); p_jumpu=exp(-lamdau*h); p_jumpu=e4[.,1].>p_jumpu; jumpu=jumpu.*p_jumpu; jumpd=-NUd*ln(e4[.,4]); p_jumpd=exp(-lamdad*h); p_jumpd=e4[.,3].>p_jumpd; jumpd=jumpd.*p_jumpd; jump=jumpu-jumpd; Pt=start1; Vt=start2; Mt=start3;// change start value for diffusion process i=0; do while i < TN; i=i+1; XX=Pt+kapa_r*(Mt-Pt)*h+sqrt(Vt)*e1[i]+jump[i]; V1=Vt+(kapa_v*(v_bar-Vt))*h+sig2*sqrt(vt)*e2[i]; if V1 < 0; V1=Vt+(kapa_v*(v_bar-Vt))*h; endif; M1=Mt+(kapa_theta*(theta_bar-Mt))*h+sig3*sqrt(Mt)*e3[i]; if M1 < 0; M1=Mt+(kapa_theta*(theta_bar-Mt))*h; endif; if i%(h^-1)==0; dat=dat|XX; v_fact=v_fact|v1; theta_fact=theta_fact|M1; endif; Pt=XX; Vt=V1; Mt=m1; endo; retp(dat,v_fact,theta_fact); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * Data statistics * * *************************************************************************************************/ proc dis_stat(yt); local avg,std,skew,kurt,stat; avg=meanc(yt); std=stdc(yt); skew=meanc(((yt-meanc(yt))/stdc(yt))^3); kurt=meanc(((yt-meanc(yt))/stdc(yt))^4); stat=avg|std|skew|kurt; retp(stat); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * estimate: estimate the parameters * ************************************************************************************************* * Inputs: dat1 - time series * * indicator - model indicator * * Output: est - estimated parameters * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ /* xt: is the lagged xt model_ind; model indicator = 1 CIR = 2 SV = 3 SVJ : = 4 CHEN = 5 CHEN_JUMP = 6 SM */ proc(1)=estimate(xt,model_ind); local b_start,b,f,g,retcode; b={}; if model_ind==1; coset; __output=0; _co_Algorithm=3; _co_LineSearch=4; _co_C = { 1 0 0, 0 1 0, 0 0 1}; _co_D = { 0, 0, 0}; _co_IneqProc = &ineq_one; b_start=0.1701|0.05|0.0249; {b,f,g,retcode} = co(&GMM_cir,b_start); elseif model_ind==2; coset; __output=0; _co_Algorithm=3; _co_LineSearch=4; _co_DirTol=1e-4; _co_C = {1 0 0 0 0 0, 0 1 0 0 0 0, 0 0 1 0 0 0, 0 0 0 1 0 0, 0 0 0 0 1 0, 0 0 0 0 0 1, 0 0 0 0 0 -1}; _co_D = {0, 0, 0, 0, 0, -1, 0}; _co_IneqProc = &ineq_two; b_start=0.27|0.05|3.5|0.00006|0.02|-0.5; {b,f,g,retcode} = co(&GMM_SV,b_start); elseif model_ind==3; coset; __output=0; _co_Algorithm=3; _co_LineSearch=1; _co_DirTol=1e-4; _co_TrustRadius=1; _co_IneqProc = &ineq_two; _co_A = {0 0 0 0 0 0 1 0 0 0, 0 0 0 0 0 0 0 1 0 0, 0 0 0 0 0 0 0 0 1 0, 0 0 0 0 0 0 0 0 0 1}; _co_B = {5.4979, 0.0006, 2.121, 0.002};// add equality constraint for this one, becoz H fails, fix lamda in SVJ estimation _co_C = { 1 0 0 0 0 0 0 0 0 0, 0 1 0 0 0 0 0 0 0 0, 0 0 1 0 0 0 0 0 0 0, 0 0 0 1 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0, 0 0 0 0 0 1 0 0 0 0, 0 0 0 0 0 -1 0 0 0 0}; _co_D = {0,0,0,0,0, -1, 0}; b_start=0.230503|0.052496|2.292159|0.000019|0.0100042|-0.105|7.899973|0.001056|1.600042|0.001922; {b,f,g,retcode} = co(&gmm_SVJ,b_start); //print "f=" f; elseif model_ind==4; coset; __output=0; _co_Algorithm=3; _co_LineSearch=1; _co_Trust=1; _co_IneqProc = &ineq_three; b_start=0.2|3|0.0003|0.01|0.3|0.06|0.04; _co_A = {0 0 0 0 1 0 0 }; _co_B = {0.29};// add equality constraint for this one, becoz H fails _co_C = {1 0 0 0 0 0 0, 0 1 0 0 0 0 0, 0 0 1 0 0 0 0, 0 0 0 1 0 0 0, 0 0 0 0 1 0 0, 0 0 0 0 0 1 0, 0 0 0 0 0 0 1}; _co_D = {0, 0, 0, 0, 0, 0, 0}; {b,f,g,retcode} = co(&gmm_C,b_start); //print "f=" f; elseif model_ind==5; coset; __output=0; _co_Algorithm=3; _co_LineSearch=1; _co_Trust=1; _co_IneqProc = &ineq_three; _co_DirTol=1e-4; b_start=0.363825|2.119171|0.000465|0.00555|0.287034|0.060095|0.184442|8.106309|0.000205|2.828375|0.000471; _co_C = {1 0 0 0 0 0 0 0 0 0 0, 0 1 0 0 0 0 0 0 0 0 0, 0 0 1 0 0 0 0 0 0 0 0, 0 0 0 1 0 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0 0, 0 0 0 0 0 1 0 0 0 0 0, 0 0 0 0 0 0 1 0 0 0 0, 0 0 0 0 0 0 0 1 0 0 0, 0 0 0 0 0 0 0 0 1 0 0, 0 0 0 0 0 0 0 0 0 1 0, 0 0 0 0 0 0 0 0 0 0 1}; _co_D = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; {b,f,g,retcode} = co(&gmm_CJ,b_start); elseif model_ind==6; coset; __output=0; _co_A = {1 0 0 0 0}; _co_B = {0.5};// add equality constraint for this one, becoz H fails _co_C = {1 0 0 0 0, 0 1 0 0 0, 0 0 1 0 0, 0 0 0 1 0, 0 0 0 0 1}; _co_D = {0, 0, 0, 0, 0}; _co_IneqProc = &ineq_two; b_start=1.8|0.03|0.3|0.06|0.04; {b,f,g,retcode} = co(&GMM_SM,b_start); //print " f=" f; endif; retp(b); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * CS_stat: Calculates The Test Statistics * ***************************************************************************************************************** * Inputs: xxt - time series R The first R observations that used for estimation * * tao - a vector indicating the different step ahead con. int. to be examined * * S - number of sample paths to be simulated N - number of sample paths to be simulated for SV * * mod_ind1 - model specification 1 * mod_ind2 - model specification 2 * * h - discretization interval * * u_bar - confidence interval * * Output: * vt - Test statistics will have rows(tao) results v1 First part of Test statistics v2 Second part of Test statistics * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc (1)= CS_stat(xxt,R,N,tao,S,mod_ind1,mod_ind2,h,u_bar); local T,tt,xp,p_true,b,f,g,retcode,vt,p_sim1,p_sim2,i,sup_vt,b1,b2,see,v1,v2; see=12343;//dont change seed, to use the same random error,but different start value of X(t) T=rows(xxt); vt=zeros(rows(tao),1); v1=zeros(rows(tao),1); v2=zeros(rows(tao),1); i=0; do while i u_bar[1,1]).*(xp .< u_bar[2,1]); xt=xxt[1:tt,.];//set xt p_sim1=Con_den(xt,tao[i],S,N,mod_ind1,h,u_bar,see);// prob(u1 u_bar[1,1]).*(xp .< u_bar[2,1]); xt=xxt[1:tt,.]; p_sim1=Con_den2(xt,tao[i],S,N,model_ind1,h,u_bar,see,svSim1,mvSim1);// prob(u1 u_bar[1,1]).*(xp .< u_bar[2,1]); xt=xxt[1:tt,.]; p_sim1=Con_den3(xt,tao[i],S,N,model_ind1,h,u_bar,see);// prob(u1 u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==2; {x_sim2,svSim}=DGP_SV(N,h,b[1],b[2],b[3],b[4],b[5],b[6],see);//simulate N times to get mvSim for j(1,N,1); {x_sim,sv}=DGP_SVsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; see=see+3; endfor; p_sim=p_sim/N; elseif model_ind==3; {x_sim2,svSim}=DGP_SVJ(N,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see);//simulate N times to get mvSim for j(1,N,1); {x_sim,sv}=DGP_SVJsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==4; {x_sim2,svSim,mvSim}=DGP_chen(N,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],see);//simulate N times to get mvSim for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chensim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==5; {x_sim2,svSim,mvSim}=DGP_chenj(N,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],b[11],see);//simulate N times to get mvSim for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chenjsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],b[11],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==6; {x_sim2,mvSim}=DGP_SM(N,h,b[1],b[2],b[3],b[4],b[5],see);//simulate N times to get mvSim for j(1,N,1); {x_sim,mv}=DGP_SMsim(tao,h,b[1],b[2],b[3],b[4],b[5],see,xtt,mvSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; endif; // calculate the probability of prob(u1 u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==2; for j(1,N,1); {x_sim,sv}=DGP_SVsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==3; for j(1,N,1); {x_sim,sv}=DGP_SVJsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==4; for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chensim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==5; for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chenjsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],b[11],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==6; {x_sim2,mvSim}=DGP_SM(N,h,b[1],b[2],b[3],b[4],b[5],see);//simulate N times to get mvSim for j(1,N,1); {x_sim,mv}=DGP_SMsim(tao,h,b[1],b[2],b[3],b[4],b[5],see,xtt,mvSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; endif; // calculate the probability of prob(u1 u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==2; {x_sim,sv}=DGP_SVsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],see,xtt,b[4]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==3; {x_sim,sv}=DGP_SVJsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see,xtt,b[4]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endif; // calculate the probability of prob(u1 u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==2; for j(1,N,1); {x_sim,sv}=DGP_SVsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==3; for j(1,N,1); {x_sim,sv}=DGP_SVJsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==4; for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chensim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==5; for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chenjsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],b[11],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==6; {x_sim2,mvSim}=DGP_SM(N,h,b[1],b[2],b[3],b[4],b[5],see);//simulate N times to get mvSim for j(1,N,1); {x_sim,mv}=DGP_SMsim(tao,h,b[1],b[2],b[3],b[4],b[5],see,xtt,mvSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; endif; // calculate the probability of prob(u10; undraw2=round((N-lval)*rndus(1,1,see)); xbl=xbl|dat1[undraw2+1:undraw2+tt,.]; endif; retp(xbl); endp; /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ * CS_statBoot2: Calculates The Test Statistics : Simulate latent ONCE, which need a latent variable series as input for the simulation process. * ***************************************************************************************************************** * Inputs: xxt - time series R The first R observations that used for estimation * * tao - a vector indicating the different step ahead con. int. to be examined * * S - number of sample paths to be simulated N - number of sample paths to be simulated for SV * * model_ind1 - model specification 1 * model_ind2 - model specification 2 * * h - discretization interval * * u_bar - confidence interval svSim1,mvSim1,svSim2,mvSim2 - simulated series for latent variables option - =1: using the same parameters =2: estimate new parameters * * Output: * vt - Test statistics will have rows(tao) results * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ proc (1)= CS_statBoot2(xxt,R,N,tao,S,model_ind1,model_ind2,h,u_bar,svSim1,mvSim1,svSim2,mvSim2,option); local T,tt,xp,p_true,b,f,g,retcode,vt,p_sim1,p_sim2,i,sup_vt,b1,b2,x_sim1, x_sim2,see,v1,v2,center1,center2,bb1,bb2,K; see=12343;//dont change seed, to use the same random error,but different start value of X(t) T=rows(xxt); //test vt=zeros(rows(tao),1); v1=zeros(rows(tao),1); v2=zeros(rows(tao),1); i=0; do while i u_bar[1,1]).*(xp .< u_bar[2,1]); xt=xxt[1:tt,.]; if option==1; p_sim1=Con_den2(xt,tao[i],S,N,model_ind1,h,u_bar,see,svSim1,mvSim1);// prob(u1 u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==2; for j(1,N,1); {x_sim,sv}=DGP_SVsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==3; for j(1,N,1); {x_sim,sv}=DGP_SVJsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see,xtt,svSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; elseif model_ind==4; for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chensim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==5; for j(1,N,1); for jj(1,N,1); {x_sim,sv,mv}=DGP_chenjsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],b[11],see,xtt,svSim[j],mvSim[jj]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; endfor; p_sim=p_sim/(N^2); elseif model_ind==6; {x_sim2,mvSim}=DGP_SM(N,h,b[1],b[2],b[3],b[4],b[5],see);//simulate N times to get mvSim for j(1,N,1); {x_sim,mv}=DGP_SMsim(tao,h,b[1],b[2],b[3],b[4],b[5],see,xtt,mvSim[j]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=p_sim +((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endfor; p_sim=p_sim/N; endif; // calculate the probability of prob(u1 u_bar[1,1]).*(xp .< u_bar[2,1]); xt=xxt[1:tt,.]; if option==1; p_sim1=Con_den3(xt,tao[i],S,N,model_ind1,h,u_bar,see);// prob(u1 u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==2; {x_sim,sv}=DGP_SVsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],see,xtt,b[4]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; elseif model_ind==3; {x_sim,sv}=DGP_SVJsim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],see,xtt,b[4]);// we need to put 2 starts here, one for Xt, one for Vt p_sim=((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))'; endif; // calculate the probability of prob(u1