Index out of range (572)

Dear professor,

I am running a Gauss program called break_fhseq and I get usually the sampe message "Index out of range (572)".
Please find attached the files (data file: Pakistan_data and Gausse files : break_fhseq + break_seq).
I am so grateful and many thaks,
N.K
*********************************break_fhseq*****************************************

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

aptech

1,773


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.

aptech

1,773


0



Dear Prof.

Exactly, the number of rows is equal for z and y1. "bigt" have to be equal to the number of rows in the original dataset (it means 36).
Many thanks,
N.K


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.

aptech

1,773


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

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.

Exactly, the number of rows is equal for z and y1. "bigt" have to be equal to the number of rows in the original dataset (it means 36).
Many thanks,
N.K
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

You must login to post answers.

Have a Specific Question?

Get a real answer from a real person

Need Support?

Get help from our friendly experts.