BRCODE.SRC

@October 11, 2004 version@
@Copyrigth, Pierre Perron (1998, 1999, 2004).@
proc(0)=pbreak(bigt,y,z,q,m,h,eps1,robust,prewhit,hetomega,
hetq,doglobal,dotest,dospflp1,doorder,dosequa,dorepart,estimbic,estimlwz,
estimseq,estimrep,estimfix,x,eps,maxi,fixb,betaini,printd,hetdat,hetvar,fixn);
local datevec,bigvec,global,i,supfl,ndat,maic,mbic,mlwz,ssrzero,nbreak;
local dateseq,ftest,wftest,cv,repartda,zz,siglev,j,nbr,datese,cvm,reparv,ii;


print "The options chosen are:";
print "h = " h;
print "eps1 = " eps1;
print "hetdat = " hetdat;
print "hetvar = " hetvar;
print "hetomega = " hetomega;
print "hetq = " hetq;
print "robust = " robust   "(prewhit = " prewhit ")";
print "The maximum number of breaks is: " m;
print "********************************************************";

siglev={10, 5, 2.5, 1};

if doglobal == 1;

/* The following calls the procedure to obtain the break dates and the
associated sum of squared residuals for all numbers of breaks between
1 and m.*/

Print "Output from the global optimization";
Print "********************************************************";

if p == 0;
{global,datevec,bigvec}=dating(y,z,h,m,q,bigt);
else;
print "This is a partial structural change model and the following";
print "specifictaions were used:";
print "The number of regressors with fixed coefficients, p: " p;
print "Whether (1) or not (0) the initial values of beta are fixed: "
fixb;
print "If so, at the following values: " betaini;
print "The convergence criterion is: " eps;
print "Whether (1) or not (0) the iterations are printed: " printd;
print "--------------------------";
{global,datevec,bigvec}=nldat(y,z,x,h,m,p,q,bigt,fixb,
eps,maxi,betaini,printd);
endif;

i=1;
do while i <= m;
print "The model with" i "breaks has SSR : " global[i,1];
print "The dates of the breaks are:  " datevec[1:i,i];
i=i+1;
endo;

endif;

if dotest == 1;

/* the following calls the procedure to estimate and print various
tests for the number of breaks. It also returns the UDmax and WDmax
tests.*/

print "********************************************************";
Print "Output from the testing procedures";
print "********************************************************";
print "a) supF tests against a fixed number of breaks";
print "--------------------------------------------------------------";

ftest=zeros(m,1);
wftest=zeros(m,1);

i=1;
do while i <= m;
ftest[i,1]=pftest(y,z,i,q,bigt,datevec,prewhit,robust,x,p,hetdat,hetvar);
print "The supF test for 0 versus" i "breaks (scaled by q) is:" ftest[i,1];
i=i+1;
endo;

print "-------------------------";
j=1;
do while j <= 4;
cv=getcv1(j,eps1);

print "The critical values at the " siglev[j,1] "% level are (for k=1 to " m "):";
print  cv[q,1:m];
j=j+1;
endo;



print "--------------------------------------------------------------";
print "b) Dmax tests against an unknown number of breaks";
print "--------------------------------------------------------------";
print "The UDmax test is: " maxc(ftest);
j=1;
do while j<= 4;
cvm=getdmax(j,eps1);
print "(the critical value at the " siglev[j,1] "% level is: " cvm[q,1] ")";
j=j+1;
endo;
print "********************************************************";


j=1;
do while j <= 4;
cv=getcv1(j,eps1);

i=1;
do while i <= m;
wftest[i,1]=cv[q,1]*ftest[i,1]/cv[q,i];
i=i+1;
endo;

print "---------------------";
print "The WDmax test at the " siglev[j,1] "% level is: " maxc(wftest);
cvm=getdmax(j,eps1);
print "(The critical value is: " cvm[q,2] ")";
j=j+1;
endo;

endif;

if dospflp1 == 1;

/* The  following calls the procedure that estimates the supF(l+1|l)
tests where the first l breaks are taken from the global minimization*/

print "********************************************************";
print "supF(l+1|l) tests using global otimizers under the null";
print "--------------------------------------------------------------";
i=1;
do while i <= m-1;

{supfl,ndat}=spflp1(bigvec,datevec[1:i,i],i+1,y,z,h,prewhit,robust,
x,p,hetdat,hetvar);

print "The supF(" i+1 "|" i ") test is : " supfl;
print "It corresponds to a new break at: " ndat;
i=i+1;
endo;
print "********************************************************";


j=1;
do while j <=4;
cv=getcv2(j,eps1);
print "The critical values of supF(i+1|i) at the " siglev[j,1] "% level are (for i=1 to " m ") are: ";
print cv[q,1:m];
j=j+1;
endo;

endif;

if doorder == 1;

/* The following calls the procedure to estimate the number of
breaks using the information criteria BIC and LWZ (Liu, Wu and
Zidek). It return the two orders selected. */

print "********************************************************";
Print "Output from the application of Information criteria";
print "--------------------------------------------------------------";

if p == 0;
zz=z;
else;
zz=z~x;
endif;
ssrzero=ssrnul(y,zz);
{mbic,mlwz}=order(ssrzero,global,bigt,m,q);

endif;

if dosequa == 1;

/*The following section calls the sequential procedure that estimate
each break one at a time. It stops when the supF(l+1|l) test is not
significant. It returns the number of breaks found and the break dates.
Note that it can be used independently of the other procedures, i.e.
global minimizers need not be obtained since it used a method to
compute the breaks in O(T) operations.*/

nbreak=zeros(4,1);
dateseq=zeros(4,m);

j =1;
do while j <= 4;
print "********************************************************";
Print "Output from the sequential procedure at significance level " siglev[j,1] "%";
print "--------------------------------------------------------------";

{nbr,datese}=sequa(m,j,q,h,bigt,robust,prewhit,z,y,
x,p,hetdat,hetvar,eps1);

print "----------------------------------------------------";
print "The sequential procedure estimated the number of breaks at:" nbr;
if nbreak >= 1;
print "The break dates are: " datese;
else;
endif;
nbreak[j,1]=nbr;
if nbr /= 0;
dateseq[j,1:nbreak[j,1]]=datese';
endif;
j=j+1;
endo;

endif;

if dorepart == 1;

