Panel unit root test with structural breaks - calculating break dates

Hello,

I have just purchased GAUSS and am a relative novice at how to use it, so have a question which I hope someone can help with. The code I have doesn't seem to calculate the break dates within this test as it is supposed to. In fact, I get this message:

"the procedure to get critical values for the break dates has reached the upper bound on the number of iterations. This may be due to incorrect specifications of the upper or lower bound in the procedure cvg."

But I am unsure how to correct the upper/lower bound in the code. The below is in a src file - I'm not sure if this is important or not but thought I'd mention it just in case.

If anyone can help I would be extremely grateful.

The code is as follows:

@November 19, 1999 version@
@Copyrigth, Pierre Perron (1998, 1999).@
proc(9)=pbreak(bigt,y,z,q,m,h,eps1,robust,prewhit,hetomega,
        hetq,doglobal,dotest,dospflp1,doorder,dosequa,dorepart,estimbic,estimlwz,
        estimseq,estimrep,estimfix,fixb,x,q,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[1,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;
                else;
                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;
                else;
                endif;
            endif;
        endo;
        
    endif;
    if estimfix == 1;
        
        print "********************************************************";
        
        Print "Output from the estimation of the model with" fixb "breaks";
        print "--------------------------------------------------------------";
        
        call estim(fixn,q,z,y,datevec[.,fixb],robust,prewhit,
            hetomega,hetq,x,p);
        
    endif;
    
    retp(datevec,nbreak,mbic,mlwz,supfl,dateseq,ftest,wftest,reparv);
    
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])-1;
    bound[.,2]=int(bound[.,2])+1;
    bound[.,3]=int(bound[.,3])-1;
    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 g;
    
    if x <= 0;
        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;
        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);
    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=200.0;
        lwb=-200.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;

Your Answer


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.

Try GAUSS for 14 days for FREE

See what GAUSS can do for your data

© Aptech Systems, Inc. All rights reserved.

Privacy Policy