new;

/* this code is written for the paper on DCC April 2007*/
library pgraph;
format /m1 /rd 8,4;

//rndseed 5454545;



// DGP.1 [CCC-Garch(1,1)]  with conditional mean of Y normalized to zero */
// first get the iid z vectors


/* Competing Test Statistics */


//=================================================//
//  DATA ENTRY
//=================================================//


//T = 100;
//load dgp[T,2] = c:\gauss6.0\GARCH\dgp1.txt;

//T =250;
//load dgp[T,2] = c:\gauss6.0\GARCH\dgp2.txt;

T =500;
//bekkpar = 0.2|0.1|0.1|0.2|0.6|0.2|0.2|0.6|0.3|0.1|0.1|0.3;

//dgp = simFBEKK(T,2,bekkpar);
//corrx(dgp);
load dgp[T,2] = c:\gauss6.0\GARCH\dgp3.txt;

//T = 1000;
//load dgp[T,2] = c:\gauss6.0\GARCH\dgp4.txt;

//corrx(dgp);
//stop;

//par = 0.01| 0.94| 0.05| 0.5| 0.5| 0.2;
//corrmat = {1 0.5, 0.5 1};
//T = 500;
//N = 2; 

//{data,Ht} = simCCC(T,N,corrmat,par);

//dgp = data;

//corrx(dgp);

//dgp1 = (dgp-meanc(dgp)')./stdc(dgp)';




           //=========================================================//
                            /* Competing Test Statistics */
           //=========================================================//

//==================================================================================================//
/* Test statistics of Bera and Kim (2002) */

proc BK(Z);  // here Z is T by 2 matrix of standardized residuals
local a1t, a2t, rhohat, z1, z2, BKSTAT, pval;
rhohat  = meanc(Z[.,1].*Z[.,2]);
a1t = (1/sqrt(1-rhohat^2))*(Z[.,1] - rhohat*Z[.,2]);
a2t = (1/sqrt(1-rhohat^2))*(Z[.,2] - rhohat*Z[.,1]);
BKSTAT = inv(4*T*(1+4*rhohat^2 +rhohat^4))*((sumc((a1t^2).*(a2t^2) - 1 - 2*rhohat^2)))^2;
pval = cdfchic(BKSTAT,1);
print " Test statistics of Bera and Kim (2002) ";
print "statistic"  " "  "P-Value"  " "   "RhoHat";
retp(BKSTAT~pval~rhohat);
endp;

//==================================================================================================//

/* Test statistic of Tse (2000) */
/* Here the input, Sh, to this statistic is the estimate of the score evaluated  
   under the null of constant conditional correlation  */

proc Tse(Sh);
local TseSTAT, pval;
TseSTAT = ones(T,1)'Sh*invpd(Sh'Sh)*Sh'ones(T,1);
pval = cdfchic(TseSTAT,1);
print " Test statistic of Tse (2000) "; 
print " statistic "  "  "  "P-Value" ;
retp(TseSTAT~pval);
endp;

//==================================================================================================//

/* Test statistic of Engle and Sheppard (2000) */
// here Z - is the vector of standardized residuals
// p - is the lag order

proc ES(Z,p);
local ES_STAT, At, AAt, rhohat, B, i, B1, gamhat, sig2, pval;
if T <= p;
errorlog " T MUST BE STRICTLY GREATER THAN p";
endif;
B = zeros(T-p,p);
rhohat  = meanc(Z[.,1].*Z[.,2]);
At = inv(rhohat)*(Z[.,1].*Z[.,2]);
AAt = trimr(At,p,0);
for i(1,p,1);
    B[.,i] = trimr(At,p-i,i);