@The following procedure construct the so-called repartition estimates
of the breaks obtained by the sequential method (see Bai (1995), Estimating
Breaks one at a time, Econometric Theory, 13, 315-352. It allows
estimates that have the same asymptotic distribution as those obtained
by global minimization. Otherwise, the output from the procedure "estim"
below do not delivers asymptotically correct confidence intervals for the
break dates.@

reparv=zeros(4,m);

j=1;
do while j <= 4;

print "********************************************************";
print "Output from the repartition procedure for the " siglev[j,1] "% significance level";

if nbreak[j,1] == 0;
print "********************************************************";
print "The sequential procedure found no break and ";
print "the repartition procedure is skipped.";
print "********************************************************";
else;
{repartda}=preparti(y,z,nbreak[j,1],dateseq[j,1:nbreak[j,1]]',h,x,p);

print "----------------------------------------";
print "The updated break dates are :" repartda;

reparv[j,1:nbreak[j,1]]=repartda';

endif;
j=j+1;
endo;
endif;

/* the following calls a procedure which estimates the model using
as the number of breaks the value specified in the first argument.
It also calls a procedure to calculate standard errors
in a way that depends on the specification for robust
and returns asymptotic confidence intervals for the break dates.*/

if estimbic == 1;

print "********************************************************";
Print "Output from the estimation of the model selected by BIC";
print "--------------------------------------------------------------";

call estim(mbic,q,z,y,datevec[.,mbic],robust,prewhit,
hetomega,hetq,x,p);

endif;


if estimlwz == 1;

print "********************************************************";
Print "Output from the estimation of the model selected by LWZ";
print "--------------------------------------------------------------";

call estim(mlwz,q,z,y,datevec[.,mlwz],robust,prewhit,
hetomega,hetq,x,p);

endif;


if estimseq == 1;
print "********************************************************";
ii=0;
j=1;
do while j <= 4;

if ii == 0;
if nbreak[j,1] /= 0;
Print "Output from the estimation of the model selected by the
 sequential method at significance level " siglev[j,1] "%";
print "--------------------------------------------------------------";

call estim(nbreak[j,1],q,z,y,dateseq[j,1:nbreak[j,1]]',robust,prewhit,
hetomega,hetq,x,p);
endif;
endif;

j=j+1;
if j <= 4;
if nbreak[j,1] == nbreak[j-1,1];
if nbreak[j,1] == 0;
print "For the " siglev[j,1] "% level, the model is the same as for the " siglev[j-1,1] "% level.";
print "The estimation is not repeated.";
print "----------------------------------------------------------------";
ii=1;
else;
if dateseq[j,1:nbreak[j,1]] == dateseq[j-1,1:nbreak[j-1,1]];
print "For the " siglev[j,1] "% level, the model is the same as for the " siglev[j-1,1] "% level.";
print "The estimation is not repeated.";
print "----------------------------------------------------------------";
ii=1;
endif;
endif;
else;
ii=0;

endif;
endif;
endo;

endif;

if estimrep == 1;
ii=0;
print "********************************************************";
j=1;
do while j <= 4;

if ii == 0;
if nbreak[j,1] == 0;
print "********************************";
print "The sequential procedure at the significance level " siglev[j,1] "%  found no break and ";
print "the repartition procedure was skipped.";
print "********************************************************";
else;

print "Output from the estimation of the model selected by the";
print "repartition method from the sequential procedure at the significance level " siglev[j,1] "%";
print "--------------------------------------------------------------";

call estim(nbreak[j,1],q,z,y,reparv[j,1:nbreak[j,1]]',robust,prewhit,
hetomega,hetq,x,p);

endif;
endif;

j=j+1;
if j <= 4;
if nbreak[j,1] == nbreak[j-1,1];
if nbreak[j,1] == 0;
print "For the " siglev[j,1] "% level, the model is the same as for the " siglev[j-1,1] "% level.";
print "The estimation is not repeated.";
print "----------------------------------------------------------------";
ii=1;
else;
if dateseq[j,1:nbreak[j,1]] == dateseq[j-1,1:nbreak[j-1,1]];
print "For the " siglev[j,1] "% level, the model is the same as for the " siglev[j-1,1] "% level.";
print "The estimation is not repeated.";
print "----------------------------------------------------------------";
ii=1;
endif;
endif;
else;
ii=0;

endif;
endif;
endo;

endif;


if estimfix == 1;

print "********************************************************";

Print "Output from the estimation of the model with" fixn "breaks";
print "--------------------------------------------------------------";

call estim(fixn,q,z,y,datevec[.,fixn],robust,prewhit,
hetomega,hetq,x,p);

endif;


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(0)=estim(m,q,z,y,b,robust,prewhit,hetomega,hetq,x,p);
local dum,nt,vdel,zbar,i,d,bound,reg;

@ This procedure estimate by OLS the model given the obtained break
dates. It also computes and report confidence intervals for the break
dates. The method used depends on the specification for robust
@

if m == 0;

print "There are no breaks in this model and estimation is skipped";
else;

nt=rows(z);
d=(m+1)*q+p;                       @total number of parameters@
vdel=zeros(d,d);

@Construct the Z_bar matrix. The diagonal partition of Z at the
estimated break dates.@

{zbar}=pzbar(z,m,b);

@Estimation and printing@

_olsres=1;
__con=0;

if p == 0;
reg=zbar;
else;
reg=x~zbar;
endif;

call ols(0,y,reg);

vdel=pvdel(y,z,m,q,bigt,b,prewhit,robust,x,p,1,hetdat,hetvar);

print "--------------------------------------------------------------";
print "Corrected standard errors for the coefficients";
print "--------------------------------------------------------------";
i=1;
do while i <= d;
print "The corrected standard error for coefficient" i "is:" sqrt(vdel[i,i]);
i=i+1;
endo;
if robust == 0 and hetdat == 1 and hetvar == 0;
print "In the case robust == 0, hetdat == 1 and hetvar == 0, the 'corrected' are the same";
print "as that of the printout except for a different small sample correction.";
endif;


@confidence interval for the break dates@


bound=interval(y,z,zbar,b,q,m,robust,prewhit,hetomega,hetq,x,p);
print "--------------------------------------------------------------";
print "Confidence intervals for the break dates";
print "--------------------------------------------------------------";
i=1;
do while i <= m;
print "The 95% C.I. for the" i "th break is: " bound[i,1] bound[i,2];
print "The 90% C.I. for the" i "th break is: " bound[i,3] bound[i,4];
i=i+1;
endo;
print "********************************************************";


endif;
endp;

@********************************************************************@

proc 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 correct(reg,res,prewhit);
local jh,imat,d,nt,i,b,bmat,vstar,hac,vmat;
@main procedures which activates the computation of robust standard
errors.@

d=cols(reg);
nt=rows(reg);
b=zeros(d,1);
bmat=zeros(d,d);
vstar=zeros(nt-1,d);
vmat=zeros(nt,d);

@First construct the matriz z_t*u_t.@

i=1;
do while i <= d;
vmat[.,i]=reg[.,i].*res;
i=i+1;
endo;

@Procedure that applies prewhitenning to the matrix vmat by filtering
with a VAR(1). If prewhit=0, it is skipped.@


if prewhit == 1;

/* estimating the VAR(1) */
i=1;
do while  i <= d;
b=olsqr(vmat[2:nt,i],vmat[1:nt-1,.]);
bmat[.,i]=b;
vstar[.,i]=vmat[2:nt,i]-vmat[1:nt-1,.]*b;
i=i+1;
endo;

/* call the kernel on the residuals */

jh=jhatpr(vstar);

/* recolor */

hac=inv(eye(d)-bmat)*jh*(inv(eye(d)-bmat))';

else;
hac=jhatpr(vmat);
endif;
retp(hac);
endp;

@********************************************************************@

proc jhatpr(vmat);
local d,nt,jhat,j,st;
@Procedure to compute the long-run covariance matrix of vmat.@

nt=rows(vmat);
d=cols(vmat);
jhat=zeros(d,d);

/* calling the automatic bandwidth selection*/
{st}=bandw(vmat);

@lag 0 covariance@
jhat=vmat'vmat;

@forward sum@
j=1;
do while j <= nt-1;
jhat=jhat+kern(j/st)*vmat[j+1:nt,.]'vmat[1:nt-j,.];
j=j+1;
endo;

@bacward sum@
j=1;
do while j <= nt-1;
jhat=jhat+kern(j/st)*vmat[1:nt-j,.]'vmat[j+1:nt,.];
j=j+1;
endo;

@small sample correction@
jhat=jhat/(nt-d);
retp(jhat);
endp;

@********************************************************************@

proc bandw(vhat);
local nt,i,b,sig,a2,a2n,a2d,st,d;
@procedure that compute the automatic bandwidth based on AR(1)
approximation for each vector of the matrix vhat. Each are given
equal weigths of 1.@

nt=rows(vhat);
d=cols(vhat);
a2n=0;
a2d=0;
i=1;
do while i <= d;
b=olsqr(vhat[2:nt,i],vhat[1:nt-1,i]);
sig=(vhat[2:nt,i]-b*vhat[1:nt-1,i])'(vhat[2:nt,i]-b*vhat[1:nt-1,i]);
sig=sig/(nt-1);
a2n=a2n+4*b*b*sig*sig/(1-b)^8;
a2d=a2d+sig*sig/(1-b)^4;
i=i+1;
endo;
a2=a2n/a2d;
st=1.3221*(a2*nt)^.2;
retp(st);
endp;

@********************************************************************@

proc kern(x);
local del,del1,ker;
@procedure to evaluate the quadratic kernel at some value x.@

del=6*pi*x/5;
ker=3*(sin(del)/del-cos(del))/(del*del);
retp(ker);
endp;

@*******************************************************************@

proc 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 interval(y,z,zbar,b,q,m,robust,prewhit,hetomega,hetq,x,p);
local nt,inter,res,qmat,delta,omega,delv,a,bound,i;
local dbdel,dd,bf,qmat1,phi1s,phi2s,eta,cvf,omega1;

@Procedure that compute confidence intervals for the break dates
based on the "shrinking shifts asymptotic framework.@

cvf=zeros(4,1);
nt=rows(z);
if p == 0;
delta=olsqr(y,zbar);
res=y-zbar*delta;
else;
dbdel=olsqr(y,zbar~x);
res=y-(zbar~x)*dbdel;
delta=dbdel[1:(m+1)*q,1];
endif;
bound=zeros(m,4);

bf=zeros(m+2,1);
bf[1,1]=0;
bf[2:m+1,1]=b[1:m,1];
bf[m+2,1]=nt;

i=1;
do while i <= m;
delv=delta[i*q+1:(i+1)*q,1]-delta[(i-1)*q+1:i*q,1];
if robust == 0;
if hetq == 1;

qmat=z[bf[i,1]+1:bf[i+1,1],.]'z[bf[i,1]+1:bf[i+1,1],.]/(bf[i+1,1]-bf[i,1]);
qmat1=z[bf[i+1,1]+1:bf[i+2,1],.]'z[bf[i+1,1]+1:bf[i+2,1],.]
/(bf[i+2,1]-bf[i+1,1]);

else;

qmat=z'z/nt;
qmat1=qmat;
endif;

if hetomega == 1;

phi1s=res[bf[i,1]+1:bf[i+1,1],1]'res[bf[i,1]+1:bf[i+1,1],1]
/(bf[i+1,1]-bf[i,1]);
phi2s=res[bf[i+1,1]+1:bf[i+2,1],1]'res[bf[i+1,1]+1:bf[i+2,1],1]
/(bf[i+2,1]-bf[i+1,1]);

else;

phi1s=res'res/nt;
phi2s=phi1s;

endif;

eta=delv'qmat1*delv/(delv'qmat*delv);

cvf=cvg(eta,phi1s,phi2s);
/*
print "------------------------------------------------------------";
print "The critical values for the construction of the confidence";
print "intervals for the " i "th break are (2.5%,5%,95%,97.5%):";
print cvf';
*/

a=(delv'qmat*delv)/phi1s;

bound[i,1]=b[i,1]-cvf[4,1]/a;
bound[i,2]=b[i,1]-cvf[1,1]/a;
bound[i,3]=b[i,1]-cvf[3,1]/a;
bound[i,4]=b[i,1]-cvf[2,1]/a;

else;

if hetq == 1;

qmat=z[bf[i,1]+1:bf[i+1,1],.]'z[bf[i,1]+1:bf[i+1,1],.]/(bf[i+1,1]-bf[i,1]);
qmat1=z[bf[i+1,1]+1:bf[i+2,1],.]'z[bf[i+1,1]+1:bf[i+2,1],.]
/(bf[i+2,1]-bf[i+1,1]);

else;

qmat=z'z/nt;
qmat1=qmat;
endif;

if hetomega == 1;

omega=correct(z[bf[i,1]+1:bf[i+1,1],.],res[bf[i,1]+1:bf[i+1,1],1],prewhit);
omega1=correct(z[bf[i+1,1]+1:bf[i+2,1],.],res[bf[i+1,1]+1:bf[i+2,1],1],
prewhit);

else;

omega=correct(z,res,prewhit);
omega1=omega;

endif;

phi1s=delv'omega*delv/(delv'qmat*delv);
phi2s=delv'omega1*delv/(delv'qmat1*delv);

