*********************************break_test2********************************************
proc(0)=pbreak(bigt,y,z,q,m,h,eps1,p,dotesting,fixb,x,q,eps,maxi,fixb,betaini,printd);

local datevec,bigvec,global,supfl,ndat,maic,mbic,mlwz,ssrzero,nbreak,i,bic,lwz,supff,cv2,k,deltaz,fz,lz,z1,y1,w1,j,wbar,yhat,uhatnull,uhatalt,a,b,c,serrors;
local dateseq,ftest,wftest,cv,repartda,zz,siglev,nbr,datese,cvm,reparv,ii,cvall,size,biclwz,fscv,udmax,wdmax,delta,zbar,theta,yhatnull,lrvar,sc1,sc2,sc3,sdelta1,sdelta2,sdelta3;
local deltak1,deltak2,con91,con92,con93,con94,con1,con2,lrvarnull,scnull,sdeltanull,serrorsnull,llz,ffz,u1,u2,u3,zu1,zu2,zu3;
local rsquare1,rsquare2,rsquare3,yd1,yd2,yd3;

print “The options chosen are:”;
print “h = ” h;
print “eps1 = ” eps1;
print “The maximum number of breaks is: ” m;
print “********************************************************”;
if dotesting==1;
lrvar=zeros(m,1);
print “Output from test procedure”;
print “__________________________________”;
print ” “;
print “The following options are used:”;
print ” “;
print “bigt: ” bigt;
supff=zeros(m,1);
k=1;
do while k<=m;
{supff[k,1],lrvar[k,1]}=wald(y,z,k,q,bigt,p);
k=k+1;
endo;
udmax=maxc(supff);
ftest=supff|udmax;
z1=z[4:bigt-2,.];
deltaz=z1-z[3:bigt-3,.];
lz=z[3:bigt-3,.]-z[2:bigt-4,.];
fz=z[5:bigt-1,.]-z[4:bigt-2,.];
llz=z[2:bigt-4,.]-z[1:bigt-5,.];
ffz=z[6:bigt,.]-z[5:bigt-1,.];
w1=deltaz~lz~fz~llz~ffz;
j=bigt-5;
z1=ones(j,1)~z1;
y1=y[4:bigt-2,1];
{global,datevec,bigvec}=nldat(y1,z1,w1,h,m,p,q,j,fixb,eps,maxi,betaini,printd);
zz=z1~w1;
delta=olsqr(y1,zz);

ssrzero=(y1-zz*delta)’(y1-zz*delta);
{bic,lwz}=order(ssrzero,global,bigt,m,q);
zbar=pzbar(z1,2,datevec[1:2,2]);
wbar=zbar~w1;
theta=olsqr(y1,wbar);
yhat=zeros(j,1);
yhat[1:datevec[1,2],1]=wbar[1:datevec[1,2],.]*theta;
yhat[datevec[1,2]+1:datevec[2,2],1]=wbar[datevec[1,2]+1:datevec[2,2],.]*theta;
yhat[datevec[2,2]+1:j,1]=wbar[datevec[2,2]+1:j,.]*theta;
yhatnull=zz*delta;
uhatalt=y1-wbar*theta;
uhatnull=y1-zz*delta;
print “the null estimated coefficients are” delta;
print “the alt.estimated coefficients are” theta;

theta=olsqr(y1,zbar);
yd1=y1[1:datevec[1,2],1]-meanc(y1[1:datevec[1,2],1]);
u1=y1[1:datevec[1,2],1]-z1[1:datevec[1,2],1]*theta[1,1]-z1[1:datevec[1,2],2]*theta[2,1];

yd2=y1[datevec[1,2]+1:datevec[2,2],1]-meanc(y1[datevec[1,2]+1:datevec[2,2],1]);
u2=y1[datevec[1,2]+1:datevec[2,2],1]-z1[datevec[1,2]+1:datevec[2,2],1]*theta[3,1]-z1[datevec[1,2]+1:datevec[2,2],2]*theta[4,1];

