gausset; library cmlmt,pgraph; #include cmlmt.sdf /* K1 - Lag number in autocorrelation test */ /* K2 - Lag number in heteroscedasticity test */ 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; /* UNIVARIATE TESTS */ /* Normality tests take 1 columns */ /* Autocorrelation tests take rows(matK1) columns */ /* Heteroscedasticity tests take rows(matK2) columns */ Result=zeros(1,1+rows(matK1)+rows(matK2)); /* The data is loaded */ load z[273,1]=EURUSD.txt; z=z[1:252,1]; #include ds.sdf struct DS z0; z0=reshape(dsCreate,3,1); z0[1].dataMatrix=z[3:rows(z)]; z0[2].dataMatrix=z[2:(rows(z)-1)]; z0[3].dataMatrix=z[1:(rows(z)-2)]; /* Initial values of parameters are given here ** varphi01 varphi02 varphi11 varphi12 ** variances: sigma1 sigma2 ** mixing coefficient: alpha */ struct PV x0; x0 = pvCreate; x0 = pvPacki(x0,0.16|-0.008|0.92|0,"varphi",1); x0 = pvPacki(x0,0.07|0.016,"sigma",2); x0 = pvPacki(x0,0.75,"alpha",3); /* If =1, a search is performed for the values of parameter alpha */ p_search=0; initialvalues = pvGetParVector(x0); AA1=((initialvalues[3]~initialvalues[4])|(1~0)); I_A1=eye(2*2)-AA1.*.AA1; Gamma1=inv(I_A1)*(initialvalues[5]|0|0|0); Gamma2=inv(I_A1)*(initialvalues[6]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); screen on; print "Eigenvalues of I_A1, Inv(I_A1)"; print eigh(I_A1)~eigh(inv(I_A1)); print "Eigenvalues of Ga1 and Ga2"; print eigh(Ga1)~eigh(Ga2); screen off; xapu=initialvalues; label_1: if p_search==1; grid_of_p=seqa(0.011,0.01,48); ll_values=zeros(48,1); xapu=pvGetParvector(x0); for i (1,48,1); x0=pvPutParVector(x0,(xapu[1:6]|grid_of_p[i])); ll_values[i,1]=meanc(fct2(x0,z0)); endfor; graphset; ylabel("Log-likelihood of the model"); xlabel("Parameter alpha"); _ptek="GMAR22_LL_wrt_p.tkf"; title("LL wrt parameter alpha for EM"); xy(grid_of_p,ll_values); ret = tkf2eps("GMAR22_LL_wrt_p.tkf", "GMAR22_LL_wrt_p.eps"); p_new=grid_of_p[maxindc(ll_values),1]; screen on; print "Search found: " p_new; screen off; x0=pvPutParVector(x0,(xapu[1:6]|p_new)); endif; /* Setting variable values for ESTIMATION */ struct cmlmtControl ctl; ctl = cmlmtControlCreate; ctl.title = "GMAR model: 2 regimes, 2 lags in each estimated using the cml-routine"; ctl.Algorithm = 3; /* = 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.PrintIters=1; ctl.maxIters=10000; ctl.Active={1,1,1,1,1,1,1,1,1}; ctl.IneqProc=&ineqp; proc ineqp (struct PV x0, struct DS z0); local v,w,AA1,I_A1,Gamma1,Gamma2,Ga1,Ga2,muu1,muu2; v = pvUnpack(x0,1); w = pvUnpack(x0,2); AA1=((v[3]~v[4])|(1~0)); I_A1=eye(2*2)-AA1.*.AA1; Gamma1=inv(I_A1)*(w[1]|0|0|0); Gamma2=inv(I_A1)*(w[2]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); retp(minc(eigh(Ga1)|eigh(Ga2)) - 1e-4); endp; ctl.Bounds = {-100 100, -100 100, -2 2, -0.99 0.99, 0 100, 0 100, 0.0001 0.9999}; ctl.DirTol = 1e-4; struct cmlmtResults out1; out1 = CMLmt(&ll,x0,z0,ctl); x0=pvPutParVector(x0,(pvGetParvector(out1.Par))); /* Log-Likelihood function for GMAR22_models */ proc ll(struct PV x0, struct DS z0, ind); local a, s, p, f, T, h11, h12, AA1,Gamma1,Gamma2; local Ga1,Ga2,muu1,muu2,h21,h22,zzz,tauhat; struct modelResults mm; a = pvUnpack(x0,1); s = pvUnpack(x0,2); p = pvUnpack(x0,3); T=rows(z0[1].dataMatrix); h11=(2*pi*s[1])^-0.5 *exp(-(z0[1].dataMatrix-a[1]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[1]/2); h12=(2*pi*s[2])^-0.5 *exp(-(z0[1].dataMatrix-a[2]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[2]/2); AA1=((a[3]~a[4])|(1~0)); Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0); Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); muu1=ones(2,1)*a[1]/(1-a[3]-a[4]); muu2=ones(2,1)*a[2]/(1-a[3]-a[4]); h21=zeros(T,1); h22=h21; for i (1,T,1); zzz=z0[2].dataMatrix[i]|z0[3].dataMatrix[i]; h21[i,1]=det(Ga1)^-0.5*exp(-(zzz-muu1)'*invpd(Ga1)*(zzz-muu1)/2); h22[i,1]=det(Ga2)^-0.5*exp(-(zzz-muu2)'*invpd(Ga2)*(zzz-muu2)/2); endfor; tauhat=p*h21./(p*h21+(1-p)*h22); f=tauhat.*h11[1:T] + (ones(T,1)-tauhat).*h12; if ind[1]; mm.Function = ln(f); endif; retp(mm); endp; /* This procedure computes the gradient for the log-likelihood function */ /* This procedure calculates the gradient for multivar. quantile residuals */ /* This procedure calculates the gradient for joint quantile residuals */ /* Returns a (T+nT+T) x 1 -matrix */ proc fct1(struct PV x0, struct DS z0); local a, s, p, T, h11, h12, AA1,Gamma1,Gamma2; local Ga1,Ga2,muu1,muu2,h21,h22,zzz,tauhat; a = pvUnpack(x0,1); s = pvUnpack(x0,2); p = pvUnpack(x0,3); T=rows(z0[1].dataMatrix); h11=(2*pi*s[1])^-0.5 *exp(-(z0[1].dataMatrix-a[1]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[1]/2); h12=(2*pi*s[2])^-0.5 *exp(-(z0[1].dataMatrix-a[2]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[2]/2); AA1=((a[3]~a[4])|(1~0)); Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0); Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); muu1=ones(2,1)*a[1]/(1-a[3]-a[4]); muu2=ones(2,1)*a[2]/(1-a[3]-a[4]); h21=zeros(T,1); h22=h21; for i (1,T,1); zzz=z0[2].dataMatrix[i]|z0[3].dataMatrix[i]; h21[i,1]=det(Ga1)^-0.5*exp(-(zzz-muu1)'*invpd(Ga1)*(zzz-muu1)/2); h22[i,1]=det(Ga2)^-0.5*exp(-(zzz-muu2)'*invpd(Ga2)*(zzz-muu2)/2); endfor; tauhat=p*h21./(p*h21+(1-p)*h22); local d,y1,y2,n,mu1,mu2; n=1; d=zeros(T,n); y1=zeros(T,1); y2=zeros(T,n); mu1=a[1].*ones(T,1)+a[3].*z0[2].dataMatrix+a[4].*z0[3].dataMatrix; mu2=a[2].*ones(T,1)+a[3].*z0[2].dataMatrix+a[4].*z0[3].dataMatrix; d=z0[1].dataMatrix; for t (1,T,1); y1[t,1]=ln(g(d[t,.]',mu1[t]',s[1],mu2[t]',s[2],tauhat[t])); y2[t,.]=resid(&g,&cumulg,d[t,.]',mu1[t],s[1],mu2[t],s[2],tauhat[t])'; endfor; retp(y1|vecr(y2')); endp; proc fct2(struct PV x0, struct DS z0); local a, s, p, T, h11, h12, AA1,Gamma1,Gamma2; local Ga1,Ga2,muu1,muu2,h21,h22,zzz,tauhat; a = pvUnpack(x0,1); s = pvUnpack(x0,2); p = pvUnpack(x0,3); T=rows(z0[1].dataMatrix); h11=(2*pi*s[1])^-0.5 *exp(-(z0[1].dataMatrix-a[1]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[1]/2); h12=(2*pi*s[2])^-0.5 *exp(-(z0[1].dataMatrix-a[2]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[2]/2); AA1=((a[3]~a[4])|(1~0)); Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0); Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); muu1=ones(2,1)*a[1]/(1-a[3]-a[4]); muu2=ones(2,1)*a[2]/(1-a[3]-a[4]); h21=zeros(T,1); h22=h21; for i (1,T,1); zzz=z0[2].dataMatrix[i]|z0[3].dataMatrix[i]; h21[i,1]=det(Ga1)^-0.5*exp(-(zzz-muu1)'*invpd(Ga1)*(zzz-muu1)/2); h22[i,1]=det(Ga2)^-0.5*exp(-(zzz-muu2)'*invpd(Ga2)*(zzz-muu2)/2); endfor; tauhat=p*h21./(p*h21+(1-p)*h22); local d,y1,n,mu1,mu2; n=1; d=zeros(T,n); y1=zeros(T,1); mu1=a[1].*ones(T,1)+a[3].*z0[2].dataMatrix+a[4].*z0[3].dataMatrix; mu2=a[2].*ones(T,1)+a[3].*z0[2].dataMatrix+a[4].*z0[3].dataMatrix; d=z0[1].dataMatrix; for t (1,T,1); y1[t,1]=ln(g(d[t,.]',mu1[t]',s[1],mu2[t]',s[2],tauhat[t])); endfor; retp(y1); endp; proc datasimu_GMAR22(struct PV x, Tboots); local a,s,p,AA1,Gamma1,Gamma2,Ga1,Ga2,muu1,muu2; a = pvUnpack(x,1); s = pvUnpack(x,2); p = pvUnpack(x,3); AA1=((a[3]~a[4])|(1~0)); Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0); Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); muu1=ones(2,1)*a[1]/(1-a[3]-a[4]); muu2=ones(2,1)*a[2]/(1-a[3]-a[4]); local epsilon, boolean, data, ey1,ey2,bool; epsilon=rndn(Tboots+200,1); data=zeros(Tboots+200,1); bool=zeros(Tboots+200,1); boolean=RNDBernoulli(p); if boolean==1; data[1:2]=muu1+chol(Ga1)'*epsilon[1:2]; else; data[1:2]=muu2+chol(Ga2)'*epsilon[1:2]; endif; ey1=p*det(Ga1)^-0.5*exp(-(data[1:2]-muu1)'*invpd(Ga1)*(data[1:2]-muu1)/2); ey2=(1-p)*det(Ga2)^-0.5*exp(-(data[1:2]-muu2)'*invpd(Ga2)*(data[1:2]-muu2)/2); boolean=RNDBernoulli(ey1/(ey1+ey2)); bool[1,1]=boolean; for k (3,Tboots+200,1); if boolean==1; data[k]=a[1]+a[3]*data[k-1]+a[4]*data[k-2]+sqrt(s[1])*epsilon[k]; else; data[k]=a[2]+a[3]*data[k-1]+a[4]*data[k-2]+sqrt(s[2])*epsilon[k]; endif; ey1=p*det(Ga1)^-0.5*exp(-(data[(k-1):k]-muu1)'*invpd(Ga1)*(data[(k-1):k]-muu1)/2); ey2=(1-p)*det(Ga2)^-0.5*exp(-(data[(k-1):k]-muu2)'*invpd(Ga2)*(data[(k-1):k]-muu2)/2); boolean=RNDBernoulli(ey1/(ey1+ey2)); bool[k,1]=boolean; endfor; retp(data[201:Tboots+200]|bool[201:Tboots+200]); endp; proc RNDBernoulli(p); local bool; bool=rndu(1,1); if p<0; retp (error(0)); elseif p>1; retp (error(0)); else; if bool<=p; retp(1); else; retp(2); endif; endif; endp; /* CALCULATION OF QUANTILE RESIDUALS */ T=rows(z0[1].dataMatrix); lqr=fct1(x0,z0); mr=lqr[(T+1):2*T,.]; /* NUMERICAL DERIVATIVES */ #include optim.sdf deriv = gradMT(&fct1,x0,z0); s = deriv[1:T,.]; mg = deriv[(T+1):2*T,.]; x0hessian = hessMT(&fct2,x0,z0); /* BOOTSRAPPED QUANTILE RESIDUALS */ struct DS z0boots; z0boots=reshape(dsCreate,3,1); apuri=datasimu_GMAR22(x0,Tboots+2); z0boots[1].dataMatrix=apuri[3:(Tboots+2)]; z0boots[2].dataMatrix=apuri[2:(Tboots+1)]; z0boots[3].dataMatrix=apuri[1:(Tboots)]; lqrboots=fct1(x0,z0boots); mrboots=lqrboots[(Tboots+1):2*Tboots,.]; /* THEIR NUMERICAL DERIVATIVES */ dim=rows(initialvalues); n=cols(apuri); sboots=zeros(Tboots,dim); mgboots=zeros(Tboots,n*dim); derivboots = gradMT(&fct1,x0,z0boots); sboots=derivboots[1:Tboots,.]; mgboots=derivboots[(Tboots+1):(2*Tboots),.]; /****************************************************************************/ /* UNIVARIATE NORMALITY TESTS */ /****************************************************************************/ Result[1,1]= ntestingboots(mr,mrboots,sboots,mgboots); /****************************************************************************/ /* UNIVARIATE AUTOCORRELATION TEST */ /****************************************************************************/ Result[1,2:1+rows(matK1)] =atestingboots(mr,mrboots,sboots,mgboots,matK1); /***************************************************************************/ /* UNIVARIATE HETEROSCEDASTICITY TEST */ /****************************************************************************/ Result[1,2+rows(matK1):1+rows(matK1)+rows(matK2)] =htestingboots(mr,mrboots,sboots,mgboots,matK2); let pr=1; outwidth 80; output file = estim_EURUSD_GMAR22_constrained.asc on; screen off; call cmlmtPrt(out1); print ""; print "Initial values used in estimation"; print pvGetParNames(x0)'; print initialvalues'; print ""; print "Sample size used for estimation " ftos(T,"%*.*lf",1,0); print ""; print "The value of the likelihood at minimum " ftos(out1.Fct,"%*.*lf",1,3); print ""; print "The AIC value " ftos(-2*(out1.Fct)+2*(rows(initialvalues)-1),"%*.*lf",1,3); print ""; print "The BIC value " ftos(-2*(out1.Fct)+ln(T)*(rows(initialvalues)-1),"%*.*lf",1,3); xy0=pvGetParVector(x0); h11=(2*pi*xy0[5])^-0.5 *exp(-(z0[1].dataMatrix-xy0[1]-xy0[3]*z0[2].dataMatrix -xy0[4]*z0[3].dataMatrix)^2/xy0[5]/2); h12=(2*pi*xy0[6])^-0.5 *exp(-(z0[1].dataMatrix-xy0[2]-xy0[3]*z0[2].dataMatrix -xy0[4]*z0[3].dataMatrix)^2/xy0[6]/2); AA1=((xy0[3]~xy0[4])|(1~0)); Gamma1=inv(eye(2*2)-AA1.*.AA1)*(xy0[5]|0|0|0); Gamma2=inv(eye(2*2)-AA1.*.AA1)*(xy0[6]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); muu1=ones(2,1)*xy0[1]/(1-xy0[3]-xy0[4]); muu2=ones(2,1)*xy0[2]/(1-xy0[3]-xy0[4]); print "Parameters of the stationary distribution"; print "Constant_1 and Myy_1:" xy0[1]~muu1[1]; print "Covariance matrix: " Ga1; print "Constant_2 and Myy_2:" xy0[2]~muu2[1]; print "Covariance matrix: " Ga2; print ""; nprtestingeb(Result[1,1]',pr); print ""; aprtestingeb(Result[1,2:1+rows(matK1)]',pr,matK1); print ""; hprtestingeb(Result[1,2+rows(matK1):1+rows(matK1)+rows(matK2)]',pr,matK2); print ""; apuria=seqa(1989,1/12,rows(z0[1].dataMatrix)); graphset; ylabel("Original Series"); xlabel("Time"); _pdate=""; if cols(z0[1].dataMatrix)>1; _pltype = seqa(1,1,cols(z0[1].dataMatrix)); endif; xtics(1989,2010,4,1); _ptek="GMAR22Series.tkf"; title("Time series of the data"); xy(apuria,z0[1].dataMatrix); ret = tkf2eps("GMAR22Series.tkf", "GMAR22Series.eps"); apuria=seqa(0,1,rows(z0boots[1].dataMatrix)); graphset; ylabel("Simulated series used in tests"); xlabel("Time"); if cols(z0boots[1].dataMatrix)>1; _pltype = seqa(1,1,cols(z0boots[1].dataMatrix)); endif; _ptek="GMAR22Simuseries.tkf"; title("Time series of the simulated data"); xy(apuria,z0boots[1].dataMatrix); ret = tkf2eps("GMAR22Simuseries.tkf", "GMAR22Simuseries.eps"); if rows(z0boots[1].dataMatrix)>2000; apuria=seqa(0,1,2000); graphset; ylabel("Simulated series used in tests"); xlabel("Time"); if cols(z0boots[1].dataMatrix)>1; _pltype = seqa(1,1,cols(z0boots[1].dataMatrix)); endif; _ptek="GMAR22Simu2000series.tkf"; title("2000 first obs. of the simulated data"); xy(apuria,z0boots[1].dataMatrix[1:2000]); ret = tkf2eps("GMAR22Simu2000series.tkf", "GMAR22Simu2000series.eps"); endif; bin=0.10; min_z=minc(z0[1].dataMatrix); max_z=maxc(z0[1].dataMatrix); Pearson=(seqa(1,1,ceil(1/bin)))/ceil(1/bin)*(max_z-min_z); Pearson=min_z*ones(rows(Pearson),1)+Pearson; apuria=seqa(1,1,rows(Pearson)*1000)/1000/rows(Pearson)*(max_z-min_z); apuria=min_z*ones(rows(apuria),1)+apuria; apuria=min_z|apuria|max_z+0.0001; luvut=counts(z0[1].dataMatrix,Pearson); luvut=luvut./rows(z0[1].dataMatrix); graphset; ylabel("Proportion"); xlabel("Original series"); xtics(floor(min_z),ceil(max_z),(1/rows(Pearson)),1); ytics(0,ceil(maxc(100.*luvut))/100,0.1,1); _pltype={6}; _pcolor={15}; _pxpmax=2; _pscreen=1; _ptek="GMAR22Histo.tkf"; axmargin(1,1,.5,.855); title("Histogram of the data"); xy(apuria,0|(luvut).*.ones(1000,1)|0); ret = tkf2eps("GMAR22Histo.tkf", "GMAR22Histo.eps"); apuria=seqa(1989,1/12,rows(mr)); graphset; ylabel("QR Series"); xlabel("Time"); if cols(mr)>1; _pltype = seqa(1,1,cols(mr)); endif; xtics(1989,2010,4,1); _pcolor={15}; _ptek="GMAR22QRSeries.tkf"; title("Time series of the QRs"); xy(apuria,mr); ret = tkf2eps("GMAR22QRSeries.tkf", "GMAR22QRSeries.eps"); bin=0.10; min_mr=minc(mr); max_mr=maxc(mr); Pearson=(seqa(1,1,ceil(1/bin)))/ceil(1/bin)*(max_mr-min_mr); Pearson=min_mr*ones(rows(Pearson),1)+Pearson; apuria=seqa(1,1,1000*rows(Pearson))/1000/rows(Pearson)*(max_mr-min_mr); apuria=min_mr*ones(rows(apuria),1)+apuria; apuria=min_mr|apuria|max_mr+0.0001; luvut=counts(mr,Pearson); luvut=luvut./rows(mr); graphset; ylabel("Proportion"); xlabel("Quantile Residuals"); xtics(floor(min_mr),ceil(max_mr),(1/rows(Pearson)),1); ytics(0,ceil(maxc(100.*luvut))/100,0.1,1); _pltype={6}; _pcolor={9}; _pxpmax=2; _pscreen=1; _ptek="GMAR22QRHisto.tkf"; axmargin(1,1,.5,.855); title("Histogram of the Quantile residuals"); xy(apuria,0|(luvut).*.ones(1000,1)|0); ret = tkf2eps("GMAR22QRHisto.tkf", "GMAR22QRHisto.eps"); order_mr=sortc(mr,1); normals_1=cdfni(seqa(1,1,rows(mr))/(rows(mr)+1)); graphset; ylabel("Quantile Residuals"); xlabel("Theoretical Quantiles"); min_val=floor(minc(order_mr[1]|normals_1[1])); max_val=ceil(maxc(order_mr[rows(mr)]|normals_1[rows(normals_1)])); line45=seqa(0,1,rows(mr))*(max_val-min_val)/(rows(mr)-1)+min_val; ytics(min_val,max_val,1,0.5); xtics(min_val,max_val,1,0.5); _pltype={2,6}; _pstype={1,0}; _pcolor={9}; _plwidth={20,0}; _pxpmax=0; _pscreen=1; _ptek="GMAR22QR_QQ.tkf"; axmargin(1,1,.5,.855); title("Q-Q Plot of Quantile residuals"); xy(normals_1,order_mr~normals_1); ret = tkf2eps("GMAR22QR_QQ.tkf","GMAR22QR_QQ.eps"); bin=1.06*stdc(mr)*rows(mr)^(-1/5); R_scale=maxc(abs(mr)); R_domain=seqa(-R_scale,0.01,2*R_scale*100+1); R_density=ones(rows(R_domain),1); for i (1,rows(R_domain),1); R_density[i]=meanc(pdfn((ones(rows(mr),1)*R_domain[i]-mr)/bin))/bin; endfor; R_normal=pdfn((R_domain-meanc(mr))/stdc(mr)); graphset; ylabel(""); xlabel(""); xtics(floor(-R_scale),ceil(R_scale),1,1); ytics(0,ceil(maxc(100.*R_density))/100,0.1,1); _pltype={6,6}; _pcolor={15,4}; _pxpmax=2; _plegctl= { 2, 3, 1.7, 4.5 }; /* legend size and locations */ _plegstr= "Kernel estimator \0"\ "Scaled normal distribution"; _pscreen=1; _ptek="GMAR22QRdensity.tkf"; axmargin(1,1,.5,.855); title("N(0,1) Kernel density of the Quantile residuals"); xy(R_domain,R_density~R_normal); ret = tkf2eps("GMAR22QRdensity.tkf", "GMAR22QRdensity.eps"); a = pvUnpack(x0,1); s = pvUnpack(x0,2); p = pvUnpack(x0,3); T=rows(z0[1].dataMatrix); h11=(2*pi*s[1])^-0.5 *exp(-(z0[1].dataMatrix-a[1]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[1]/2); h12=(2*pi*s[2])^-0.5 *exp(-(z0[1].dataMatrix-a[2]-a[3]*z0[2].dataMatrix -a[4]*z0[3].dataMatrix)^2/s[2]/2); AA1=(a[3]~a[4])|(1~0); Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0); Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0); Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]); Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]); muu1=ones(2,1)*a[1]/(1-a[3]-a[4]); muu2=ones(2,1)*a[2]/(1-a[3]-a[4]); h21=zeros(T,1); h22=h21; for i (1,T,1); zzz=z0[2].dataMatrix[i]|z0[3].dataMatrix[i]; h21[i,1]=det(Ga1)^-0.5*exp(-(zzz-muu1)'*invpd(Ga1)*(zzz-muu1)/2); h22[i,1]=det(Ga2)^-0.5*exp(-(zzz-muu2)'*invpd(Ga2)*(zzz-muu2)/2); endfor; tau1=p[1]*h21./(p[1]*h21+(1-p[1])*h22); tau2=1-tau1; for i (1,rows(tau2),1); if tau2[i,1]<=0; tau1[i,1]=tau1[i,1]+tau2[i,1]-1e-8; tau2[i,1]=1e-8; endif; endfor; zhat1=(h11.*tau1)./(h11.*tau1 + h12.*tau2); zhat2=(1-zhat1); for i (1,rows(zhat2),1); if zhat2[i,1]<0; zhat1[i,1]=zhat1[i,1]+zhat2[i,1]; zhat2[i,1]=0; endif; endfor; apuria=seqa(1989,1/12,rows(tau1)); graphset; ylabel("Mixing proportions"); xlabel("Time"); _pltype ={1,3}; _ptek="GMAR22MixingSeriesf.tkf"; ytics(0,1,0.1,1); xtics(1989,2010,3,1); title("Mixing propostions: Forecast"); xy(apuria,tau1~(ones(T,1)-tau2)); ret = tkf2eps("GMAR22MixingSeriesf.tkf", "GMAR22MixingSeriesf.eps"); apuria=seqa(1989,1/12,rows(zhat1)); graphset; ylabel("Mixing proportions"); xlabel("Time"); _pltype ={1,3}; _ptek="GMAR22MixingSeriesi.tkf"; ytics(0,1,0.1,1); xtics(1989,2010,3,1); title("Mixing propostions: Inference"); xy(apuria,zhat1~(ones(T,1)-zhat2)); ret = tkf2eps("GMAR22MixingSeriesi.tkf", "GMAR22MixingSeriesi.eps"); m1t=a[1]+a[3]*z0[2].dataMatrix+a[4]*z0[3].dataMatrix; m2t=a[2]+a[3]*z0[2].dataMatrix+a[4]*z0[3].dataMatrix; conde=zeros(T,1); condv=conde; for i (1,T,1); conde[i]=tau1[i]*m1t[i]+tau2[i]*m2t[i]; condv[i]=tau1[i]*s[1]+tau2[i]*s[2] + tau1[i]*(m1t[i]-conde[i])^2 + tau2[i]*(m2t[i]-conde[i])^2; endfor; graphset; ylabel("Conditional expectation"); xlabel("Time"); _pltype ={6}; _ptek="GMAR22CondE.tkf"; xtics(1989,2010,3,1); title("Conditional expectation"); xy(apuria,conde); ret = tkf2eps("GMAR22CondE.tkf", "GMAR22CondE.eps"); graphset; ylabel("Conditional variance"); xlabel("Time"); _pltype ={6}; _ptek="GMAR22CondV.tkf"; xtics(1989,2010,3,1); title("Conditional variance"); xy(apuria,condv); ret = tkf2eps("GMAR22CondV.tkf", "GMAR22CondV.eps"); graphset; ylabel("Conditional variance"); xlabel("Conditional expectation"); _plctrl = {-1}; _pstype={1}; /*_pltype={0}; _plwidth={2}; _pxpmax=0; _pscreen=1;*/ _ptek="GMAR22CondEV.tkf"; title("Observations"); xy(conde,condv); ret = tkf2eps("GMAR22CondEV.tkf", "GMAR22CondEV.eps"); graphset; fonts("simplex complex microb simgrma"); ylabel(""); xlabel(""); _pcolor={0, 0}; _pltype ={6, 3}; _ptek="GMAR22All.tkf"; _plegctl= { 3, 5, 300, 400 }; /* legend size and locations */ _plegstr= "y]t[\0"\ /* 4 lines legend text */ "Scaled \204a\201]1,t["; title(""); xtics(1989,2010,3,1); /*tau1 and zhat1 need to be (0,1)*/ min_z=minc(z0[1].dataMatrix); max_z=maxc(z0[1].dataMatrix); diff_z=max_z-min_z; tau1_scaled=tau1*diff_z+min_z*ones(rows(tau1),1); zhat1_scaled=zhat1*diff_z+min_z*ones(rows(tau1),1); xy(apuria,z0[1].dataMatrix~tau1_scaled); ret = tkf2eps("GMAR22All.tkf", "GMAR22All.eps"); graphset; makewind(4.5,3.4275,0,0,0); makewind(4.5,3.4275,4.5,0,0); setwind(1); fonts("simplex complex microb simgrma"); ylabel(""); xlabel(""); _pcolor={0, 0}; _pltype ={6, 3}; _ptek="GMAR22All_2.tkf"; title(""); xtics(1989,2010,3,1); /*tau1 and zhat1 need to be (0,1)*/ min_z=minc(z0[1].dataMatrix); max_z=maxc(z0[1].dataMatrix); diff_z=max_z-min_z; tau1_scaled=tau1*diff_z+min_z*ones(rows(tau1),1); zhat1_scaled=zhat1*diff_z+min_z*ones(rows(tau1),1); xy(apuria,z0[1].dataMatrix~tau1_scaled); nextwind; graphset; bin=1.06*stdc(z0[1].dataMatrix)*rows(z0[1].dataMatrix)^(-1/5); screen on; print bin; screen off; R_scale=maxc(abs(min_z|max_z)); R_domain=seqa(-R_scale,0.01,2*R_scale*100+1); R_density=ones(rows(R_domain),1); for i (1,rows(R_domain),1); R_density[i]=meanc(pdfn((ones(rows(z0[1].dataMatrix),1)*R_domain[i]-z0[1].dataMatrix)/bin))/bin; endfor; mu1=-0.321; mu2=1.175; sigma1=0.221; sigma2=1.214; alpha=0.586; R_mixture=alpha*(2*pi*sigma1)^(-0.5)*exp(-0.5*(R_domain-mu1)^2/(sigma1)) +(1-alpha)*(2*pi*sigma2)^(-0.5)*exp(-0.5*(R_domain-mu2)^2/(sigma2)); ylabel(""); xlabel(""); xtics(floor(-R_scale),ceil(R_scale),1,1); ytics(0,ceil(maxc(100.*(R_density|R_mixture)))/100,0.1,1); _pltype={6,1}; _pcolor={0,0}; _pdate=""; _pxpmax=2; _pscreen=1; _ptek="GMAR22All_2.tkf"; title(""); xy(R_domain,R_density~R_mixture); endwind; ret = tkf2eps("GMAR22All_2.tkf", "GMAR22All_2.eps"); graphset; makewind(4.5,3.4275,0,0,0); makewind(4.5,3.4275,4.5,0,0); setwind(1); fonts("simplex complex microb simgrma"); ylabel(""); xlabel(""); _pcolor={0, 0}; _pltype ={6, 3}; _ptek="GMAR22All_3.tkf"; _plegctl= { 3, 5, 300, 400 }; /* legend size and locations */ _plegstr= "y]t["; /* 4 lines legend text */ title(""); xtics(1989,2010,3,1); /*tau1 and zhat1 need to be (0,1)*/ min_z=minc(z0[1].dataMatrix); max_z=maxc(z0[1].dataMatrix); diff_z=max_z-min_z; tau1_scaled=tau1*diff_z+min_z*ones(rows(tau1),1); zhat1_scaled=zhat1*diff_z+min_z*ones(rows(tau1),1); xy(apuria,z0[1].dataMatrix); nextwind; graphset; bin=1.06*stdc(z0[1].dataMatrix)*rows(z0[1].dataMatrix)^(-1/5); screen on; print bin; screen off; R_scale=maxc(abs(min_z|max_z)); R_domain=seqa(-R_scale,0.01,2*R_scale*100+1); R_density=ones(rows(R_domain),1); for i (1,rows(R_domain),1); R_density[i]=meanc(pdfn((ones(rows(z0[1].dataMatrix),1)*R_domain[i]-z0[1].dataMatrix)/bin))/bin; endfor; mu1=-0.32; mu2=1.18; sigma1=0.22; sigma2=1.21; alpha=0.59; R_mixture=alpha*(2*pi*sigma1)^(-0.5)*exp(-0.5*(R_domain-mu1)^2/(sigma1)) +(1-alpha)*(2*pi*sigma2)^(-0.5)*exp(-0.5*(R_domain-mu2)^2/(sigma2)); ylabel(""); xlabel(""); xtics(floor(-R_scale),ceil(R_scale),1,1); ytics(0,ceil(maxc(100.*(R_density|R_mixture)))/100,0.1,1); _pltype={6,1}; _pcolor={0,0}; _pdate=""; _pxpmax=2; _plegctl= { 3, 5, 300, 400 }; /* legend size and locations */ _plegstr= "Kernel estimate for observations"; _pscreen=1; _ptek="GMAR22All_3.tkf"; title(""); xy(R_domain,R_density); endwind; ret = tkf2eps("GMAR22All_3.tkf", "GMAR22All_3.eps"); graphset; fonts("simplex complex microb simgrma"); x = seqa(min_z-2,.01,(max_z+4-min_z)*100); /* generate data */ y1 = (1/sqrt(2*pi*Ga1[1,1]))*exp(-(x-muu1[1,1]).*(x-muu1[1,1])./Ga1[1,1]/2); y2 = (1/sqrt(2*pi*Ga2[1,1]))*exp(-(x-muu2[1,1]).*(x-muu2[1,1])./Ga2[1,1]/2); _pdate = ""; /* date is not printed */ _pcolor= {0, 0, 0}; _pltype = { 6, 3, 3}; /* dotted, solid lines */ _plegctl= { 3, 5, 2000, 2000 }; /* legend size and locations */ _plegstr= "\204a\201*N(\204m\201]1[,\204g\201]1,0[)+(1-\204a\201)*N(\204m\201]2[,\204g\201]2,0[) \0"\ "N(\204m\201]1[,\204g\201]1,0[)\0"\ "N(\204m\201]2[,\204g\201]2,0[) \0"; xtics(-3,4,1,0.5); _ptek="GMAR22Mixturedensity.tkf"; title(""); /* main title */ xy(x,p[1]*y1+(1-p[1])*y2~p[1]*y1~(1-p[1])*y2); ret = tkf2eps("GMAR22Mixturedensity.tkf", "GMAR22Mixturedensity.eps"); fonts("simplex"); graphset; x = seqa(min_z-2,.01,(max_z+4-min_z)*100); h21=zeros(rows(x),1); h22=h21; tau1=zeros(rows(x),rows(x)); zz=x; zzz=x~x; for i (1,rows(x),1); zzz[.,1]=x[i,1].*ones(rows(x),1); for j (1,rows(x),1); h21[j,1]=det(Ga1)^-0.5*exp(-(zzz[j,.]-muu1')*invpd(Ga1)*(zzz[j,.]-muu1')'/2); h22[j,1]=det(Ga2)^-0.5*exp(-(zzz[j,.]-muu2')*invpd(Ga2)*(zzz[j,.]-muu2')'/2); endfor; tau1[.,i]=p[1]*h21./(p[1]*h21+(1-p[1])*h22); endfor; tau=tau1[.,1]; for i (2,floor(cols(tau1)/50),1); tau=tau~tau1[.,i]; endfor; tau=tau~tau1[.,cols(tau1)]; ytics(0,1,1,0.5); _pdate = ""; /* date is not printed */ _ptek="GMAR22Mixturefunction.tkf"; title("Mixture function"); /* main title */ xy(x,tau); ret = tkf2eps("GMAR22Mixturefunction.tkf", "GMAR22Mixturefunction.eps"); graphset; _pdate = ""; /* date is not printed */ _ptek="GMAR22Mixturefunction_2.tkf"; xtics(-3,-1,1,0.5); title("Mixture function"); /* main title */ xy(x,tau); ret = tkf2eps("GMAR22Mixturefunction_2.tkf", "GMAR22Mixturefunction_2.eps"); graphset; /* reset global variables */ x = seqa(min_z-2,.1,(max_z+1-min_z)*10); h21=zeros(rows(x),1); h22=h21; tau1=zeros(rows(x),rows(x)); zz=x; zzz=x~x; for i (1,rows(x),1); zzz[.,1]=x[i,1].*ones(rows(x),1); for j (1,rows(x),1); h21[j,1]=det(Ga1)^-0.5*exp(-(zzz[j,.]-muu1')*invpd(Ga1)*(zzz[j,.]-muu1')'/2); h22[j,1]=det(Ga2)^-0.5*exp(-(zzz[j,.]-muu2')*invpd(Ga2)*(zzz[j,.]-muu2')'/2); endfor; tau1[.,i]=p[1]*h21./(p[1]*h21+(1-p[1])*h22); endfor; _pzclr = { 1, 2, 3, 4 }; /* set Z level colors */ _ptek="3dmixturesGMAR22.tkf"; view(-5,-20,10); title(""); /* set main title */ xlabel("y_{t-2}"); /* set X axis label */ ylabel("y_{t-1}"); /* set Y axis label */ scale3d(miss(0,0),miss(0,0),0|1.0); /* scale Z axis */ surface(zz',zz,tau1); ret = tkf2eps("3dmixturesGMAR22.tkf", "3dmixturesGMAR22.eps"); KK1=10; T=rows(mr); Tboots=rows(mrboots); k=cols(sboots); covboots=invpd((sboots'*sboots)/Tboots); GG2=zeros(KK1,k); RR=zeros(KK1,1); guboots=zeros(Tboots-KK1,KK1); help=0; for t1 (1,(Tboots-KK1),1); help=((mrboots[t1+KK1])*mrboots[t1+KK1-1]); for i (2,KK1,1); help=help~((mrboots[t1+KK1])*mrboots[t1+KK1-i]); endfor; guboots[t1,.]=help; endfor; Sigma=(sboots[KK1+1:Tboots,.]~guboots)'*(sboots[KK1+1:Tboots,.]~guboots)/(Tboots-KK1); for i1 (1,KK1,1); GG2[i1,.]=(meanc((mrboots[1+i1:Tboots]*~mgboots[1:Tboots-i1,.]) +(mrboots[1:Tboots-i1]*~mgboots[1+i1:Tboots,.])))'; RR[i1]=meanc(mr[1+i1:T]*~mr[1:T-i1]); endfor; O2=(GG2*covboots~eye(KK1))*Sigma*(GG2*covboots~eye(KK1))'; autoc=0; for lk (1,rows(RR),1); autoc=autoc|(RR[lk]./sqrt(O2[lk,lk])); endfor; autoc=autoc[2:rows(autoc)]; screen on; print "ACF of Quantile residuals"; print RR'; print "Scaled ACF of Quantile residuals"; print autoc'; screen off; lags=seqa(0,1/100,100*rows(autoc)); autocu=autoc[1]|zeros(99,1); for kkm (2,rows(autoc),1); autocu=autocu|autoc[kkm]|zeros(99,1); endfor; autoc=autocu; autoc[1]=0; ylar=maxc(abs(autoc))+0.02; critic=0.995; bounds=(cdfni(critic).*ones(rows(autoc),1))./sqrt(T-KK1); ylar=maxc(abs(autoc)|bounds)+0.02; autoc1=autoc; bounds1=bounds; graphset; begwind; makewind(9,6.855,0,0,0); makewind(9,6.855,0,0,1); setwind(1); graphset; ylabel("Autocorrelation"); xlabel("Lags"); title("99% Critical bounds"); xtics(0,rows(lags)/100,1,1); ytics(-ylar,ylar,(ylar/3),1); _pltype={6}; _pcolor={9}; _plwidth=15; _pxpmax=0; _pypmax=2; _ptek="GMAR22AcorQR.tkf"; axmargin(1,1,.5,.855); xy(lags,autoc); nextwind; graphset; axmargin(1,1,.5,.855); ylabel("Autocorrelation"); xlabel("Lags"); title("99% Critical bounds"); xtics(0,rows(lags)/100,1,1); ytics(-ylar,ylar,(ylar/3),1); _pltype={1,1}; _pcolor={9,9}; _pxpmax=0; _pypmax=2; _ptek="GMAR22AcorQR.tkf"; xy(lags,bounds~(-bounds)); endwind; ret = tkf2eps("GMAR22AcorQR.tkf", "GMAR22AcorQR.eps"); KK2=10; T=rows(mr); Tboots=rows(mrboots); k=cols(sboots); covboots=invpd((sboots'*sboots)/Tboots); GG3=zeros(KK2,k); RRR=zeros(KK2,1); guboots=zeros(Tboots-KK2,KK2); for i2 (1,KK2,1); GG3[i2,.]=2.*meanc( (mrboots[KK2+1-i2:Tboots-i2]^2-1)*~mrboots[KK2+1:Tboots]*~mgboots[KK2+1:Tboots,.]+ ((mrboots[KK2+1:Tboots]^2-1)*~mrboots[KK2+1-i2:Tboots-i2])*~mgboots[KK2+1-i2:Tboots-i2,.])'; RRR[i2,1]=meanc((mr[KK2+1:T]^2-1)*~(mr[KK2+1-i2:T-i2]^2-1)); guboots[.,i2]=(mrboots[KK2+1:Tboots]^2-1)*~(mrboots[KK2+1-i2:Tboots-i2]^2-1); endfor; Sigma=(sboots[KK2+1:Tboots,.]~guboots)'*(sboots[KK2+1:Tboots,.]~guboots)/(Tboots-KK2); O3=(GG3*covboots~eye(KK2))*Sigma*(GG3*covboots~eye(KK2))'; autoc=0; for lk (1,rows(RRR),1); autoc=autoc|(RRR[lk]/sqrt(O3[lk,lk])); endfor; autoc=autoc[2:rows(autoc)]; screen on; print "ACF of Squared Quantile residuals"; print RRR'; print "Scaled ACF of Squared Quantile residuals"; print autoc'; screen off; lags=seqa(0,1/100,100*rows(autoc)); autocu=autoc[1]|zeros(99,1); for kkm (2,rows(autoc),1); autocu=autocu|autoc[kkm]|zeros(99,1); endfor; autoc=autocu; autoc[1]=0; critic=0.995; bounds=(cdfni(critic).*ones(rows(autoc),1))./sqrt(T-KK2); ylar=maxc(abs(autoc)|bounds)+0.02; autoc2=autoc; bounds2=bounds; graphset; begwind; makewind(9,6.855,0,0,0); makewind(9,6.855,0,0,1); setwind(1); graphset; ylabel("Autocorrelation"); xlabel("Lags"); title("99% Critical bounds"); xtics(0,rows(lags)/100,1,1); ytics(-ylar,ylar,(ylar/3),1); _pltype={6}; _pcolor={9}; _plwidth=15; _pxpmax=0; _pypmax=2; _ptek="GMAR22AcorQR2.tkf"; axmargin(1,1,.5,.855); xy(lags,autoc); nextwind; graphset; axmargin(1,1,.5,.855); ylabel("Autocorrelation"); xlabel("Lags"); title("99% Critical bounds"); xtics(0,rows(lags)/100,1,1); ytics(-ylar,ylar,(ylar/3),1); _pltype={1,1}; _pcolor={9,9}; _pxpmax=0; _pypmax=2; _ptek="GMAR22AcorQR2.tkf"; xy(lags,bounds~(-bounds)); endwind; ret = tkf2eps("GMAR22AcorQR2.tkf", "GMAR22AcorQR2.eps"); graphset; begwind; window(2,2,0); setwind(1); ylabel("QR Series"); xlabel("Time"); xtics(1989,2010,4,1); _pcolor={15}; _ptek="GMAR22diagnostics.tkf"; title("Time series of QRs"); xy(apuria,mr); nextwind; graphset; ylabel("Observed Quantiles"); xlabel("Theoretical Quantiles"); min_val=floor(minc(order_mr[1]|normals_1[1])); max_val=ceil(maxc(order_mr[rows(mr)]|normals_1[rows(normals_1)])); line45=seqa(0,1,rows(mr))*(max_val-min_val)/(rows(mr)-1)+min_val; ytics(min_val,max_val,1,0.5); xtics(min_val,max_val,1,0.5); _pltype={2,6}; _pstype={1,0}; _pcolor={0,0}; _plwidth={20,0}; _pxpmax=0; _pscreen=1; _ptek="GMAR22diagnostics.tkf"; title("Q-Q Plot of Quantile residuals"); xy(normals_1,order_mr~normals_1); nextwind; graphset; ylabel("Scaled autocovariance"); xlabel("Lags"); title("99% Critical bounds"); xtics(0,rows(lags)/100,1,1); ytics(-ylar,ylar,(ylar/3),1); _pltype={6,2,2}; _pcolor={0,0,0}; _plwidth=15; _pxpmax=0; _pypmax=2; _ptek="GMAR22diagnostics.tkf"; xy(lags,autoc1~bounds1~(-bounds1)); nextwind; graphset; ylabel("Scaled autocovariance"); xlabel("Lags"); title("99% Critical bounds"); xtics(0,rows(lags)/100,1,1); ytics(-ylar,ylar,(ylar/3),1); _pltype={6,2,2}; _pcolor={0,0,0}; _plwidth=15; _pxpmax=0; _pypmax=2; _ptek="GMAR22diagnostics.tkf"; xy(lags,autoc2~bounds2~(-bounds2)); endwind; ret = tkf2eps("GMAR22diagnostics.tkf", "GMAR22diagnostics.eps"); apurib=seqa(0,0.002,250); apuria=ones(rows(apurib),1); apuric=ones(T,1); struct PV uusib; uusib = pvCreate; uusib=x0; uus=pvGetParVector(x0); for i (1,rows(apurib),1); uusib=pvPutParVector(uusib,(uus[1:6,1]|apurib[i,1])); apuric=fct2(uusib,z0); apuria[i,1]=meanc(apuric[1:T]); endfor; graphset; ylabel("Log-likelihood of the model"); xlabel("Parameter p"); if cols(apuria)>1; _pltype = seqa(1,1,cols(apuria)); endif; _ptek="GMAR22_LogL_wrt_p.tkf"; title("Log-L wrt parameter p"); xy(apurib,apuria); ret = tkf2eps("GMAR22_LogL_wrt_p.tkf", "GMAR22_LogL_wrt_p.eps"); /* CODES FOR COMPUTING QUANTILE RESIDUAL TESTS */ /* NORMALITY TEST */ proc ntestingboots(r,rboots,sboots,gboots); local Res; Res=zeros(1,1); Res[1,1]=cdfchic(ntestboots(r,rboots,sboots,gboots),3); retp(Res); endp; proc ntestboots(r,rboots,sboots,gboots); local G1,O1,T1,T,mv,k,Sigma,Tboots,cov; T=rows(r); Tboots=rows(rboots); k=cols(sboots); cov=invpd((sboots'*sboots)/Tboots); G1=(2.*meanc(rboots*~gboots))'|(3.*meanc((rboots^2)*~gboots))' |(4.*meanc((rboots^3)*~gboots))'; Sigma=(sboots~(rboots^2-1)~(rboots^3)~(rboots^4-3))' *(sboots~(rboots^2-1)~(rboots^3)~(rboots^4-3))/Tboots; O1=(G1*cov~eye(3))*Sigma*(G1*cov~eye(3))'; mv=(meanc(r^2)-1)|meanc(r^3)|(meanc(r^4)-3); T1=T*mv'*(invpd(O1))*mv; retp(T1); endp; proc (0) = nprtestingeb(tp,pr); print ""; print "*******************************************************************"; print "LM Normality test bootstrapped, r^2-1 included "; print " H estim. value Chi 3 / 3 "; print tp[1,.]'~pr; print ""; print "*******************************************************************"; endp; /* AUTOCORRELATION TEST */ proc atestingboots(r,rboots,sboots,gboots,matK1); local Res,K1,LK1; K1=1; LK1=rows(matK1); Res=zeros(1,LK1); for ak1 (1,LK1,1); K1=matK1[ak1]; Res[1,ak1]=cdfchic(atestboots(r,rboots,sboots,gboots,K1),K1); endfor; retp(Res); endp; proc atestboots(r,rboots,sboots,gboots,K1); local O2,G2,T2,T,RR,guboots,k,Sigma,help,Tboots,cov; T=rows(r); Tboots=rows(rboots); k=cols(sboots); cov=invpd((sboots'*sboots)/Tboots); G2=zeros(K1,k); RR=zeros(K1,1); guboots=zeros(Tboots-K1,K1); help=0; for t1 (1,(Tboots-K1),1); help=((rboots[t1+K1])*rboots[t1+K1-1]); for i (2,K1,1); help=help~((rboots[t1+K1])*rboots[t1+K1-i]); endfor; guboots[t1,.]=help; endfor; Sigma=(sboots[K1+1:Tboots,.]~guboots)'*(sboots[K1+1:Tboots,.]~guboots)/(Tboots-K1); for i1 (1,K1,1); G2[i1,.]=(meanc((rboots[1+i1:Tboots]*~gboots[1:Tboots-i1,.]) +(rboots[1:Tboots-i1]*~gboots[1+i1:Tboots,.])))'; RR[i1]=meanc(r[1+i1:T]*~r[1:T-i1]); endfor; O2=(G2*cov~eye(K1))*Sigma*(G2*cov~eye(K1))'; T2=(T-K1)*(RR')*(invpd(O2))*RR; retp(T2); endp; proc (0) = aprtestingeb(tp,pr,matK1); local LK1,label; LK1=rows(matK1); print ""; print "*******************************************************************"; print "LM AUTOCORRELATION TESTS with bootstrapped omega"; label="LM Autocorrelation tests"; call pTEST(label,matK1,tp[1:LK1,.],pr); print ""; print "*******************************************************************"; endp; /* CONDITIONAL HETEROSKEDASTICITY TEST */ proc htestingboots(r,rboots,sboots,gboots,matK2); local Res,K2,LK2; K2=1; LK2=rows(matK2); Res=zeros(1,LK2); for ak2 (1,LK2,1); K2=matK2[ak2]; Res[1,ak2]=cdfchic(htestboots(r,rboots,sboots,gboots,K2),K2); endfor; retp(Res); endp; proc htestboots(r,rboots,sboots,gboots,K2); local O3,G3,T3,T,RRR,k,Sigma,guboots,Tboots,cov; T=rows(r); Tboots=rows(rboots); k=cols(sboots); cov=invpd((sboots'*sboots)/Tboots); G3=zeros(K2,k); RRR=zeros(K2,1); guboots=zeros(Tboots-K2,K2); for i2 (1,K2,1); G3[i2,.]=2.*meanc( (rboots[K2+1-i2:Tboots-i2]^2-1)*~rboots[K2+1:Tboots]*~gboots[K2+1:Tboots,.]+ ((rboots[K2+1:Tboots]^2-1)*~rboots[K2+1-i2:Tboots-i2])*~gboots[K2+1-i2:Tboots-i2,.])'; RRR[i2,1]=meanc((r[K2+1:T]^2-1)*~(r[K2+1-i2:T-i2]^2-1)); guboots[.,i2]=(rboots[K2+1:Tboots]^2-1)*~(rboots[K2+1-i2:Tboots-i2]^2-1); endfor; Sigma=(sboots[K2+1:Tboots,.]~guboots)'*(sboots[K2+1:Tboots,.]~guboots)/(Tboots-K2); O3=(G3*cov~eye(K2))*Sigma*(G3*cov~eye(K2))'; T3=T*(RRR')*(invpd(O3))*RRR; retp(T3); endp; proc (0) = hprtestingeb(tp,pr,matK2); local LK2,label; LK2=rows(matK2); print ""; print "*******************************************************************"; print "LM HETEROSCEDASTICITY TEST with bootstrapped omega"; label="LM Heteroscedasticity test"; call pTEST(label,matK2,tp[1:LK2,.],pr); print ""; print "*******************************************************************"; endp; proc (0) = pTEST(label,matK,tp,pr); local LK; LK=rows(matK); for i(1,rows(matK),3); print label; if LK>=3; print " Lags " ftos(matK[i],"%*.*lf",1,0) " Lags " ftos(matK[i+1],"%*.*lf",1,0) " Lags " ftos(matK[i+2],"%*.*lf",1,0) " Chi lags"; print tp[i:i+2,.]'~pr; print ""; endif; if LK==2; print " Lags " ftos(matK[i],"%*.*lf",1,0) " Lags " ftos(matK[i+1],"%*.*lf",1,0) " Chi lags"; print tp[i:i+1,.]'~pr; print ""; endif; if LK==1; print " Lags " ftos(matK[i],"%*.*lf",1,0) " Chi lags"; print tp[i,.]'~pr; print ""; endif; LK=LK-3; endfor; endp; 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; proc g(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; proc cumulg(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;