eta=delv'qmat1*delv/(delv'qmat*delv);

cvf=cvg(eta,phi1s,phi2s);
/*
print "------------------------------------------------------------";
print "The critical values for the construction of the confidence";
print "intervals for the " i "th break are (2.5%,5%,95%,97.5%):";
print cvf';
*/

a=(delv'qmat*delv)^2/(delv'omega*delv);

bound[i,1]=b[i,1]-cvf[4,1]/a;
bound[i,2]=b[i,1]-cvf[1,1]/a;
bound[i,3]=b[i,1]-cvf[3,1]/a;
bound[i,4]=b[i,1]-cvf[2,1]/a;
endif;
i=i+1;
endo;

bound[.,1]=int(bound[.,1]);
bound[.,2]=int(bound[.,2])+1;
bound[.,3]=int(bound[.,3]);
bound[.,4]=int(bound[.,4])+1;

retp(bound);
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 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 pvdel(y,z,i,q,bigt,b,prewhit,robust,x,p,withb,hetdat,hetvar);
local zbar,delv,vdel,hac,res,j,sig,sigmq,lambda,q22,reg;
local vv,ev,q21,q11,aaa,bbb,xbar,aa,bb,cc,dd,mat,ff,gam,ie,w,wbar,ww,gg;

@procedure that compute the covariance matrix of the estimates delta.@
ev=ones(i+1,1);
{zbar}=pzbar(z,i,b);

if p == 0;
   delv=olsqr(y,zbar);
   res=y-zbar*delv;
   reg=zbar;
