1 Star 0 Fork 0

ymhong / Computer Codes

Create your Gitee Account
Explore and code with more than 8 million developers,Free private repositories !:)
Sign up
Clone or Download
ccc-dcc 25.98 KB
Copy Edit Web IDE Raw Blame History
ymhong authored 2021-04-23 09:37 . paper on DCC April 2007
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973
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;
R
1
https://gitee.com/ymhong/code.git
git@gitee.com:ymhong/code.git
ymhong
code
Computer Codes
master

Search