endfor;
B1 = ones(T-p,1)~B;
//det(B1'B1);
//endif;
//ols(0,AAt,B1);
gamhat = invpd(B1'B1)*B1'AAt;
//gamhat;
sig2 = inv(T-p)*((AAt'AAt)- gamhat'B1'AAt);
//sig2;
ES_STAT = inv(sig2)*(gamhat'B1'B1*gamhat);
pval = cdfchic(ES_STAT,p+1);
print " Test statistic of Engle and Sheppard (2000) ";
print "statistic"  " "  "P-value"   "   "   "lagnum";
retp(ES_STAT~pval~p);
endp;


//================================================================//




et=hsec;

library cml,pgraph;
#include cml.ext;
#include pgraph.ext;
cmlset;
graphset;

N = cols(dgp);
Y0 = dgp;
Y01 = dgp[.,1];
Y02 = dgp[.,2];




//================================================================//
// UNIVARIATE LOGLIKELIHOOD 
//================================================================//
    

    proc ulikelhd(B,data);
    local e,hh,j,logL;
    e=zeros(T,1);
    hh=e;

    e[1]=data[1];
    hh[1]=stdc(data)^2;
    for j(2,T,1);
        e[j]=data[j];
        hh[j]=B[1]+B[2]*hh[j-1]+B[3]*e[j-1]^2;
   endfor;

    logL = -0.918938533204673 - 0.5*ln(hh) - 0.5*(e^2)./hh;
    retp(logL);
endp;


/* getting estimated residual, standardized residual, and cond variance */

proc (5) = ezh(B,data);
    local e,hh,j,logL,e2,z,z2;
    e=zeros(T,1); hh=e;

    e[1]=data[1];
    hh[1]=stdc(data)^2;

    for j(2,T,1);
        e[j]=data[j];
        hh[j]=B[1]+B[2]*hh[j-1]+B[3]*e[j-1]^2;
   endfor;
    
    e2=e.*e;
    z=e./sqrt(hh);
    z2=e2./hh;
    retp(e,e2,z,z2,hh);
endp;

//================================================================//
// ESTIMATING UNIVARIATE GARCH(1,1) by CML
//================================================================//


// Starting values

B00=0.05|0.80|0.15;
B01=0.5|0.5|0.2;

@ Constraints @
_cml_Bounds ={1e-8 1e256, 0 1e256, 0 1e256};

/*
  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(3,2).*_ww_;
  _cml_Bounds[2:3,1] = zeros(2,1);

_cml_Bounds[1,1] = 0.001;
_cml_c = zeros(1,3);
_cml_c[1,2:3] = -ones(1,2);
_cml_d = -0.99;
*/

_cml_Algorithm = 4;
_cml_DirTol = 1e-5;
_cml_MaxIters = 100;
format /rd 13,8;


 {B1,f1,g1,hh1,retcode}=cmlPrt(cml(Y01,0,&ulikelhd,B00));
 
 {B2,f2,g2,hh2,retcode}=cmlPrt(cml(Y02,0,&ulikelhd,B00));

  "";

 B = B1|B2;

 {e1,e12,z1,z12,hh11} = ezh(B1,Y01);

 {e2,e22,z2,z22,hh21} = ezh(B2,Y02);

 hmat = hh11~hh21; //The T by N matrix of conditional variances from the univariate models

 ustdres = z1~z2;  // univariate standardized residuals

RR = corrx(ustdres); // sample correlation matrix from the univariate standardized residuals;

RR;

"";
B;



//==================================================================================================//



/*
proc ccivech(rho);
    local sz, fz, corrz;
    sz = {-1,1};  // shift amount
    fz = rho*ones(dim,1);  // replacement vector
    corrz = shiftr(ones(dim,dim),sz,fz);
retp(corrz);
endp;

*/

/*
//==================================================================================================//

library cml,pgraph;
#include cml.ext;
#include pgraph.ext;
cmlset;
graphset;

B00=B00|0.05|0.80|0.15;

@ Constraints @
//_cml_Bounds ={1e-8 1e256, 0 1e256, 0 1e256, 1e-8 1e256, 0 1e256, 0 1e256};


  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(6,2).*_ww_;
  _cml_Bounds[2:3,1] = zeros(2,1);

_cml_Bounds[1,1] = 0.01;
_cml_Bounds[4,1] = 0.01;

_cml_c = zeros(2,6);
_cml_c[1,2:3] = -ones(1,2);
_cml_c[2,5:6] = -ones(1,2);
_cml_d = -0.99*ones(2,1);


_cml_Algorithm = 4;
_cml_DirTol = 1e-5;
_cml_MaxIters = 50;
format /rd 13,8;


proc likelihd(B,d);
    local zz, hh, DD, j, logL, ccmat, Rinv, Rdetlog;
    logL = zeros(T,1);
    zz = zeros(T,N);
    hh = zz;
    DD = zeros(T,N);
    //zz[1,.] = ustdres[1,.];
    zz[1,.] = Y0[1,.];
    hh[1,.] = (stdc(Y0).^2)';
   // hh[1,.] = hmat[1,.];
    DD[1,.]  = sqrt(hh[1,.]);
    for j(2,T,1);
        hh[j,1] = B[1] + B[2]*hh[j-1,1] + B[3]*Y0[j-1,1]^2;
        hh[j,2] = B[4] + B[5]*hh[j-1,2] + B[6]*Y0[j-1,2]^2;
        DD[j,.]  = sqrt(hh[j,.]);
        zz[j,.] = Y0[j,.]./DD[j,.];
    endfor;
   
  ccmat = corrx(zz);

  Rinv = inv(ccmat);
  Rdetlog = log(det(ccmat));

   for j(1,T,1);
        logL[j] = -0.5*(N*log(2*pi) + sumc(log(hh[j,.])') + Rdetlog + zz[j.,]*Rinv*zz[j,.]');
        //logL[j] = logL[j];
    endfor;
retp(logL);
endp;



//{BB,f,g,hh,retcode}=cmlPrt(cml(Y0,0,&likelihd,B00));

"";

//B~BB;

//==================================================================================================//

// the procedure infomat(BB,d) computes the T by N score matrix 
// BB - is the vector of MLE
// d - is the data matrix to be passed through the likelihood function

proc infomat(BB,d);
    local kk, score, lf, lf1, scovmat;
    kk = rows(BB);
    score = zeros(T,kk);
    lf = likelihd(BB,d);
    for i(1,kk,1);
        BB[i] = BB[i] + 0.00001;
        lf1 = likelihd(BB,d);
        score[.,i] = (lf1 - lf)/0.00001;
        BB[i] = BB[i] - 0.00001;
    endfor;
    scovmat = score;    
    retp(scovmat);
endp;
"";

//Sh = infomat(BB,Y0);
//Tse(Sh);

*/
//==================================================================================================//


//=================================================================================//
//  ESTIMATING Engle's DCC-MVGARCH(1,1) MODEL
//=================================================================================//
/*
library cml,pgraph;
#include cml.ext;
#include pgraph.ext;
cmlset;
graphset;

_cml_Bounds ={0 1e256, 0 1e256};


_cml_Algorithm = 4;
_cml_DirTol = 1e-5;
_cml_MaxIters = 50;
format /rd 13,8;


stdresid = ustdres;

// restricted (partial) likelihood 

proc plikelhd(dccpar,stdresid);
    local a, b, Qbar, m, Qt, Rt, liklhd, j;
    a = dccpar[1];
    b = dccpar[2];
    Qbar = corrx(stdresid);
    m = 1;
    Qt = zeros(N,N*(T+m));
    Rt = zeros(N,N*(T+m));
    Qt[.,1:N] = Qbar;
    Rt[.,1:N] = Qbar;
    liklhd = zeros(T+m,1);
    stdresid = zeros(m,N)|stdresid;
    for j(2,T+m,1);
        Qt[.,(N*j)-1:(N*j)] = Qbar*(1-a-b) + a*stdresid[j-1,.]'stdresid[j-1,.] + b*Qt[.,(N*(j-1))-1:(N*(j-1))];
        Rt[.,(N*j)-1:(N*j)] = Qt[.,(N*j)-1:(N*j)] ./ ( sqrt(diag(Qt[.,(N*j)-1:(N*j)]))* ( sqrt(diag(Qt[.,(N*j)-1:(N*j)])))');
        liklhd[j] = -0.5*(  log(det(Rt[.,(N*j)-1:(N*j)])) + stdresid[j,.]*inv(Rt[.,(N*j)-1:(N*j)])*stdresid[j,.]');
    endfor;
        liklhd = liklhd[2:T+m];
    retp(liklhd);
endp;


// estimating the time-varying DCC parameters  
dccstart = 0.01|0.97;
{dccpar,fdcc,gdcc,hhdcc,retcode}=cmlPrt(cml(stdresid,0,&plikelhd,dccstart));
*/


/*
//==================================================================================================//

// full concentrated likelihood for the DCC Model: this procedure is NOT used in maximization

proc(4) = fconlike(dccpar,stdresid);
    local a, b, Qbar, m, Qt, Rt, liklhd, j, likeval;
    a = dccpar[1];
    b = dccpar[2];
    Qbar = corrx(stdresid);
    m = 1;
    Qt = zeros(N,N*(T+m));
    Rt = zeros(N,N*(T+m));
    Qt[.,1:N] = Qbar;
    Rt[.,1:N] = Qbar;
    liklhd = zeros(T+m,1);
    stdresid = zeros(m,N)|stdresid;
    hmat = zeros(m,N)|hmat;
    for j(2,T+m,1);
        Qt[.,(N*j)-1:(N*j)] = Qbar*(1-a-b) + a*stdresid[j-1,.]'stdresid[j-1,.] + b*Qt[.,(N*(j-1))-1:(N*(j-1))];
        Rt[.,(N*j)-1:(N*j)] = Qt[.,(N*j)-1:(N*j)] ./ ( sqrt(diag(Qt[.,(N*j)-1:(N*j)]))* ( sqrt(diag(Qt[.,(N*j)-1:(N*j)])))');
        liklhd[j] = -0.5*(N*log(2*pi) + 2*sumc(log(hmat[j,.])') +  log(det(Rt[.,(N*j)-1:(N*j)])) + stdresid[j,.]*inv(Rt[.,(N*j)-1:(N*j)])*stdresid[j,.]');
    endfor;
        Rt = Rt[.,3:T+m];
        Qt = Qt[.,3:T+m];
        liklhd = liklhd[2:T+m];
        likeval = sumc(liklhd);
    retp(Rt,Qt,liklhd,likeval);
endp;

{Rt,Qt,liklhd,likeval} = fconlike(dccpar,stdresid);
*/

//==================================================================================================//

// Procedure for simulating data from a DCC MODEL: ONLY WORKS FOR SPECIFIC DIMENSIONS !!!!
// Corrmat - the unconditional correlation matrix
// garchpar - 6 by 1 parameter vector for the 2-dimensional GARCH(1,1) process


proc(2) = simDCC(T,N,Corrmat,dccpar,garchpar);
    local a, b, Qbar, m, Qt, Rt, j;
    local rawdata, a01, alpha1, beta1, a02, alpha2, beta2, H, Ht, data, parvec1, parvec2;
    local stdresid;
    a = dccpar[1];
    b = dccpar[2];
    Qbar = Corrmat;
    T = T+500;
    m = 1;
    rawdata = rndn(T+m,N);
    stdresid = rndn(T+m,N);
    H = zeros(T+m,N);

    Qt = zeros(N,N*(T+m));
    Rt = zeros(N,N*(T+m));
    Ht = zeros(N,N*(T+m));

    Qt[.,1:N] = Qbar;
    Rt[.,1:N] = Qbar;
   
    for j(2,T+m,1);
        Qt[.,(N*j)-1:(N*j)] = Qbar*(1-a-b) + a*stdresid[j-1,.]'stdresid[j-1,.] + b*Qt[.,(N*(j-1))-1:(N*(j-1))];
        Rt[.,(N*j)-1:(N*j)] = Qt[.,(N*j)-1:(N*j)] ./ ( sqrt(diag(Qt[.,(N*j)-1:(N*j)]))* ( sqrt(diag(Qt[.,(N*j)-1:(N*j)])))');
        stdresid[j,.] = rawdata[j,.]*sqrt( Rt[.,(N*j)-1:(N*j)] );
    endfor;

        parvec1 = garchpar[1:3];
        a01 = parvec1[1]; alpha1=parvec1[2]; beta1=parvec1[3];

        parvec2 = garchpar[4:6];
        a02 = parvec2[1]; alpha2=parvec2[2]; beta2=parvec2[3];
        
        e1 = stdresid[.,1];
        e2 = stdresid[.,2];
        H[.,1] = garchsim(T+m,e1,a01,alpha1,beta1);
        H[.,2] = garchsim(T+m,e2,a02,alpha2,beta2);

        data = stdresid .* sqrt(H);
        for j(2,T+m,1);
            Ht[.,(N*j)-1:(N*j)] = sqrt(diagrv(zeros(N,N),H[j,.]'))* Rt[.,(N*j)-1:(N*j)]*sqrt(diagrv(zeros(N,N),H[j,.]'));
        endfor;
        
        Rt = Rt[.,N*(500+m)+1:N*(T+m)];
        Qt = Qt[.,N*(500+m)+1:N*(T+m)];
        Ht = Ht[.,N*(500+m)+1:N*(T+m)];

        data = data[500+m+1:T+m,.];
    retp(data,Ht);
endp;

//==================================================================================================//

proc garchsim(T,e,a0,a,b);
  local disc, h2,  i;
    /* first, check for validity of inputs */

    if (( a+b >= 1) or (a0 <= 0) );
    errorlog "Unable to generate GARCH(1,1) process with positive variance";
    end;
    endif;
    disc = 500;                         /* increase this if desired; do not
                                        change the value of ini below */
    h2 = zeros(T+disc,1);               /* initialize output vectors */
    h2[1] = a0*inv(1-a-b);             /* for first value, h2:= unconditional variance */
    e = sqrt(h2[1])*rndn(T+disc,1);
    for i(2,T+disc,1);
        h2[i] = a0 + a*h2[i-1] + b*e[i-1]^2;
        //e[i]  = sqrt(h2[i]) * rndn(1,1);
    endfor;

    h2 = h2[disc+1:T+disc];             /* discard disc startup values */
    retp(h2);
endp; 

//==================================================================================================//

// Procedure for simulating data from a CCC MODEL: ONLY WORKS FOR SPECIFIC DIMENSIONS !!!!
// Corrmat - the unconditional correlation matrix
// garchpar - 6 by 1 parameter vector for the 2-dimensional GARCH(1,1) process

proc(2) = simCCC(T,N,Corrmat,garchpar);
    local m, j, uncond1, uncond2, R, e1, e2, stdresid;
    local rawdata, a01, alpha1, beta1, a02, alpha2, beta2, H, Ht, data, parvec1, parvec2;

   // T = T + 500;
    m = 1;
    H = zeros(T+m,N);
    Ht = zeros(N,N*(T+m));
    parvec1 = garchpar[1:3];
    a01 = parvec1[1]; alpha1 = parvec1[2]; beta1 = parvec1[3];

    parvec2 = garchpar[4:6];
    a02 = parvec2[1]; alpha2 = parvec2[2]; beta2 = parvec2[3];
    
    rawdata = rndn(T+m,N);
    R = CorrMat;
    stdresid = rawdata*chol(R);
    
   // stdresid = rawdata*sqrt(R); // this overestimates rho!!

    uncond1 = sqrt(a01*inv(1-alpha1-beta1));
    uncond2 = sqrt(a02*inv(1-alpha2-beta2));
    
   // e1 = stdresid[.,1];
   // e2 = stdresid[.,2];
    H[.,1] = uncond1^2*ones(T+m,1);
    H[.,2] = uncond2^2*ones(T+m,1);

    e1 = uncond1*ones(T+m,1);
    e2 = uncond2*ones(T+m,1);

    for j(2,T+m,1);
        H[j,1] = a01 + alpha1*H[j-1,1] + beta1*e1[j-1]^2;
        e1[j]  = sqrt(H[j,1]) * stdresid[j,1];

        H[j,2] = a02 + alpha2*H[j-1,2] + beta2*e2[j-1]^2;
        e2[j]  = sqrt(H[j,2]) * stdresid[j,2];
    endfor;
    data = e1~e2;

 //   H[.,1] = garchsim(T+m,e1,a01,alpha1,beta1);
 //   H[.,2] = garchsim(T+m,e2,a02,alpha2,beta2);

   // data = stdresid .* sqrt(H);
    for j(2,T+m,1);
       Ht[.,(N*j)-1:(N*j)] = sqrt(diagrv(zeros(N,N),H[j,.]'))*R*sqrt(diagrv(zeros(N,N),H[j,.]'));
       
    endfor;
        
 //  Ht = Ht[.,N*(500+m)+1:N*(T+m)];

 //  data = data[500+m+1:T+m,.];

    Ht = Ht[.,N*m+1:N*(T+m)];

    data = data[m+1:T+m,.];

    retp(data,Ht);
endp;



//==================================================================================================//
// Procedure for simulating data from a BEKK(1,1) MODEL: ONLY WORKS FOR SPECIFIC DIMENSIONS !!!!
// Corrmat - the unconditional correlation matrix
// garchpar - 6 by 1 parameter vector for the 2-dimensional GARCH(1,1) process







//==================================================================================================//

// Procedure for simulating data from a VEC(1,1) MODEL: ONLY WORKS FOR SPECIFIC DIMENSIONS !!!!
// Corrmat - the unconditional correlation matrix
// garchpar - 6 by 1 parameter vector for the 2-dimensional GARCH(1,1) process







/*
//==================================================================================================//

// Estimating CCC GARCH Using Tsay's method in Chapter 10

proc clikelhd(B,data);
   local zz, hh, DD, j;
   local likecc;
    likecc = zeros(T,1); 
    zz = zeros(T,N);
    hh = zz;
    DD = zeros(T,N);
    zz[1,.] = Y0[1,.];
    hh[1,.] = (stdc(Y0).^2)';
    DD[1,.]  = sqrt(hh[1,.]);
    for j(2,T,1);
        hh[j,1] = B[1] + B[2]*hh[j-1,1] + B[3]*Y0[j-1,1]^2;
        hh[j,2] = B[4] + B[5]*hh[j-1,2] + B[6]*Y0[j-1,2]^2;
        DD[j,.]  = sqrt(hh[j,.]);
        zz[j,.] = Y0[j,.]./DD[j,.];
    endfor;
        for j(1,T,1);
        likecc[j] = -0.5*( sumc(log(hh[j,.])') + log(1-B[7]^2)  );
        likecc[j] = likecc[j] - 0.5*inv(1-B[7]^2)*( sumc((zz[j,.] .^2 )') - 2*B[7]*zz[j,1]*zz[j,2] );
        endfor;
    retp(likecc);
endp;


library cml,pgraph;
#include cml.ext;
#include pgraph.ext;
cmlset;
graphset;
B00=0.05|0.80|0.15|0.5|0.5|0.15|0.7;

@ Constraints @
//_cml_Bounds ={1e-8 1e256, 0 1e256, 0 1e256, 1e-8 1e256, 0 1e256, 0 1e256, 0 1e256};


  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(7,2).*_ww_;
  _cml_Bounds[2:3,1] = zeros(2,1);

_cml_Bounds[1,1] = 0.01;
_cml_Bounds[4,1] = 0.01;

_cml_c = zeros(4,7);
_cml_c[1,2:3] = -ones(1,2);
_cml_c[2,5:6] = -ones(1,2);
_cml_c[3,7]   = -ones(1,1);
_cml_c[4,7]   =  ones(1,1);
_cml_d = -0.99*ones(4,1);



_cml_Algorithm = 4;
_cml_DirTol = 1e-5;
_cml_MaxIters = 50;
format /rd 13,8;

{CC,f,g,hh,retcode}=cmlPrt(cml(Y0,0,&clikelhd,B00));
*/


//==================================================================================================//
/*

// Code for estimating Bivariate VC-MGARCH(1,1) of Tse and Tsui (2002) 
// using a variant of what is in Tsay's chapter 10

proc ttlikehd(B,data);
    local theta1, theta2, phi1, phi2, rhot;
    local zz, hh, DD, j;
    local likecc;
    likecc = zeros(T,1); 

    rhot = zeros(T,1);
    phi1 = zeros(T,1);
    phi2 = zeros(T,1);
    zz = zeros(T,N);
    hh = zz;
    DD = zeros(T,N);
    zz[1:2,.] = Y0[1:2,.];
    hh[1:2,.] = (stdc(Y0).^2)'|(stdc(Y0).^2)';
    DD[1:2,.]  = sqrt(hh[1:2,.]);
    for j(3,T,1);
        hh[j,1] = B[1] + B[2]*hh[j-1,1] + B[3]*Y0[j-1,1]^2;
        hh[j,2] = B[4] + B[5]*hh[j-1,2] + B[6]*Y0[j-1,2]^2;
        DD[j,.]  = sqrt(hh[j,.]);
        zz[j,.] = Y0[j,.]./DD[j,.];
    
       // theta1 = (B[8]^2)/(1+B[8]^2 + B[9]^2 );
       // theta2 = (B[9]^2)/(1+B[8]^2 + B[9]^2 );


        phi1[j] = zz[j-1,1]*zz[j-1,2] + zz[j-2,1]*zz[j-2,2];
        phi2[j] = sqrt(zz[j-1,1]^2 + zz[j-2,1]^2)*sqrt(zz[j-1,2]^2 + zz[j-2,2]^2);

       // rhot[j] = B[7] + theta1*rhot[j-1] + theta2*(phi1[j]/phi2[j]);

        rhot[j] = B[7] + ((B[8]^2)/(1+B[8]^2 + B[9]^2 ))*rhot[j-1];
        rhot[j] = rhot[j] + ((B[9]^2)/(1+B[8]^2 + B[9]^2 ))*(phi1[j]/phi2[j]);
    endfor;
        for j(1,T,1);
        likecc[j] = -0.5*( sumc(log(hh[j,.])') + log(1-rhot[j]^2) );
        likecc[j] = likecc[j] - 0.5*inv(1-rhot[j]^2)*( sumc((zz[j,.] .^2 )') - 2*rhot[j]*zz[j,1]*zz[j,2] );
        endfor;
    retp(likecc);
endp;
        

library cml,pgraph;
#include cml.ext;
#include pgraph.ext;
cmlset;
graphset;
B00=0.05|0.80|0.15|0.5|0.5|0.15|0.01|0|0;

@ Constraints @
//_cml_Bounds ={1e-8 1e256, 0 1e256, 0 1e256, 1e-8 1e256, 0 1e256, 0 1e256, 0 1e256};


  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(9,2).*_ww_;
  _cml_Bounds[2:3,1] = zeros(2,1);
  _cml_Bounds[5:6,1] = zeros(2,1);

_cml_Bounds[1,1] = 0.01;
_cml_Bounds[4,1] = 0.01;

_cml_c = zeros(4,9);
_cml_c[1,2:3] = -ones(1,2);
_cml_c[2,5:6] = -ones(1,2);
_cml_c[3,7]   = -ones(1,1);
_cml_c[4,7]   =  ones(1,1);
_cml_d = -0.99*ones(4,1);



_cml_Algorithm = 4;
_cml_DirTol = 1e-5;
_cml_MaxIters = 50;
format /rd 13,8;

//{tt,f,g,hh,retcode}=cmlPrt(cml(Y0,0,&ttlikehd,B00));

*/

//==================================================================================================//

// this procedure simulates time-varying correlation model of Tse and Tsui (2002)

proc(2) = simVCMG(T,N,Corrmat,dccpar,garchpar);
    local a, b, Qbar, m, Qt, Rt, j;
    local rawdata, a01, alpha1, beta1, a02, alpha2, beta2, H, Ht, data, parvec1, parvec2;
    local c, stdresid, e1, e2 ;
    a = dccpar[1];
    b = dccpar[2];
    c = dccpar[3];
    Qbar = Corrmat;
    T = T+500;
    m = 1;
    rawdata = rndn(T+m,N);
    stdresid = rndn(T+m,N);
    H = zeros(T+m,N);

    Qt = zeros(N,N*(T+m));
    Rt = zeros(N,N*(T+m));
    Ht = zeros(N,N*(T+m));

    Qt[.,1:N] = Qbar;
    Rt[.,1:N] = Qbar;
   
    for j(2,T+m,1);
        Qt[.,(N*j)-1:(N*j)] = Qbar*(1-a-b) + a*stdresid[j-1,.]'stdresid[j-1,.] + b*Qt[.,(N*(j-1))-1:(N*(j-1))];
        Rt[.,(N*j)-1:(N*j)] = Qt[.,(N*j)-1:(N*j)] ./ ( sqrt(diag(Qt[.,(N*j)-1:(N*j)]))* ( sqrt(diag(Qt[.,(N*j)-1:(N*j)])))');
        stdresid[j,.] = rawdata[j,.]*sqrt( Rt[.,(N*j)-1:(N*j)] );
    endfor;

        parvec1 = garchpar[1:3];
        a01 = parvec1[1]; alpha1=parvec1[2]; beta1=parvec1[3];

        parvec2 = garchpar[4:6];
        a02 = parvec2[1]; alpha2=parvec2[2]; beta2=parvec2[3];
        
        e1 = stdresid[.,1];
        e2 = stdresid[.,2];
        H[.,1] = garchsim(T+m,e1,a01,alpha1,beta1);
        H[.,2] = garchsim(T+m,e2,a02,alpha2,beta2);

        data = stdresid .* sqrt(H);
        for j(2,T+m,1);
            Ht[.,(N*j)-1:(N*j)] = sqrt(diagrv(zeros(N,N),H[j,.]'))* Rt[.,(N*j)-1:(N*j)]*sqrt(diagrv(zeros(N,N),H[j,.]'));
        endfor;
        
        Rt = Rt[.,N*(500+m)+1:N*(T+m)];
        Qt = Qt[.,N*(500+m)+1:N*(T+m)];
        Ht = Ht[.,N*(500+m)+1:N*(T+m)];

        data = data[500+m+1:T+m,.];
    retp(data,Ht);
endp;



//==================================================================================================//

// this procedure simulates data from a Bivariate Full BEKK model
// bekkpar - is a [ N(N+1)/2 + N^2 + N^2 ] by 1 parameter vector obtained by vecr
// H_{t} = C*C' + A*(e_{t-1}* e_{t-1}')* A' + B*H_{t-1}*B'

proc simFBEKK(T,N,bekkpar);
    local nn, C, A, B, mat, Ht, data, Cvec, Uncond, const;
    T = T + 500;
    nn = 0.5*N*(N+1);
    C = bekkpar[1:nn];
    A = bekkpar[nn+1:N^2+nn];
    B = bekkpar[nn+N^2+1:nn+2*N^2];
    C = xpnd(C);
    C = lowmat(C);
    Cvec = vec(C*C');
    //Cvec = (C.*.C)'vec(eye(N));
    A = reshape(A,N,N);
    B = reshape(B,N,N);
    mat = eye(N^2);
    mat = mat - A .*.A - B.*.B;
    Ht = zeros(N,N*T);
   
    Uncond = inv(mat)*Cvec;
    Uncond  = reshape(Uncond,N,N);
    const = C*C';
    data = rndn(T,N)*chol(Uncond);
    Ht[.,1:2] = const;
    for i(2,T,1);
        Ht[.,(N*i)-1:N*i] = const + A*(data[i-1,.]'data[i-1,.])*A' + B*Ht[.,(N*(i-1))-1:N*(i-1)]*B';
        data[i,.] = rndn(1,N)*chol(Ht[.,(N*i)-1:N*i]);
    endfor;
    retp(data[500+1:T,.]);

endp;

//==================================================================================================//

// 

proc bekklike(bekkpar,e);
    local nn, C, A, B, mat, Ht, data, Cvec, Uncond, const;
    local llike, i;
    nn = 0.5*N*(N+1);
    C = bekkpar[1:nn];
    A = bekkpar[nn+1:N^2+nn];
    B = bekkpar[nn+N^2+1:nn+2*N^2];
    C = xpnd(C);
    C = lowmat(C);
    const = C*C';
    Cvec = vec(C*C');
    A = reshape(A,N,N);
    B = reshape(B,N,N);
    Uncond = corrvc(e);
    Ht = zeros(N,N*(T+1));
    Ht[.,1:2] = Uncond;
    e = sqrt(diag(Uncond))'|e;
    llike = zeros(T+1,1);
    for i(2,T+1,1);
        Ht[.,(N*i)-1:N*i] = const + A*(e[i-1,.]'e[i-1,.])*A' + B*Ht[.,(N*(i-1))-1:N*(i-1)]*B';
        llike[i] = -0.5*(N*log(2*pi) + log(det(Ht[.,(N*i)-1:N*i])) + e[i,.]*inv(Ht[.,(N*i)-1:N*i])*e[i,.]');
    endfor;
    retp(llike[2:T+1]);
endp;

//==================================================================================================//


// this procedure gives the MLE for a scalar BEKK. These MLE's are used as starting values for
// estimating the Complete BEKK model

sbekkpar = scaleBEK(dgp);

proc sbeklike(sbekkpar,e);
    local nn, C, A, B, mat, Ht, data, Uncond, const;
    local llike, i;
    nn = 0.5*N*(N+1);
    C = sbekkpar[1:nn];
    A = sbekkpar[nn+1];  // A is a scalar
    B = sbekkpar[nn+2]; // A is a scalar
    C = xpnd(C);
    C = lowmat(C);
    const = C*C';
    Uncond = corrvc(e);
    Ht = zeros(N,N*(T+1));
    Ht[.,1:2] = Uncond;
    e = sqrt(diag(Uncond))'|e;
    llike = zeros(T+1,1);
    for i(2,T+1,1);
        Ht[.,(N*i)-1:N*i] = const + A*(e[i-1,.]'e[i-1,.])*A' + B*Ht[.,(N*(i-1))-1:N*(i-1)]*B';
        llike[i] = -0.5*(N*log(2*pi) + log(det(Ht[.,(N*i)-1:N*i])) + e[i,.]*inv(Ht[.,(N*i)-1:N*i])*e[i,.]');
    endfor;
    retp(llike[2:T+1]);
endp;


// this procedure yields the starting values for the estimation of the 
// Scalar BEKK model using the observed data
    B00=0.05|0.80|0.15;
    


proc scaleBEK(data);
    local startg, i, nn, A, B, C, Cstar, CC, beginSB;
    nn = 0.5*N*(N+1);
    startg = zeros(N,3);  // this stores the univariate GARCH(1,1) parameters
    for i(1,N,1);
    //di = ;
    {B1,f1,g1,hh1,retcode}=cmlPrt(cml(dgp[.,i],0,&ulikelhd,B00));
    startg[i,.] = B1';
    endfor;
    A = meanc(startg[.,2]);
    B = meanc(startg[.,3]);
    C = corrvc(data);
    Cstar = C*(1-A^2-B^2);
    CC = chol(Cstar)';
    beginSB = vech(CC)|A|B;
    retp(beginSB);
endp;



    
// now use the starting values ``beginSB" to maximize the likelihood function
// of the scalar BEKK model
    {B1,f1,g1,hh1,retcode}=cmlPrt(cml(dgp,0,&sbeklike,sbekkpar));
    scalepar = B1;

// this procedure converts the output, the MLE, from sbeklike to a form 
// that is conformable for use as the starting values in the Complete BEKK
// estimation

proc convert(scalepar);
    local nn, Cstar, A, B, bekkpar;
    nn = 0.5*N*(N+1);
    Cstar = scalepar[1:nn];
    A = eye(N)*scalepar[nn+1];
    B = eye(N)*scalepar[nn+2];
    A = reshape(A,N^2,1);
    B = reshape(B,N^2,1);
    bekkpar = Cstar|A|B;
    retp(bekkpar);
endp;
   
 //   {B1,f1,g1,hh1,retcode}=cmlPrt(cml(e,0,&bekklike,bekkpar));






end;