/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*  1. estimate model using first R obs
   2. simulate x_sim(t+tao), for the out-of-sample P
   3. calculate conditional distribution  prob(u1<x_sim(t+tao)|x(t)<u2)= sumc (I(u1<x_sim(t+tao)|x(t)<u2))/N  
   4. CS test
   5. Bootstrap Critical Value

                             *
*************************************************************************************************
                                                             *
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
new; cls;
library co,pgraph;
#include co.ext
clearg xt;
see = 39378;
/*
xt: is the lagged xt
model_ind; model indicator
    = 1 SM
    = 2 SV
    = 3 SVJ : need to change the jumps
    = 4 CHEN
    = 5 CHEN_JUMP
    =  FEEDBACk // we don't want this model    
    = 6 CIR
*/
//output file = D:\research\continuous\code\BCS\BCSresult_euro_whole.txt reset; 
outwidth 255; output on; format /MA1 /LD 6,4;
load data[] = D:\research\continuous\code\BCS\1monthEuro.txt;//weekly data


/*  1978-1987: 368-894
    1988-1997: 895-1419
    1998-2008: 1420-1961 
    1971-1980: 1-525  */

xt=data/100;




/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
* CS out-of-sample specification test                                                          
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */


tao=1|2|4|12;  // steps of ahead simulation
                             
//S=2*(rows(xt)-tao[rows(tao)]); //S=1,2,5*T
S=1000; //simulation times
N=1;//Simulation for SV, here we use V_bar as start value, so the simulation time is 1

hh=1/2; // descretlize each interval
u_bar=(meanc(xt)-0.5*stdc(xt))|(meanc(xt)+0.5*stdc(xt));//the interval u_bar


//**************** the BCS specification test ****************************
T=rows(data); 
R=T;// whole sample estimation


y1=time;

b1=estimate(xt[1:R,1],1);
b2=estimate(xt[1:R,1],2);
b3=estimate(xt[1:R,1],3);
b4=estimate(xt[1:R,1],4);
b5=estimate(xt[1:R,1],5);
b6=estimate(xt[1:R,1],6);//delete FEED, put CIR here 


parameter=b1|b2|b3|b4|b5|b6;




// BCS test, we estimate each model in BCS_stat

lval=20;

for mod_ind(1,6,1);
    {vt,sup_vt}= BCS_stat(xt,tao,S,N,mod_ind,xt,hh,u_bar,see,parameter);   
    {b_sup_vt,cv95,cv90,cv85,cv80}=BCS_boot(xt,tao,S,N, mod_ind,lval,100,vt,hh,u_bar,see,parameter);
    print " this is results for model=" mod_ind;
    print " The BCS stat for each model is :" sup_vt;
    print " The bootstrap critical cv95,cv90,cv85,cv80 are :" cv95'~cv90'~cv85'~cv80';
endfor;


y2=time;
print "time is :" y2-y1;
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 SM
    = 2 SV
    = 3 SVJ : need to change the jumps
    = 4 CHEN
    = 5 FEEDBACk
    = 6 CHEN_JUMP
    = 7 CIR
*/
proc(1)=estimate(xt,model_ind);
local b_start,b,f,g,retcode,cov0,SD_fnl;
b={};
    if model_ind==1;
        coset;
        __output=0;
        _co_Diagnostic = 1;
        _co_A = {0 0 1 0 0};
        _co_B = {0.300};// 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);
        // get std of parameters
        cov0=inv(hessp(&GMM_SM,b)); 
        SD_fnl =sqrt(diag(cov0)); @Standard errors of the estimated coefficients@
        print " std are=" SD_fnl;

print " f=" f;

    elseif model_ind==2;
        coset;
        __output=0;
        _co_Algorithm=3;
        _co_LineSearch=4;
        _co_DirTol=1e-3;
        _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.3|0.06|3.5|0.00006|0.02|-0.5;
        {b,f,g,retcode} = co(&GMM_SV,b_start);

        // get std of parameters
        cov0=inv(hessp(&GMM_SV,b)); 
        SD_fnl =sqrt(diag(cov0)); @Standard errors of the estimated coefficients@
        print " std are=" SD_fnl;

    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.0011,
                3.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 1 0 0 0 0 0 };
        //_co_B = {2.9051};// 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;//CIR
        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.1801|0.05|0.02549;            
        {b,f,g,retcode} = co(&GMM_cir,b_start);


  endif;
