/* 23.05.2011 10:00 ** estimates the MAR(3,4,1)+GARCH(1,1) model with 3 regimes ** calculates quantile residuals ** calculates test statistics using them ** estimator of omega is done using a simulated sample */ library cmlmt,LK_QRtest_univariate; #include cmlmt.sdf /* Globals are set to default values */ gausset; screen on; /* 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 */ load z[347,1]=DEM1.dat; /* 347 obsesrvations are used for estimation */ maxdat=-maxc(z); mindat=minc(z); /* These are needed in constraint estimaton! */ struct DS z0; z0=reshape(dsCreate,6,1); z0[1].dataMatrix=z[6:347]; z0[2].dataMatrix=z[5:346]; z0[3].dataMatrix=z[4:345]; z0[4].dataMatrix=z[3:344]; z0[5].dataMatrix=z[2:343]; z0[6].dataMatrix=z[1:342]; /* Initial values of parameters are given here ** "mu1"|"c1"|"mu2"|"c2"|"mu3"|"s1"|"s2"|"s3"|"se"|"b1"|"b2"|"b3"|"b4"| ** "alfa"|"beta" */ struct PV x0; x0 = pvCreate; x0 = pvPacki(x0,0.205|6.929|0.570|12.158|0.628,"mean",1); x0 = pvPacki(x0,0.010|0.004|0.950|6.957,"s",2); x0 = pvPacki(x0,1.005|-0.08|0.188|-0.171,"b",3); x0 = pvPacki(x0,0.451|0.299,"garch",4); initialvalues = pvGetParVector(x0); /* Setting variable values for ESTIMATION */ struct cmlmtControl ctl; ctl = cmlmtControlCreate; ctl.title ="MAR(3,4,1)-GARCH(1,1)-model__estimated"; ctl.maxIters = 10000; ctl.Algorithm=4; /* = 0, Modifed BFGS ** = 1, BFGS (Broyden, Fletcher, Goldfarb, Shanno) (default) ** = 2, DFP (Davidon, Fletcher, Powell) ** = 3, NEWTON (Newton-Raphson) ** = 4, BHHH */ ctl.linesearch=1; /* = 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: "mu1"|"c1"|"mu2"|"c2"|"mu3"|"s1"|"s2"|"s3"|"se" ** |"b1"|"b2"|"b3"|"b4"|"alfa"|"beta" ** The standard errors are limited to be strictly positive (4 rows) ** Thresholds are limited from below and above by the data (4 rows) ** c1 is the smaller threshold value (1 row) */ /* ctl.C ={0 0 0 0 0 1 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 1 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0, 0 -1 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 1 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 -1 0 1 0 0 0 0 0 0 0 0 0 0 0}; ctl.D = {1e-5,1e-5,1e-5,1e-5,-maxdat,-maxdat,mindat,mindat,1e-5}; */ struct cmlmtResults out1; out1 = CMLmt(&ll,x0,z0,ctl); /* CALCULATION OF QUANTILE RESIDUALS */ T=rows(z0[1].dataMatrix); apu1=fct1(out1.Par,z0); r=apu1[(T+1):2*T,.]; /* NUMERICAL DERIVATIVES & QUANTILE RESIDUALS OF SIMULATED DATA */ struct DS z0boots; z0boots=reshape(dsCreate,6,1); apu2=datasimuMAR431_GARCH11(out1.Par,Tboots+5); z0boots[1].dataMatrix=apu2[6:(Tboots+5)]; z0boots[2].dataMatrix=apu2[5:(Tboots+4)]; z0boots[3].dataMatrix=apu2[4:(Tboots+3)]; z0boots[4].dataMatrix=apu2[3:(Tboots+2)]; z0boots[5].dataMatrix=apu2[2:(Tboots+1)]; z0boots[6].dataMatrix=apu2[1:Tboots]; apu3 = fct1(out1.Par,z0boots); rboots = apu3[(T+1):2*T,.]; /*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[(T+1):2*T,.]; sboots = apu4[1:T,.]; /* TESTS ** The p - values are printed to this file */ outwidth 80; output file = DEM_withMAR341GARCH11.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.src for the proc "prtesting" */ prtesting(r,rboots,sboots,gboots,matK1,matK2); print ""; format /mb1 /ros 16,8; outwidth 256; output off; screen on; /****************************************************************************/ /* THE NEEDED PROCEDURE DEFINITIONS USED BY THE MAIN FILE */ /****************************************************************************/ /* Log-Likelihood function for a MAR341GARCH11 model */ proc ll(struct PV x0, struct DS z0, ind); local mean, s, b, garch; struct modelResults mm; mean = pvUnpack(x0,1); s = pvUnpack(x0,2); b = pvUnpack(x0,3); garch = pvUnpack(x0,4); local u1,u2,u3,lag1,f,ar1,pii1,pii2,T,ss,ss1,ss2,ss3; T=rows(z0[1].dataMatrix); ar1=b[1].*z0[2].dataMatrix + b[2].*z0[3].dataMatrix+ b[3].*z0[4].dataMatrix + b[4].*z0[5].dataMatrix; lag1=b[1]*z0[3].dataMatrix +b[2]*z0[4].dataMatrix+ b[3]*z0[5].dataMatrix+b[4]*z0[6].dataMatrix; u1=z0[2].dataMatrix-lag1-mean[1]; u2=z0[2].dataMatrix-lag1-mean[3]; u3=z0[2].dataMatrix-lag1-mean[5]; pii1=cdfn((z0[2].dataMatrix-mean[2])/sqrt(s[4])); pii2=cdfn((z0[2].dataMatrix-mean[4])/sqrt(s[4])); /* Computation of the standard errors */ ss=zeros(T+1,3); /* Starting value for the standard errors caluculated using the data */ ss[1,1]=meanc(((z0[1].dataMatrix-mean[1]-ar1).*(1-pii1)+ (z0[1].dataMatrix-mean[3]-ar1).*(pii1-pii2)+ (z0[1].dataMatrix-mean[5]-ar1).*pii2)^2); ss[1,2]=ss[1,1]; ss[1,3]=ss[1,1]; for i(1,T,1); ss[i+1,1]=s[1]+garch[1]*u1[i]^2+garch[2]*ss[i,1]; ss[i+1,2]=s[2]+garch[1]*u2[i]^2+garch[2]*ss[i,2]; ss[i+1,3]=s[3]+garch[1]*u3[i]^2+garch[2]*ss[i,3]; endfor; /* ss1 - ss3 are the standard error vectors ** Initial value is removed from the vector of standard errors */ ss1=ss[2:T+1,1]; ss2=ss[2:T+1,2]; ss3=ss[2:T+1,3]; f=((pdfn((z0[1].dataMatrix-mean[1]-ar1)./sqrt(ss1)))./sqrt(ss1)).*(1-pii1) + ((pdfn((z0[1].dataMatrix-mean[3]-ar1)./sqrt(ss2)))./sqrt(ss2)).*(pii1-pii2) + ((pdfn((z0[1].dataMatrix-mean[5]-ar1)./sqrt(ss3)))./sqrt(ss3)).*pii2; mm.Function = ln(f); retp(mm); endp; /* This procedure calculates the gradient for the likelihood function ** and the gradient for quantile residuals ** Returns a (2*T) x rows(parameters) -matrix */ proc fct1(struct PV x0, struct DS z0); local mean,s,b,garch; mean = pvUnPack(x0,1); s = pvUnPack(x0,2); b = pvUnPack(x0,3); garch = pvUnPack(x0,4); local u1,u2,u3,lag1,f,ar1,pii1,pii2,T; local ss,ss1,ss2,ss3,d,yy; ar1=b[1]*z0[2].dataMatrix+b[2]*z0[3].dataMatrix+ b[3]*z0[4].dataMatrix+b[4]*z0[5].dataMatrix; lag1=b[1]*z0[3].dataMatrix+b[2]*z0[4].dataMatrix+ b[3]*z0[5].dataMatrix+b[4]*z0[6].dataMatrix; u1=z0[2].dataMatrix-lag1-mean[1]; u2=z0[2].dataMatrix-lag1-mean[3]; u3=z0[2].dataMatrix-lag1-mean[5]; pii1=cdfn((z0[2].dataMatrix-mean[2])/sqrt(s[4])); pii2=cdfn((z0[2].dataMatrix-mean[4])/sqrt(s[4])); /* Computation of the standard errors */ T=rows(z0[1].dataMatrix); ss=zeros(T+1,3); /* Starting value for the standard errors caluculated using the data */ ss[1,1]=meanc(((z0[1].dataMatrix-mean[1]-ar1).*(1-pii1)+ (z0[1].dataMatrix-mean[3]-ar1).*(pii1-pii2)+ (z0[1].dataMatrix-mean[5]-ar1).*pii2)^2); ss[1,2]=ss[1,1]; ss[1,3]=ss[1,1]; for i(1,T,1); ss[i+1,1]=s[1]+garch[1]*u1[i]^2+garch[2]*ss[i,1]; ss[i+1,2]=s[2]+garch[1]*u2[i]^2+garch[2]*ss[i,2]; ss[i+1,3]=s[3]+garch[1]*u3[i]^2+garch[2]*ss[i,3]; endfor; /* ss1 - ss3 are the standard error vectors ** Initial value is removed from the vector of standard errors */ ss1=ss[2:T+1,1]; ss2=ss[2:T+1,2]; ss3=ss[2:T+1,3]; /* log-likelihood */ f=((pdfn((z0[1].dataMatrix-mean[1]-ar1)./sqrt(ss1)))./sqrt(ss1)).*(1-pii1) + ((pdfn((z0[1].dataMatrix-mean[3]-ar1)./sqrt(ss2)))./sqrt(ss2)).*(pii1-pii2) + ((pdfn((z0[1].dataMatrix-mean[5]-ar1)./sqrt(ss3)))./sqrt(ss3)).*pii2; f=log(f); yy = cdfni( (1-pii1).*cdfn((z0[1].dataMatrix-mean[1]-ar1)./sqrt(ss1)) +(pii1-pii2).*cdfn((z0[1].dataMatrix-mean[3]-ar1)./sqrt(ss2)) +(pii2).*cdfn((z0[1].dataMatrix-mean[5]-ar1)./sqrt(ss3)) ); retp(f|yy); endp; /* ************************************************************************ */ /* Simulates a MAR341_GARCH(1,1) model with two endogenous regimes ** ** Result = datasimuMAR341GARCH11(x,T) ** returns T x 1 matrix of simulated values ** x - simulation values ** "mu1"|"c1"|"mu2"|"c2"|"mu3"|"s1"|"s2"|"s3"|"se"|"b1"|"b2"| ** "b3"|"b4"|"alfa"|"beta" ** T - number of observations to be simulated */ proc datasimuMAR431_GARCH11(struct PV x0, T); local mean,s,b,garch; mean = pvUnPack(x0,1); s = pvUnPack(x0,2); b = pvUnPack(x0,3); garch = pvUnPack(x0,4); local epsilon,eta,data,sigmas,pii1,pii2,boolean1,boolean2; epsilon=rndn(200+T,1); eta=s[4].*rndn(200+T,1); sigmas=zeros(200+T,3); data=zeros(200+T,1); data[1]=(mean[1]+mean[3]+mean[5])/3; data[2]=data[1]; data[3]=data[1]; data[4]=data[1]; data[5]=data[1]; boolean1=data[5]-eta[5]-mean[2]; boolean2=data[5]-eta[5]-mean[4]; for k (6,200+T,1); sigmas[k,1]=sqrt(s[1]^2+garch[2]*sigmas[k-1,1]^2+garch[1]* (data[k-1]-mean[1]-b[1]*data[k-2]-b[2]*data[k-3] -b[3]*data[k-4]-b[4]*data[k-5])^2); sigmas[k,2]=sqrt(s[2]^2+garch[2]*sigmas[k-1,2]^2+garch[1]* (data[k-1]-mean[3]-b[1]*data[k-2]-b[2]*data[k-3] -b[3]*data[k-4]-b[4]*data[k-5])^2); sigmas[k,3]=sqrt(s[3]^2+garch[2]*sigmas[k-1,3]^2+garch[1]* (data[k-1]-mean[5]-b[1]*data[k-2]-b[2]*data[k-3] -b[3]*data[k-4]-b[4]*data[k-5])^2); if boolean1<0; data[k]=mean[1]+b[1]*data[k-1]+b[2]*data[k-2] +b[3]*data[k-3]+b[4]*data[k-4] +sigmas[k,1]*epsilon[k]; elseif boolean2>=0; data[k]=mean[5]+b[1]*data[k-1]+b[2]*data[k-2] +b[3]*data[k-3]+b[4]*data[k-4] +sigmas[k,3]*epsilon[k]; else; data[k]=mean[3]+b[1]*data[k-1]+b[2]*data[k-2] +b[3]*data[k-3]+b[4]*data[k-4] +sigmas[k,2]*epsilon[k]; endif; boolean1=data[k]-eta[k]-mean[2]; boolean2=data[k]-eta[k]-mean[4]; endfor; data=data[201:200+T]; retp(data); endp;