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
2. Chen and Hong 2012 ET 16.40 KB
Copy Edit Web IDE Raw Blame History
ymhong authored 2021-03-05 00:34 . first
12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109
/* PRG0: Program for Markov Testing */
/* omnibus testing M(0,0) */
/* 11-3-2003 */
new;
/* seed numbers */
iseed1=2555130;
iseed2=2660427;
/* N= sample size */
/* For applican of real data, N should be changed */
N=100;
Nu=30;
n0=100;
Iv0=2000;
ep=10^(-4);
ITRBmax=100;
ITR=1;
ITRmax=500;
/* U00 and V00 are grid points of (u,v) */
/* We use normal variates for (u,v) */
U00=rndns(1,0.5*Nu,iseed1);
u00=u00~-u00;
V00=rndns(1,0.5*Nu,iseed2);
v00=v00~-v00;
/* RANGE OF PRELIMINARY BANDWIDTH Pbar*/
Pmax=16;
Pmin=4;
/* note: The total number of preliminary lags : 5 to 30 */
/* P0max is the maximum lag order to be integrated numerically.
It can be set to N-1. But for the kernel with bounded support,
it can be must less than N. This will speed up the computation */
P0max=50;
bandh=0;
/* GAUSSIAN KERNEL N(0,1) */
const2=1/sqrt(2*pi);
con1=pi;
proc GK(u);
local Y;
Y=const2*exp(-0.5*u.*u);
retp(Y);
endp;
/*
/*higher order kernel*/
proc hk(u);
local y;
y= 0.5*(3-u.^2).*pdfn(u);
retp(y);
*/
/* KERNEL FUNCTIONS */
/* BARTLETT KERNEL */
proc Wbar(u);
local Y;
Y=(1-u).*(u.<=1);
retp(Y);
endp;
/* DANIELL KERNEL */
proc Wdan(u);
local Y;
Y=sin(Con1*u)./(Con1*u);
retp(Y);
endp;
/* PARZEN KERNEL */
proc Wpar(u);
local X,Y;
X=(Con1/6)*u;
Y = (1-6.*X^2+6.*X^3).*(X.<=0.5) + (2.*(1-X)^3).*(X.>0.5).*(X.<=1.0);
retp(Y);
endp;
/* QS KERNEL */
proc Wqs(u);
local Y;
Y=sqrt(5/3)*Con1*u;
Y=(3./Y^2).*(sin(Y)./Y - cos(Y));
retp(Y);
endp;
/*E kernel*/
proc ekernel(u);
local Y;
Y=0.75.*((1-u.*u).>=0).*(1-u.*u);
retp(Y);
endp;
/*nonparametric estimation*/
proc ccfm(x,u00);
local y,t,yreal,yimag,ccfreal,ccfimag,xnew,ymi,ymr,ccfrealsort,ccfimagsort;
yreal=cos(x*u00);
yimag=sin(x*u00);
xnew=zeros(n-1,1);
xnew=x[1:n-1,.];
ccfreal=zeros(n,nu);
ccfimag=zeros(n,nu);
ymr=yreal[2:n,.];
ymi=yimag[2:n,.];
ccfreal[1,.]=cos(x[1,1]*u00);
ccfimag[1,.]=sin(x[1,1]*u00);
t=1;
do until t==nu+1;
ccfrealsort=lpoly(ymr[.,t],xnew,bandh,xnew);
ccfimagsort=lpoly(ymi[.,t],xnew,bandh,xnew);
ccfreal[2:n,t]=ccfrealsort;
ccfimag[2:n,t]=ccfimagsort;
t=t+1;
endo;
y=x~ccfreal~ccfimag;
retp(y);
endp;
/*nonparamtric conditional distribution function*/
proc cdf(x,v0);
local y,cdfhat,xnew, temp,t;
y=zeros(n-1,iv0);
y=(x[2:n,.].<=v0);
xnew=zeros(n-1,1);
xnew=x[1:n-1,.];
cdfhat=zeros(n,iv0);
cdfhat[1,.]=(x[1,1].<=v0);
t=1;
do until t==iv0+1;
cdfhat[2:n,t]=lpoly(y[.,t],xnew,bandh,xnew);
t=t+1;
endo;
retp(cdfhat);
endp;
proc Kbar(n0,p0);
local j,Y,c0;
c0=cols(p0);
Y=zeros(N0,c0);
j=1;
do until j==N0+1;
Y[j,.]=Wbar(j./p0);
j=j+1;
endo;
retp(Y);
endp;
proc Kdan(n0,p0);
local j,Y,c0;
c0=cols(p0);
Y=zeros(n0,c0);
j=1;
do until j==n0+1;
Y[j,.]=Wdan(j./p0);
j=j+1;
endo;
retp(Y);
endp;
proc Kpar(n0,p0);
local j,Y,c0;
c0=cols(p0);
Y=zeros(n0,c0);
j=1;
do until j==n0+1;
Y[j,.]=Wpar(j./p0);
j=j+1;
endo;
retp(Y);
endp;
proc Kqs(n0,p0);
local j,Y,c0;
c0=cols(p0);
Y=zeros(n0,c0);
j=1;
do until j==n0+1;
Y[j,.]=Wqs(j./p0);
j=j+1;
endo;
retp(Y);
endp;
/* JUST USE THE BARTLETT KERNEL FOR AUTOMATIC P */
/* P00 is the preliminary bandwidth vector */
P00=zeros(1,Pmax);
j=1;
do until j==Pmax+1;
P00[1,j]=Pmin+j;
j=j+1;
endo;
Bar=Kbar(P0max,P00);
Bar2=Bar^2;
itr=1;
do until ITR==ITRmax+1;
/*simiulate regime swching**/
itrb=0;
e=rndns(N+N0,1,iseed1);
s=rndns(n+n0,1,iseed1);
y=s;
t=1;
do until t==N+N0+1;
if t==1;
ru=rndu(1,1);
if (ru<0.6);
s[t]=0;
y[t]=e[1];
else;
s[t]=1;
y[t]=e[1];
endif;
else;
ru=rndu(1,1);
if (s[t-1]==0);
if(ru<0.6);
s[t]=1;
y[t]=-0.3*y[t-1]+e[t];
else;
s[t]=0;
y[t]=0.7*y[t-1]+e[t];
endif;
else;
if(ru<0.6);
s[t]=0;
y[t]=0.7*y[t-1]+e[t];
else;
s[t]=1;
y[t]=-0.3*y[t-1]+e[t];
endif;
endif;
endif;
t=t+1;
endo;
y=y[(n0+1):(n0+n),.];
/*scale y on support [0,1]*/
/***
y=1/(1+exp(-y));
***/
/*y=(y-minc(y))/(maxc(y)-minc(y));*/
step=(maxc(y)-minc(y))/iv0;
V0_sm=seqa(minc(y),step,iv0)';
do until itrb==itrbmax+1;
if itrb==0;
x=y;
bandh=stdc(x)*n^(-1/5);
/* Iv0 x 1 cdf estimates */
g0_cdf_sm=cdf(x,v0_sm);
/*normalization*/
g0_cdf_sm=(g0_cdf_sm-minc(g0_cdf_sm)')./(maxc(g0_cdf_sm)'-minc(g0_cdf_sm)');
/* Below are for BOOTSTRAP Residuals */
else;
/* Generate U[0,1] randam variables */
X=zeros(N,1);
IB=rndus(N,1,iseed2);
i=1;
do until i==N+1;
Y=1000*(g0_cdf_sm[i,.]-IB[i]);
Y=abs(Y);
Iy=indexcat(Y',minc(Y'));
X[i,1]=v0_sm[1,Iy[1]];
i=i+1;
endo;
/* FOR BOOTSTRAP: corresponding to "If B0 >= 1" */
endif;
/* SCALE {X} ON SUPPORT [0,1] */
/*X=(X-minc(X))/(maxc(X)-minc(X));*/
/***simulate MA(1)***/
/**
yy=zeros(500+n,1);
estep=zeros(2,1);
estep[1,1]=rndns(1,1,iseed1);
simindex=1;
ap=0.5;
do until simindex==500+n+1;
estep[2,1]=rndns(1,1,iseed1);
yy[simindex,1]=estep[2,1]+ap*estep[1,1];
estep[1,1]=estep[2,1];
simindex=simindex+1;
endo;
x=yy[501:(500+n),.];
***/
/*y=rndns(N,1,iseed1);
X=Y;
x=x/stdc(x);*/
stdy=stdc(x);
meany=meanc(x);
U00=rndns(1,0.5*Nu,iseed1).*stdy+meany;
u00=u00~-u00;
V00=rndns(1,0.5*Nu,iseed2).*stdy+meany;
v00=v00~-v00;
bandh=stdc(x)*N^(-1/5);
ccf=ccfm(x,u00);
x=ccf[.,1];
realccf=ccf[.,2:(Nu+1)];
imagccf=ccf[.,(Nu+2):(2*Nu+1)];
/*use true*/
/*realccf=ones(n,1)*exp(-u00.^2/2);
imagccf=zeros(n,nu);*/
imagccf2=-imagccf;
realccf2=realccf;
/* THIS IS THE ORIGINAL RANDOM VARIABLE X(t) */
/* t=j+1 to N */
proc X1(j);
local Z;
Z=X[j+1:N,.];
retp(Z);
endp;
/* t=1 to N-j */
proc X2(j);
local Z;
Z=X[1:N-j,.];
retp(Z);
endp;
/* Data CCF */
proc CX1(u,j);
local Z1,Z2,r,q;
r=rows(u);
Z1=X1(j);
/*
Z2=Z1[.,1].*u[1,.];
q=2; do until q==r+1;
Z2=Z2+Z1[.,q].*u[q,.];
q=q+1; endo;*/
z2=z1*u;
Z2=cos(Z2);
retp(Z2);
endp;
proc CX2(u,j);
local Z1,Z2,r,q;
r=rows(u);
Z1=X2(j);
/*Z2=Z1[.,1].*u[1,.];
q=2; do until q==r+1;
Z2=Z2+Z1[.,q].*u[q,.];
q=q+1; endo;*/
z2=z1*u;
Z2=cos(Z2);
retp(Z2);
endp;
proc SX1(u,j);
local Z1,Z2,r,q;
r=rows(u);
Z1=X1(j);
/*Z2=Z1[.,1].*u[1,.];
q=2; do until q==r+1;
Z2=Z2+Z1[.,q].*u[q,.];
q=q+1; endo;*/
z2=z1*u;
Z2=sin(Z2);
retp(Z2);
endp;
proc SX2(u,j);
local Z1,Z2,r,q;
r=rows(u);
Z1=X2(j);
/*Z2=Z1[.,1].*u[1,.];
q=2;do until q==r+1;
Z2=Z2+Z1[.,q].*u[q,.];
q=q+1;endo;*/
z2=z1*u;
Z2=sin(Z2);
retp(Z2);
endp;
proc Re_CCF1(u,j);
local Z;
Z=realccf[j+1:N,.];
retp(Z);
endp;
/*I add this for u00[.,q]*/
proc Re_CCFq1(u,j,q);
local Z;
Z=realccf[j+1:N,q];
retp(Z);
endp;
/* t=1 to N-j */
proc Re_CCF2(u,j);
local Z;
Z=realccf[1:N-j,.];
retp(Z);
endp;
/*for -u00*/
proc Re_CCF3(u,j);
local Z;
Z=realccf2[1:N-j,.];
retp(Z);
endp;
proc IM_CCF1(u,j);
local Z;
Z=imagccf[j+1:N,.];
retp(Z);
endp;
proc IM_CCFq1(u,j,q);
local Z;
Z=imagccf[j+1:N,q];
retp(Z);
endp;
/* t=1 to N-j */
proc IM_CCF2(u,j);
local Z;
Z=imagccf[1:N-j,.];
retp(Z);
endp;
/*for -u00*/
proc IM_CCF3(u,j);
local Z;
Z=imagccf2[1:N-j,.];
retp(Z);
endp;
/* exp(u*X(t))- Model CCF */
proc CM1(u,j);
local Z;
Z=CX1(u,j)-Re_CCF1(u,j);
retp(Z);
endp;
/*I add this for u00[.,q]*/
proc CMq1(u,j,q);
local Z;
Z=CX1(u,j)-Re_CCFq1(u,j,q);
retp(Z);
endp;
proc CM2(u,j);
local Z;
Z=CX2(u,j)-Re_CCF2(u,j);
retp(Z);
endp;
/*for -u00*/
proc CM3(u,j);
local z;
z=cx2(u,j)-Re_CCF3(u,j);
retp(Z);
endp;
proc SM1(u,j);
local Z;
Z=SX1(u,j)-Im_CCF1(u,j);
retp(Z);
endp;
/*I add this for u00[.,q]*/
proc SMq1(u,j,q);
local Z;
Z=SX1(u,j)-Im_CCFq1(u,j,q);
retp(Z);
endp;
proc SM2(u,j);
local Z;
Z=SX2(u,j)-Im_CCF2(u,j);
retp(Z);
endp;
proc SM3(u,j);
local Z;
Z=SX2(u,j)-Im_CCF3(u,j);
retp(Z);
endp;
proc Re_xxsgm(u,v,j);
local a1,b1,a2,b2,Z;
a1=CX1(u,j); a1=a1-meanc(a1)';
b1=SX1(u,j); b1=b1-meanc(b1)';
a2=CX2(v,j); a2=a2-meanc(a2)';
b2=SX2(v,j); b2=b2-meanc(b2)';
Z=meanc(a1.*a2-b1.*b2);
retp(Z);
endp;
proc Im_xxsgm(u,v,j);
local a1,b1,a2,b2,Z;
a1=CX1(u,j); a1=a1-meanc(a1)';
b1=SX1(u,j); b1=b1-meanc(b1)';
a2=CX2(v,j); a2=a2-meanc(a2)';
b2=SX2(v,j); b2=b2-meanc(b2)';
Z=meanc(a1.*b2+b1.*a2);
retp(Z);
endp;
proc Re_mmsgm(u,v,j);
local a1,a2,b1,b2,Z;
a1=CM1(u,j);a1=a1-meanc(a1)';
b1=SM1(u,j);b1=b1-meanc(b1)';
a2=CM3(v,j);a2=a2-meanc(a2)';
b2=SM3(v,j);b2=b2-meanc(b2)';
Z=meanc(a1.*a2-b1.*b2);
retp(Z);
endp;
proc Im_mmsgm(u,v,j);
local a1,a2,b1,b2,Z;
a1=CM1(u,j);a1=a1-meanc(a1)';
b1=SM1(u,j);b1=b1-meanc(b1)';
a2=CM3(v,j);a2=a2-meanc(a2)';
b2=SM3(v,j);b2=b2-meanc(b2)';
Z=meanc(a1.*b2+b1.*a2);
retp(Z);
endp;
/* Below is used to determine lag order that from data */
/* COMPUTE THE INTEGRALS */
/* WRSXX will be used to select lag orders for M(0,0) */
WRSMM=zeros(1,P0max+1);
WRSXX=zeros(1,P0max+1);
WISMM=WRSMM;
WISXX=WRSXX;
j=0;
do until j==P0max+1;
/* For lag dedetrmination */
/* Compute the integral of real Sigma_mm(u,-u).*G(u) and Sigma_xx(v,-v,j).*GK(v)*/
WRSMM[1,j+1]=meanc(Re_MMsgm(U00,-U00,j));
WRSXX[1,j+1]=meanc(Re_XXsgm(V00,-V00,j));
/*
WISMM[1,j+1]=meanc(Im_MMsgm(U00,-U00,j));
WISXX[1,j+1]=meanc(Im_XXsgm(V00,-V00,j));
*/
/* NOTE: I have Checked INTQUAD1(&WIXsgm, V1) is indeed zero */
j=j+1;
endo;
/* TEST FOR MARTINGALE DIFFERENCE */
/* TEST 1: Based on Cov( exp(u'X(t))-CCF(u), exp(ivX(t-j)) */
/* MS1: M(0,0) */
/* Smx00s(u,v,j) */
proc Smx00s(u,v,j);
local a1,a2,b1,b2, RZ,IZ;
a1=CM1(u,j);a1=a1-meanc(a1)';
b1=SM1(u,j);b1=b1-meanc(b1)';
a2=CX2(v,j);a2=a2-meanc(a2)';
b2=SX2(v,j);b2=b2-meanc(b2)';
RZ=meanc(a1.*a2-b1.*b2);
IZ=meanc(a1.*b2+b1.*a2);
retp(RZ.*RZ+IZ.*IZ);
endp;
MS1=zeros(1,P0max);
M1_1=MS1;
M1_2=MS1;
j=0;
do until j==P0Max+1;
if j==0;
M1_10 = (1-j/N)*WRSmm[1,j+1].*WRSxx[1,j+1];
else;
MS1[1,j]= (1-j/N)*meanc(Smx00s(U00,V00,j));
M1_2[1,j] = j*j*MS1[1,j];
M1_1[1,j] = (1-j/N)*WRSmm[1,j+1]*WRSxx[1,j+1];
endif;
j=j+1;
endo;
MS1=N*MS1;
M1_1=2*M1_1*Bar2+M1_10;
M1_2=2*M1_2*Bar2;
MM1=N*2*M1_2./M1_1;
P1_1=1.1447142*MM1^(1/3);
P1_2=1.4017333*MM1^0.2;
P1_3=1.3934819*MM1^0.2;
P1_4=1.4223747*MM1^0.2;
lnn=ln(n);
P1_1=(P1_1.>lnn).*P1_1+(P1_1.<=lnn).*lnn;
P1_2=(P1_2.>lnn).*P1_2+(P1_2.<=lnn).*lnn;
P1_3=(P1_3.>lnn).*P1_3+(P1_3.<=lnn).*lnn;
P1_4=(P1_4.>lnn).*P1_4+(P1_4.<=lnn).*lnn;
/* THIS IS 1 by 4*Pmax vector */
P0=P1_1~P1_2~P1_3~P1_4;
/*p0=p00~p00~p00~p00;*/
/* PART II */
/* COMPUTE TEST STATISTICS */
N1=P0max;
Ka=Kbar(n1,P1_1)~Kdan(n1,P1_2)~Kpar(n1,P1_3)~Kqs(n1,P1_4);
/*ka=kbar(n1,p00)~kdan(n1,p00)~kpar(n1,p00)~kqs(n1,p00);*/
Ka2=Ka.*Ka;
Ka4=Ka2.*Ka2;
Cka=sumc(Ka2);
Cka=Cka';
Dka=sumc(Ka4);
Dka=Dka';
/* means/variance of M(0,0) */
/* For M(0,0) */
C00=zeros(P0max,1);
D00=zeros(P0max,P0max);
/* This helps to reduce computational cost */
Pdim=CEIL(maxc(P0'))+1;
proc C00j(u,v,j);
local a1,b1,a2,b2,Z1,Z2,Z;
a1=CM1(u,j); a1=a1-meanc(a1)';
b1=SM1(u,j); b1=b1-meanc(b1)';
a2=CX2(v,j); a2=a2-meanc(a2)';
b2=SX2(v,j); b2=b2-meanc(b2)';
Z1=a1.*a1+b1.*b1;
Z2=a2.*a2+b2.*b2;
Z=Z1.*Z2;
Z=meanc(Z);
retp(Z);
endp;
proc D00JL(u1,v1,u2,v2,j,L,q);
local a1,b1,a2,b2,a3,b3,a4,b4, Z0, Z,RZ, IZ, R12,I12,R34,I34, r,p;
a1=CM1(u1,j); a1=a1-meanc(a1)';
b1=SM1(u1,j); b1=b1-meanc(b1)';
a2=CMq1(u2,j,q); a2=a2-meanc(a2)';
b2=SMq1(u2,j,q); b2=b2-meanc(b2)';
a3=CX2(v1,j); a3=a3-meanc(a3)';
b3=SX2(v1,j); b3=b3-meanc(b3)';
Z0=X[j-L+1:N-L,.];
r=rows(v2);
/*Z=Z0[.,1].*v2[1,.];
p=2; do until p==r+1;
Z=Z+Z0[.,p].*v2[p,.];
p=p+1; endo;*/
z=z0*v2;
a4=cos(Z); a4=a4-meanc(a4)';
b4=sin(Z); b4=b4-meanc(b4)';
R12=a1.*a2-b1.*b2;
I12=a1.*b2+b1.*a2;
R34=a3.*a4-b3.*b4;
I34=a3.*b4+b3.*a4;
RZ=meanc(R12.*R34-I12.*I34);
IZ=meanc(R12.*I34+I12.*R34);
retp(RZ.*RZ+IZ.*IZ);
endp;
j=1;
do until j==Pdim+1;
C00[j,1]=meanc(C00j(U00,V00,j));
L=1;
do until L==j+1;
D00tmp=zeros(Nu,1);
q=1;
do until q==Nu+1;
D00tmp[q,1]=meanc(D00jL(U00,V00,U00[.,q],V00[.,q],j,L,q));
q=q+1;
endo;
D00[j,L]=meanc(D00tmp);
D00[L,j]=D00[j,L];
L=L+1;
endo;
j=j+1;
endo;
W00m=Ka2'D00*Ka2;
Wdim=rows(W00m);
W00=zeros(1,Wdim);
j=1;
do until j==Wdim+1;
W00[1,j]=W00m[j,j];
j=j+1;
endo;
/* TEST FOR MARTINGALE DIFFERENCE Mxx(0,0) */
MS1=MS1*Ka2;
MS1=(MS1-C00'Ka2)./sqrt(2*W00);
MS=MS1;
if ITRb==0;
/*original test statistics*/
statis=ms;
elseif itrb==1;
data0=MS;
else;
data0=data0|MS;
endif;
ITRb=ITRb+1;
endo;
/*compute bootstrap p-values*/
pv1=(statis-data0).<0;
pv1=meanc(pv1);
pv1=pv1';
pv1;
/*asymptotic p value*/
asympv=cdfnc(statis);
output file=c:\users\Public\Documents\p100markov_switching_100.dat on;
itr~pv1~asympv~statis;
output off;
/* INDEX for lag Pbar*/
Pbar0=zeros(Pmax,1);
j=1;
do until j==Pmax+1;
Pbar0[j,1]=Pmin+j;
j=j+1;
endo;
j=1;
do until j==1+1;
if j==1;
Ipbar=Pbar0;
Itest=ones(Pmax,1)*j;
else;
Ipbar=Ipbar|Pbar0;
Itest=Itest|ones(Pmax,1)*j;
endif;
j=j+1;
endo;
/* Results based on BARTLETT KERNEL */
j=1;
do until j==1+1;
if j==1;
data1=data0[1:itrbmax,(j-1)*4*Pmax+0*Pmax+1:(j-1)*4*Pmax+0*Pmax+Pmax];
else;
data1=data1~data0[1:itrbmax,(j-1)*4*Pmax+0*Pmax+1:(j-1)*4*Pmax+0*Pmax+Pmax];
endif;
j=j+1;
endo;
/* sort data to compute ecv */
j=1;
do until j==cols(data1)+1;
data1[.,j]=sortc(data1[.,j],1);
j=j+1;
endo;
Bar_ecv=data1[90,.]~data1[95,.]~data1[99,.];
CDim=cols(data1);
if itr==1;
BAR_rej=zeros(itrmax,3*cdim);
par_rej=zeros(itrmax,3*cdim);
bar_rej2=zeros(itrmax,cdim);
par_rej2=zeros(itrmax,cdim);
endif;
/* Test statistics values */
/* rejections based on asymptotic critical values */
statisb=statis[.,1:pmax];
statisp=statis[.,2*pmax+1:2*pmax+pmax];
bar_rej2[itr,.]=meanc((statisb-data1).<0)';
BAR_REJ[itr,.]=(statisb.>BAR_ECV[1,1:CDim])~(statisb.>BAR_ECV[1,CDim+1:2*CDim])~(statisb.>BAR_ECV[1,2*CDim+1:3*CDim]);
/* Results based on PARZEN KERNEL */
j=1;
do until j==1+1;
if j==1;
data2=data0[1:itrbmax,(j-1)*4*Pmax+2*Pmax+1:(j-1)*4*Pmax+2*Pmax+Pmax];
else;
data2=data2~data0[1:itrbmax,(j-1)*4*Pmax+2*Pmax+1:(j-1)*4*Pmax+2*Pmax+Pmax];
endif;
j=j+1;
endo;
par_rej2[itr,.]=meanc((statisp-data2).<0)';
/* sort data to compute ecv */
j=1;
do until j==cols(data2)+1;
data2[.,j]=sortc(data2[.,j],1);
j=j+1;
endo;
Par_ecv=data2[90,.]~data2[95,.]~data2[99,.];
/* Statistics Values */
pAR_REJ[itr,.]=(statisp.>pAR_ECV[1,1:CDim])~(statisp.>pAR_ECV[1,CDim+1:2*CDim])~(statisp.>pAR_ECV[1,2*CDim+1:3*CDim]);
bar_rej[itr,.];
par_rej[itr,.];
itr=itr+1;
endo;
bar_rejmean=meanc(bar_rej);
par_rejmean=meanc(par_rej);
/*output file=c:\Markov\ar5_200.out on;
bar_rejmean;
par_rejmean;
output off;*/
/*****************************
/*****************
USE THIS for DATA UNDER THE ALTERNATIVE MODELS (e.g., A1(3))
******************/
/* compute p-values */
/* INDEX for lag Pbar*/
Pbar0=zeros(Pmax,1);
j=1;
do until j==Pmax+1;
Pbar0[j,1]=Pmin+j;
j=j+1;
endo;
j=1;
do until j==1+1;
if j==1;
Ipbar=Pbar0;
Itest=ones(Pmax,1)*j;
else;
Ipbar=Ipbar|Pbar0;
Itest=Itest|ones(Pmax,1)*j;
endif;
j=j+1;
endo;
/* Results based on BARTLETT KERNEL */
j=1;
do until j==1+1;
if j==1;
data1=data0[.,1+1+(j-1)*4*Pmax+0*Pmax+1:1+1+(j-1)*4*Pmax+0*Pmax+Pmax];
else;
data1=data1~data0[.,1+1+(j-1)*4*Pmax+0*Pmax+1:1+1+(j-1)*4*Pmax+0*Pmax+Pmax];
endif;
j=j+1;
endo;
/* Test statistics values */
/* rejections based on asymptotic critical values */
pvalue1=cdfnc(data1);
BAR_ACVREJ=meanc(pvalue1.<0.10)~meanc(pvalue1.<0.05)~meanc(pvalue1.<0.01);
/* rejection using ecv */
/* Here 2 = Bartlett and Parzen, 3 = 10%, 5% and 1% levels */
CDim=cols(data1);
/*load ECV[1,2*3*CDim]=c:\diagonalk\hlp0_ecv500.out;*/
/*********
BAR_ECV=ECV[1,0*CDim+1:3*CDim];
PAR_ECV=ECV[1,3*CDim+1:2*3*CDim];
BAR_ECVREJ=meanc(data1.>BAR_ECV[1,1:CDim])~meanc(data1.>BAR_ECV[1,CDim+1:2*CDim])~meanc(data1.>BAR_ECV[1,2*CDim+1:3*CDim]);
************/
/* Results based on PARZEN KERNEL */
j=1;
do until j==1+1;
if j==1;
data2=data0[.,1+1+(j-1)*4*Pmax+2*Pmax+1:1+1+(j-1)*4*Pmax+2*Pmax+Pmax];
else;
data2=data2~data0[.,1+1+(j-1)*4*Pmax+2*Pmax+1:1+1+(j-1)*4*Pmax+2*Pmax+Pmax];
endif;
j=j+1;
endo;
/* compute ecv */
j=1;
do until j==cols(data2)+1;
data2[.,j]=sortc(data2[.,j],1);
j=j+1;
endo;
/*Par_ecv=data2[900,.]~data2[950,.]~data2[990,.];*/
/* Statistics Values */
pvalue2=cdfnc(data2);
PAR_ACVREJ=meanc(pvalue2.<0.10)~meanc(pvalue2.<0.05)~meanc(pvalue2.<0.01);
/*PAR_ECVREJ=meanc(data2.>PAR_ECV[1,1:CDim])~meanc(data2.>PAR_ECV[1,CDim+1:2*CDim])~meanc(data2.>PAR_ECV[1,2*CDim+1:3*CDim]);*/
output file=c:\Markov\ibmtest_y.out on;
Itest~Ipbar~BAR_ACVREJ~PAR_ACVREJ;
output off;
output file=c:\Markov\ibmtest_statistics_y.out on;
data1~data2;
output off;
*************************************/
proc lpoly(y,x,h,targmat);
local m,k,yhat,n,l,cont,kernarg,ki,estden,estreg;
n=rows(y);
targmat=targmat';
m=n;
k=cols(x);
yhat=zeros(m,1);
/* START OF MAIN LOOP */
cont=1;
do while cont==1;
kernarg=(x-targmat)/h;
ki=0.5*(3-kernarg.^2).*pdfn(kernarg);
estden=sumc(ki);
estreg=sumc(ki.*y);
if abs(estden) > (10^(-20));
yhat=estreg./estden;
cont=0;
else;
h=h*1.1;
cont=1;
endif;
endo;
retp(yhat);
endp;
R
1
https://gitee.com/ymhong/code.git
git@gitee.com:ymhong/code.git
ymhong
code
Computer Codes
master

Search

113223 674803ea 1850385 170725 2838fb2a 1850385