/* ** Estimates the Model 4 of HECER Discussion Paper No 247 / December 2008 ** ** Computes the residuls of the estimated M-GARCH model ** - Mixtures of normal distribution ** - Two regimes with one factor in the first regime and two ** factors in the second regime ** */ library cmlmt,LK_QRtest_multivariate; #include cmlmt.sdf gausset; screen on; /* n - dimension of the data r - number of factors */ n=4; r=2; /* matK1 - Lag numbers of the autocorrelation test ** matK2 - Lag numbers of the heteroscedasticity test ** Tboots- Size of the simulated sample used in computing the covariance ** matrix estimate in the test statistic */ let matK1=1 2 3 4 5 6 7 8 9 10 20 30; let matK2=1 2 3 4 5 6 7 8 9 10 20 30; let Tboots=20000; /* The data is loaded and modified (first 4*52 rows removed & centered) */ load fxweekr; data=fxweekr[.,5]~fxweekr[.,4]~fxweekr[.,6]~fxweekr[.,8]; data=trimr(data,4*52,0); data=data-meanc(data)'; struct DS z0; z0 = dsCreate; z0.dataMatrix=data; /* Initial values of parameters are given here ** "alpha11"|"beta11" |"alpha21"|"alpha22"|"beta21"|"beta22" |"b11"|"b21"|"b31"|"b41"|"b12"|"b22"|"b32"|"b42" |"b13"|"b23"|"b33"|"b43"|"b14"|"b24"|"b34"|"b44" |"psii1"|"psii2"|"psii3"|"psii4"|"p" */ struct PV x0; x0 = pvCreate; x0 = pvPacki(x0,0.0202|0.9789,"garch1",1); x0 = pvPacki(x0,0.2726|0.1684|0.6441|0.8239,"garch2",2); x0 = pvPacki(x0,0.1901|-8.9138|8.7480|-0.0924| 2.7383|-1.1200|-1.6036|0.0068| 0.7186|1.1959|-0.5037|-0.8058| 0.0528|-1.8266|0.5527|1.4895,"B",3); x0 = pvPacki(x0,0.1106|0.1390|0.2005|0.3879,"psii",4); x0 = pvPacki(x0,0.0588,"p",5); initialvalues = pvGetParVector(x0); /* Setting variable values for ESTIMATION */ struct cmlmtControl ctl; ctl = cmlmtControlCreate; ctl.title ="Model 4"; ctl.maxIters = 10000; ctl.Algorithm=1; /* = 0, Modifed BFGS ** = 1, BFGS (Broyden, Fletcher, Goldfarb, Shanno) (default) ** = 2, DFP (Davidon, Fletcher, Powell) ** = 3, NEWTON (Newton-Raphson) ** = 4, BHHH */ ctl.linesearch=2; /* = 0, augmented trust region method (requires constraints) ** = 1, STEPBT (quadratic and cubic curve ¯t) (default) ** = 2, Brent's method ** = 3, BHHHStep ** = 4, half ** = 5, Wolfe's condition */ ctl.CovParType=1; /* = 0, no covariance matrix is computed ** = 1, ML covariance matrix is computed (Default = 1) ** = 2, QML covariance matrix, */ ctl.GradMethod=1; /* = 0, central difference ** = 1, forward difference (default) ** = 2, backward difference */ ctl.DirTol=1e-5; /* PARAMETERS ARE: "alpha11"|"beta11" |"alpha21"|"alpha22"|"beta21"|"beta22" |"b11"|"b21"|"b31"|"b41"|"b12"|"b22"|"b32"|"b42" |"b13"|"b23"|"b33"|"b43"|"b14"|"b24"|"b34"|"b44" |"psii1"|"psii2"|"psii3"|"psii4"|"p" ** garch coefficients positive and sum <1 ** p restricted between 0.01 and 0.5 */ ctl.C = {1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0, 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0, 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0}; ctl.D = {1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,0.01,-0.5,1e-5,1e-5,1e-5,1e-5, -1+1e-5,-1+1e-5,-1+1e-5}; struct cmlmtResults out1; out1 = CMLmt(&ll,x0,z0,ctl); /* CALCULATION OF QUANTILE RESIDUALS */ T=rows(z0.dataMatrix); apu1=fct1(out1.Par,z0); k1=apu1[(T+1):(2*T),.]~apu1[(2*T+1):(3*T),.]~ apu1[(3*T+1):(4*T),.]~apu1[(4*T+1):(5*T),.]; jk=apu1[(5*T+1):(6*T),.]; /* NUMERICAL DERIVATIVES & QUANTILE RESIDUALS OF SIMULATED DATA */ struct DS z0boots; z0boots = dsCreate; z0boots.dataMatrix=datasimuModel4(out1.Par,Tboots); apu3 = fct1(out1.Par,z0boots); k1boots = apu3[(Tboots+1):(2*Tboots),.]~apu3[(2*Tboots+1):(3*Tboots),.]~ apu3[(3*Tboots+1):(4*Tboots),.]~apu3[(4*Tboots+1):(5*Tboots),.]; jkboots= apu3[(5*Tboots+1):(6*Tboots),.]; /*The numerical derivatives of the simulated sample ** gboots is the derivatives of quantile residuals ** sboots is the score of the simulated sample */ apu4 = gradMT(&fct1,out1.Par,z0boots); gboots = apu4[(Tboots+1):(2*Tboots),.]~apu4[(2*Tboots+1):(3*Tboots),.]~ apu4[(3*Tboots+1):(4*Tboots),.]~apu4[(4*Tboots+1):(5*Tboots),.]; jgboots = apu4[(5*Tboots+1):(6*Tboots),.]; sboots = apu4[1:Tboots,.]; /* TESTS ** The p - values are printed to this file */ outwidth 80; output file = FXweekr_withModel4.asc on; screen off; call cmlmtPrt(out1); print ""; print "Initial values used in estimation"; print pvGetParNames(x0)'; print initialvalues'; print "The value of the likelihood at minimum " ftos(out1.Fct*T,"%*.*lf",1,3); print ""; print ""; /* ** See LK_QRtest_multivariate.src for the proc's "prtesting" & "multiprtesting" */ print "Tests based on multivariate quantile residuals"; multiprtesting(k1,k1boots,sboots,gboots,matK1,matK2); print ""; print ""; print "Tests based on joint quantile residuals"; prtesting(jk,jkboots,sboots,jgboots,matK1,matK2); format /mb1 /ros 16,8; outwidth 256; output off; screen on; /****************************************************************************/ /* THE NEEDED PROCEDURE DEFINITIONS USED BY THE MAIN FILE */ /****************************************************************************/ /* Computes the log-likelihood of the model */ proc ll(struct PV x0, struct DS z0,ind); local garch1, garch2, B, psii, p; struct modelResults mm; garch1 = pvUnpack(x0,1); garch2 = pvUnpack(x0,2); B = pvUnpack(x0,3); psii = pvUnpack(x0,4); p = pvUnpack(x0,5); local T,B1,B11,psi,Psi1,Psi2,f,v1,v2,H1,H2; T =rows(z0[1].datamatrix); B=reshape(B,n,n)'; /*B on nxn matrix , the first n elements are put into the first row of B*/ B1=B[.,1:r]; B11=B[.,1]; psi=diagrv(eye(n),psii); Psi1=p*eye(n)+(1-p)*psi; Psi2=Psi1*inv(psi); f=zeros(T,1); /*vcx() computes a variance-covariance matrix. In this case rxr -matrix diag() creates a column vector from the diagonal of a matrix. v1 & v2 are T x r -matrices, contain the initial values. */ v1=diag(vcx(z0.dataMatrix*B11))'.*ones(T,1); v2=diag(vcx(z0.dataMatrix*B1))'.*ones(T,r); for i (2,T,1); v1[i,1] = (1-garch1[1]-garch1[2])+ garch1[2]*v1[i-1,1] + garch1[1]*(B1[.,1]'z0.dataMatrix[i-1,.]')^2; for j (1,r,1); v2[i,j] = (1-garch2[j]-garch2[2+j])+ garch2[2+j]*v2[i-1,j] + garch2[j]*(B1[.,j]'z0.dataMatrix[i-1,.]')^2; endfor; /*diagrv(a,b) is a matrix, b is a column vector */ H1 = diagrv(eye(n),(v1[i,.]~ones(1,n-r+1))'); H2 = diagrv(eye(n),(v2[i,.]~ones(1,n-r))'); f[i] = -0.5*n*ln(2*pi) + ln(det(B)) + ln(p*sqrt(det(inv(H1)*Psi1))* exp(-0.5*z0.dataMatrix[i,.]*B*Psi1*inv(H1)*B'z0.dataMatrix[i,.]') + (1-p)*sqrt(det(inv(H2)*Psi2))* exp(-0.5*z0.dataMatrix[i,.]*B*Psi2*inv(H2)*B'z0.dataMatrix[i,.]')); endfor; f=trimr(f,1,0); mm.Function = f; retp(mm); endp; /* This procedure calculates the gradient for the likelihood function ** and the gradient for quantile residuals ** Returns a (1+n*T+1) x rows(parameters) -matrix */ proc fct1(struct PV x0, struct DS z0); local n, r; n=4; r=2; local garch1, garch2, B, psii, p; struct modelResults mm; garch1 = pvUnpack(x0,1); garch2 = pvUnpack(x0,2); B = pvUnpack(x0,3); psii = pvUnpack(x0,4); p = pvUnpack(x0,5); local T,B1,B11,psi,Psi1,Psi2,f,v1,v2; T =rows(z0.datamatrix); B=reshape(B,n,n)'; /*B on nxn matrix , the first n elements are put into the first row of B*/ B1=B[.,1:r]; B11=B[.,1]; psi=diagrv(eye(n),psii); Psi1=p*eye(n)+(1-p)*psi; Psi2=Psi1*inv(psi); /*vcx() computes a variance-covariance matrix. In this case rxr -matrix diag() creates a column vector from the diagonal of a matrix. v1 & v2 are T x r -matrices, contain the initial values. */ v1=diag(vcx(z0.dataMatrix*B11))'.*ones(T,1); v2=diag(vcx(z0.dataMatrix*B1))'.*ones(T,r); for i (2,T,1); v1[i,1] = (1-garch1[1]-garch1[2])+ garch1[2]*v1[i-1,1] + garch1[1]*(B1[.,1]'z0.dataMatrix[i-1,.]')^2; for j (1,r,1); v2[i,j] = (1-garch2[j]-garch2[2+j])+ garch2[2+j]*v2[i-1,j] + garch2[j]*(B1[.,j]'z0.dataMatrix[i-1,.]')^2; endfor; endfor; local mu, sigma, nu, omega, H1t, H2t, d,y,y1,y2,y3,apu,be; mu=zeros(n,1); sigma=zeros(n,n); nu=zeros(n,1); omega=zeros(n,n); H1t=zeros(n,n); H2t=zeros(n,n); d=ones(1,n); y1=zeros(T,1); y2=zeros(T,n); y3=zeros(T,1); for tt (1,T,1); H1t=diagrv(eye(n),(v1[tt,1]~ones(1,n-1))); H2t=diagrv(eye(n),(v2[tt,1:r]~ones(1,n-r))); sigma=inv(B')*H1t*invpd(Psi1)*inv(B); omega=inv(B')*H2t*invpd(Psi2)*inv(B); d=z0.dataMatrix[tt,.]'; y1[tt,1]=ln(mixdistr(d,mu,sigma,nu,omega,p)); y2[tt,.]=(resid(&mixdistr,&cumulmixdistr,d,mu,sigma,nu,omega,p))'; y3[tt,1]=cdfn(y2[tt,1]); for i (2,n,1); y3[tt,1]=y3[tt,1]*cdfn(y2[tt,i]); endfor; apu=y3[tt,1]; be=0; for j (0,n-1,1); be=be+((-1)^j)*(ln(apu)^j)/(gamma(j+1)); endfor; y3[tt,1]=cdfni(apu*be); endfor; y=y1|vecr(y2')|y3; retp(y); endp; proc datasimuModel4(struct PV p1, T); local n, r; n=4; r=2; local garch1, garch2, B, psii, p; struct modelResults mm; garch1 = pvUnpack(x0,1); garch2 = pvUnpack(x0,2); B = pvUnpack(x0,3); psii = pvUnpack(x0,4); p = pvUnpack(x0,5); B=reshape(B,n,n)'; /*B on nxn matrix , the first n elements are put into the first row of B*/ local B1, B11, psi, psi1, psi2; B1=B[.,1:r]; B11=B[.,1]; psi=diagrv(eye(n),psii); Psi1=p*eye(n)+(1-p)*psi; Psi2=Psi1*inv(psi); local epsilon,data,bool,st,st1,st2,raja; epsilon=rndn(200+T,n); data=zeros(200+T,n); bool=rndn(200+T,1); st=zeros(200+T,1); raja=cdfni(p); for i (1,200+T,1); if bool[i]>raja; st[i]=1; endif; endfor; st1=zeros(200+T,1); st2=st; for i (1,200+T,1); if st[i]==0; st1[i]=1; endif; endfor; local v1,v2,H1t,H2t; /* Initial values for vi's */ v1=ones(200+T,1); v2=ones(200+T,r); H1t=zeros(n,n); H2t=zeros(n,n); data[1,.]=epsilon[1,.]; /*v -matrices */ for tt (2,200+T,1); v1[tt,1] = (1-garch1[1]-garch1[2]) + garch1[2]*v1[tt-1,1] + garch1[1]*(B1[.,1]'data[tt-1,.]')^2; for i (1,r,1); v2[tt,i] = (1-garch2[i]-garch2[i+2]) + garch2[i+2]*v2[tt-1,i] + garch2[i]*(B1[.,i]'data[tt-1,.]')^2; endfor; H1t=diagrv(eye(n),(v1[tt,1]~ones(1,n-1))); H2t=diagrv(eye(n),(v2[tt,1:r]~ones(1,n-r))); data[tt,.]=(inv(B')*(st1[tt]*sqrt(H1t)*sqrt(invpd(Psi1))+ st2[tt]*sqrt(H2t)*sqrt(invpd(Psi2)))*epsilon[tt,.]')'; endfor; retp(data[201:200+T,.]); endp; /* ************************************************************************ */ /* >mixdistr ** Purpose: computes the joint, marginal or conditional density function of ** mixture model ** ** Format: e = mixdistr(data,m1,s1,m2,s2,p); ** ** Input: data kx1, observation at time point t ** m1 kx1, corresponding 1. mean vector of the model ** s1 kxk, corresponding 1. covariance matrix of the model, ** must be symmetric & positive definite ** m2 kx1, corresponding 2. mean vector of the model ** s2 kxk, corresponding 2. covariance matrix of the model, ** must be symmetric & positive definite ** p 1x1, value of the corresponding mixing proportion ** of the 1. component, takes values in [0,1] ** ** Output: e 1x1 value of the joint, marginal or conditional density function ** */ proc mixdistr(data,m1,s1,m2,s2,p); local xx1,xx2, k, y; xx1=(data-m1)'*invpd(s1)*(data-m1); xx2=(data-m2)'*invpd(s2)*(data-m2); k=rows(data); y= p*det(s1)^(-0.5)*(2*pi)^(-k/2)*exp(-0.5*xx1) +(1-p)*det(s2)^(-0.5)*(2*pi)^(-k/2)*exp(-0.5*xx2); retp (y); endp; /* >cumulmixdistr ** Purpose: computes the cumulative univariate distribution ** function of mixture model ** ** Format: e = cumulmixdistr(data,m1,s1,m2,s2,p); ** ** Input: data 1x1, observation at time point t ** m1 1x1, corresponding 1. mean of the model ** s1 1x1, corresponding 1. variance of the model, ** must be positive ** m2 1x1, corresponding 2. mean of the model ** s2 1x1, corresponding 2. variance of the model, ** must be positive ** p 1x1, value of the corresponding mixing proportion ** of the 1. component, takes values in [0,1] ** ** Output: e 1x1 value of the cumulative univariate distribution function ** */ proc cumulmixdistr(data,m1,s1,m2,s2,p); local x1,x2; x1=(data-m1)/sqrt(s1); x2=(data-m2)/sqrt(s2); retp (p*cdfn(x1)+(1-p)*cdfn(x2)); endp; /* >resid ** Purpose: computes the vector of quantile residuals ** the order of conditioning is: ** 1. marginal - 1. component ** 2. conditional marginal - 2. component ** ... ** n. conditional marginal - n. component ** ** Format: e = resid(&g,&cumulg,data,mu,sigma,nu,omega,p); ** ** Input: data kx1 matrix of the observed dataset ** g scalar, pointer to a procedure returning ** 1x1 scalar, the value of the (joint) ** density function for one observation ** procedure call uses/needs parameters ** (data,mu,sigma,nu,omega,p) ** cumulg scalar, pointer to a procedure returning ** 1x1 scalar, the value of the (joint) ** cumulative density function for one observation ** procedure call uses/needs parameters ** (data,mu,sigma,nu,omega,p) ** mu kx1 vector, corresponding mean vector of the 1. component ** sigma kxk matrix, corresponding covariance matrix of the 1. component ** nu kx1 vector, corresponding mean vector of the 2. component ** omega kxk matrix, corresponding covariance matrix of the 2. component ** p 1x1, value of the corresponding mixing proportion ** of the 1. component, takes values in (0,1) ** ** Output: e kx1 vector of quantile residuals ** */ proc resid(&g,&cumulg,data,mu,sigma,nu,omega,p); local g:proc, cumulg:proc; local eps; eps=10^-10; if rows(data)<=0; retp (error(0)); elseif rows(data)==1; if cumulg(data,mu,sigma,nu,omega,p)>=1; retp (1-eps); else; if cumulg(data,mu,sigma,nu,omega,p)<=0; retp (0+eps); else; retp (cdfni(cumulg(data,mu,sigma,nu,omega,p))); endif; endif; else; local x1,x2,mu1,mu2,sigma1,sigma2,nu1,nu2,omega1,omega2,p2,e; x1=data[1,1]; x2=data[2:rows(data),1]; mu1=mu[1,1]; nu1=nu[1,1]; sigma1=sigma[1,1]; omega1=omega[1,1]; mu2=mu[2:rows(mu),1]+sigma[2:rows(sigma),1]*(x1-mu1)/sigma1; nu2=nu[2:rows(nu),1]+omega[2:rows(omega),1]*(x1-nu1)/omega1; sigma2=sigma[2:rows(sigma),2:cols(sigma)] -sigma[2:rows(sigma),1]*sigma[1,2:cols(sigma)]/sigma1; omega2=omega[2:rows(omega),2:cols(omega)] -omega[2:rows(omega),1]*omega[1,2:cols(omega)]/omega1; p2=p*g(x1,mu1,sigma1,nu1,omega1,1)/g(x1,mu1,sigma1,nu1,omega1,p); e=zeros(rows(data),1); eps=10^-10; if cumulg(x1,mu1,sigma1,nu1,omega1,p)>=1; e[1,1]=(1-eps); else; if cumulg(x1,mu1,sigma1,nu1,omega1,p)<=0; e[1,1]=(0+eps); else; e[1,1]=cdfni(cumulg(x1,mu1,sigma1,nu1,omega1,p)); endif; endif; e[2:rows(data),1]=resid(&g,&cumulg,x2,mu2,sigma2,nu2,omega2,p2); retp (e); endif; endp;