retp(b);
endp;


/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
* CS_stat: Calculates The Test Statistics                                                           		    *
*****************************************************************************************************************
* Inputs:	xt      -   	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                                	    *
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
proc (1)= CS_stat(xt,R,N,tao,S,mod_ind1,mod_ind2,h,u_bar,parameter);
local T,tt,xp,p_true,b,f,g,retcode,see,vt,p_sim1,p_sim2,i,sup_vt,b1,b2;
    see=1234;//dont change seed, to use the same random error,but different start value of X(t)
    T=rows(xt); 
    vt=zeros(rows(tao),1);
    b1=parameter[6:11,1];// default model is SV, mod_ind1=2;
    if mod_ind2==1;
        b2=parameter[1:5,1];
    elseif mod_ind2==3;
        b2=parameter[12:21,1];
    elseif mod_ind2==4;
        b2=parameter[22:28,1];
    elseif mod_ind2==5;
        b2=parameter[29:39,1];
    elseif mod_ind2==6;
        b2=parameter[40:42,1];
    endif;

    i=0;
    do while i <rows(tao); 
        i=i+1;
        tt=R;// the start of out-of-sample 
        do while tt<T-tao[i];       
            xp=xt[tt+tao[i],.];// the actual value of X(P), tt is the in-sample obs
            p_true=(xp .> u_bar[1,1]).*(xp .< u_bar[2,1]); 
            p_sim1=Con_den(xt[1:tt,.],tao[i],S,N,mod_ind1,h,u_bar,see,b1);// prob(u1<x_sim(t+tao)<u2|x(t))
            p_sim2=Con_den(xt[1:tt,.],tao[i],S,N,mod_ind2,h,u_bar,see,b2);// prob(u1<x_sim(t+tao)<u2|x(t))
            vt[i]=vt[i]+(p_sim1-p_true).^2-(p_sim2-p_true).^2; 
    //print " p_sim1~p_sim2~p_true~vt[i]==" p_sim1~p_sim2~p_true; 
            tt=tt+1;
        endo;
        vt[i]=vt[i]./sqrt(T-R-tao[i]+1); 
        
    endo;
retp(vt);
endp;


/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*get simulated conditonal density  
  
format:  prob(u1<x_sim(t+tao)<u2|x(t))= sumc (I(u1<x_sim(t+tao)<u2|x(t)))/S                                                                   *

*****************************************************************************************************************
* Inputs:	xt      -   	time series, The first R observations that used for estimation                                                             	        *
*		    tao     -   	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 or MV part 
            model_ind -     model indicator 
*		    h       -   	discretization interval, for each step,make it 1/h small intervals                                                		        *
*		    u_bar   -   	confidence interval  
            see     -       simulation seed
                                   		            *
