code
/
1. Chen and Hong 2012 Econometrica


/* 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;