else;
   delv=olsqr(y,zbar~x);
   res=y-(zbar~x)*delv;

   if withb == 0;
      reg=zbar-x*inv(x'x)*x'zbar;
   else;
      reg=x~zbar;
   endif;
endif;


if robust == 0;
@
section on testing with no serial correlation in errors
@

if p == 0;

     if hetdat == 1 and hetvar == 0;
        sig=res'res/bigt;
        vdel=sig*inv(reg'reg);
     endif;
     if hetdat == 1 and hetvar == 1;
         sig=psigmq(res,b,q,i,bigt);
          vdel=(sig.*.eye(q))*inv(reg'reg);
      endif;
     if hetdat == 0 and hetvar == 0;
        lambda=plambda(b,i,bigt);
        sig=res'res/bigt;
        vdel=sig*inv(lambda.*.(z'z));
     endif;
     if hetdat ==  0 and hetvar == 1;
        lambda=plambda(b,i,bigt);
        sig=psigmq(res,b,q,i,bigt);
        vdel=(sig.*.eye(q))*inv(lambda.*.(z'z));
     endif;

else;

if hetdat == 0;
print "the case hetdat = 0 is not allowed.";
print "vdel is returned as zeros.";
vdel=zeros(q*(i+1),q*(i+1));
retp(vdel);
endif;


if hetdat == 1 and hetvar == 0;

   sig=res'res/bigt;
   vdel=sig*inv(reg'reg);
endif;

if hetdat == 1 and hetvar == 1;


wbar=pzbar(reg,i,b);
ww=wbar'wbar;
sig=psigmq(res,b,q,i,bigt);
gg=zeros((i+1)*q+p*withb,(i+1)*q+p*withb);
ie=1;
do while ie <= i+1;
gg=gg+sig[ie,ie]*ww[(ie-1)*((i+1)*q+p*withb)+1:ie*((i+1)*q+p*withb),(ie-1)*((i+1)*q+p*withb)+1:ie*((i+1)*q+p*withb)];
ie=ie+1;
endo;
vdel=inv(reg'reg)*gg*inv(reg'reg);




endif;

endif;

else;           @robust=1@

if hetdat == 0;
print "the case hetdat = 0 is case is not allowed.";
print "vdel is returned as zeros.";
vdel=zeros(q*(i+1),q*(i+1));
retp(vdel);
endif;

if p == 0;

    if hetvar == 1;
      hac=zeros((i+1)*q,(i+1)*q);
      vdel=zeros((i+1)*q,(i+1)*q);
      hac[1:q,1:q]=b[1,1]*correct(z[1:b[1,1],.],res[1:b[1,1],1],prewhit);
       j=2;
do while j <= i;
      hac[(j-1)*q+1:j*q,(j-1)*q+1:j*q]=(b[j,1]-b[j-1,1])*correct(z[b[j-1,1]+1:b[j,1],.],res[b[j-1,1]+1:b[j,1],1],prewhit);
       j=j+1;
endo;
    hac[i*q+1:(i+1)*q,i*q+1:(i+1)*q]=(bigt-b[i,1])*correct(z[b[i,1]+1:bigt,.],res[b[i,1]+1:bigt,1],prewhit);
    vdel=inv(reg'reg)*hac*inv(reg'reg);
         else;
      hac=correct(z,res,prewhit);
      lambda=plambda(b,i,bigt);
      vdel=bigt*inv(reg'reg)*(lambda.*.hac)*inv(reg'reg);
    endif;


else;

   hac=correct(reg,res,prewhit);
   vdel=bigt*inv(reg'reg)*hac*inv(reg'reg);

endif;
endif;
retp(vdel);
endp;


@********************************************************************@


proc pftest(y,z,i,q,bigt,datevec,prewhit,robust,x,p,hetdat,hetvar);
local ftest,vdel,fstar,rsub,j,rmat,zbar,delta,dbdel;

@procedure that computes the supF test for i breaks.@

@construct the R matrix@
rsub=zeros(i,i+1);
j=1;
do while j <= i;
rsub[j,j]=-1;
rsub[j,j+1]=1;
j=j+1;
endo;
rmat=rsub.*.eye(q);

{zbar}=pzbar(z,i,datevec[1:i,i]);
if p == 0;
delta=olsqr(y,zbar);
else;
dbdel=olsqr(y,zbar~x);
delta=dbdel[1:(i+1)*q,1];
endif;

{vdel}=pvdel(y,z,i,q,bigt,datevec[1:i,i],prewhit,robust,
x,p,0,hetdat,hetvar);

fstar=delta'rmat'inv(rmat*vdel*rmat')*rmat*delta;
ftest=(bigt-(i+1)*q-p)*fstar/(bigt*i);

retp(ftest);
endp;

@********************************************************************@

proc(2)=order(ssrzero,global,bigt,m,q);
local bic,lwz,mbic,mlwz,i,glob;

@
Determination of the order using BIC and the criterion of Liu, Wu and
Zidek.
@
glob=zeros(m+1,1);
glob[1,1]=ssrzero;
glob[2:m+1,1]=global;

bic=zeros(m+1,1);
lwz=zeros(m+1,1);
i=0;
do while i <= m;
bic[i+1,1]=ln(glob[i+1,1]/bigt)+ln(bigt)*i*(q+1)/bigt;
lwz[i+1,1]=ln(glob[i+1,1]/(bigt-(i+1)*q-i))
+(i*(q+1)*.299*(ln(bigt))^2.1)/bigt;
print "Values of BIC and lwz with " i " breaks: " 
bic[i+1,1] lwz[i+1,1];
i=i+1;
endo;
mbic=minindc(bic);
mlwz=minindc(lwz);
print "The number of breaks chosen by BIC is :" mbic-1;
print "The number of breaks chosen by LWZ is :" mlwz-1;

retp(mbic-1,mlwz-1);
endp;

@********************************************************************@

proc getcv2(signif,eps1);
local cv;

@input of the critical values of the supF(l+1|l) test@
cv=zeros(10,10);

if eps1 == .05;

if signif == 1;
cv={
 8.02   9.56  10.45  11.07  11.65  12.07  12.47  12.70  13.07  13.34,
11.02  12.79  13.72  14.45  14.90  15.35  15.81  16.12  16.44  16.58,
13.43  15.26  16.38  17.07  17.52  17.91  18.35  18.61  18.92  19.19,
15.53  17.54  18.55  19.30  19.80  20.15  20.48  20.73  20.94  21.10,
17.42  19.38  20.46  21.37  21.96  22.47  22.77  23.23  23.56  23.81,
19.38  21.51  22.81  23.64  24.19  24.59  24.86  25.27  25.53  25.87,
21.23  23.41  24.51  25.07  25.75  26.30  26.74  27.06  27.46  27.70,
22.92  25.15  26.38  27.09  27.77  28.15  28.61  28.90  29.19  29.49,
24.75  26.99  28.11  29.03  29.69  30.18  30.61  30.93  31.14  31.46,
26.13  28.40  29.68  30.62  31.25  31.81  32.37  32.78  33.09  33.53
};
endif;
if signif == 2;
cv={
 9.63  11.14  12.16  12.83  13.45  14.05  14.29  14.50  14.69  14.88,
12.89  14.50  15.42  16.16  16.61  17.02  17.27  17.55  17.76  17.97,
15.37  17.15  17.97  18.72  19.23  19.59  19.94  20.31  21.05  21.20,
17.60  19.33  20.22  20.75  21.15  21.55  21.90  22.27  22.63  22.83,
19.50  21.43  22.57  23.33  23.90  24.34  24.62  25.14  25.34  25.51,
21.59  23.72  24.66  25.29  25.89  26.36  26.84  27.10  27.26  27.40,
23.50  25.17  26.34  27.19  27.96  28.25  28.64  28.84  28.97  29.14,
25.22  27.18  28.21  28.99  29.54  30.05  30.45  30.79  31.29  31.75,
27.08  29.10  30.24  30.99  31.48  32.46  32.71  32.89  33.15  33.43,
28.49  30.65  31.90  32.83  33.57  34.27  34.53  35.01  35.33  35.65
};
endif;
if signif == 3;
cv={
11.17  12.88  14.05  14.50  15.03  15.37  15.56  15.73  16.02  16.39,
14.53  16.19  17.02  17.55  17.98  18.15  18.46  18.74  18.98  19.22,
17.17  18.75  19.61  20.31  21.33  21.59  21.78  22.07  22.41  22.73,
19.35  20.76  21.60  22.27  22.84  23.44  23.74  24.14  24.36  24.54,
21.47  23.34  24.37  25.14  25.58  25.79  25.96  26.39  26.60  26.84,
23.73  25.41  26.37  27.10  27.42  28.02  28.39  28.75  29.13  29.44,
25.23  27.24  28.25  28.84  29.14  29.72  30.41  30.76  31.09  31.43,
27.21  29.01  30.09  30.79  31.80  32.50  32.81  32.86  33.20  33.60,
29.13  31.04  32.48  32.89  33.47  33.98  34.25  34.74  34.88  35.07,
30.67  32.87  34.27  35.01  35.86  36.32  36.65  36.90  37.15  37.41
};
endif;
if signif == 4;
cv={
13.58  15.03  15.62  16.39  16.60  16.90  17.04  17.27  17.32  17.61,
16.64  17.98  18.66  19.22  20.03  20.87  20.97  21.19  21.43  21.74,
19.25  21.33  22.01  22.73  23.13  23.48  23.70  23.79  23.84  24.59,
21.20  22.84  24.04  24.54  24.96  25.36  25.51  25.58  25.63  25.88,
23.99  25.58  26.32  26.84  27.39  27.86  27.90  28.32  28.38  28.39,
25.95  27.42  28.60  29.44  30.18  30.52  30.64  30.99  31.25  31.33,
28.01  29.14  30.61  31.43  32.56  32.75  32.90  33.25  33.25  33.85,
29.60  31.80  32.84  33.60  34.23  34.57  34.75  35.01  35.50  35.65,
31.66  33.47  34.60  35.07  35.49  37.08  37.12  37.23  37.47  37.68,
33.62  35.86  36.68  37.41  38.20  38.70  38.91  39.09  39.11  39.12
};
endif;
endif;


if eps1 == .10;

if signif == 1;
cv={
  7.42     9.05     9.97   10.49   10.91   11.29   11.86   12.26   12.57   12.84, 
10.37   12.19   13.20   13.79   14.37   14.68   15.07   15.42   15.81   16.09,
12.77   14.54   15.64   16.46   16.94   17.35   17.68   17.93   18.35   18.55, 
14.81   16.70   17.84   18.51   19.13   19.50   19.93   20.15   20.46   20.67,  
16.65   18.61   19.74   20.46   21.04   21.56   21.96   22.46   22.72   22.96,   
18.65   20.63   22.03   22.90   23.57   24.08   24.38   24.73   25.10   25.29,   
20.34   22.55   23.84   24.59   24.97   25.48   26.18   26.48   26.86   26.97,   
22.01   24.24   25.49   26.31   26.98   27.55   27.92   28.16   28.64   28.89,  
23.79   26.14   27.34   28.16   28.83   29.33   29.86   30.23   30.46   30.74,    
25.29   27.59   28.75   29.71   30.35   30.99   31.41   31.82   32.25   32.61 
};
endif;

if signif == 2;  
cv={
  9.10   10.55   11.36   12.35   12.97   13.45   13.88   14.12   14.45   14.51,   
12.25   13.83   14.73   15.46   16.13   16.55   16.82   17.07   17.34   17.58,
14.60   16.53   17.43   17.98   18.61   19.02   19.25   19.61   19.94   20.35,  
16.76   18.56   19.53   20.24   20.72   21.13   21.55   21.83   22.08   22.40,  
18.68   20.57   21.60   22.55   23.00   23.63   24.13   24.48   24.82   25.14,   
20.76   23.01   24.14   24.77   25.48   25.89   26.25   26.77   26.96   27.14,  
22.62   24.64   25.57   26.54   27.04   27.51   28.14   28.44   28.74   28.87,   
24.34   26.42   27.66   28.25   28.99   29.34   29.86   30.29   30.50   30.68,   
26.20   28.23   29.44   30.31   30.77   31.35   31.91   32.60   32.71   32.86,   
27.64   29.78   31.02   31.90   32.71   33.32   33.95   34.29   34.52   34.81
};
endif;
if signif == 3;   
cv={
10.56   12.37   13.46   14.13   14.51   14.88   15.37   15.47   15.62   15.79,  
13.86   15.51   16.55   17.07   17.58   17.98   18.19   18.55   18.92   19.02,   
16.55   17.99   19.06   19.65   20.35   21.40   21.57   21.76   22.07   22.53,   
18.62   20.30   21.18   21.86   22.40   22.83   23.42   23.63   23.77   24.14,   
20.59   22.57   23.66   24.50   25.14   25.46   25.77   25.87   26.02   26.34,   
23.05   24.79   25.91   26.80   27.14   27.42   27.85   28.10   28.55   28.89,  
24.65   26.56   27.53   28.51   28.87   29.08   29.43   29.85   30.35   30.68,    
26.50   28.29   29.36   30.34   30.68   31.82   32.42   32.64   32.82   33.08,   
28.25   30.31   31.41   32.60   32.86   33.39   33.79   34.00   34.35   34.75,   
29.80   31.90   33.34   34.31   34.81   35.65   36.23   36.36   36.65   36.72 
};
endif;
if signif == 4;
cv={
13.00   14.51   15.44   15.73   16.39   16.60   16.78   16.90   16.99   17.04,   
16.19   17.58   18.31   18.98   19.63   20.09   20.30   20.87   20.97   21.13,   
18.72   20.35   21.60   22.35   22.96   23.37   23.53   23.71   23.79   23.84,   
20.75   22.40   23.55   24.13   24.54   24.96   25.11   25.50   25.56   25.58,   
23.12   25.14   25.79   26.32   26.60   26.96   27.39   27.51   27.75   27.75,   
25.50   27.14   27.92   28.75   29.44   30.12   30.18   30.29   30.52   30.64,   
27.19   28.87   29.51   30.43   31.38   32.56   32.62   32.87   32.90   33.25,   
29.01   30.68   32.52   32.86   33.27   34.10   34.26   34.38   34.57   34.72,   
30.81   32.86   33.92   34.60   35.07   35.66   37.08   37.12   37.22   37.23,   
32.80   34.81   36.32   36.65   37.15   38.20   38.60   38.70   38.80   39.09 
};
endif;
endif;  

if eps1 == .15;

if signif == 1;
cv={
 7.04      8.51     9.41   10.04   10.58   11.03   11.43   11.75   12.01   12.20, 
 9.81    11.40   12.29   12.90   13.47   13.98   14.36   14.70   15.11   15.28,   
12.08   13.91   14.96   15.68   16.35   16.81   17.24   17.51   17.87   18.12,   
14.26   16.11   17.31   18.00   18.45   18.84   19.22   19.61   19.92   20.07,   
16.14   18.14   19.10   19.84   20.50   20.96   21.42   21.68   21.95   22.28,   
17.97   20.01   21.16   22.08   22.64   23.02   23.35   23.70   24.10   24.37,   
19.70   21.79   22.87   24.06   24.68   25.10   25.66   25.97   26.29   26.50,   
21.41   23.62   24.74   25.63   26.39   26.73   27.29   27.56   28.06   28.46,   
23.06   25.54   26.68   27.60   28.25   28.79   29.19   29.52   29.94   30.43,   
24.65   26.92   28.26   29.18   29.88   30.40   30.90   31.40   31.75   32.03
};
endif;
if signif == 2;
cv={   
  8.58   10.13   11.14   11.83   12.25   12.66   13.08   13.35   13.75   13.89,   
11.47   12.95   14.03   14.85   15.29   15.80   16.16   16.44   16.77   16.84,   
13.98   15.72   16.83   17.61   18.14   18.74   19.09   19.41   19.68   19.77,   
16.19   18.11   18.93   19.64   20.19   20.54   21.21   21.42   21.72   21.97,  
18.23   19.91   20.99   21.71   22.37   22.77   23.15   23.42   24.04   24.42,   
20.08   22.11   23.04   23.77   24.43   24.75   24.96   25.22   25.61   25.93,    
21.87   24.17   25.13   26.03   26.65   27.06   27.37   27.90   28.18   28.36,
23.70   25.75   26.81   27.65   28.48   28.80   29.08   29.30   29.50   29.69,     
25.65   27.66   28.91   29.67   30.52   30.96   31.48   31.77   31.94   32.33,   
27.03   29.24   30.45   31.45   32.12   32.50   32.84   33.12   33.22   33.85 
};
endif;
if signif == 3;
cv={  
10.18   11.86   12.66   13.40   13.89   14.32   14.73   14.89   15.22   15.29,   
12.96   14.92   15.81   16.51   16.84   17.18   17.61   17.84   18.32   18.76,   
15.76   17.70   18.87   19.42   19.77   20.45   20.57   20.82   21.51   22.00,   
18.13   19.70   20.66   21.46   21.97   22.52   22.79   22.82   23.03   23.13,  
19.95   21.72   22.81   23.47   24.42   24.83   25.28   25.59   25.98   26.29,    
22.15   23.79   24.76   25.22   25.93   26.58   26.99   27.11   27.40   27.76,
24.20   26.03   27.06   27.91   28.36   28.72   29.17   29.43   29.66   30.00,     
25.77   27.72   28.80   29.33   29.69   30.02   30.46   30.74   30.90   31.07,    
27.69   29.67   31.00   31.78   32.33   33.06   33.51   33.68   34.16   34.58,   
29.27   31.47   32.54   33.15   33.85   34.32   34.45   34.76   34.94   35.15 
};
endif;
if signif == 4;
cv={  
12.29   13.89   14.80   15.28   15.76   16.27   16.63   16.77   16.81   17.01,   
15.37   16.84   17.72   18.67   19.17   19.46   19.74   19.93   20.12   20.53,   
18.26   19.77   20.75   21.98   22.46   22.69   22.93   23.11   23.12   23.15,   
20.23   21.97   22.80   23.06   23.76   24.55   24.85   25.11   25.53   25.57,   
22.40   24.42   25.53   26.17   26.53   26.77   26.96   27.10   27.35   27.37,   
24.45   25.93   27.09   27.56   28.20   29.61   29.62   30.27   30.45   30.56,   
26.71   28.36   29.30   29.86   30.52   30.89   30.95   31.03   31.11   31.17,   
28.51   29.69   30.65   31.03   31.87   32.42   32.67   33.00   33.11   33.45,   
30.62   32.33   33.51   34.28   34.94   35.71   36.03   36.34   36.48   36.49,   
32.16   33.85   34.58   35.14   36.15   36.76   36.92   37.37   37.87   37.96
};
endif;
endif;   


if eps1 == .20;

if signif == 1;
cv={
  6.72     8.13     9.07     9.66   10.17   10.59   10.95   11.28   11.64   11.89,   
  9.37   10.92   11.90   12.50   12.89   13.38   13.84   14.15   14.41   14.66,  
11.59   13.43   14.43   15.16   15.72   16.24   16.69   16.95   17.32   17.42,    
13.72   15.59   16.67   17.53   18.17   18.52   18.84   19.12   19.43   19.67, 
15.51   17.59   18.76   19.43   20.02   20.53   20.91   21.21   21.59   21.70,    
17.39   19.49   20.65   21.37   22.07   22.57   22.90   23.12   23.38   23.63,    
19.11   21.24   22.42   23.20   24.13   24.68   25.00   25.31   25.76   26.03,   
20.86   23.09   24.30   25.14   25.76   26.27   26.59   27.06   27.41   27.58,   
22.38   24.80   26.10   26.88   27.47   28.05   28.40   28.79   29.16   29.51,   
23.95   26.33   27.50   28.50   29.13   29.52   30.07   30.43   30.87   31.17
};
endif;
if signif == 2;
cv={   
  8.22     9.71   10.66   11.34   11.93   12.30   12.68   12.92   13.21   13.61,   
10.98   12.55   13.46   14.22   14.78   15.37   15.81   16.13   16.44   16.69,   
13.47   15.25   16.36   17.08   17.51   18.08   18.44   18.89   19.01   19.35,   
15.67   17.61   18.54   19.21   19.80   20.22   20.53   21.06   21.31   21.55,  
17.66   19.50   20.63   21.40   21.72   22.19   22.72   23.01   23.24   23.67,  
19.55   21.44   22.64   23.19   23.75   24.28   24.46   24.75   24.96   25.02,    
21.33   23.31   24.75   25.38   26.10   26.47   26.87   27.15   27.37   27.74,   
23.19   25.23   26.39   27.19   27.63   28.09   28.49   28.70   28.83   29.02,   
24.91   26.92   28.10   28.93   29.64   30.29   30.87   31.09   31.39   31.67,   
26.38   28.56   29.62   30.48   31.23   31.96   32.20   32.38   32.72   32.90
};
endif;   
if signif == 3;
cv={
  9.77   11.34   12.31   12.99   13.61   13.87   14.25   14.37   14.73   14.86, 
12.59   14.22   15.39   16.14   16.69   17.00   17.18   17.53   17.65   17.83,    
15.28   17.08   18.10   18.91   19.35   19.70   20.00   20.21   20.53   20.72,    
17.67   19.22   20.25   21.19   21.55   21.88   22.18   22.52   22.77   22.82,   
19.51   21.42   22.28   23.04   23.67   24.20   24.47   24.79   24.94   25.28,   
21.47   23.21   24.28   24.76   25.02   25.70   26.07   26.43   26.73   26.95,  
23.36   25.47   26.47   27.20   27.74   28.21   28.40   28.63   29.09   29.29,   
25.26   27.19   28.10   28.70   29.02   29.41   29.62   29.91   30.11   30.46,    
26.96   28.98   30.34   31.13   31.67   31.89   32.26   32.84   33.14   33.51,
28.62   30.50   31.97   32.39   32.90   33.20   33.90   34.33   34.53   34.76 
};
endif;
if signif == 4;
cv={     
11.94   13.61   14.31   14.80   15.26   15.76   15.87   16.23   16.33   16.63,   
14.92   16.69   17.41   17.72   18.27   19.06   19.17   19.23   19.54   19.74,   
17.60   19.35   20.02   20.64   21.23   21.98   22.19   22.54   22.90   22.93,   
19.82   21.55   22.27   22.80   23.06   23.76   23.97   24.55   24.78   24.85,   
21.75   23.67   24.60   25.18   25.76   26.29   26.42   26.53   26.65   26.67,   
23.80   25.02   26.24   26.77   27.27   27.76   28.12   28.48   28.56   28.80,   
26.16   27.74   28.50   29.17   29.66   30.52   30.66   30.89   30.93   30.95,   
27.71   29.02   29.71   30.20   30.78   31.03   31.80   32.42   32.42   32.47,   
29.67   31.67   32.52   33.28   33.81   34.81   35.22   35.54   35.71   36.03,   
31.38   32.90   34.12   34.68   35.00   36.15   36.76   36.92   37.14   37.37 
};
endif;
endif;  


if eps1 == .25;

if signif == 1;
cv={
  6.35     7.79     8.70     9.22     9.71   10.06   10.45   10.89   11.16   11.30,  
  8.96   10.50   11.47   12.13   12.56   12.94   13.29   13.76   14.03   14.22,    
11.17   12.96   13.96   14.58   15.13   15.54   15.93   16.47   16.79   16.96, 
13.22   15.16   16.14   16.94   17.52   17.97   18.34   18.67   18.84   19.04,   
14.98   16.98   18.12   18.87   19.47   19.90   20.47   20.74   21.00   21.44,    
16.77   18.88   20.03   20.83   21.41   21.83   22.28   22.58   22.83   23.04,    
18.45   20.69   21.81   22.73   23.49   24.19   24.60   24.87   25.08   25.60,   
20.15   22.51   23.56   24.42   25.11   25.61   25.95   26.43   26.59   26.90,  
21.69   24.08   25.45   26.19   26.79   27.33   27.78   28.23   28.54   28.99,   
23.29   25.72   26.97   27.69   28.55   29.13   29.48   29.90   30.30   30.75
};
endif;
if signif == 2;
cv={    
  7.86     9.29   10.12   10.93   11.37   11.82   12.20   12.65   12.79   13.09,
10.55   12.19   12.97   13.84   14.32   14.92   15.28   15.48   15.87   16.34, 
13.04   14.65   15.60   16.51   17.08   17.39   17.76   18.08   18.32   18.72, 
15.19   17.00   18.10   18.72   19.14   19.63   20.10   20.50   20.98   21.23,    
17.12   18.94   20.02   20.81   21.45   21.72   22.10   22.69   22.98   23.15,  
18.97   20.89   21.92   22.66   23.09   23.42   23.96   24.28   24.46   24.75,  
20.75   22.78   24.24   24.93   25.66   26.03   26.28   26.56   26.87   27.21,     
22.56   24.54   25.71   26.50   27.01   27.51   27.74   28.09   28.48   28.70, 
24.18   26.28   27.42   28.27   29.03   29.67   30.34   30.79   30.93   31.13, 
25.77   27.75   29.18   30.02   30.83   31.40   31.92   32.20   32.38   32.72
};
endif;
if signif == 3;
cv={          
  9.32   10.94   11.86   12.66   13.09   13.51   13.85   14.16   14.37   14.70,   
12.21   13.85   14.94   15.48   16.34   16.55   16.80   16.82   17.06   17.34,   
14.66   16.56   17.40   18.12   18.72   19.01   19.40   19.73   20.02   20.50,   
17.04   18.73   19.64   20.52   21.23   21.71   21.95   22.24   22.56   22.79,   
18.96   20.83   21.75   22.69   23.15   23.82   24.20   24.43   24.77   24.83,   
20.93   22.70   23.49   24.35   24.75   25.02   25.58   25.83   26.30   26.68,  
22.85   24.93   26.05   26.66   27.21   27.75   27.99   28.36   28.61   29.09, 
24.56   26.51   27.52   28.10   28.70   29.01   29.46   29.69   29.93   30.11,    
26.31   28.28   29.67   30.81   31.13   31.73   32.26   32.84   33.14   33.28,   
27.80   30.07   31.40   32.20   32.72   33.00   33.20   34.02   34.37   34.68
};
endif;
if signif == 4;
cv={   
 11.44   13.09   14.02   14.63   14.89   15.29   15.76   16.13   16.17   16.23,   
 14.34   16.34   16.81   17.18   17.61   17.83   17.85   18.32   18.67   19.06,   
 17.08   18.72   19.58   20.45   20.72   21.27   21.98   22.46   22.54   22.57,   
 19.22   21.23   22.07   22.61   22.89   23.17   23.77   23.97   24.55   24.78,   
 21.51   23.15   24.36   24.82   25.18   25.76   25.98   26.42   26.43   26.53,   
 23.12   24.75   25.73   26.58   26.99   27.44   27.56   28.00   28.12   28.56,   
 25.67   27.21   28.21   28.80   29.43   29.86   30.38   30.55   30.71   30.89,   
 27.10   28.70   29.53   30.02   30.74   31.01   31.80   32.42   32.42   32.47,   
 29.12   31.13   32.52   33.25   33.62   34.65   34.81   35.22   35.54   35.71,   
 30.86   32.72   33.54   34.58   34.94   35.58   36.15   36.91   36.92   37.14 
};
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,robust,prewhit,z,y,x,p,hetdat,hetvar,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],prewhit,robust,x,p,hetdat,hetvar);

if ftest < cv[q,1];
     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],prewhit,robust,0,p,hetdat,hetvar);
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],prewhit,robust,
x[dv[is,1]+1:dv[is+1,1],.],p,hetdat,hetvar);
endif;

else;
ftestv[is,1]=0.0;
endif;
is=is+1;
endo;
maxf=maxc(ftestv[1:nseg,1]);

if maxf < cv[q,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(2)=spflp1(bigvec,dt,nseg,y,z,h,prewhit,robust,x,p,hetdat,hetvar);
local length,news,dv,is,ssr,ds,ssr0,ftestv,maxf,newd,bigt,in;

@procedure that computes the supF(l+1|l) test l=nseg-1. The l breaks
under the null are taken from the global minimization (in dt).@

ssr=zeros(nseg,1);
ftestv=zeros(nseg,1);
bigt=rows(z);

ftestv=zeros(nseg,1);
dv=zeros(nseg+1,1);
dv[2:nseg,1]=dt;
dv[nseg+1,1]=bigt;
ds=zeros(nseg,1);

in=0;
is=1;
do while is <= nseg;
length=dv[is+1,1]-dv[is,1];
if length >= 2*h;

if p == 0;
{ssr[is,1],ds[is,1]}=parti(dv[is,1]+1,dv[is,1]+h,
          dv[is+1,1]-h,dv[is+1,1],bigvec,bigt);
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]-dv[is,1],prewhit,robust,
0,p,hetdat,hetvar);
else;
{ssr[is,1],ds[is,1]}=onebp(y,z,x,h,dv[is,1]+1,dv[is+1,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]-dv[is,1],prewhit,robust,
x[dv[is,1]+1:dv[is+1,1],.],p,hetdat,hetvar);
endif;


else;
in=in+1;
ftestv[is,1]=0.0;
endif;
is=is+1;
endo;
if in == nseg;
print "Given the location of the breaks from the global optimization";
print "with " nseg-1 "breaks there was no more place to insert ";
print "an additional breaks that satisfy the minimal length requirement.";
endif;
maxf=maxc(ftestv[1:nseg,1]);
news=maxindc(ftestv[1:nseg,1]);
newd=ds[news,1];
retp(maxf,newd);
endp;

@********************************************************************@

proc getcv1(signif,eps1);
local cv;



@input of the critical values of the supF(k) test@

if eps1 == .05;
cv=zeros(10,10);

if signif == 1;
cv={
 8.02   7.87   7.07   6.61   6.14   5.74   5.40   5.09   4.81,
11.02  10.48   9.61   8.99   8.50   8.06   7.66   7.32   7.01,
13.43  12.73  11.76  11.04  10.49  10.02   9.59   9.21   8.86,
15.53  14.65  13.63  12.91  12.33  11.79  11.34  10.93  10.55,
17.42  16.45  15.44  14.69  14.05  13.51  13.02  12.59  12.18,
19.38  18.15  17.17  16.39  15.74  15.18  14.63  14.18  13.74,
21.23  19.93  18.75  17.98  17.28  16.69  16.16  15.69  15.24,
22.92  21.56  20.43  19.58  18.84  18.21  17.69  17.19  16.70,
24.75  23.15  21.98  21.12  20.37  19.72  19.13  18.58  18.09,
26.13  24.70  23.48  22.57  21.83  21.16  20.57  20.03  19.55};

endif;
if signif == 2;
cv={
 9.63   8.78   7.85   7.21   6.69   6.23   5.86   5.51   5.20,
12.89  11.60  10.46   9.71   9.12   8.65   8.19   7.79   7.46,
15.37  13.84  12.64  11.83  11.15  10.61  10.14   9.71   9.32,
17.60  15.84  14.63  13.71  12.99  12.42  11.91  11.49  11.04,
19.50  17.60  16.40  15.52  14.79  14.19  13.63  13.16  12.70,
21.59  19.61  18.23  17.27  16.50  15.86  15.29  14.77  14.30,
23.50  21.30  19.83  18.91  18.10  17.43  16.83  16.28  15.79,
25.22  23.03  21.48  20.46  19.66  18.97  18.37  17.80  17.30,
27.08  24.55  23.16  22.08  21.22  20.49  19.90  19.29  18.79,
28.49  26.17  24.59  23.59  22.71  21.93  21.34  20.74  20.17};

endif;
if signif == 3;
cv={
11.17   9.81   8.52   7.79   7.22   6.70   6.27   5.92   5.56,
14.53  12.64  11.20  10.29   9.69   9.10   8.64   8.18   7.80,
17.17  14.91  13.44  12.49  11.75  11.13  10.62  10.14   9.72,
19.35  16.85  15.44  14.43  13.64  13.01  12.46  11.94  11.49,
21.47  18.75  17.26  16.13  15.40  14.75  14.19  13.66  13.17,
23.73  20.80  19.15  18.07  17.21  16.49  15.84  15.29  14.78,
25.23  22.54  20.85  19.68  18.79  18.03  17.38  16.79  16.31,
27.21  24.20  22.41  21.29  20.39  19.63  18.98  18.34  17.78,
29.13  25.92  24.14  22.97  21.98  21.28  20.59  19.98  19.39,
30.67  27.52  25.69  24.47  23.45  22.71  21.95  21.34  20.79};

endif;

if signif == 4;
cv={
13.58  10.95   9.37   8.50   7.85   7.21   6.75   6.33   5.98,
16.64  13.78  12.06  11.00  10.28   9.65   9.11   8.66   8.22,
19.25  16.27  14.48  13.40  12.56  11.80  11.22  10.67  10.19,
21.20  18.21  16.43  15.21  14.45  13.70  13.04  12.48  12.02,
23.99  20.18  18.19  17.09  16.14  15.34  14.81  14.26  13.72,
25.95  22.18  20.29  18.93  17.97  17.20  16.54  15.94  15.35,
28.01  24.07  21.89  20.68  19.68  18.81  18.10  17.49  16.96,
29.60  25.66  23.44  22.22  21.22  20.40  19.66  19.03  18.46,
31.66  27.42  25.13  24.01  23.06  22.18  21.35  20.63  19.94,
33.62  29.14  26.90  25.58  24.44  23.49  22.75  22.09  21.47};
else;

endif;
endif;


if eps1 == .10;

cv=zeros(10,8);
if signif == 1;
cv={
  7.42    6.93    6.09    5.44    4.85    4.32    3.83    3.22   ,
10.37    9.43    8.48    7.68    7.02    6.37    5.77    4.98 ,
12.77   11.61   10.53    9.69    8.94    8.21    7.49    6.57 ,
14.81   13.56   12.36   11.43   10.61    9.86    9.04    8.01  ,  
16.65   15.32   14.06   13.10   12.20   11.40   10.54    9.40 ,    
18.65   17.01   15.75   14.70   13.78   12.92   11.98   10.80 , 
20.34   18.71   17.26   16.19   15.26   14.35   13.40   12.13 ,    
22.01   20.32   18.90   17.75   16.79   15.82   14.80   13.45 ,  
23.79   21.88   20.43   19.28   18.22   17.24   16.19   14.77 ,
25.29   23.33   21.89   20.71   19.63   18.59   17.50   16.00 
};
endif;
if signif == 2;
cv={                 
   9.10    7.92    6.84    6.03    5.37    4.80    4.23    3.58  ,
 12.25   10.58    9.29    8.37    7.62    6.90    6.21    5.41 ,
 14.60   12.82   11.46   10.41    9.59    8.80    8.01    7.03  ,
 16.76   14.72   13.30   12.25   11.29   10.42    9.58    8.46 ,   
 18.68   16.50   15.07   13.93   13.00   12.10   11.16    9.96 ,  
 20.76   18.32   16.81   15.67   14.65   13.68   12.63   11.34 ,  
 22.62   20.04   18.45   17.19   16.14   15.11   14.09   12.71 ,    
 24.34   21.69   20.01   18.74   17.66   16.65   15.54   14.07 ,  
 26.20   23.36   21.63   20.32   19.19   18.09   16.89   15.40 ,  
 27.64   24.87   23.11   21.79   20.58   19.47   18.29   16.70 
};
endif;
if signif == 3;
cv={           
10.56    8.90    7.55    6.64    5.88    5.22    4.61    3.90  ,   
13.86   11.63   10.14    9.05    8.17    7.40    6.63    5.73  ,  
16.55   13.90   12.35   11.12   10.19    9.28    8.43    7.40  ,  
18.62   15.88   14.22   12.96   11.94   11.05   10.06    8.93  ,   
20.59   17.71   16.02   14.68   13.67   12.71   11.68   10.42  ,      
23.05   19.69   17.82   16.47   15.31   14.24   13.20   11.89  ,  
24.65   21.34   19.41   18.13   16.90   15.84   14.67   13.25  ,   
26.50   22.98   20.95   19.69   18.52   17.35   16.15   14.67 ,   
28.25   24.73   22.68   21.29   20.01   18.76   17.56   16.00  ,      
29.80   26.37   24.27   22.71   21.42   20.21   18.94   17.33   
};
endif;
if signif == 4;
cv={       
 13.00   10.14    8.42    7.31    6.48    5.74    5.05    4.28  ,  
 16.19   12.90   11.12    9.87    8.84    8.01    7.18    6.18  ,   
 18.72   15.38   13.38   11.97   10.93    9.94    8.99    7.85  ,    
 20.75   17.24   15.30   13.93   12.78   11.67   10.64    9.47  ,   
 23.12   18.93   16.91   15.61   14.42   13.31   12.30   11.00  ,    
 25.50   21.15   19.04   17.48   16.19   15.11   13.88   12.55  ,     
 27.19   22.97   20.68   19.14   17.81   16.59   15.43   13.92 ,   
 29.01   24.51   22.40   20.68   19.41   18.08   16.83   15.30  ,    
 30.81   26.30   23.95   22.33   20.88   19.56   18.35   16.79  ,  
 32.80   28.24   25.63   23.83   22.32   21.04   19.73   18.10  
};
endif;
endif;  


if eps1 == .15;

cv=zeros(10,5);
if signif == 1;
cv={
 7.04    6.28    5.21    4.41    3.47   ,
 9.81    8.63    7.54    6.51    5.27  ,
12.08   10.75    9.51    8.29    6.90  ,     
14.26   12.60   11.21    9.97    8.37  , 
16.14   14.37   12.90   11.50    9.79 ,  
17.97   16.02   14.45   13.00   11.19 ,     
19.70   17.67   16.04   14.55   12.59 ,   
21.41   19.16   17.47   15.88   13.89 ,  
23.06   20.82   19.07   17.38   15.23 ,  
24.65   22.26   20.42   18.73   16.54  
};
endif;
if signif == 2;
cv={                                                 
  8.58    7.22    5.96    4.99    3.91 ,
11.47    9.75    8.36    7.19    5.85  , 
13.98   11.99   10.39    9.05    7.46 ,  
16.19   13.77   12.17   10.79    9.09 ,    
18.23   15.62   13.93   12.38   10.52  ,   
20.08   17.37   15.58   13.90   11.94 ,   
21.87   18.98   17.23   15.55   13.40 ,
23.70   20.62   18.69   16.96   14.77 ,   
25.65   22.35   20.18   18.40   16.11 , 
27.03   23.80   21.62   19.79   17.44 
};
endif;
if signif == 3;
cv={                                          
10.18    8.14    6.72    5.51    4.34  ,   
12.96   10.75    9.15    7.81    6.38  ,
15.76   13.13   11.23    9.72    8.03  , 
18.13   14.99   13.06   11.55    9.66 ,     
19.95   16.92   14.98   13.25   11.21  , 
22.15   18.62   16.50   14.68   12.63  ,   
24.20   20.40   18.25   16.41   14.18  , 
25.77   21.97   19.71   17.91   15.52  ,   
27.69   23.68   21.28   19.29   16.88  ,   
29.27   24.99   22.74   20.81   18.26  
};
endif;
if signif == 4; 
cv={
 12.29    9.36    7.60    6.19    4.91  ,               
 15.37   12.15   10.27    8.65    7.00 ,               
 18.26   14.45   12.16   10.56    8.71 ,               
 20.23   16.55   14.26   12.42   10.53 ,               
 22.40   18.37   16.16   14.25   12.14  ,              
 24.45   20.06   17.57   15.73   13.44  ,              
 26.71   21.87   19.42   17.44   15.02  ,              
 28.51   23.58   20.96   19.00   16.56  ,              
 30.62   25.32   22.72   20.38   17.87  ,              
 32.16   26.82   24.41   22.09   19.27  
};
endif;
endif;              


if eps1 == .20;

cv=zeros(10,3);
if signif == 1;
cv={
  6.72    5.59    4.37, 
 9.37    7.91    6.43, 
 11.59    9.93    8.21,    
13.72   11.70    9.90, 
15.51   13.46   11.50,      
17.39   15.05   12.91,    
19.11   16.67   14.46,    
 20.86   18.16   15.88,  
 22.38   19.71   17.30,   
 23.95   21.13   18.65
};
endif;
if signif == 2;
cv={                       
  8.22    6.53    5.08,      
 10.98    8.98    7.13,  
13.47   11.09    9.12,    
 15.67   12.94   10.78,    
17.66   14.69   12.45,   
 19.55   16.35   13.91,  
 21.33   18.14   15.55,   
23.19   19.58   17.10,   
24.91   21.23   18.58,   
 26.38   22.62   19.91
};
endif;
if signif == 3;
cv={
  9.77    7.49    5.73,
12.59   10.00    7.92,     
 15.28   12.25    9.91,    
17.67   14.11   11.66,   
 19.51   15.96   13.49,    
 21.47   17.66   14.97,     
 23.36   19.41   16.56,  
 25.26   20.94   18.03,   
 26.96   22.69   19.51,  
28.62   24.04   20.96 
};
endif;                           
if signif == 4;
cv={
 11.94    8.77    6.58,                   
 14.92   11.30    8.95,                  
 17.60   13.40   10.91,               
 19.82   15.74   12.99 ,                  
 21.75   17.21   14.60 ,                   
 23.80   19.25   16.29 ,                
 26.16   21.03   17.81 ,                  
 27.71   22.71   19.37 ,                 
 29.67   24.43   20.74 ,                    
 31.38   25.73   22.34                
};
endif;
endif;

if eps1 == .25;

cv=zeros(10,2);
if signif ==  1;
cv={
  6.35    4.88,
 8.96    7.06, 
11.17    9.01,  
13.22   10.74,  
14.98   12.39,  
16.77   13.96,   
18.45   15.53, 
20.15   16.91,  
21.69   18.42, 
23.29   19.84 
};
endif;
if signif == 2;
cv={
  7.86    5.80,
10.55    8.17,  
13.04   10.16, 
15.19   11.91,  
17.12   13.65, 
18.97   15.38,    
20.75   16.97,   
22.56   18.43,  
24.18   19.93,   
25.77   21.34 
};
endif;
if signif == 3;
cv={                    
  9.32    6.69,
 12.21    9.16,   
14.66   11.22,   
17.04   13.00,  
18.96   14.86, 
20.93   16.53,  
22.85   18.25, 
24.56   19.68, 
26.31   21.38,   
27.80   22.79 
};
endif;
if signif == 4;
cv={                        
 11.44    7.92,
 14.34   10.30,                        
 17.08   12.55,                        
 19.22   14.65,                        
 21.51   16.18,                        
 23.12   18.10,                        
 25.67   19.91,                        
 27.10   21.41,                        
 29.12   23.23,                        
 30.86   24.51
};
endif;
endif;                        



retp(cv);
endp;


@********************************************************************@

proc(1)=preparti(y,z,nbreak,dateseq,h,x,p);
local bigt,q,dv,is,length,vssr,vssrev,dum,ds,dr;

@procedure that construct modified sequential estimates using the
repartition method of Bai (1995), Estimating Breaks one at a time. @

bigt=rows(z);
q=cols(z);

dv=zeros(nbreak+2,1);

dv[1,1]=0;                  @the first date of a segment is set to 0@
dv[2:nbreak+1,1]=dateseq;   @the segments are determined by the dates
                             obtained from the sequential method@
dv[nbreak+2,1]=bigt;        @set the last date of the segment to T@
ds=zeros(nbreak,1);
dr=zeros(nbreak,1);

is=1;
do while is <= nbreak;

length=dv[is+2,1]-dv[is,1];

if p == 0;
{vssr}=ssr(1,y[dv[is,1]+1:dv[is+2,1],1],
z[dv[is,1]+1:dv[is+2,1],.],h,length);
{vssrev}=ssr(1,rev(y[dv[is,1]+1:dv[is+2,1],1]),
rev(z[dv[is,1]+1:dv[is+2,1],.]),h,length);
{dum,ds[is,1]}=partione(h,length-h,length,vssr,vssrev);
dr[is,1]=ds[is,1]+dv[is,1];
else;
{dum,ds[is,1]}=onebp(y,z,x,h,dv[is,1]+1,dv[is+2,1]);
dr[is,1]=ds[is,1];
endif;


is=is+1;
endo;
retp(dr);
endp;

proc(1)= funcg(x,bet,alph,b,deld,gam);
local aa,xb,g;


if x <= 0;

xb=bet*sqrt(abs(x));

if abs(xb) <= 30;

g=-sqrt(-x/(2*pi))*exp(x/8)
-(bet/alph)*exp(-alph*x)*cdfn(-bet*sqrt(abs(x)))
+((2*bet*bet/alph)-2-x/2)*cdfn(-sqrt(abs(x))/2);

else;
aa=ln(bet/alph)-alph*x-xb^2/2-ln(sqrt(2*pi))-ln(xb);
g=-sqrt(-x/(2*pi))*exp(x/8)-exp(aa)*cdfn(-sqrt(abs(x))/2)
+((2*bet*bet/alph)-2-x/2)*cdfn(-sqrt(abs(x))/2);
endif;

else;

xb=deld*sqrt(x);


if abs(xb) <= 30;

g=1+(b/sqrt(2*pi))*sqrt(x)*exp(-b*b*x/8)
+(b*deld/gam)*exp(gam*x)*cdfn(-deld*sqrt(x))
+(2-b*b*x/2-2*deld*deld/gam)*cdfn(-b*sqrt(x)/2);

else;

aa=ln((b*deld/gam))+gam*x-xb^2/2-ln(sqrt(2*pi))-ln(xb);
g=1+(b/sqrt(2*pi))*sqrt(x)*exp(-b*b*x/8)
+exp(aa)
+(2-b*b*x/2-2*deld*deld/gam)*cdfn(-b*sqrt(x)/2);

endif;

endif;

retp(g);
endp;

proc(1)= cvg(eta,phi1s,phi2s);
local a,gam,b,deld,alph,bet;
local cct,isig,sig,upb,lwb,crit,xx,pval,cvec;

cvec=zeros(4,1);
a=phi1s/phi2s;
gam=((phi2s/phi1s)+1)*eta/2;
b=sqrt(phi1s*eta/phi2s);
deld=sqrt(phi2s*eta/phi1s)+b/2;
alph=a*(1+a)/2;
bet=(1+2*a)/2;
sig={0.025 0.05 0.95 0.975};

isig=1;
do while isig <= 4;

upb=2000.0;
lwb=-2000.0;
crit=999999.0;
cct=1;

do while abs(crit) >= .000001;
cct=cct+1;
if cct > 100;
print "the procedure to get critical values for the break dates has";
print "reached the upper bound on the number of iterations. This may";
print "be due to incorrect specifications of the upper or lower bound";
print "in the procedure cvg. The resulting confidence interval for this";
print "break date is incorrect.";
retp(cvec);
else;

xx=lwb+(upb-lwb)/2;
pval=funcg(xx,bet,alph,b,deld,gam);
crit=pval-sig[isig];
if crit <= 0;
lwb=xx;
else;
upb=xx;
endif;
endif;
endo;
cvec[isig,1]=xx;

isig=isig+1;
endo;
retp(cvec);
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;

proc getdmax(signif,eps1);
local cv;

@input of the critical values of the supF(l+1|l) test@
cv=zeros(10,2);
if eps1 == .05;
if signif == 1;
cv={
8.78   9.14,
11.69  12.33,
14.05  14.76,
16.17  16.95,
17.94  18.85,
19.92  20.89,
21.79  22.81,
23.53  24.55,
25.19  26.40,
26.66  27.79
};
endif;
if signif == 2;
cv={
10.17  10.91,
13.27  14.19,
15.80  16.82,
17.88  19.07,
19.74  20.95,
21.90  23.27,
23.77  25.02,
25.51  26.83,
27.28  28.78,
28.75  30.16
};
endif;
if signif == 3;
cv={
11.52  12.53,
14.69  16.04,
17.36  18.79,
19.51  20.89,
21.57  23.04,
23.83  25.22,
25.46  26.92,
27.32  28.98,
29.20  30.82,
30.84  32.46
};
endif;
if signif == 4;
cv={
13.74  15.02,
16.79  18.11,
19.38  20.81,
21.25  22.81,
24.00  25.46,
26.07  27.63,
28.02  29.57,
29.60  31.32,
31.72  33.32,
33.86  35.47
};
endif;
endif;

if eps1 == .10;

if signif == 1;
cv={
8.05  	8.63,
10.86	11.71,
13.26	14.14,
15.23	16.27,
17.06	18.14,
19.06	20.22,
20.76	22.03,
22.42	23.71,
24.24	25.66,
25.64	27.05
};
endif;
if signif == 2;
cv={
9.52 	10.39,
12.59	13.66,
14.85	16.07,
17.00      18.38,
18.91	20.3,
21.01	22.55,
22.8	        24.34,
24.56	26.1,
26.48	27.99,
27.82	29.46
};
endif;
if signif == 3;
cv={
10.83	12.06,
14.15	15.33,
16.64	18.04,
18.75	20.3,
20.68	22.22,
23.25	24.66,
24.75	26.47,
26.54	28.24,
28.33	30.02,
29.9	         31.58
};
endif;
if signif == 4;
cv={
13.07	14.53,
16.19	17.8,
18.75	20.42,
20.75	22.35,
23.16	24.81,
25.55	27.28,
27.23	28.87,
29.01	30.62,
30.81	32.74,
32.82	34.51
};
endif;
endif;

if eps1 == .15;
if signif == 1;
cv={
7.46	     8.2,
10.16  11.15,
12.4	   13.58,
14.58   15.88,
16.49  17.8,
18.23  19.66,
20	  21.46,
21.7	  23.31,
23.38  24.99,
24.9	 26.62
};
endif;
if signif == 2;
cv={
8.88	  9.91,
11.7	  12.81,
14.23  15.59,
16.37  17.83,
18.42  19.96,
20.3	  21.86,
22.04  23.81,
23.87  25.63,
25.81  27.53,
27.23  29.06
};
endif;
if signif == 3;
cv={
10.39  11.67,
13.18  14.58,
15.87  17.41,
18.24  19.82,
20.1	  21.76,
22.27  23.97,
24.26  26.1,
25.88  27.8,
27.78  29.78,
29.36  31.47
};
endif;
if signif == 4;
cv={
12.37  13.83,
15.41  17.01,
18.26  19.86,
20.39  21.95,
22.49   24.5,
24.55  26.68,
26.75  28.76,
28.51  30.4,
30.62  32.71,
32.17   34.25
};
endif;
endif;

if eps1 == .20;
if signif == 1;
cv={
6.96	 7.67,
9.66  10.46,
11.84  12.79,
13.94  15.05,
15.74  16.8,
17.62  18.76,
19.3	  20.56,
21.09  22.45,
22.55  24,
24.17  25.6
};
endif;
if signif == 2;
cv={
8.43	 9.27,
11.16  12.15,
13.66  14.73,
15.79  17.04,
17.76  19.11,
19.69  21.04,
21.46  22.76,
23.28  24.68,
25.04  26.4,
26.51  28.02
};
endif;
if signif == 3;
cv={
9.94	 10.93,
12.68  13.87,
15.31  16.65,
17.73  19.12,
19.59  20.84,
21.56   23.09,
23.4	  25.03,
25.31  26.91,
27.02  28.54,
28.67  30.31
};
endif;
if signif == 4;
cv={
12.02   13.16,
14.92  16.52,
17.6	  18.89,
19.9	  21.27,
21.75  23.39,
23.8	  25.17,
26.16  27.71,
27.71  29.3,
29.67  31.67,
31.4	   32.99
};
endif;
endif;

if eps1 == .25;
if signif == 1;
cv={
6.55	 7.09,
9.16	  9.8,
11.31  12.01,
13.36  14.16,
15.12  15.93,
16.94  17.92,
18.6	  19.61,
20.3	  21.38,
21.81  22.81,
23.43  24.43
};
endif;
if signif == 2;
cv={
8.01  8.69,
10.67  11.49,
13.15  13.99,
15.28  16.13,
17.14  18.11,
19.1	  20.02,
20.84  21.81,
22.62  23.6,
24.28  25.4,
25.8	  27.01
};
endif;
if signif == 3;
cv={
9.37  10.24,
12.25  13.02,
14.67  15.64,
17.13  18.17,
18.97  19.92,
20.97  22.08,
22.89  24.38,
24.61  25.76,
26.38  27.55,
27.91  29.13
};
endif;
if signif == 4;
cv={
11.5	 12.27,
14.34 15.41,
17.08   18.03,
19.22  20.34,
21.51  22.39,
23.12  24.27,
25.67  26.77,
27.12  28.29,
29.12  30.57,
30.87  32.2
};
endif;
endif;




retp(cv);
endp;