* Output:   p_sim      -   	generalized residual ,    			            *
*           	                             	    *
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
proc (1)= Con_den(xt,tao,S,N,model_ind,h,u_bar,see,b);
local x_sim,x_sim2,p_sim,pp_sim,i,ii,j,jj,mv,xtt,mvSim,sv,svSim;
        p_sim=0; // prob(u1<x_sim<u2) for each simulation, may be averaged value for SV
        pp_sim=0;// average [ prob(u1<x_sim<u2) ]       
        xtt=xt[rows(xt),.];// the initial value for simulation
        ii=0;

        do while ii<S;ii=ii+1;
            x_sim={}; x_sim2={}; mvSim={}; mv={}; 
            see=see+3;// change seed to change each path;
           
            if model_ind==6;       
                x_sim=DGP_CIRsim(tao,h,b[1],b[2],b[3],see,xtt); // simulate based on Xt(i), simulat S times;
           
            elseif model_ind==1;                   
                    {x_sim,mv}=DGP_SMsim(tao,h,b[1],b[2],b[3],b[4],b[5],see,xtt,b[4]);// we need to put 2 starts here, one for Xt, one for Vt
            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  
                         
            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
    
            elseif model_ind==4;                  
                    {x_sim,sv,mv}=DGP_chensim(tao,h,b[1],b[2],b[3],b[4],b[5],b[6],b[7],see,xtt,b[3],b[6]);// we need to put 2 starts here, one for Xt, one for Vt

            elseif model_ind==5;   
                    {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,b[3],b[6]);// we need to put 2 starts here, one for Xt, one for Vt
                     
            endif;
  //print "x_sim=" x_sim;               
            p_sim=((x_sim[tao] > u_bar[1,1]).*(x_sim[tao]< u_bar[2,1]))';          


            // calculate the probability of  prob(u1<x_sim(t+tao)<u2|x(t)) for S simulations
            pp_sim=pp_sim +p_sim;     
                 
        endo; 
  
pp_sim=pp_sim/S;  

retp(pp_sim);
endp;



/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
* boot1: Generates the block bootstrap sample                                                   *
*************************************************************************************************
* Inputs:	dat1    -   time series                                                             *
*		    lval    -   block boot length                                                       *
* Output:   xb1     -   bootstrap sample                                                        *
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ 
proc (1) = boot1(dat1,lval,see);
local N,num_uns,undraw1,xbl,ib;
    N=rows(dat1);
    num_uns=N/lval;
    /* draw uniforms U[0,T-Lval+1] */
    undraw1=round((N-lval)*rndus(num_uns,1,see));
    xbl={};
    ib=1;
    do while ib<=num_uns;
    xbl=xbl|dat1[undraw1[ib]+1:undraw1[ib]+lval,.];
    ib=ib+1;
    endo;
retp(xbl);
endp;


/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
* BCS_stat: Calculates The Test Statistics                                                           		    *
*****************************************************************************************************************
* Inputs:	xt      -   	time series                                                             	        *
*		    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 Volatility, stochastic mean                
                                            *
*		    mod_ind2 -   	model specification                          	        *
*        	xt_v    -   	time series for calculating v 		                                     	        *
*		    h       -   	discretization interval                                                		        *
*		    u_bar   -   	confidence interval    
            see       -     seed    
            parameter   estimated out of program
                                         		            *
