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
1. Chen and Hong 2012 Econometrica 7.37 KB
Copy Edit Web IDE Raw Blame History
ymhong authored 2021-03-05 00:34 . first
/* PROGRAM TO COMPUTE TEST STATISTICS OF smooth structural change */
/*September, 2006*/
new;
/* NOTE: NUMBER OF ITERATIONS */
ITR=1;
ITRmax=19;
/* N= sample size */
/* For application of real data, N should be changed */
load datam[708,15]=D:\structurechange\data\monthly47.txt;
n=rows(datam)-1;
y=datam[2:n+1,15];
Bmax=9999;
/*bandwidth*/
h=1/sqrt(12)*n^(-1/5);
proc k0(u);
local Y;
Y=(1/2).*(abs(u).<=1.0);
retp(Y);
endp;
/*seed number*/
iseed1=2555130;
iseed2=2660427;
bound=ceil(n*h)+1;
/*second component of Chow's mean*/
j=-bound;
ca22=0;
ha22=0;
ca12=0;
cb12=0;
ha12=0;
hb12=0;
cb22=0;
hb22=0;
ca13=0;
ca23=0;
do until j==2*bound+1;
kj=j/(h*N);
/*the part related to kernel
Note: 2 due to symmetry
Need to add k(0)^2*/
ca12=ca12+(1-abs(j)/N)*(k0(kj)^2);
/*boundary correction*/
proc FSfun1(u);
local Y;
y=k0(2*u+kj);
retp(y);
endp;
Let U1={1,
-1};
_INTREC=0;
_INTORD=40;
ca13k=intquad1(&fsfun1,u1);
/**2 is from 2h, the other 2 is from symmetry!!*/
ca13=ca13+h*(1-abs(j)/n)*k0(kj)*ca13k;
proc FSfun2(u);
local Y;
y=k0(u).*k0(u+kj);
retp(y);
endp;
Let U5={1,
-1};
_INTREC=0;
_INTORD=40;
if j>=1;
cb12k=intquad1(&fsfun2,u5);
cb12=cb12+(1-j/N)*(2*k0(kj)-cb12k)^2;
hb12=hb12+(1-j/N)*cb12k^2;
endif;
j=j+1;
endo;
itr=1;
cpvalue1=zeros(itrmax,1);
hpvalue1=cpvalue1;
cpvalue2=cpvalue1;
hpvalue2=cpvalue1;
do until ITR==ITRmax+1;
if itr==15;
x=datam[1:n,1]~datam[1:n,8];
elseif itr==16;
x=datam[1:n,1]~datam[1:n,3];
elseif itr==17;
x=datam[1:n,1]~datam[1:n,3]~datam[1:n,8];
elseif itr==18;
x=datam[1:n,1]~datam[1:n,8]~datam[1:n,11]~datam[1:n,12];
elseif itr==19;
x=datam[1:n,1:3]~datam[1:n,5:8]~datam[1:n,10:14];
else;
x=datam[1:n,itr];
endif;
intercept=ones(N,1);
data=intercept~x;
/*dimension of regressor*/
d=cols(data);
/*ols estimation*/
xx=data'*data;
aols=invpd(xx)*data'*y;
yhatols=data*aols;
resid=y-data*aols;
/*sum of square residuals of OLS*/
resols=(y-data*aols).^2;
varols=meanc(resols);
ssrols=sumc(resols);
/*nonparametric estimation*/
/*nonparametric estimation*/
xleft=rev(data[2:bound,.]);
yleft=rev(y[2:bound,.]);
xright=rev(data[(n-bound+1):n-1,.]);
yright=rev(y[(n-bound+1):n-1,.]);
datasym=xleft|data|xright;
ysym=yleft|y|yright;
apoly=lpolydb(datasym,ysym,h);
yhatpoly=sumc((data.*apoly)');
respoly=y-yhatpoly;
respoly=respoly-meanc(respoly);
respolys=(y-yhatpoly).^2;
/*nonparametric SSR*/
varpoly=meanc(respolys);
ssrpoly=sumc(respolys);
/*chow's test*/
chow=h^(1/2)*(ssrols-ssrpoly);
/*Hausman's test*/
hau=h^(1/2)*sumc((yhatols-yhatpoly).^2);
/*sample average of XtXt'*/
m=data'*data/n;
invm=invpd(m);
sigma=data'*(data.*respolys)/n;
ca11=h^(-1/2)*(2*k0(0)+h)*d*varpoly;
ca21=h^(-1/2)*(2*k0(0)+h)*sumc(diag(sigma*invm));
/*ca1 is Chow's mean assuming Xt is independent*/
ca1=ca11-(ca12+ca13)/N*h^(-3/2)*varpoly*d;
cb1=4/(N*h)*varpoly^2*d*cb12;
ha1=(ca12+ca13)/N*h^(-3/2)*varpoly*d;
hb1=4/(N*h)*varpoly^2*d*hb12;
ca2=ca21-(ca12+ca13)/N*h^(-3/2)*sumc(diag(sigma*invm));
cb2=4/(N*h)*sumc(diag(invm*sigma*invm*sigma))*cb12;
ha2=(ca12+ca13)/N*h^(-3/2)*sumc(diag(sigma*invm));
hb2=4/(N*h)*sumc(diag(invm*sigma*invm*sigma))*hb12;
/*standarized Chow's test*/
/*chow1: chow~*/
chow1=(chow-ca1)./sqrt(cb1);
/*standarized Hausman's test*/
hau1=(hau-ha1)./sqrt(hb1);
chow2=(chow-ca2)./sqrt(cb2);
hau2=(hau-ha2)./sqrt(hb2);
cboot1=zeros(bmax,1);
hboot1=cboot1;
cboot2=cboot1;
hboot2=hboot1;
b=1;
do until b==bmax+1;
/*wild bootstrap*/
a=(1+sqrt(5))/2;
bi=rbinom(n,1,a/sqrt(5),1);
eb=(1-a)*respoly.*bi+a*respoly.*(1-bi);
yb=yhatols+eb;
aolsb=invpd(xx)*data'*yb;
yhatolsb=data*aolsb;
/*sum of square residuals of OLS*/
resolsb=(yb-data*aolsb).^2;
varolsb=meanc(resolsb);
ssrolsb=sumc(resolsb);
/*nonparametric estimation*/
bound=ceil(n*h)+1;
xleft=rev(data[2:bound,.]);
ybleft=rev(yb[2:bound,.]);
xright=rev(data[(n-bound+1):n-1,.]);
ybright=rev(yb[(n-bound+1):n-1,.]);
datasym=xleft|data|xright;
ybsym=ybleft|yb|ybright;
apolyb=lpolydb(datasym,ybsym,h);
yhatpolyb=sumc((data.*apolyb)');
respolyb=yb-yhatpolyb;
respolyb=respolyb-meanc(respolyb);
respolysb=(yb-yhatpolyb).^2;
/*nonparametric SSR*/
varpolyb=meanc(respolysb);
ssrpolyb=sumc(respolysb);
/*chow's test*/
/*chow's test*/
chowb=h^(1/2)*(ssrolsb-ssrpolyb);
/*Hausman's test*/
haub=h^(1/2)*sumc((yhatolsb-yhatpolyb).^2);
/*sample average of XtXt'*/
m=data'*data/n;
invm=invpd(m);
sigmab=data'*(data.*respolysb)/n;
ca11b=h^(-1/2)*(2*k0(0)+h)*d*varpolyb;
ca21b=h^(-1/2)*(2*k0(0)+h)*sumc(diag(sigmab*invm));
/*ca1 is Chow's mean assuming Xt is independent*/
ca1b=ca11b-(ca12+ca13)/N*h^(-3/2)*varpolyb*d;
cb1b=4/(N*h)*varpolyb^2*d*cb12;
ha1b=(ca12+ca13)/N*h^(-3/2)*varpolyb*d;
hb1b=4/(N*h)*varpolyb^2*d*hb12;
ca2b=ca21b-(ca12+ca13)/N*h^(-3/2)*sumc(diag(sigmab*invm));
cb2b=4/(N*h)*sumc(diag(invm*sigmab*invm*sigmab))*cb12;
ha2b=(ca12+ca13)/N*h^(-3/2)*sumc(diag(sigmab*invm));
hb2b=4/(N*h)*sumc(diag(invm*sigmab*invm*sigmab))*hb12;
/*standarized Chow's test*/
/*chow1: chow~*/
chow1b=(chowb-ca1b)./sqrt(cb1b);
/*standarized Hausman's test*/
hau1b=(haub-ha1b)./sqrt(hb1b);
chow2b=(chowb-ca2b)./sqrt(cb2b);
hau2b=(haub-ha2b)./sqrt(hb2b);
cboot1[b,.]=chow1b;
cboot2[b,.]=chow2b;
/*Hausman's test*/
hboot1[b,.]=hau1b;
hboot2[b,.]=hau2b;
b=b+1;
endo;
cpvalue1[itr,.]=meanc(chow1.<cboot1);
cpvalue2[itr,.]=meanc(chow2.<cboot2);
hpvalue1[itr,.]=meanc(hau1.<hboot1);
hpvalue2[itr,.]=meanc(hau2.<hboot2);
if ITR==1;
data0=N~ITR~chow1~chow2~hau1~hau2;
else;
data0=data0|(N~ITR~chow1~chow2~hau1~hau2);
endif;
ITR;
ITR=ITR+1;
endo;
means=meanc(data0[.,3:6]);
vars=stdc(data0[.,3:6]).^2;
output file=z:\structuralchange\empirical\smooth1947month.out on;
n;
"nonparametric residuals";
"bootstrap: b=9999";
"chow";
cpvalue1;
"chow robust";
cpvalue2;
"Hausman";
hpvalue1;
"hausman robust";
hpvalue2;
output off;
proc lpolydb(x,y,h);
local dx,nx,kk,z,in,tc,tm,tmx,tgrid,jj,w2,xw,xwx,beta,ixwx,w;
/*dimension of x*/
nx = rows(x) ;
dx=cols(x);
z=ones(nx,dx);
tgrid = zeros(nx,1) ;
in=1;
do until in==nx+1;
if in<bound;
tgrid[in,1]=-(bound-(in-1));
else;
tgrid[in,1]=in-bound+1;
endif;
in=in+1;
endo;
in=1;
do until in==nx+1;
tc = (tgrid - tgrid[in,1])/n ;
/*weighting function*/
w=k0(tc/h);
tm = ones(nx,1) ;
tmx=tm*~x;
/*generate weighting matrix*/
w2 = (w)*ones(1,dx);
xw = w2.*tmx ;
xwx = tmx' * xw ;
ixwx = inv(xwx) ;
beta = ixwx * xw' * y ;
z[in,.] = beta[1:dx,1]';
in=in+1;
endo;
z=z[bound:(bound+n-1),.];
retp(z);
endp;
proc dbinom(x,p,n);
local xx,d;
if p==0;
d = (x.==0);
elseif p==1;
d = (x.==n);
else;
xx = x.*(x.==TRUNC(x)).*(x.>=0).*(x.<=n);
d = EXP(LNFACT(n)-LNFACT(xx)-LNFACT(n-xx)+xx*LN(p)+(n-xx)*LN(1-p));
endif;
retp(d.*(x.==TRUNC(x)).*(x.>=0).*(x.<=n));
endp;
/********************************************************************/
proc pbinom(x,p,n);
local xx,d;
xx = TRUNC(x).*(x.>=0);
if MAXC(MAXC(xx)) ge 0;
d = CUMSUMC(dbinom(SEQA(0,1,MAXC(MAXC(xx))+1),p,n));
d = RESHAPE(d[xx+1],ROWS(x),COLS(x));
else;
d = ZEROS(ROWS(x),COLS(x));
endif;
retp(d.*(x.>=0));
endp;
/********************************************************************/
proc ibinom(h,p,n);
local x,qi;
x = SEQA(0,1,n+1);
qi = SUMC((VECR(h).>pbinom(x,p,n)')')+1;
retp(RESHAPE(x[qi],ROWS(h),COLS(h)));
endp;
/********************************************************************/
proc rbinom(r,c,p,n);
retp(RESHAPE(SUMC((RNDU(r*c,1).>(0|pbinom(SEQA(0,1,n),p,n))')')-1,r,c));
endp;
R
1
https://gitee.com/ymhong/code.git
git@gitee.com:ymhong/code.git
ymhong
code
Computer Codes
master

Search