yd3=y1[datevec[2,2]+1:rows(z1),1]-meanc(y1[datevec[2,2]+1:rows(z1),1]);
u3=y1[datevec[2,2]+1:j,1]-z1[datevec[2,2]+1:j,1]*theta[5,1]-z1[datevec[2,2]+1:j,2]*theta[6,1];

 

rsquare1=1-(u1′*u1)/(yd1′*yd1);
rsquare2=1-(u2′*u2)/(yd2′*yd2);
rsquare3=1-(u3′*u3)/(yd3′*yd3);

print “the ftests are” ftest;
print “the BIC order is” bic;
print “the LWZ order is” lwz;
print “the break dates are” datevec;
print “the rsquare in regime 1 is” rsquare1;
print “the rsquare in regime 2 is” rsquare2;
print “the rsquare in regime 3 is” rsquare3;

 

sc1=((lrvar[2,1]/(bigt-5))*inv(datevec[1,2]/(bigt-5)-((bigt-5)^(-3/2)*z1[1:datevec[1,2],2]’*ones(datevec[1,2],1))^2*inv(z1[1:datevec[1,2],2]’*z1[1:datevec[1,2],2])))^0.5;
sc2=((lrvar[2,1]/(bigt-5))*inv((datevec[2,2]-datevec[1,2])/(bigt-5)-((bigt-5)^(-3/2)*z1[datevec[1,2]+1:datevec[2,2],2]’*ones(datevec[2,2]-datevec[1,2],1))^2*inv(z1[datevec[1,2]+1:datevec[2,2],2]’*z1[datevec[1,2]+1:datevec[2,2],2])))^0.5;
sc3=((lrvar[2,1]/(bigt-5))*inv((bigt-5-datevec[2,2])/(bigt-5)-((bigt-5)^(-3/2)*z1[datevec[2,2]+1:bigt-5,2]’*ones(bigt-5-datevec[2,2],1))^2*inv(z1[datevec[2,2]+1:bigt-5,2]’*z1[datevec[2,2]+1:bigt-5,2])))^0.5;
sdelta1=(lrvar[2,1]*inv((z1[1:datevec[1,2],2]-meanc(z1[1:datevec[1,2],2]))’*(z1[1:datevec[1,2],2]-meanc(z1[1:datevec[1,2],2]))))^0.5;
sdelta2=(lrvar[2,1]*inv((z1[datevec[1,2]+1:datevec[2,2],2]-meanc(z1[datevec[1,2]+1:datevec[2,2],2]))’*(z1[datevec[1,2]+1:datevec[2,2],2]-meanc(z1[datevec[1,2]+1:datevec[2,2],2]))))^0.5;
sdelta3=(lrvar[2,1]*inv((z1[datevec[2,2]+1:bigt-5,2]-meanc(z1[datevec[2,2]+1:bigt-5,2]))’*(z1[datevec[2,2]+1:bigt-5,2]-meanc(z1[datevec[2,2]+1:bigt-5,2]))))^0.5;
lrvarnull=lrvar0(y,z,q,bigt,p);
scnull=lrvarnull*(j)^(-1)*inv(1-((j)^(-1.5)*z1[.,2]‘*z1[.,2])^2*inv(j^(-2)*z1[.,2]‘*z1[.,2]));
scnull=scnull^(0.5);
sdeltanull=(lrvarnull*inv((z1[1:bigt-5,2]-meanc(z1[1:bigt-5,2]))’*(z1[1:bigt-5,2]-meanc(z1[1:bigt-5,2]))))^0.5;

serrors=sc1|sc2|sc3|sdelta1|sdelta2|sdelta3;
serrorsnull=scnull|sdeltanull;
print “the standard errors are” serrors;
print “the standard errors under the null are” serrorsnull;

endif;
endp;

 

