/*
** Lesson 14.1: Unit Roots Test with Structural Break
** Data Series: Personal Consumption
*/
use gpe2;

load z[67,3] = gpe\usyc87.txt;
year = z[2:67,1];
y = z[2:67,2];
c = z[2:67,3];

/* select one variable to work on */
x = y;

/* difference the data if needed */
diff = 0;
j = 1;
do until j > diff;
	x = x-lagn(x,1);
	j = j+1;
endo;

x1 = lagn(x,1);
dx = x-x1;
data = packr(year~x~x1~dx); @ data truncated @
year = data[.,1];
x = data[.,2];
x1 = data[.,3];
dx = data[.,4];
t=rows(x);
trend = seqa(1,1,t);

call reset;
_rstat = 1;
_dlags = 1; @ augmented terms if needed @
_names = {"dx","du","trend","x1"};
_print = -1;

/* Model III */

tstat = miss(zeros(t,1),0);
tb = _dlags+1;
do until tb > t-1;
    du = (trend.>tb);
    dt = (trend-tb).*du;
    call estimate(dx,dt~trend~x1);
    k=rows(__b);
    tstat[tb] = __b[k-1]./sqrt(__vb[k-1,k-1]);
    /*    
    print trend[tb]~year[tb]~(__b[k-1]./sqrt(__vb[k-1,k-1]));
    */
    tb = tb+1;
endo;
print tstat;

library pgraph;
graphset;
call xy(trend,abs(tstat));

_print = 1;
tb = minindc(tstat);
du = (trend.>tb);
dt = (trend-tb).*du;
call estimate(dx,dt~trend~x1);

end;