* Output:   	vt      -   	pre sup pre sum vt to use in bootstrap procedure      				            *
*           	sup_vt  -   	Test statistics will have rows(tao) results                                	    *
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
proc (2)= BCS_stat(xt,tao,S,N,mod_ind2,xt_v,h,u_bar,see,parameter);
local v,p_true,b,i,tt,x_sim,p_sim,jj,v_t,vt,sup_vt,v_ind,xp,p_sim1,xtt,sv,mv,T,sup_v;
    see=123456;//dont change seed, to use the same random error,but different start value of X(t)
    T=rows(xt); 
    v_t=zeros(rows(tao),6); //save vt in eq(11)
    sup_vt={};//eq (10)

    
    v=sortc(xt_v,1);
	v=v[0.1*rows(v)]|v[0.3*rows(v)]|v[0.5*rows(v)]|v[0.7*rows(v)]|v[0.9*rows(v)]|v[rows(v)];
       
    /*estimate the model using whole sample obs 
    if mod_ind2==1;
        b=estimate(xt,1);
    elseif mod_ind2==2;
        b=estimate(xt,2);
    elseif mod_ind2==3;
        b=estimate(xt,3);
    elseif mod_ind2==4;
        b=estimate(xt,4);
    elseif mod_ind2==5;
        b=estimate(xt,5);
     elseif mod_ind2==6;
        b=estimate(xt,6);
    endif;
    */

    if mod_ind2==1;
        b=parameter[1:5,1];
    elseif mod_ind2==2;
        b=parameter[6:11,1];
    elseif mod_ind2==3;
        b=parameter[12:21,1];
    elseif mod_ind2==4;
        b=parameter[22:28,1];
    elseif mod_ind2==5;
        b=parameter[29:39,1];
    elseif mod_ind2==6;
        b=parameter[40:42,1];//CIR
    endif;

    i=0;
    do while i <rows(tao); 
        i=i+1;
        vt={};// eq (11)
        tt=1;// the start of simulation, start from 1.
        do while tt<=T-tao[i];       
            xp=xt[tt+tao[i],.];// the actual value of X(P), tt is the in-sample obs
            p_true=(xp .> u_bar[1,1]).*(xp .< u_bar[2,1]); 
            p_sim1=Con_den(xt[1:tt,.],tao[i],S,N,mod_ind2,h,u_bar,see,b);// prob(u1<x_sim(t+tao)<u2|x(t))           
            vt=vt|(p_sim1-p_true); 
           // print " p_sim1~p_true~vt[i]==" p_sim1~p_true; 
            tt=tt+1;
        endo;

       // get {x(t)<v}
        v_ind=(xt[1:T-tao[i]] .< v[1])~(xt[1:T-tao[i]] .< v[2]);
        v_ind=v_ind~(xt[1:T-tao[i]] .< v[3])~(xt[1:T-tao[i]] .< v[4]);
        v_ind=v_ind~(xt[1:T-tao[i]] .< v[5])~(xt[1:T-tao[i]] .< v[6]);
        
        v_t[i,.]=vt'*v_ind/sqrt(T-tao[i]);

        sup_vt=sup_vt~(maxc(abs(v_t[i,.]')));       
    //print "sup=" sup_vt;    
    endo;

retp(v_t,sup_vt);
endp;

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
* BCS_boot: Calculates The Statistics                                                           *
*************************************************************************************************
* Inputs:	xt      -   time series                                                             *
*		    tao     -   a vector indicating the different step ahead con. int. to be examined   *
*		    s       -   number of sample paths to be simulated                                  *
*		    mod_ind -   model specification 1 for CIR and 2 for log ou                          *
*           vt      -   pre sup pre sum vt to use in bootstrap will have rows(tao)*18 cols      *
*           bL      -   block boot lengt                                                        *
*           bs      -   number of bootstrap replication  
            h       -   discretization                                       *
*		    u_bar   -   confidence interval       
            see     -   seed                                         		*
            parameter   estimated out of program
* Output:   b_sup_vt-   boot strap distribution, sorted 100 rows and rows(tao)*3 cols           *
*           cv95    -   95% critical value                                                      *
*           cv90    -   90% critical value                                                      *
*           cv80    -   80% critical value                                                      *
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/  
proc (5)=BCS_boot(xt,tao,S,N, mod_ind,bL,bs,vt,h,u_bar,see,parameter);
local xb,ii,b_vt,sup_vt,t_tao_b,vt_store,b_sup_vt,jj,cv95,cv90,cv85,cv80;
        
    ii=0;b_sup_vt={};
    do while ii < bs;ii=ii+1;
        xb=boot1(xt,bL,see);
        {b_vt,sup_vt}=BCS_stat(xb,tao,S,N,mod_ind,xb,h,u_bar,see,parameter);
        b_vt=abs(b_vt-vt);//recenter
        b_sup_vt=b_sup_vt|maxc(b_vt')';
        see=see+100;//change seed
    endo;

//print " bootstrap sup_vt for tao=1,2,4,12, are: " b_sup_vt;
    jj=0;
    do while jj< cols(b_sup_vt);
        jj=jj+1;
        b_sup_vt[.,jj]=sortc(b_sup_vt[.,jj],1);
    endo;

    cv95=b_sup_vt[round(0.95*rows(b_sup_vt)),.];
    cv90=b_sup_vt[round(0.90*rows(b_sup_vt)),.];
    cv85=b_sup_vt[round(0.85*rows(b_sup_vt)),.];
    cv80=b_sup_vt[round(0.80*rows(b_sup_vt)),.];
retp(b_sup_vt,cv95,cv90,cv85,cv80);
endp;