@*******************************************************************@
proc(2)=wald(y,z,g,q,bigt,p);
local ftest,zbar,gammahat,deltahat,ssr0,ssrk,wbar,w,e,eb,eb1,eb2,ef,ef1,ef2,u,ub,uf,uhat,uhatnull,deltaz,lz,fz,llz,ffz,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,global,datevec,bigvec;

y1=y[4:bigt-2,1];
z1=z[4:bigt-2,.];
deltaz=z1-z[3:bigt-3,.];
lz=z[3:bigt-3,.]-z[2:bigt-4,.];
fz=z[5:bigt-1,.]-z1;
llz=z[2:bigt-4,.]-z[1:bigt-5,.];
ffz=z[6:bigt,.]-z[5:bigt-1,.];
w1=deltaz~lz~fz~llz~ffz;
j=bigt-5;
z1=z1~ones(j,1);
if g==1;
fsample=zeros(j,1);
for b(eps1*j,(1-eps1)*j,1);
zbar=pzbar(z1,1,b);
wbar=zbar~w1;
w=z1~w1;
gammahat=olsqr(y1,wbar);
deltahat=olsqr(y1,w);
ssr0=(y1-w*deltahat)’*(y1-w*deltahat);
ssrk=(y1-wbar*gammahat)’*(y1-wbar*gammahat);
uhat=y1-wbar*gammahat;
uhatnull=y1-w*deltahat;
u=uhat;
e = u;
te = bigt-5;
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*j:(1-eps1)*j,1];
ftest=maxc(fsample);
else;

if p==5;
{global,datevec,bigvec}=nldat(y1,z1,w1,h,m,p,q,j,fixb,eps,maxi,betaini,printd);
endif;
{zbar}=pzbar(z1,g,datevec[1:g,g]);
wbar=zbar~w1;
w=z1~w1;
gammahat=olsqr(y1,wbar);
deltahat=olsqr(y1,w);
ssrk=(y1-wbar*gammahat)’*(y1-wbar*gammahat);
ssr0=(y1-w*deltahat)’*(y1-w*deltahat);
uhat=y1-wbar*gammahat;
uhatnull=y1-w*deltahat;
u=uhat;
e = u;
te = bigt-5;
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,omega[1,1]);
endp;

proc(1)=lrvar0(y,z,q,bigt,p);
local ftest,zbar,gammahat,deltahat,ssr0,ssrk,wbar,w,e,eb,eb1,eb2,ef,ef1,ef2,u,ub,uf,uhat,deltaz,lz,fz,llz,ffz;
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,global,datevec,bigvec;

z1=z[4:bigt-2,.];
deltaz=z1-z[3:bigt-3,.];
lz=z[3:bigt-3,.]-z[2:bigt-4,.];
fz=z[5:bigt-1,.]-z[4:bigt-2,.];
llz=z[2:bigt-4,.]-z[1:bigt-5,.];
ffz=z[6:bigt,.]-z[5:bigt-1,.];
w1=deltaz~lz~fz~llz~ffz;
j=bigt-5;
z1=ones(j,1)~z1;
y1=y[4:bigt-2,1];
w=z1~w1;
deltahat=olsqr(y1,w);
ssr0=(y1-w*deltahat)’*(y1-w*deltahat);
uhat=y1-w*deltahat;
u=uhat;
e = u;
te = bigt-5;
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 = u’*u;
lam = 0;
j = 1;
do while j<=te-1;
lam = lam + (u[1:te-j,.]‘*u[1+j:te,.])*kern[j];
j=j+1;
endo;
delta = sig + lam;
omega = sig +lam+lam’ ;
omega[1,1]=omega[1,1]/te;

retp(omega[1,1]);
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(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+1)*q+i+p)/bigt;
lwz[i+1,1]=ln(glob[i+1,1]/(bigt-(i+1)*q-i-p))
+(((i+1)*q+i+p)*.299*(ln(bigt))^2.1)/bigt;
i=i+1;
endo;
mbic=minindc(bic);
mlwz=minindc(lwz);

retp(mbic-1,mlwz-1);
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(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;