Dear professor,
/*sequential procedure to select number of breaks*/
new;
format /ld 6,3;
cls;
load mat[36, 6] = "Pakistan_data.txt";
Year = mat[.,1];
CO = mat[.,2]; GI = mat[.,3]; LY = mat[.,4]; EC = mat[.,5];
GL = mat[.,6];
bigt = rows(mat);
y = CO;
/* Change Y and ZZ according to what you want to test */
ZZ = GI~LY~EC~GL; /* do not include (or add) a column of ones in ZZ matrix */
z = ZZ; /* Here, I select the first column of ZZ as a regressor to test the code as it is,
change the definition of z according to your needs */
x = 0;
T = rows(y);
y1 = y[4:T-2,1];
z1 = z[4:T-2,.];
deltaz = z1-z[3:T-3,.];
lz = z[3:T-3,.]-z[2:T-4,.];
llz = z[2:T-4,.]-z[1:T-5,.];
fz = z[5:T-1,.]-z1;
ffz = z[6:T,.]-z[5:T-1,.];
x = deltaz~lz~llz~fz~ffz;
y = y1;
z = z1;
q = 4; @ number of regressors whose coefficients are allowed to change: intercept and slope @
p = 2; @ number of first differenced regressors used in DOLS estimation @
m = 3; @ maximum number of structural changes allowed @
eps1 = .15; @ Value of the trimming (in percentage) @
h = int(eps1*(bigt)); @ minimal length of a segment (h >= q). @
/* the following are options if p > 0.
----------------------------------- */
fixb = 0; @ set to 1 if use fixed initial values for beta @
betaini = 0; @ if fixb = 1, load the initial value of beta. @
maxi = 20; @ maximum number of iterations for the nonlinear
procedure to obtain global minimizers. @
printd = 0; @ set to 1 if want the output from the iterations to be printed. @
eps = 0.001; @ criterion for the convergence. @
/*--------------------------------- */
doseq = 1;
call pbreak(bigt,y,z,deltaz,q,m,h,eps1,p,fixb,x,q,eps,maxi,fixb,betaini,printd,doseq);
#include break_seq.src @ set the path to where you store the file brcode.src @
end;
*************************************break_seq*****************************************
proc(0)=pbreak(bigt,y,z,deltaz,q,m,h,eps1,p,fixb,x,q,eps,maxi,fixb,betaini,printd,doseq);
local datevec,bigvec,global,supfl,ndat,maic,mbic,mlwz,ssrzero,nbreak,i;
local dateseq,ftest,wftest,cv,repartda,zz,siglev,j,nbr,datese,cvm,reparv,ii,cvall,size,seqtests;
print "The options chosen are:";
print "h = " h;
print "eps1 = " eps1;
print "The maximum number of breaks is: " m;
print "********************************************************";
if doseq==1;
seqtests=zeros(m+1,1);
print "The following options are used:";
print "h = " h;
print "bigt: " bigt;
{nbr,datese}=sequa(m,2,q,h,bigt,z,y,
x,p,eps1);
print"the seqtest selects the number of breaks at" nbr;
endif;
endp;
@*******************************************************************@
proc(1)=pftest(y,z,g,q,bigt,datevec,x,p);
local ftest,zbar,gammahat,deltahat,ssr0,ssrk,wbar,w,e,eb,eb1,eb2,ef,ef1,ef2,u,ub,uf,uhat,uhatnull,fsample,b;
local te,ae,ae1,ae2,se1,se2,ee,se,ad,a2,eband,jb,jband,kern,sig,lam,j,delta,omega,a,uu,ai,y1,z1,w1;
if g==1;
fsample=zeros(bigt,1);
for b(eps1*bigt,(1-eps1)*bigt,1);
zbar=pzbar(z,1,b);
wbar=zbar~x;
w=z~x;
gammahat=olsqr(y,wbar);
deltahat=olsqr(y,w);
ssr0=(y-w*deltahat)'*(y-w*deltahat);
ssrk=(y-wbar*gammahat)'*(y-wbar*gammahat);
uhat=y-wbar*gammahat;
uhatnull=y-w*deltahat;
u=uhat;
e = u;
te = bigt;
eb = e[1:te-1,.];
ef = e[2:te,.];
ae=eb'*ef/(eb'*eb);
ee=zeros(te-1,1);
ee[.,1]=ef[.,1]-eb[.,1]*(ae);
se=sqrt(inv(te-1)*ee[.,1]'*ee[.,1]);
ad=((se)^2)/(1-ae)^4;
a2=(4*(ae*se/(1-ae)^4)^2/ad);
eband = 1.3221*((a2*te)^.2);
jb = seqa(1,1,te-1)/eband;
jband = jb*1.2*pi;
kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
sig = uhatnull'*uhatnull;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (uhatnull[1:te-j,.]'*uhatnull[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam' ;
omega[1,1]=omega[1,1]/te;
fsample[b,1]=(ssr0-ssrk)/(g*omega[1,1]);
endfor;
fsample=fsample[eps1*bigt:(1-eps1)*bigt,1];
ftest=maxc(fsample);
else;
{zbar}=pzbar(z,g,datevec[1:g,g]);
wbar=zbar~x;
w=z~x;
gammahat=olsqr(y,wbar);
deltahat=olsqr(y,w);
ssrk=(y-wbar*gammahat)'*(y-wbar*gammahat);
ssr0=(y-w*deltahat)'*(y-w*deltahat);
uhat=y-wbar*gammahat;
uhatnull=y-w*deltahat;
u=uhat;
e = u;
te = bigt;
eb = e[1:te-1,.];
ef = e[2:te,.];
ae=eb'*ef/(eb'*eb);
ee=zeros(te-1,1);
ee[.,1]=ef[.,1]-eb[.,1]*(ae);
se=sqrt(inv(te-1)*ee[.,1]'*ee[.,1]);
ad=((se)^2)/(1-ae)^4;
a2=(4*(ae*se/(1-ae)^4)^2/ad);
eband = 1.3221*((a2*te)^.2);
jb = seqa(1,1,te-1)/eband;
jband = jb*1.2*pi;
kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
sig = uhatnull'*uhatnull;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (uhatnull[1:te-j,.]'*uhatnull[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam' ;
omega[1,1]=omega[1,1]/te;
ftest=(ssr0-ssrk)/(g*omega[1,1]);
endif;
retp(ftest);
endp;
proc(3)=dating(y,z,h,m,q,bigt);
@This is the main procedure which calculates the break points that
globally minimizes the SSR. It returns optimal dates and associated SSR
for all numbers of breaks less than or equal to m.@
local datevec,optdat,optssr,dvec,i,ssrmin,datx,j1,ib,jlast;
local global,vecssr,jb,xx,bigvec;
datevec=zeros(m,m);
optdat=zeros(bigt,m);
optssr=zeros(bigt,m);
dvec=zeros(bigt,1);
global=zeros(m,1);
bigvec=zeros(bigt*(bigt+1)/2,1);
@
segment to construct the vector of SSR of dimension T*(T+1)/2
@
i=1;
do while i <= bigt-h+1;
{vecssr}=ssr(i,y,z,h,bigt);
bigvec[(i-1)*bigt+i-(i-1)*i/2:i*bigt-(i-1)*i/2,1]=vecssr[i:bigt,1];
i=i+1;
endo;
@
Section that applies the dynamic programming algorithm to look for the
optimal breaks
The first case is when m = 1. Here the dynamic programming algorithm is not
needed and one call to the parti(.) procedure is enough.
@
if m == 1;
{ssrmin,datx}=parti(1,h,bigt-h,bigt,bigvec,bigt);
datevec[1,1]=datx;
global[1,1]=ssrmin;
else;
@ when m > 1, a dynamic programming algorithm is used.
The first step is to obtain the optimal one break partitions for all
possible ending dates from 2h to T-mh+1.
The optimal dates are stored in a vector optdat.
The associated ssr are stored in a vector optssr.
@
j1=2*h; @First loop. Looking for the@
do while j1 <= bigt; @optimal one break partitions@
{ssrmin,datx}=parti(1,h,j1-h,j1,bigvec,bigt);@for break dates between h and@
optssr[j1,1]=ssrmin; @T-h. j1 is the last date of the@
optdat[j1,1]=datx; @segments.@
j1=j1+1;
endo;
global[1,1]=optssr[bigt,1];
datevec[1,1]=optdat[bigt,1];
@
When this is done the algorithm looks for optimal 2,3,... breaks
partitions. The index used is ib.
@
ib=2;
do while ib <= m;
if ib == m; @if we have reached the highest number of breaks
considered, only one segment is considered, that
which ends at the last date of the sample.@
jlast=bigt;
jb=ib*h; /* date of the break */
do while jb <=jlast-h;
dvec[jb,1] = optssr[jb,ib-1]+bigvec[(jb+1)*bigt-jb*(jb+1)/2,1];
jb=jb+1;
endo;
optssr[jlast,ib]=minc(dvec[ib*h:jlast-h,1]);
optdat[jlast,ib]=(ib*h-1)+minindc(dvec[ib*h:jlast-h,1]);
else;
@if we have not reached the highest number of breaks
considered, we need to loop over the last date of
the segment, between (ib+1)*h and T.@
jlast=(ib+1)*h;
do while jlast <= bigt;
jb=ib*h; /* date of the break */
do while jb <=jlast-h;
dvec[jb,1] = optssr[jb,ib-1]+bigvec[jb*bigt-jb*(jb-1)/2+jlast-jb,1];
jb=jb+1;
endo;
optssr[jlast,ib]=minc(dvec[ib*h:jlast-h,1]);
optdat[jlast,ib]=(ib*h-1)+minindc(dvec[ib*h:jlast-h,1]);
jlast=jlast+1;
endo;
endif;
@At each time we loop the results with ib breaks are retrieved
and printed@
datevec[ib,ib]=optdat[bigt,ib];
i=1;
do while i <= ib-1;
xx=ib-i;
datevec[xx,ib]=optdat[datevec[xx+1,ib],xx];
i=i+1;
endo;
global[ib,1]=optssr[bigt,ib];
ib=ib+1;
endo;
endif; /*closing the if for the case m >1*/
retp(global,datevec,bigvec);
endp;
@********************************************************************@
proc(1)=ssr(start,y,z,h,last);
local vecssr,delta1,delta2,inv1,inv2,invz,res,v,f,r;
@
This procedure computes recursive residuals from a data set that
starts at date "start" and ends at date "last". It returns a
vector of sum of squared residuals (SSR) of length last-start+1 (stored for
convenience in a vector of length T).
start: starting entry of the sample used.
last: ending date of the last segment considered
y: dependent variable
z: matrix of regressors of dimension q
h: minimal length of a segment
@
/* initialize the vectors */
vecssr=zeros(last,1);
/* initialize the recursion with the first h data points */
inv1=inv(z[start:start+h-1,.]'z[start:start+h-1,.]);
delta1=inv1*(z[start:start+h-1,.]'y[start:start+h-1,1]);
res=y[start:start+h-1,1]-z[start:start+h-1,.]*delta1;
vecssr[start+h-1,1]=res'res;
/* loop to construct the recursive residuals and update the SSR */
r=start+h;
do while r <= last;
v=y[r,1]-z[r,.]*delta1;
invz=inv1*z[r,.]';
f=1+z[r,.]*invz;
delta2=delta1+(invz*v)/f;
inv2=inv1-(invz*invz')/f;
inv1=inv2;
delta1=delta2;
vecssr[r,1]=vecssr[r-1,1]+v*v/f;
r=r+1;
endo;
retp(vecssr);
endp;
@********************************************************************@
proc(2)=parti(start,b1,b2,last,bigvec,bigt);
local k,dvec,j,ssrmin,dx,ini,jj;
@
procedure to obtain an optimal one break partitions for a segment that
starts at date start and ends at date last. It return the optimal break
and the associated SSR.
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
@
dvec=zeros(bigt,1);
ini=(start-1)*bigt-(start-2)*(start-1)/2+1;
j=b1;
do while j<=b2;
jj=j-start;
k=j*bigt-(j-1)*j/2+last-j;
dvec[j,1]=bigvec[ini+jj,1]+bigvec[k,1];
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@********************************************************************@
proc(1)=ssrnul(y,zz);
local delta,ssrn;
@Computes the SSR under the null of no break@
delta=olsqr(y,zz);
ssrn=(y-zz*delta)'(y-zz*delta);
retp(ssrn);
endp;
@********************************************************************@
proc(1)=pzbar(zz,m,bb);
local nt,q1,zb,i;
@procedure to construct the diagonal partition of z with m break at
dates b.@
nt=rows(zz);
q1=cols(zz);
zb=zeros(nt,(m+1)*q1);
zb[1:bb[1,1],1:q1]=zz[1:bb[1,1],.];
i=5;
do while i <= m;
zb[bb[i-1,1]+1:bb[i,1],(i-1)*q1+1:i*q1]=zz[bb[i-1,1]+1:bb[i,1],.];
i=i+1;
endo;
zb[bb[m,1]+1:nt,m*q1+1:(m+1)*q1]=zz[bb[m,1]+1:nt,.];
retp(zb);
endp;
@********************************************************************@
proc plambda(b,m,bigt);
local lambda,k;
@procedure that construct a diagonal matrix of dimension m+1 with ith
entry (T_i - T_i-1)/T.@
lambda=zeros(m+1,m+1);
lambda[1,1]=b[1,1]/bigt;
k=2;
do while k <= m;
lambda[k,k]=(b[k,1]-b[k-1,1])/bigt;
k=k+1;
endo;
lambda[m+1,m+1]=(bigt-b[m,1])/bigt;
retp(lambda);
endp;
@********************************************************************@
proc psigmq(res,b,q,i,nt);
local sigmat,kk,sigmq;
@procedure that computes a diagonal matrix of dimension i+1 with ith
entry the estimate of the variance of the residuals for segment i.
@
sigmat=zeros(i+1,i+1);
sigmat[1,1]=res[1:b[1,1],1]'res[1:b[1,1],1]/b[1,1];
kk = 2;
do while kk <= i;
sigmat[kk,kk]=res[b[kk-1,1]+1:b[kk,1],1]'
res[b[kk-1,1]+1:b[kk,1],1]/(b[kk,1]-b[kk-1,1]);
kk=kk+1;
endo;
sigmat[i+1,i+1]=res[b[i,1]+1:nt,1]'res[b[i,1]+1:nt,1]/(nt-b[i,1]);
retp(sigmat);
endp;
@********************************************************************@
proc getcv2(signif,eps1);
local cv;
if eps1 == .15;
if signif == 2;
cv={12.11 13.79 15.25 16.38 17.03 17.70 17.98 18.25 18.33 18.53};
endif;
endif;
retp(cv);
endp;
@********************************************************************@
proc(2)=partione(b1,b2,last,vssr,vssrev);
local dvec,j,ssrmin,dx;
@
procedure to obtain an optimal one break partitions for a segment that
starts at date start and ends at date last. It return the optimal break
and the associated SSR.
Procedure used with the sequential method when the T*(T+1)/2 vector
of SSR is not computed.
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
@
dvec=zeros(last,1);
j=b1;
do while j<=b2;
dvec[j,1]=vssr[j,1]+vssrev[last-j,1];
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@********************************************************************@
proc(2)=sequa(m,signif,q,h,bigt,z,y,x,p,eps1);
local cv,ssrmin,datx,dv,nbreak,ds,ftest,nseg,is,maxf,newseg;
local ftestv,dv2,length,vssr,vssrev,dum,dbreak;
@This procedure applies the sequential procedure to obtain the number of
breaks and the break dates.
The current version only allow pure structural changes. This will
be generalized.
Here m is the maximum number of breaks allowed.@
dv=zeros(m+2,1);
dv2=zeros(m+2,1);
ftestv=zeros(m+1,1);
@input of the critical values of the supF(l+1|l) test@
cv=getcv2(signif,eps1);
dv[1,1]=0; @the first date of a segment is set to 0@
@Obtain the first break and test if significant.@
if p == 0;
{vssrev}=ssr(1,rev(y),rev(z),h,bigt);
{vssr}=ssr(1,y,z,h,bigt);
{ssrmin,datx}=partione(h,bigt-h,bigt,vssr,vssrev);
else;
{ssrmin,datx}=onebp(y,z,x,h,1,bigt);
endif;
dv[2,1]=datx;
ftest=pftest(y,z,1,q,bigt,dv[2,1],x,p);
if ftest < 12.11;
nbreak=0;
dv[2,1]=0;
retp(nbreak,0);
else;
print "The first break found is at: " datx;
nbreak=1;
nseg=2;
dv[nseg+1,1]=bigt; @ set the last date of the segment to T @
endif;
do while nseg <= m;
ds=zeros(nseg+1,1);
ftestv=zeros(nseg+1,1);
is=1;
do while is <= nseg;
length=dv[is+1,1]-dv[is,1];
if length >= 2*h;
if p == 0;
{vssr}=ssr(1,y[dv[is,1]+1:dv[is+1,1],1],
z[dv[is,1]+1:dv[is+1,1],.],h,length);
{vssrev}=ssr(1,rev(y[dv[is,1]+1:dv[is+1,1],1]),
rev(z[dv[is,1]+1:dv[is+1,1],.]),h,length);
{dum,ds[is,1]}=partione(h,length-h,length,vssr,vssrev);
ftestv[is,1]=pftest(y[dv[is,1]+1:dv[is+1,1],1],z[dv[is,1]+1:dv[is+1,1],.],
1,q,length,ds[is,1],0,p);
else;
{dum,ds[is,1]}=onebp(y,z,x,h,dv[is,1]+1,dv[is+1,1]);
ds[is,1]=ds[is,1]-dv[is,1];
ftestv[is,1]=pftest(y[dv[is,1]+1:dv[is+1,1],1],z[dv[is,1]+1:dv[is+1,1],.],
1,q,length,ds[is,1],x[dv[is,1]+1:dv[is+1,1],.],p);
endif;
else;
ftestv[is,1]=0.0;
endif;
is=is+1;
endo;
maxf=maxc(ftestv[1:nseg,1]);
if maxf < cv[q-1,nseg];
retp(nbreak,dv[2:nbreak+1,1]);
else;
newseg=maxindc(ftestv[1:nseg,1]);
print "The next break found is at: " ds[newseg,1]+dv[newseg,1];
dv[nseg+2,1]=ds[newseg,1]+dv[newseg,1];
nbreak=nbreak+1;
dv2=sortc(dv[2:nseg+2,1],1);
dv[1,1]=0;
dv[2:nseg+2,1]=dv2;
endif;
nseg=nseg+1;
endo;
print "The sequential procedure has reached the upper limit";
retp(nbreak,dv[2:nbreak+1,1]);
endp;
@********************************************************************@
proc(3)=nldat(y,z,x,h,m,p,q,bigt,fixb,eps,maxi,betaini,printd);
local qq,zz,xbar,zbar,teta,delta1,beta1,ssr1,i,mi,length,globnl;
local datenl,teta1,ssrn,global,datevec,bigvec;
global=zeros(m,1);
globnl=zeros(m,1);
datevec=zeros(m,m);
datenl=zeros(m,m);
mi=1;
do while mi <= m;
if printd == 1;
print "Model with " mi "breaks";
else;
endif;
if fixb == 0;
@******************************************************************@
/*Segment to initialize the beta vector. We apply treat the model
as one of pure structural change and apply the dynamic programming
algorithm to get the break dates.
*/
qq=p+q;
zz=x~z;
{globnl,datenl,bigvec}=dating(y,zz,h,mi,qq,bigt);
xbar=pzbar(x,mi,datenl[1:mi,mi]);
zbar=pzbar(z,mi,datenl[1:mi,mi]);
teta=olsqr(y,zbar~xbar);
delta1=teta[1:q*(mi+1),1];
beta1=olsqr(y-zbar*delta1,x);
/* Calculate SRR. */
ssr1=(y-x*beta1-zbar*delta1)'(y-x*beta1-zbar*delta1);
if printd == 1;
print "The iterations are initialized with: ";
print "Break dates: " datenl[1:mi,mi]';
print "Delta1: " delta1;
print "Beta1: " beta1;
print "SSR1: " ssr1;
else;
endif;
@******************************************************************@
else;
@******************************************************************@
/*
if fix = 1, the initial value of beta is fixed. We initialize
SSR at some negative number to start the iterations*/
beta1=betaini;
ssr1=-5;
endif;
@*****************************************************************@
/* Starting the iterations. */
if printd == 1;
print "Output from the iterations";
else;
endif;
length=99999999.0;
i=1;
do while length > eps;
{globnl,datenl,bigvec}=dating(y-x*beta1,z,h,mi,q,bigt);
zbar=pzbar(z,mi,datenl[1:mi,mi]);
teta1=olsqr(y,x~zbar);
beta1=teta1[1:p,1];
delta1=teta1[p+1:p+q*(mi+1),1];
ssrn=(y-(x~zbar)*teta1)'(y-(x~zbar)*teta1);
/* Calculate overall SRR and check if significantly smaller.
*/
length=abs(ssrn-ssr1);
if printd == 1;
print "Iteration " i;
print "Break dates: " datenl[1:mi,mi]';
print "Delta1: " delta1;
print "Beta1: " beta1;
print "SSRN: " ssrn;
print "length" length;
else;
endif;
if i >= maxi;
print "The number of iterations has reached the upper limit";
goto itmax;
else;
i=i+1;
ssr1=ssrn;
global[mi,1]=ssrn;
datevec[1:mi,mi]=datenl[1:mi,mi];
endif;
endo;
mi=mi+1;
endo;
itmax:
retp(global,datevec,bigvec);
endp;
proc(2)=onebp(y,z,x,h,start,last);
local ssrind,i,zb,bb,rr,bdat;
/*procedure that computes an optimal one break partition in the case
of a partial structural by searching over all possible breaks. */
ssrind=9999999999;
i=h;
do while i <= last-start+1-h;
zb=pzbar(z[start:last,.],1,i);
bb=olsqr(y[start:last,1],x[start:last,.]~zb);
rr=(y[start:last,1]-(x[start:last,.]~zb)*bb)'
(y[start:last,.]-(x[start:last,.]~zb)*bb);
if rr < ssrind;
ssrind=rr;
bdat=i;
else;
endif;
i=i+1;
endo;
retp(ssrind,bdat+start-1);
endp;
************************************Pakistan_data**************************************
1975 -1,118973 -1,052683 9,597891 4,650136 3,289109
1976 -1,167021 -1,041287 9,61668 4,618381 3,287532
1977 -1,132484 -1,03002 9,624223 4,638634 3,306165
1978 -1,093896 -1,018877 9,670957 4,738784 3,307904
1979 -1,046148 -1,007858 9,677892 4,878953 3,312336
1980 -0,948549 -0,996959 9,746017 4,854814 3,296726
1981 -0,906538 -0,986177 9,794051 4,942109 3,277744
1982 -0,850456 -0,987967 9,830249 5,023859 3,269124
1983 -0,802308 -0,989757 9,868803 5,101533 3,274657
1984 -0,767842 -0,991553 9,89128 5,179186 3,269642
1985 -0,698645 -0,993353 9,93762 5,264911 3,295545
1986 -0,678145 -0,995153 9,964535 5,34484 3,305095
1987 -0,625194 -0,996959 10,00071 5,414522 3,317013
1988 -0,567751 -1,035637 10,04785 5,53168 3,321395
1989 -0,547182 -1,061316 10,07078 5,570929 3,473463
1990 -0,454928 -1,055553 10,08903 5,622542 3,485259
1991 -0,485021 -1,000578 10,11303 5,694621 3,511054
1992 -0,445589 -0,948471 10,16219 5,812414 3,55146
1993 -0,401418 -0,898942 10,15455 5,816492 3,583626
1994 -0,342549 -0,895263 10,16616 5,847647 3,604149
1995 -0,371346 -0,891598 10,19 5,884014 3,622598
1996 -0,284364 -0,916291 10,21283 5,888835 3,673784
1997 -0,305582 -0,916291 10,19891 5,899259 3,712659
1998 -0,298924 -0,916291 10,20006 5,846497 3,794483
1999 -0,295532 -0,916291 10,21192 5,879936 3,798949
2000 -0,260982 -0,903868 10,22953 5,923037 3,836792
2001 -0,268019 -0,891598 10,22504 5,936442 3,872089
2002 -0,239944 -1,090927 10,23267 5,95101 3,935485
2003 -0,222749 -1,290257 10,25588 6,010738 3,897321
2004 -0,145186 -1,275111 10,30286 6,052151 3,904995
2005 -0,131891 -1,26019 10,35263 6,122985 3,931617
2006 -0,088021 -1,245489 10,39115 6,173071 3,94753
2007 -0,039678 -1,231001 10,425 6,161466 3,964407
2008 0,002418 -1,197991 10,41942 6,077841 3,981168
2009 0,04378 -1,178087 10,43367 6,039584 3,997817
2010 0,083498 -1,156593 10,45491 5,986284 4,014193
6 Answers
0
After changing the commas to periods in Pakistan.txt, the first big problem is that your code assigns bigt to equal the number of rows of z right after z is created, which is 36. However, the code strips off 5 of the observations so that z only has 31 rows, but the code is trying to access 36 rows.
You need to either change bigt to equal the number of rows of z after it is decreased to 31, or you need to you need to index into z with bigt - 5 if bigt needs to remain the same size. I don't know which is correct for your analysis, but either way will allow this part of the code to run.
0
Many thanks for your response!
would you like please to indicate me which line should I change (your modification to the code).
I am grateful,
N.K
0
Around lines 27 to 37 or so of break_fseq, some of the first observations are trimmed of some of the variables. If you assign bigt to be equal to the number of rows of z or y1 (they are the same at this point, I think), this will solve your first problem.
However, before I did that I would make sure that I understood whether or not big was supposed to be equal to total number of observations in the original dataset or the number after trimming the observations.
0
Dear Prof.
0
Do you have an example dataset for this code that works successfully? If so, that would make it much, much easier to sort out.
0
Please find the original code for an example which the code run.
Many Thanks,
************************************break_fhseq****************************************
/*sequential procedure to select number of breaks*/
new;
format /ld 6,3;
cls;
vls = reshape(error(0),9,1);
vls[4] = 9999.99;
mat = xlsreadm("austria","b2:g142",1,vls);
gpdi=mat[.,6];
gnsave=mat[.,5];
x=0;
bigt=rows(mat);
y=gpdi;
z=gnsave;
T=rows(y);
y1=y[4:T-2,1];
z1=z[4:T-2,1];
deltaz=z1-z[3:T-3,1];
lz=z[3:T-3,1]-z[2:T-4,1];
llz=z[2:T-4,1]-z[1:T-5,1];
fz=z[5:T-1,1]-z1;
ffz=z[6:T,1]-z[5:T-1,1];
x=deltaz~lz~llz~fz~ffz;
y=y1;
z=z1;
bigt=rows(y);
z=z~ones(bigt,1);
q=2; @number of regressors whose coefficients are allowed to change: intercept and slope@
p=5; @number of first differenced regressors used in DOLS estimation@
m=5; @maximum number of structural changes allowed@
eps1=.15; @Value of the trimming (in percentage)@
h=int(eps1*(bigt)); @minimal length of a segment (h >= q).@
/* the following are options if p > 0.
----------------------------------- */
fixb=0; @set to 1 if use fixed initial values for beta@
betaini={0,0,0}; @if fixb=1, load the initial value of beta.@
maxi=20; @maximum number of iterations for the nonlinear
procedure to obtain global minimizers.@
printd=0; @set to 1 if want the output from the iterations
to be printed.@
eps=0.001; @criterion for the convergence.@
/*--------------------------------- */
doseq=1;
call pbreak(bigt,y,z,deltaz,q,m,h,eps1,p,fixb,x,q,eps,maxi,fixb,betaini,printd,doseq);
#include break_seq.src @set the path to where you store the file brcode.src@
end;
************************************break_seq*********************************************
proc(0)=pbreak(bigt,y,z,deltaz,q,m,h,eps1,p,fixb,x,q,eps,maxi,fixb,betaini,printd,doseq);
local datevec,bigvec,global,supfl,ndat,maic,mbic,mlwz,ssrzero,nbreak,i;
local dateseq,ftest,wftest,cv,repartda,zz,siglev,j,nbr,datese,cvm,reparv,ii,cvall,size,seqtests;
print "The options chosen are:";
print "h = " h;
print "eps1 = " eps1;
print "The maximum number of breaks is: " m;
print "********************************************************";
if doseq==1;
seqtests=zeros(m+1,1);
print "The following options are used:";
print "h = " h;
print "bigt: " bigt;
{nbr,datese}=sequa(m,2,q,h,bigt,z,y,
x,p,eps1);
print"the seqtest selects the number of breaks at" nbr;
endif;
endp;
@*******************************************************************@
proc(1)=pftest(y,z,g,q,bigt,datevec,x,p);
local ftest,zbar,gammahat,deltahat,ssr0,ssrk,wbar,w,e,eb,eb1,eb2,ef,ef1,ef2,u,ub,uf,uhat,uhatnull,fsample,b;
local te,ae,ae1,ae2,se1,se2,ee,se,ad,a2,eband,jb,jband,kern,sig,lam,j,delta,omega,a,uu,ai,y1,z1,w1;
if g==1;
fsample=zeros(bigt,1);
for b(eps1*bigt,(1-eps1)*bigt,1);
zbar=pzbar(z,1,b);
wbar=zbar~x;
w=z~x;
gammahat=olsqr(y,wbar);
deltahat=olsqr(y,w);
ssr0=(y-w*deltahat)'*(y-w*deltahat);
ssrk=(y-wbar*gammahat)'*(y-wbar*gammahat);
uhat=y-wbar*gammahat;
uhatnull=y-w*deltahat;
u=uhat;
e = u;
te = bigt;
eb = e[1:te-1,.];
ef = e[2:te,.];
ae=eb'*ef/(eb'*eb);
ee=zeros(te-1,1);
ee[.,1]=ef[.,1]-eb[.,1]*(ae);
se=sqrt(inv(te-1)*ee[.,1]'*ee[.,1]);
ad=((se)^2)/(1-ae)^4;
a2=(4*(ae*se/(1-ae)^4)^2/ad);
eband = 1.3221*((a2*te)^.2);
jb = seqa(1,1,te-1)/eband;
jband = jb*1.2*pi;
kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
sig = uhatnull'*uhatnull;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (uhatnull[1:te-j,.]'*uhatnull[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam' ;
omega[1,1]=omega[1,1]/te;
fsample[b,1]=(ssr0-ssrk)/(g*omega[1,1]);
endfor;
fsample=fsample[eps1*bigt:(1-eps1)*bigt,1];
ftest=maxc(fsample);
else;
{zbar}=pzbar(z,g,datevec[1:g,g]);
wbar=zbar~x;
w=z~x;
gammahat=olsqr(y,wbar);
deltahat=olsqr(y,w);
ssrk=(y-wbar*gammahat)'*(y-wbar*gammahat);
ssr0=(y-w*deltahat)'*(y-w*deltahat);
uhat=y-wbar*gammahat;
uhatnull=y-w*deltahat;
u=uhat;
e = u;
te = bigt;
eb = e[1:te-1,.];
ef = e[2:te,.];
ae=eb'*ef/(eb'*eb);
ee=zeros(te-1,1);
ee[.,1]=ef[.,1]-eb[.,1]*(ae);
se=sqrt(inv(te-1)*ee[.,1]'*ee[.,1]);
ad=((se)^2)/(1-ae)^4;
a2=(4*(ae*se/(1-ae)^4)^2/ad);
eband = 1.3221*((a2*te)^.2);
jb = seqa(1,1,te-1)/eband;
jband = jb*1.2*pi;
kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
sig = uhatnull'*uhatnull;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (uhatnull[1:te-j,.]'*uhatnull[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam' ;
omega[1,1]=omega[1,1]/te;
ftest=(ssr0-ssrk)/(g*omega[1,1]);
endif;
retp(ftest);
endp;
proc(3)=dating(y,z,h,m,q,bigt);
@This is the main procedure which calculates the break points that
globally minimizes the SSR. It returns optimal dates and associated SSR
for all numbers of breaks less than or equal to m.@
local datevec,optdat,optssr,dvec,i,ssrmin,datx,j1,ib,jlast;
local global,vecssr,jb,xx,bigvec;
datevec=zeros(m,m);
optdat=zeros(bigt,m);
optssr=zeros(bigt,m);
dvec=zeros(bigt,1);
global=zeros(m,1);
bigvec=zeros(bigt*(bigt+1)/2,1);
@
segment to construct the vector of SSR of dimension T*(T+1)/2
@
i=1;
do while i <= bigt-h+1;
{vecssr}=ssr(i,y,z,h,bigt);
bigvec[(i-1)*bigt+i-(i-1)*i/2:i*bigt-(i-1)*i/2,1]=vecssr[i:bigt,1];
i=i+1;
endo;
@
Section that applies the dynamic programming algorithm to look for the
optimal breaks
The first case is when m = 1. Here the dynamic programming algorithm is not
needed and one call to the parti(.) procedure is enough.
@
if m == 1;
{ssrmin,datx}=parti(1,h,bigt-h,bigt,bigvec,bigt);
datevec[1,1]=datx;
global[1,1]=ssrmin;
else;
@ when m > 1, a dynamic programming algorithm is used.
The first step is to obtain the optimal one break partitions for all
possible ending dates from 2h to T-mh+1.
The optimal dates are stored in a vector optdat.
The associated ssr are stored in a vector optssr.
@
j1=2*h; @First loop. Looking for the@
do while j1 <= bigt; @optimal one break partitions@
{ssrmin,datx}=parti(1,h,j1-h,j1,bigvec,bigt);@for break dates between h and@
optssr[j1,1]=ssrmin; @T-h. j1 is the last date of the@
optdat[j1,1]=datx; @segments.@
j1=j1+1;
endo;
global[1,1]=optssr[bigt,1];
datevec[1,1]=optdat[bigt,1];
@
When this is done the algorithm looks for optimal 2,3,... breaks
partitions. The index used is ib.
@
ib=2;
do while ib <= m;
if ib == m; @if we have reached the highest number of breaks
considered, only one segment is considered, that
which ends at the last date of the sample.@
jlast=bigt;
jb=ib*h; /* date of the break */
do while jb <=jlast-h;
dvec[jb,1] = optssr[jb,ib-1]+bigvec[(jb+1)*bigt-jb*(jb+1)/2,1];
jb=jb+1;
endo;
optssr[jlast,ib]=minc(dvec[ib*h:jlast-h,1]);
optdat[jlast,ib]=(ib*h-1)+minindc(dvec[ib*h:jlast-h,1]);
else;
@if we have not reached the highest number of breaks
considered, we need to loop over the last date of
the segment, between (ib+1)*h and T.@
jlast=(ib+1)*h;
do while jlast <= bigt;
jb=ib*h; /* date of the break */
do while jb <=jlast-h;
dvec[jb,1] = optssr[jb,ib-1]+bigvec[jb*bigt-jb*(jb-1)/2+jlast-jb,1];
jb=jb+1;
endo;
optssr[jlast,ib]=minc(dvec[ib*h:jlast-h,1]);
optdat[jlast,ib]=(ib*h-1)+minindc(dvec[ib*h:jlast-h,1]);
jlast=jlast+1;
endo;
endif;
@At each time we loop the results with ib breaks are retrieved
and printed@
datevec[ib,ib]=optdat[bigt,ib];
i=1;
do while i <= ib-1;
xx=ib-i;
datevec[xx,ib]=optdat[datevec[xx+1,ib],xx];
i=i+1;
endo;
global[ib,1]=optssr[bigt,ib];
ib=ib+1;
endo;
endif; /*closing the if for the case m >1*/
retp(global,datevec,bigvec);
endp;
@********************************************************************@
proc(1)=ssr(start,y,z,h,last);
local vecssr,delta1,delta2,inv1,inv2,invz,res,v,f,r;
@
This procedure computes recursive residuals from a data set that
starts at date "start" and ends at date "last". It returns a
vector of sum of squared residuals (SSR) of length last-start+1 (stored for
convenience in a vector of length T).
start: starting entry of the sample used.
last: ending date of the last segment considered
y: dependent variable
z: matrix of regressors of dimension q
h: minimal length of a segment
@
/* initialize the vectors */
vecssr=zeros(last,1);
/* initialize the recursion with the first h data points */
inv1=inv(z[start:start+h-1,.]'z[start:start+h-1,.]);
delta1=inv1*(z[start:start+h-1,.]'y[start:start+h-1,1]);
res=y[start:start+h-1,1]-z[start:start+h-1,.]*delta1;
vecssr[start+h-1,1]=res'res;
/* loop to construct the recursive residuals and update the SSR */
r=start+h;
do while r <= last;
v=y[r,1]-z[r,.]*delta1;
invz=inv1*z[r,.]';
f=1+z[r,.]*invz;
delta2=delta1+(invz*v)/f;
inv2=inv1-(invz*invz')/f;
inv1=inv2;
delta1=delta2;
vecssr[r,1]=vecssr[r-1,1]+v*v/f;
r=r+1;
endo;
retp(vecssr);
endp;
@********************************************************************@
proc(2)=parti(start,b1,b2,last,bigvec,bigt);
local k,dvec,j,ssrmin,dx,ini,jj;
@
procedure to obtain an optimal one break partitions for a segment that
starts at date start and ends at date last. It return the optimal break
and the associated SSR.
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
@
dvec=zeros(bigt,1);
ini=(start-1)*bigt-(start-2)*(start-1)/2+1;
j=b1;
do while j<=b2;
jj=j-start;
k=j*bigt-(j-1)*j/2+last-j;
dvec[j,1]=bigvec[ini+jj,1]+bigvec[k,1];
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@********************************************************************@
proc(1)=ssrnul(y,zz);
local delta,ssrn;
@Computes the SSR under the null of no break@
delta=olsqr(y,zz);
ssrn=(y-zz*delta)'(y-zz*delta);
retp(ssrn);
endp;
@********************************************************************@
proc(1)=pzbar(zz,m,bb);
local nt,q1,zb,i;
@procedure to construct the diagonal partition of z with m break at
dates b.@
nt=rows(zz);
q1=cols(zz);
zb=zeros(nt,(m+1)*q1);
zb[1:bb[1,1],1:q1]=zz[1:bb[1,1],.];
i=2;
do while i <= m;
zb[bb[i-1,1]+1:bb[i,1],(i-1)*q1+1:i*q1]=zz[bb[i-1,1]+1:bb[i,1],.];
i=i+1;
endo;
zb[bb[m,1]+1:nt,m*q1+1:(m+1)*q1]=zz[bb[m,1]+1:nt,.];
retp(zb);
endp;
@********************************************************************@
proc plambda(b,m,bigt);
local lambda,k;
@procedure that construct a diagonal matrix of dimension m+1 with ith
entry (T_i - T_i-1)/T.@
lambda=zeros(m+1,m+1);
lambda[1,1]=b[1,1]/bigt;
k=2;
do while k <= m;
lambda[k,k]=(b[k,1]-b[k-1,1])/bigt;
k=k+1;
endo;
lambda[m+1,m+1]=(bigt-b[m,1])/bigt;
retp(lambda);
endp;
@********************************************************************@
proc psigmq(res,b,q,i,nt);
local sigmat,kk,sigmq;
@procedure that computes a diagonal matrix of dimension i+1 with ith
entry the estimate of the variance of the residuals for segment i.
@
sigmat=zeros(i+1,i+1);
sigmat[1,1]=res[1:b[1,1],1]'res[1:b[1,1],1]/b[1,1];
kk = 2;
do while kk <= i;
sigmat[kk,kk]=res[b[kk-1,1]+1:b[kk,1],1]'
res[b[kk-1,1]+1:b[kk,1],1]/(b[kk,1]-b[kk-1,1]);
kk=kk+1;
endo;
sigmat[i+1,i+1]=res[b[i,1]+1:nt,1]'res[b[i,1]+1:nt,1]/(nt-b[i,1]);
retp(sigmat);
endp;
@********************************************************************@
proc getcv2(signif,eps1);
local cv;
if eps1 == .15;
if signif == 2;
cv={12.11 13.79 15.25 16.38 17.03 17.70 17.98 18.25 18.33 18.53};
endif;
endif;
retp(cv);
endp;
@********************************************************************@
proc(2)=partione(b1,b2,last,vssr,vssrev);
local dvec,j,ssrmin,dx;
@
procedure to obtain an optimal one break partitions for a segment that
starts at date start and ends at date last. It return the optimal break
and the associated SSR.
Procedure used with the sequential method when the T*(T+1)/2 vector
of SSR is not computed.
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
@
dvec=zeros(last,1);
j=b1;
do while j<=b2;
dvec[j,1]=vssr[j,1]+vssrev[last-j,1];
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@********************************************************************@
proc(2)=sequa(m,signif,q,h,bigt,z,y,x,p,eps1);
local cv,ssrmin,datx,dv,nbreak,ds,ftest,nseg,is,maxf,newseg;
local ftestv,dv2,length,vssr,vssrev,dum,dbreak;
@This procedure applies the sequential procedure to obtain the number of
breaks and the break dates.
The current version only allow pure structural changes. This will
be generalized.
Here m is the maximum number of breaks allowed.@
dv=zeros(m+2,1);
dv2=zeros(m+2,1);
ftestv=zeros(m+1,1);
@input of the critical values of the supF(l+1|l) test@
cv=getcv2(signif,eps1);
dv[1,1]=0; @the first date of a segment is set to 0@
@Obtain the first break and test if significant.@
if p == 0;
{vssrev}=ssr(1,rev(y),rev(z),h,bigt);
{vssr}=ssr(1,y,z,h,bigt);
{ssrmin,datx}=partione(h,bigt-h,bigt,vssr,vssrev);
else;
{ssrmin,datx}=onebp(y,z,x,h,1,bigt);
endif;
dv[2,1]=datx;
ftest=pftest(y,z,1,q,bigt,dv[2,1],x,p);
if ftest < 12.11;
nbreak=0;
dv[2,1]=0;
retp(nbreak,0);
else;
print "The first break found is at: " datx;
nbreak=1;
nseg=2;
dv[nseg+1,1]=bigt; @set the last date of the segment to T@
endif;
do while nseg <= m;
ds=zeros(nseg+1,1);
ftestv=zeros(nseg+1,1);
is=1;
do while is <= nseg;
length=dv[is+1,1]-dv[is,1];
if length >= 2*h;
if p == 0;
{vssr}=ssr(1,y[dv[is,1]+1:dv[is+1,1],1],
z[dv[is,1]+1:dv[is+1,1],.],h,length);
{vssrev}=ssr(1,rev(y[dv[is,1]+1:dv[is+1,1],1]),
rev(z[dv[is,1]+1:dv[is+1,1],.]),h,length);
{dum,ds[is,1]}=partione(h,length-h,length,vssr,vssrev);
ftestv[is,1]=pftest(y[dv[is,1]+1:dv[is+1,1],1],z[dv[is,1]+1:dv[is+1,1],.],
1,q,length,ds[is,1],0,p);
else;
{dum,ds[is,1]}=onebp(y,z,x,h,dv[is,1]+1,dv[is+1,1]);
ds[is,1]=ds[is,1]-dv[is,1];
ftestv[is,1]=pftest(y[dv[is,1]+1:dv[is+1,1],1],z[dv[is,1]+1:dv[is+1,1],.],
1,q,length,ds[is,1],x[dv[is,1]+1:dv[is+1,1],.],p);
endif;
else;
ftestv[is,1]=0.0;
endif;
is=is+1;
endo;
maxf=maxc(ftestv[1:nseg,1]);
if maxf < cv[q-1,nseg];
retp(nbreak,dv[2:nbreak+1,1]);
else;
newseg=maxindc(ftestv[1:nseg,1]);
print "The next break found is at: " ds[newseg,1]+dv[newseg,1];
dv[nseg+2,1]=ds[newseg,1]+dv[newseg,1];
nbreak=nbreak+1;
dv2=sortc(dv[2:nseg+2,1],1);
dv[1,1]=0;
dv[2:nseg+2,1]=dv2;
endif;
nseg=nseg+1;
endo;
print "The sequential procedure has reached the upper limit";
retp(nbreak,dv[2:nbreak+1,1]);
endp;
@********************************************************************@
proc(3)=nldat(y,z,x,h,m,p,q,bigt,fixb,eps,maxi,betaini,printd);
local qq,zz,xbar,zbar,teta,delta1,beta1,ssr1,i,mi,length,globnl;
local datenl,teta1,ssrn,global,datevec,bigvec;
global=zeros(m,1);
globnl=zeros(m,1);
datevec=zeros(m,m);
datenl=zeros(m,m);
mi=1;
do while mi <= m;
if printd == 1;
print "Model with " mi "breaks";
else;
endif;
if fixb == 0;
@******************************************************************@
/*Segment to initialize the beta vector. We apply treat the model
as one of pure structural change and apply the dynamic programming
algorithm to get the break dates.
*/
qq=p+q;
zz=x~z;
{globnl,datenl,bigvec}=dating(y,zz,h,mi,qq,bigt);
xbar=pzbar(x,mi,datenl[1:mi,mi]);
zbar=pzbar(z,mi,datenl[1:mi,mi]);
teta=olsqr(y,zbar~xbar);
delta1=teta[1:q*(mi+1),1];
beta1=olsqr(y-zbar*delta1,x);
/* Calculate SRR.
*/
ssr1=(y-x*beta1-zbar*delta1)'(y-x*beta1-zbar*delta1);
if printd == 1;
print "The iterations are initialized with: ";
print "Break dates: " datenl[1:mi,mi]';
print "Delta1: " delta1;
print "Beta1: " beta1;
print "SSR1: " ssr1;
else;
endif;
@******************************************************************@
else;
@******************************************************************@
/*
if fix = 1, the initial value of beta is fixed. We initialize
SSR at some negative number to start the iterations*/
beta1=betaini;
ssr1=-5;
endif;
@*****************************************************************@
/* Starting the iterations. */
if printd == 1;
print "Output from the iterations";
else;
endif;
length=99999999.0;
i=1;
do while length > eps;
{globnl,datenl,bigvec}=dating(y-x*beta1,z,h,mi,q,bigt);
zbar=pzbar(z,mi,datenl[1:mi,mi]);
teta1=olsqr(y,x~zbar);
beta1=teta1[1:p,1];
delta1=teta1[p+1:p+q*(mi+1),1];
ssrn=(y-(x~zbar)*teta1)'(y-(x~zbar)*teta1);
/* Calculate overall SRR and check if significantly smaller.
*/
length=abs(ssrn-ssr1);
if printd == 1;
print "Iteration " i;
print "Break dates: " datenl[1:mi,mi]';
print "Delta1: " delta1;
print "Beta1: " beta1;
print "SSRN: " ssrn;
print "length" length;
else;
endif;
if i >= maxi;
print "The number of iterations has reached the upper limit";
goto itmax;
else;
i=i+1;
ssr1=ssrn;
global[mi,1]=ssrn;
datevec[1:mi,mi]=datenl[1:mi,mi];
endif;
endo;
mi=mi+1;
endo;
itmax:
retp(global,datevec,bigvec);
endp;
proc(2)=onebp(y,z,x,h,start,last);
local ssrind,i,zb,bb,rr,bdat;
/*procedure that computes an optimal one break partition in the case
of a partial structural by searching over all possible breaks. */
ssrind=9999999999;
i=h;
do while i <= last-start+1-h;
zb=pzbar(z[start:last,.],1,i);
bb=olsqr(y[start:last,1],x[start:last,.]~zb);
rr=(y[start:last,1]-(x[start:last,.]~zb)*bb)'
(y[start:last,.]-(x[start:last,.]~zb)*bb);
if rr < ssrind;
ssrind=rr;
bdat=i;
else;
endif;
i=i+1;
endo;
retp(ssrind,bdat+start-1);
endp;
**************************************austria***************************************
Q3 1992 | 101,473 | 135,098 | 309,542 | 548,755 | 0,251005 | 0,24619 | 544,605 |
Q4 1992 | 102,343 | 135,795 | 328,093 | 559,386 | 0,230521 | 0,242757 | 556,998 |
Q1 1993 | 108,671 | 90,512 | 292,205 | 514,714 | 0,221167 | 0,175849 | 510,877 |
Q2 1993 | 108,893 | 126,607 | 316,955 | 544,092 | 0,217324 | 0,232694 | 539,95 |
Q3 1993 | 108,777 | 140,481 | 320,846 | 568,962 | 0,2449 | 0,246908 | 565,604 |
Q4 1993 | 109,769 | 138,068 | 337,902 | 577,66 | 0,225027 | 0,239013 | 575,243 |
Q1 1994 | 114,858 | 102,399 | 311,433 | 547,736 | 0,221722 | 0,18695 | 543,945 |
Q2 1994 | 115,073 | 133,867 | 336,489 | 569,695 | 0,207362 | 0,23498 | 565,132 |
Q3 1994 | 114,616 | 145,262 | 340,515 | 591,787 | 0,230921 | 0,245463 | 588,643 |
Q4 1994 | 114,943 | 150,913 | 356,258 | 615,486 | 0,234425 | 0,245193 | 609,541 |
Q1 1995 | 121,035 | 103,681 | 326,099 | 573,317 | 0,220093 | 0,180844 | 560,532 |
Q2 1995 | 120,826 | 137,972 | 342,406 | 598,439 | 0,225933 | 0,230553 | 588,823 |
Q3 1995 | 120,818 | 146,286 | 346,951 | 612,718 | 0,236567 | 0,238749 | 606,237 |
Q4 1995 | 122,583 | 146,199 | 362,657 | 630,81 | 0,230767 | 0,231764 | 618,219 |
Q1 1996 | 123,597 | 108,507 | 343,218 | 593,159 | 0,213002 | 0,182931 | 586,716 |
Q2 1996 | 124,357 | 142,841 | 366,553 | 621,653 | 0,210315 | 0,229776 | 614,531 |
Q3 1996 | 124,418 | 151,396 | 361,589 | 638,093 | 0,238345 | 0,237263 | 634,71 |
Q4 1996 | 125,191 | 150,574 | 381,426 | 649,706 | 0,220237 | 0,231757 | 641,302 |
Q1 1997 | 119,705 | 114,075 | 348,547 | 601,033 | 0,220921 | 0,189798 | 590,582 |
Q2 1997 | 120,234 | 144,327 | 363,578 | 627,099 | 0,228492 | 0,23015 | 619,285 |
Q3 1997 | 120,674 | 153,416 | 371,14 | 651,413 | 0,245004 | 0,235513 | 641,721 |
Q4 1997 | 121,585 | 156,345 | 391,328 | 668,045 | 0,232218 | 0,234034 | 658,096 |
Q1 1998 | 124,12 | 120,551 | 353,956 | 630,917 | 0,242252 | 0,191073 | 620,493 |
Q2 1998 | 124,846 | 149,853 | 373,254 | 658,368 | 0,243432 | 0,227613 | 647,687 |
Q3 1998 | 124,598 | 158,422 | 381,212 | 670,142 | 0,24522 | 0,236401 | 660,664 |
Q4 1998 | 124,667 | 163,342 | 394,513 | 687,837 | 0,245199 | 0,237472 | 676,776 |
Q1 1999 | 123,6 | 114,3 | 342,7 | 645,893 | 0,278054 | 0,176964 | 611,9 |
Your Answer
6 Answers
After changing the commas to periods in Pakistan.txt, the first big problem is that your code assigns bigt to equal the number of rows of z right after z is created, which is 36. However, the code strips off 5 of the observations so that z only has 31 rows, but the code is trying to access 36 rows.
You need to either change bigt to equal the number of rows of z after it is decreased to 31, or you need to you need to index into z with bigt - 5 if bigt needs to remain the same size. I don't know which is correct for your analysis, but either way will allow this part of the code to run.
Many thanks for your response!
would you like please to indicate me which line should I change (your modification to the code).
I am grateful,
N.K
Around lines 27 to 37 or so of break_fseq, some of the first observations are trimmed of some of the variables. If you assign bigt to be equal to the number of rows of z or y1 (they are the same at this point, I think), this will solve your first problem.
However, before I did that I would make sure that I understood whether or not big was supposed to be equal to total number of observations in the original dataset or the number after trimming the observations.
Dear Prof.
Do you have an example dataset for this code that works successfully? If so, that would make it much, much easier to sort out.
Please find the original code for an example which the code run.
Many Thanks,
************************************break_fhseq****************************************
/*sequential procedure to select number of breaks*/
new;
format /ld 6,3;
cls;
vls = reshape(error(0),9,1);
vls[4] = 9999.99;
mat = xlsreadm("austria","b2:g142",1,vls);
gpdi=mat[.,6];
gnsave=mat[.,5];
x=0;
bigt=rows(mat);
y=gpdi;
z=gnsave;
T=rows(y);
y1=y[4:T-2,1];
z1=z[4:T-2,1];
deltaz=z1-z[3:T-3,1];
lz=z[3:T-3,1]-z[2:T-4,1];
llz=z[2:T-4,1]-z[1:T-5,1];
fz=z[5:T-1,1]-z1;
ffz=z[6:T,1]-z[5:T-1,1];
x=deltaz~lz~llz~fz~ffz;
y=y1;
z=z1;
bigt=rows(y);
z=z~ones(bigt,1);
q=2; @number of regressors whose coefficients are allowed to change: intercept and slope@
p=5; @number of first differenced regressors used in DOLS estimation@
m=5; @maximum number of structural changes allowed@
eps1=.15; @Value of the trimming (in percentage)@
h=int(eps1*(bigt)); @minimal length of a segment (h >= q).@
/* the following are options if p > 0.
----------------------------------- */
fixb=0; @set to 1 if use fixed initial values for beta@
betaini={0,0,0}; @if fixb=1, load the initial value of beta.@
maxi=20; @maximum number of iterations for the nonlinear
procedure to obtain global minimizers.@
printd=0; @set to 1 if want the output from the iterations
to be printed.@
eps=0.001; @criterion for the convergence.@
/*--------------------------------- */
doseq=1;
call pbreak(bigt,y,z,deltaz,q,m,h,eps1,p,fixb,x,q,eps,maxi,fixb,betaini,printd,doseq);
#include break_seq.src @set the path to where you store the file brcode.src@
end;
************************************break_seq*********************************************
proc(0)=pbreak(bigt,y,z,deltaz,q,m,h,eps1,p,fixb,x,q,eps,maxi,fixb,betaini,printd,doseq);
local datevec,bigvec,global,supfl,ndat,maic,mbic,mlwz,ssrzero,nbreak,i;
local dateseq,ftest,wftest,cv,repartda,zz,siglev,j,nbr,datese,cvm,reparv,ii,cvall,size,seqtests;
print "The options chosen are:";
print "h = " h;
print "eps1 = " eps1;
print "The maximum number of breaks is: " m;
print "********************************************************";
if doseq==1;
seqtests=zeros(m+1,1);
print "The following options are used:";
print "h = " h;
print "bigt: " bigt;
{nbr,datese}=sequa(m,2,q,h,bigt,z,y,
x,p,eps1);
print"the seqtest selects the number of breaks at" nbr;
endif;
endp;
@*******************************************************************@
proc(1)=pftest(y,z,g,q,bigt,datevec,x,p);
local ftest,zbar,gammahat,deltahat,ssr0,ssrk,wbar,w,e,eb,eb1,eb2,ef,ef1,ef2,u,ub,uf,uhat,uhatnull,fsample,b;
local te,ae,ae1,ae2,se1,se2,ee,se,ad,a2,eband,jb,jband,kern,sig,lam,j,delta,omega,a,uu,ai,y1,z1,w1;
if g==1;
fsample=zeros(bigt,1);
for b(eps1*bigt,(1-eps1)*bigt,1);
zbar=pzbar(z,1,b);
wbar=zbar~x;
w=z~x;
gammahat=olsqr(y,wbar);
deltahat=olsqr(y,w);
ssr0=(y-w*deltahat)'*(y-w*deltahat);
ssrk=(y-wbar*gammahat)'*(y-wbar*gammahat);
uhat=y-wbar*gammahat;
uhatnull=y-w*deltahat;
u=uhat;
e = u;
te = bigt;
eb = e[1:te-1,.];
ef = e[2:te,.];
ae=eb'*ef/(eb'*eb);
ee=zeros(te-1,1);
ee[.,1]=ef[.,1]-eb[.,1]*(ae);
se=sqrt(inv(te-1)*ee[.,1]'*ee[.,1]);
ad=((se)^2)/(1-ae)^4;
a2=(4*(ae*se/(1-ae)^4)^2/ad);
eband = 1.3221*((a2*te)^.2);
jb = seqa(1,1,te-1)/eband;
jband = jb*1.2*pi;
kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
sig = uhatnull'*uhatnull;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (uhatnull[1:te-j,.]'*uhatnull[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam' ;
omega[1,1]=omega[1,1]/te;
fsample[b,1]=(ssr0-ssrk)/(g*omega[1,1]);
endfor;
fsample=fsample[eps1*bigt:(1-eps1)*bigt,1];
ftest=maxc(fsample);
else;
{zbar}=pzbar(z,g,datevec[1:g,g]);
wbar=zbar~x;
w=z~x;
gammahat=olsqr(y,wbar);
deltahat=olsqr(y,w);
ssrk=(y-wbar*gammahat)'*(y-wbar*gammahat);
ssr0=(y-w*deltahat)'*(y-w*deltahat);
uhat=y-wbar*gammahat;
uhatnull=y-w*deltahat;
u=uhat;
e = u;
te = bigt;
eb = e[1:te-1,.];
ef = e[2:te,.];
ae=eb'*ef/(eb'*eb);
ee=zeros(te-1,1);
ee[.,1]=ef[.,1]-eb[.,1]*(ae);
se=sqrt(inv(te-1)*ee[.,1]'*ee[.,1]);
ad=((se)^2)/(1-ae)^4;
a2=(4*(ae*se/(1-ae)^4)^2/ad);
eband = 1.3221*((a2*te)^.2);
jb = seqa(1,1,te-1)/eband;
jband = jb*1.2*pi;
kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
sig = uhatnull'*uhatnull;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (uhatnull[1:te-j,.]'*uhatnull[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam' ;
omega[1,1]=omega[1,1]/te;
ftest=(ssr0-ssrk)/(g*omega[1,1]);
endif;
retp(ftest);
endp;
proc(3)=dating(y,z,h,m,q,bigt);
@This is the main procedure which calculates the break points that
globally minimizes the SSR. It returns optimal dates and associated SSR
for all numbers of breaks less than or equal to m.@
local datevec,optdat,optssr,dvec,i,ssrmin,datx,j1,ib,jlast;
local global,vecssr,jb,xx,bigvec;
datevec=zeros(m,m);
optdat=zeros(bigt,m);
optssr=zeros(bigt,m);
dvec=zeros(bigt,1);
global=zeros(m,1);
bigvec=zeros(bigt*(bigt+1)/2,1);
@
segment to construct the vector of SSR of dimension T*(T+1)/2
@
i=1;
do while i <= bigt-h+1;
{vecssr}=ssr(i,y,z,h,bigt);
bigvec[(i-1)*bigt+i-(i-1)*i/2:i*bigt-(i-1)*i/2,1]=vecssr[i:bigt,1];
i=i+1;
endo;
@
Section that applies the dynamic programming algorithm to look for the
optimal breaks
The first case is when m = 1. Here the dynamic programming algorithm is not
needed and one call to the parti(.) procedure is enough.
@
if m == 1;
{ssrmin,datx}=parti(1,h,bigt-h,bigt,bigvec,bigt);
datevec[1,1]=datx;
global[1,1]=ssrmin;
else;
@ when m > 1, a dynamic programming algorithm is used.
The first step is to obtain the optimal one break partitions for all
possible ending dates from 2h to T-mh+1.
The optimal dates are stored in a vector optdat.
The associated ssr are stored in a vector optssr.
@
j1=2*h; @First loop. Looking for the@
do while j1 <= bigt; @optimal one break partitions@
{ssrmin,datx}=parti(1,h,j1-h,j1,bigvec,bigt);@for break dates between h and@
optssr[j1,1]=ssrmin; @T-h. j1 is the last date of the@
optdat[j1,1]=datx; @segments.@
j1=j1+1;
endo;
global[1,1]=optssr[bigt,1];
datevec[1,1]=optdat[bigt,1];
@
When this is done the algorithm looks for optimal 2,3,... breaks
partitions. The index used is ib.
@
ib=2;
do while ib <= m;
if ib == m; @if we have reached the highest number of breaks
considered, only one segment is considered, that
which ends at the last date of the sample.@
jlast=bigt;
jb=ib*h; /* date of the break */
do while jb <=jlast-h;
dvec[jb,1] = optssr[jb,ib-1]+bigvec[(jb+1)*bigt-jb*(jb+1)/2,1];
jb=jb+1;
endo;
optssr[jlast,ib]=minc(dvec[ib*h:jlast-h,1]);
optdat[jlast,ib]=(ib*h-1)+minindc(dvec[ib*h:jlast-h,1]);
else;
@if we have not reached the highest number of breaks
considered, we need to loop over the last date of
the segment, between (ib+1)*h and T.@
jlast=(ib+1)*h;
do while jlast <= bigt;
jb=ib*h; /* date of the break */
do while jb <=jlast-h;
dvec[jb,1] = optssr[jb,ib-1]+bigvec[jb*bigt-jb*(jb-1)/2+jlast-jb,1];
jb=jb+1;
endo;
optssr[jlast,ib]=minc(dvec[ib*h:jlast-h,1]);
optdat[jlast,ib]=(ib*h-1)+minindc(dvec[ib*h:jlast-h,1]);
jlast=jlast+1;
endo;
endif;
@At each time we loop the results with ib breaks are retrieved
and printed@
datevec[ib,ib]=optdat[bigt,ib];
i=1;
do while i <= ib-1;
xx=ib-i;
datevec[xx,ib]=optdat[datevec[xx+1,ib],xx];
i=i+1;
endo;
global[ib,1]=optssr[bigt,ib];
ib=ib+1;
endo;
endif; /*closing the if for the case m >1*/
retp(global,datevec,bigvec);
endp;
@********************************************************************@
proc(1)=ssr(start,y,z,h,last);
local vecssr,delta1,delta2,inv1,inv2,invz,res,v,f,r;
@
This procedure computes recursive residuals from a data set that
starts at date "start" and ends at date "last". It returns a
vector of sum of squared residuals (SSR) of length last-start+1 (stored for
convenience in a vector of length T).
start: starting entry of the sample used.
last: ending date of the last segment considered
y: dependent variable
z: matrix of regressors of dimension q
h: minimal length of a segment
@
/* initialize the vectors */
vecssr=zeros(last,1);
/* initialize the recursion with the first h data points */
inv1=inv(z[start:start+h-1,.]'z[start:start+h-1,.]);
delta1=inv1*(z[start:start+h-1,.]'y[start:start+h-1,1]);
res=y[start:start+h-1,1]-z[start:start+h-1,.]*delta1;
vecssr[start+h-1,1]=res'res;
/* loop to construct the recursive residuals and update the SSR */
r=start+h;
do while r <= last;
v=y[r,1]-z[r,.]*delta1;
invz=inv1*z[r,.]';
f=1+z[r,.]*invz;
delta2=delta1+(invz*v)/f;
inv2=inv1-(invz*invz')/f;
inv1=inv2;
delta1=delta2;
vecssr[r,1]=vecssr[r-1,1]+v*v/f;
r=r+1;
endo;
retp(vecssr);
endp;
@********************************************************************@
proc(2)=parti(start,b1,b2,last,bigvec,bigt);
local k,dvec,j,ssrmin,dx,ini,jj;
@
procedure to obtain an optimal one break partitions for a segment that
starts at date start and ends at date last. It return the optimal break
and the associated SSR.
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
@
dvec=zeros(bigt,1);
ini=(start-1)*bigt-(start-2)*(start-1)/2+1;
j=b1;
do while j<=b2;
jj=j-start;
k=j*bigt-(j-1)*j/2+last-j;
dvec[j,1]=bigvec[ini+jj,1]+bigvec[k,1];
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@********************************************************************@
proc(1)=ssrnul(y,zz);
local delta,ssrn;
@Computes the SSR under the null of no break@
delta=olsqr(y,zz);
ssrn=(y-zz*delta)'(y-zz*delta);
retp(ssrn);
endp;
@********************************************************************@
proc(1)=pzbar(zz,m,bb);
local nt,q1,zb,i;
@procedure to construct the diagonal partition of z with m break at
dates b.@
nt=rows(zz);
q1=cols(zz);
zb=zeros(nt,(m+1)*q1);
zb[1:bb[1,1],1:q1]=zz[1:bb[1,1],.];
i=2;
do while i <= m;
zb[bb[i-1,1]+1:bb[i,1],(i-1)*q1+1:i*q1]=zz[bb[i-1,1]+1:bb[i,1],.];
i=i+1;
endo;
zb[bb[m,1]+1:nt,m*q1+1:(m+1)*q1]=zz[bb[m,1]+1:nt,.];
retp(zb);
endp;
@********************************************************************@
proc plambda(b,m,bigt);
local lambda,k;
@procedure that construct a diagonal matrix of dimension m+1 with ith
entry (T_i - T_i-1)/T.@
lambda=zeros(m+1,m+1);
lambda[1,1]=b[1,1]/bigt;
k=2;
do while k <= m;
lambda[k,k]=(b[k,1]-b[k-1,1])/bigt;
k=k+1;
endo;
lambda[m+1,m+1]=(bigt-b[m,1])/bigt;
retp(lambda);
endp;
@********************************************************************@
proc psigmq(res,b,q,i,nt);
local sigmat,kk,sigmq;
@procedure that computes a diagonal matrix of dimension i+1 with ith
entry the estimate of the variance of the residuals for segment i.
@
sigmat=zeros(i+1,i+1);
sigmat[1,1]=res[1:b[1,1],1]'res[1:b[1,1],1]/b[1,1];
kk = 2;
do while kk <= i;
sigmat[kk,kk]=res[b[kk-1,1]+1:b[kk,1],1]'
res[b[kk-1,1]+1:b[kk,1],1]/(b[kk,1]-b[kk-1,1]);
kk=kk+1;
endo;
sigmat[i+1,i+1]=res[b[i,1]+1:nt,1]'res[b[i,1]+1:nt,1]/(nt-b[i,1]);
retp(sigmat);
endp;
@********************************************************************@
proc getcv2(signif,eps1);
local cv;
if eps1 == .15;
if signif == 2;
cv={12.11 13.79 15.25 16.38 17.03 17.70 17.98 18.25 18.33 18.53};
endif;
endif;
retp(cv);
endp;
@********************************************************************@
proc(2)=partione(b1,b2,last,vssr,vssrev);
local dvec,j,ssrmin,dx;
@
procedure to obtain an optimal one break partitions for a segment that
starts at date start and ends at date last. It return the optimal break
and the associated SSR.
Procedure used with the sequential method when the T*(T+1)/2 vector
of SSR is not computed.
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
@
dvec=zeros(last,1);
j=b1;
do while j<=b2;
dvec[j,1]=vssr[j,1]+vssrev[last-j,1];
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@********************************************************************@
proc(2)=sequa(m,signif,q,h,bigt,z,y,x,p,eps1);
local cv,ssrmin,datx,dv,nbreak,ds,ftest,nseg,is,maxf,newseg;
local ftestv,dv2,length,vssr,vssrev,dum,dbreak;
@This procedure applies the sequential procedure to obtain the number of
breaks and the break dates.
The current version only allow pure structural changes. This will
be generalized.
Here m is the maximum number of breaks allowed.@
dv=zeros(m+2,1);
dv2=zeros(m+2,1);
ftestv=zeros(m+1,1);
@input of the critical values of the supF(l+1|l) test@
cv=getcv2(signif,eps1);
dv[1,1]=0; @the first date of a segment is set to 0@
@Obtain the first break and test if significant.@
if p == 0;
{vssrev}=ssr(1,rev(y),rev(z),h,bigt);
{vssr}=ssr(1,y,z,h,bigt);
{ssrmin,datx}=partione(h,bigt-h,bigt,vssr,vssrev);
else;
{ssrmin,datx}=onebp(y,z,x,h,1,bigt);
endif;
dv[2,1]=datx;
ftest=pftest(y,z,1,q,bigt,dv[2,1],x,p);
if ftest < 12.11;
nbreak=0;
dv[2,1]=0;
retp(nbreak,0);
else;
print "The first break found is at: " datx;
nbreak=1;
nseg=2;
dv[nseg+1,1]=bigt; @set the last date of the segment to T@
endif;
do while nseg <= m;
ds=zeros(nseg+1,1);
ftestv=zeros(nseg+1,1);
is=1;
do while is <= nseg;
length=dv[is+1,1]-dv[is,1];
if length >= 2*h;
if p == 0;
{vssr}=ssr(1,y[dv[is,1]+1:dv[is+1,1],1],
z[dv[is,1]+1:dv[is+1,1],.],h,length);
{vssrev}=ssr(1,rev(y[dv[is,1]+1:dv[is+1,1],1]),
rev(z[dv[is,1]+1:dv[is+1,1],.]),h,length);
{dum,ds[is,1]}=partione(h,length-h,length,vssr,vssrev);
ftestv[is,1]=pftest(y[dv[is,1]+1:dv[is+1,1],1],z[dv[is,1]+1:dv[is+1,1],.],
1,q,length,ds[is,1],0,p);
else;
{dum,ds[is,1]}=onebp(y,z,x,h,dv[is,1]+1,dv[is+1,1]);
ds[is,1]=ds[is,1]-dv[is,1];
ftestv[is,1]=pftest(y[dv[is,1]+1:dv[is+1,1],1],z[dv[is,1]+1:dv[is+1,1],.],
1,q,length,ds[is,1],x[dv[is,1]+1:dv[is+1,1],.],p);
endif;
else;
ftestv[is,1]=0.0;
endif;
is=is+1;
endo;
maxf=maxc(ftestv[1:nseg,1]);
if maxf < cv[q-1,nseg];
retp(nbreak,dv[2:nbreak+1,1]);
else;
newseg=maxindc(ftestv[1:nseg,1]);
print "The next break found is at: " ds[newseg,1]+dv[newseg,1];
dv[nseg+2,1]=ds[newseg,1]+dv[newseg,1];
nbreak=nbreak+1;
dv2=sortc(dv[2:nseg+2,1],1);
dv[1,1]=0;
dv[2:nseg+2,1]=dv2;
endif;
nseg=nseg+1;
endo;
print "The sequential procedure has reached the upper limit";
retp(nbreak,dv[2:nbreak+1,1]);
endp;
@********************************************************************@
proc(3)=nldat(y,z,x,h,m,p,q,bigt,fixb,eps,maxi,betaini,printd);
local qq,zz,xbar,zbar,teta,delta1,beta1,ssr1,i,mi,length,globnl;
local datenl,teta1,ssrn,global,datevec,bigvec;
global=zeros(m,1);
globnl=zeros(m,1);
datevec=zeros(m,m);
datenl=zeros(m,m);
mi=1;
do while mi <= m;
if printd == 1;
print "Model with " mi "breaks";
else;
endif;
if fixb == 0;
@******************************************************************@
/*Segment to initialize the beta vector. We apply treat the model
as one of pure structural change and apply the dynamic programming
algorithm to get the break dates.
*/
qq=p+q;
zz=x~z;
{globnl,datenl,bigvec}=dating(y,zz,h,mi,qq,bigt);
xbar=pzbar(x,mi,datenl[1:mi,mi]);
zbar=pzbar(z,mi,datenl[1:mi,mi]);
teta=olsqr(y,zbar~xbar);
delta1=teta[1:q*(mi+1),1];
beta1=olsqr(y-zbar*delta1,x);
/* Calculate SRR.
*/
ssr1=(y-x*beta1-zbar*delta1)'(y-x*beta1-zbar*delta1);
if printd == 1;
print "The iterations are initialized with: ";
print "Break dates: " datenl[1:mi,mi]';
print "Delta1: " delta1;
print "Beta1: " beta1;
print "SSR1: " ssr1;
else;
endif;
@******************************************************************@
else;
@******************************************************************@
/*
if fix = 1, the initial value of beta is fixed. We initialize
SSR at some negative number to start the iterations*/
beta1=betaini;
ssr1=-5;
endif;
@*****************************************************************@
/* Starting the iterations. */
if printd == 1;
print "Output from the iterations";
else;
endif;
length=99999999.0;
i=1;
do while length > eps;
{globnl,datenl,bigvec}=dating(y-x*beta1,z,h,mi,q,bigt);
zbar=pzbar(z,mi,datenl[1:mi,mi]);
teta1=olsqr(y,x~zbar);
beta1=teta1[1:p,1];
delta1=teta1[p+1:p+q*(mi+1),1];
ssrn=(y-(x~zbar)*teta1)'(y-(x~zbar)*teta1);
/* Calculate overall SRR and check if significantly smaller.
*/
length=abs(ssrn-ssr1);
if printd == 1;
print "Iteration " i;
print "Break dates: " datenl[1:mi,mi]';
print "Delta1: " delta1;
print "Beta1: " beta1;
print "SSRN: " ssrn;
print "length" length;
else;
endif;
if i >= maxi;
print "The number of iterations has reached the upper limit";
goto itmax;
else;
i=i+1;
ssr1=ssrn;
global[mi,1]=ssrn;
datevec[1:mi,mi]=datenl[1:mi,mi];
endif;
endo;
mi=mi+1;
endo;
itmax:
retp(global,datevec,bigvec);
endp;
proc(2)=onebp(y,z,x,h,start,last);
local ssrind,i,zb,bb,rr,bdat;
/*procedure that computes an optimal one break partition in the case
of a partial structural by searching over all possible breaks. */
ssrind=9999999999;
i=h;
do while i <= last-start+1-h;
zb=pzbar(z[start:last,.],1,i);
bb=olsqr(y[start:last,1],x[start:last,.]~zb);
rr=(y[start:last,1]-(x[start:last,.]~zb)*bb)'
(y[start:last,.]-(x[start:last,.]~zb)*bb);
if rr < ssrind;
ssrind=rr;
bdat=i;
else;
endif;
i=i+1;
endo;
retp(ssrind,bdat+start-1);
endp;
**************************************austria***************************************
Q3 1992 | 101,473 | 135,098 | 309,542 | 548,755 | 0,251005 | 0,24619 | 544,605 |
Q4 1992 | 102,343 | 135,795 | 328,093 | 559,386 | 0,230521 | 0,242757 | 556,998 |
Q1 1993 | 108,671 | 90,512 | 292,205 | 514,714 | 0,221167 | 0,175849 | 510,877 |
Q2 1993 | 108,893 | 126,607 | 316,955 | 544,092 | 0,217324 | 0,232694 | 539,95 |
Q3 1993 | 108,777 | 140,481 | 320,846 | 568,962 | 0,2449 | 0,246908 | 565,604 |
Q4 1993 | 109,769 | 138,068 | 337,902 | 577,66 | 0,225027 | 0,239013 | 575,243 |
Q1 1994 | 114,858 | 102,399 | 311,433 | 547,736 | 0,221722 | 0,18695 | 543,945 |
Q2 1994 | 115,073 | 133,867 | 336,489 | 569,695 | 0,207362 | 0,23498 | 565,132 |
Q3 1994 | 114,616 | 145,262 | 340,515 | 591,787 | 0,230921 | 0,245463 | 588,643 |
Q4 1994 | 114,943 | 150,913 | 356,258 | 615,486 | 0,234425 | 0,245193 | 609,541 |
Q1 1995 | 121,035 | 103,681 | 326,099 | 573,317 | 0,220093 | 0,180844 | 560,532 |
Q2 1995 | 120,826 | 137,972 | 342,406 | 598,439 | 0,225933 | 0,230553 | 588,823 |
Q3 1995 | 120,818 | 146,286 | 346,951 | 612,718 | 0,236567 | 0,238749 | 606,237 |
Q4 1995 | 122,583 | 146,199 | 362,657 | 630,81 | 0,230767 | 0,231764 | 618,219 |
Q1 1996 | 123,597 | 108,507 | 343,218 | 593,159 | 0,213002 | 0,182931 | 586,716 |
Q2 1996 | 124,357 | 142,841 | 366,553 | 621,653 | 0,210315 | 0,229776 | 614,531 |
Q3 1996 | 124,418 | 151,396 | 361,589 | 638,093 | 0,238345 | 0,237263 | 634,71 |
Q4 1996 | 125,191 | 150,574 | 381,426 | 649,706 | 0,220237 | 0,231757 | 641,302 |
Q1 1997 | 119,705 | 114,075 | 348,547 | 601,033 | 0,220921 | 0,189798 | 590,582 |
Q2 1997 | 120,234 | 144,327 | 363,578 | 627,099 | 0,228492 | 0,23015 | 619,285 |
Q3 1997 | 120,674 | 153,416 | 371,14 | 651,413 | 0,245004 | 0,235513 | 641,721 |
Q4 1997 | 121,585 | 156,345 | 391,328 | 668,045 | 0,232218 | 0,234034 | 658,096 |
Q1 1998 | 124,12 | 120,551 | 353,956 | 630,917 | 0,242252 | 0,191073 | 620,493 |
Q2 1998 | 124,846 | 149,853 | 373,254 | 658,368 | 0,243432 | 0,227613 | 647,687 |
Q3 1998 | 124,598 | 158,422 | 381,212 | 670,142 | 0,24522 | 0,236401 | 660,664 |
Q4 1998 | 124,667 | 163,342 | 394,513 | 687,837 | 0,245199 | 0,237472 | 676,776 |
Q1 1999 | 123,6 | 114,3 | 342,7 | 645,893 | 0,278054 | 0,176964 | 611,9 |