hi, i downloaded some codes written by others and changed some parameters, however it did not work. could you help me find the errors.
load data[61,3]=market.txt;
@Set your T by n matrix of dependent variables@
y=data[2:61,.];
@Set your T by q matrix of regressors@
z=ones(60,1)~data[1:60,2];
@Specify the number of observations@
T=rows(y);
@Specify the matrix _S (here the first and second regressors are used in the
first equation and the first and third regressors are used in the second equation)@
_S=eye(2);
@set the maximum number of breaks allowed. If use the WDmax test set m=5, since the critical values reported
correspond to this case@
m=5;
@If you have restrictions, specify the matrix R (here only the constants (first
regressor) is allowed to change in each equation).
if you do not have restrictions, set R=eye(cols(_S)*(m+1))@
R=eye(cols(_S)*(m+1));
@set the trimming parameter that specifies the minimal length of a segment as a proportion of the sample size.@
trm=0.2;
@set brv=1 when allowing breaks in the covariance matrix of the errors; otherwise, set brv=0@
brv=0;
@set brbeta=1 when allowing breaks in regression coefficients; otherwise, set brbeta=0.@
brbeta=1;
@set vauto=1 when applying a correction for serial correlation in the errors
(in which case Andrews' (1991) method is used to construct the robust
covariance matrix.@
vauto=0;
@set to prewhit=1 if you want to apply an AR(1) pre-whitening when estimating the
robust covariance matrix (see Andrews and Monahan (1992).@
prewhit=0;
@Set hetq=1 if you want to allow the distribution of the regressors to change across
regimes (this is used only when constructing confidence intervals for the estimates
of the break dates. For the construction of the test hetq=1 always.@
hetq=0;
/*Call the main procedure*/
call mainp(m,cols(_S),z,y,cols(y),trm,T,brv,brbeta,vauto,hetq,prewhit);i didn't change the main procedure. it looks like following:
/*Main procedure*/
proc(0)=mainp(m,sumq,z,y,n,trm,bigT,brv,brbeta,vauto,hetq,prewhit);
@
***Input: m: number of breaks allowed
sumq: number of regression coefficients in each regime
z: regressors as a matrix, with the jth row containing
observations at time j
y, dependent variable as a matrix, with the jth row containing
observations at time j
n: the dimension of the dependent variable
trm: trimming
bigT: sample size
brv: equal to 1 if covariance matrix is allowed to change
brbeta: equal to 1 if regression coefficients is allowed to change
vauto: equal to 1 if the error is allowed to be autocorrelated
hetq: equal to 1 if distribution of the regressors is allowed to change
prewhit: equal to 1 if use prewhittening
***Global: _S, matrix of constants selecting variables appearing in each
regression.
R: matrix of constant imposing restrictions across regimes
@
local i,lrtest,maty,matz,global,datevec,bigvec,h,beta,vv,bound,jk,rbeta,rvv,dx,bound2;
local dxi,rbetai,rlrtest,cvs,clevel,wd,seq,cq,rftest;
print "------------------------------------------------------------------------";
print "The basic specifications for testing and estimation:";
print "-------------------------------------------";
print "(1) M=" m;
print "(2) Trimming=" trm;
print "(3) T=" bigt;
if brv==1;
print "(4) The ovariance matrix of errors is allowed to change";
print " Normality is assumed when testing changes in the covariance matrix";
else;
print "(4) The covariance matrix of errors is not allowed to change";
endif;
print "(5) The number of coefficients (beta) in each regime is: " sumq;
print "-------------------------------------------";
print "other specificiations:";
print "-------------------------------------------";
print "(1) the code reports results for an unrestricted model (apart";
print " from the basic specifications and for a restricted model";
print " with restrictions specified by the user)";
if vauto==1;
print "(2) the error is allowed to be autocorrelated";
else;
print "(2) the error is serially uncorrelated";
endif;
if hetq==1;
print "(3) the distribution of the regressors is allowed to change";
else;
print "(3) the distributions of the regressors is NOT allowed to change";
endif;
if prewhit==1;
print "(4) Pre-whittening with VAR(1) when constructing confidence intervals";
else;
print "(4) No pre-whittening when constructing confidence intervals";
endif;
@print "***********************************************************************";@
print "------------------------------------------------------------------------";
h=round(trm*bigt); @trimming@
@transform the data@
{maty,matz}=transf(y,z);
print " ";
print "***********************************************************************";
print "* OUTPUT FROM THE UNRESTRICTED MODEL: TESTING AND ESTIMATION *";
print "***********************************************************************";
@
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 "Output from the testing procedures:";
print "------------------------------------------------------------------------";
lrtest=zeros(m,1);
cvs=zeros(m,4);
clevel=10~5~2.5~1; @significance level@
/* The variance-covariance matrix of the errors is not allowed to change*/
if brv==0 and brbeta==1;
print "a) SupLR test for structural changes in the regression coefficients";
print "(Number of changing parameters" sumq ")";
print "-----------------";
dxi=zeros(m,m);
if vauto==0; @no serial correlation@
i=1;
do while i <= m;
{dxi[1:i,i],rbetai}=est(maty,matz,n,i,bigt,trm,eye(cols(matz)*(i+1)),1,0);
lrtest[i]=suplr(maty,matz,i,n,dxi[1:i,i],rbetai,bigt,1,0);
print "The SupLR test for 0 versus" i "break" lrtest[i,1];
i=i+1;
endo;
print "-----------------";
i=1;
do while i <= m;
cvs[i,.]=cvlr(i,sumq,trm);
i=i+1;
endo;
i=1;
do while i<=4;
print "The critical values at the" clevel[1,i] "% level are (for 1 to" m "breaks)";
print cvs[1:m,i]';
i=i+1;
endo;
@The WD test@
print "-------------------------------------------";
print "b) The WDmax test for up to" m "breaks";
@print "-------------------------------------------";@
WD=zeros(1,4);
i=1;
do while i<=m;
wd=(maxc(WD|(lrtest[i]*(cvs[1,.]./cvs[i,.]))))';
i=i+1;
endo;
cvs=cvwd(sumq,trm);
i=1;
do while i<=4;
print "-----------------";
print "The WDmax test at the" clevel[1,i] "%level is" Wd[i];
print "(The critical value is" cvs[i] ")";
i=i+1;
endo;
print "-------------------------------------------";
print "c) The output from the seq(l+1|l) test";
i=1;
do while i<=m-1;
seq=seqlr(maty,matz,n,dxi[1:i,i],i+1,bigt,h,brbeta,brv);
cvs=supseq(i,sumq,trm);
print "-----------------";
print "The Seq(" i+1 "|" i ") test is : " seq;
print "(The critical values are" cvs ")";
i=i+1;
endo;
elseif vauto==1; @Allowing serial correlation in the errors@
i=1;
do while i <= m;
{dxi[1:i,i],rbetai}=est(maty,matz,n,i,bigt,trm,eye(cols(matz)*(i+1)),1,0);
lrtest[i]=supft(maty,matz,i,n,dxi[1:i,i],rbetai,bigt,brbeta,brv,prewhit);
print "The Supf test for 0 versus" i "break" lrtest[i,1];
i=i+1;
endo;
print "-----------------";
i=1;
do while i <= m;
cvs[i,.]=cvlr(i,sumq,trm);
i=i+1;
endo;
i=1;
do while i<=4;
print "The critical values at the" clevel[1,i] "% level are (for 1 to" m "breaks)";
print cvs[1:m,i]';
i=i+1;
endo;
@The WD test@
print "-------------------------------------------";
print "b) The WDmax test for up to" m "breaks";
@print "-------------------------------------------";@
WD=zeros(1,4);
i=1;
do while i<=m;
wd=(maxc(WD|(lrtest[i]*(cvs[1,.]./cvs[i,.]))))';
i=i+1;
endo;
cvs=cvwd(sumq,trm);
i=1;
do while i<=4;
print "-----------------";
print "The WDmax test at the" clevel[1,i] "%level is" Wd[i];
print "(The critical value is" cvs[i] ")";
i=i+1;
endo;
print "-------------------------------------------";
print "c) The output from the seq(l+1|l) test";
i=1;
do while i<=m-1;
@the sequential test@
seq=seqf(maty,matz,n,dxi[1:i,i],i+1,bigt,h,brbeta,brv,prewhit);
cvs=supseq(i,sumq,trm);
print "-----------------";
print "The Seq(" i+1 "|" i ") test is : " seq;
print "(The critical values are" cvs ")";
i=i+1;
endo;
endif;
print "-------------------------------------------";
print "------------------------------------------------------------------------";
print "Output from the estimation procedures allowing" m "breaks";
{rbeta,rvv}=r_estim(maty,matz,bigt,n,m,dxi[.,m],eye(cols(matz)*(i+1)),brbeta,brv);
bound2=interval2(maty,matz,m,n,bigt,dxi[.,m],rbeta,rvv,hetq,vauto,brv,brbeta,prewhit);
print "--------------------------------------------";
print " The estimated breaks are:" dxi[.,m]';
print "--------------------------------------------";
print "Confidence intervals for the break dates are:";
i=1;
do while i <= m;
print "The 95% C.I. for the" i "th break is: " bound2[i,1] bound2[i,2];
print "The 90% C.I. for the" i "th break is: " bound2[i,3] bound2[i,4];
i=i+1;
endo;
print "--------------------------------------------";
print " The estimated coefficients are:" ;
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rbeta[(jk-1)*(rows(rbeta)/(m+1))+1:jk*(rows(rbeta)/(m+1))]';
jk=jk+1;
endo;
print "---------------------------------------------";
print "The estimated covariance matrices of the errors are";
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rvv[(jk-1)*n+1:jk*n,.];
jk=jk+1;
endo;
print "---------------------------------------------";
endif;
if brv==1 and brbeta==0;
lrtest=zeros(m,1);
print "-------------------------------------------";
print "a) SupLR test for pure structural change in covariance matrix (sigma)";
print "Number of changing parameters" (n+1)*n/2;
print "-----------------";
if vauto==1;
print "NOTE! The code does not allow to test for structural change in the";
print "covariance matrix with serially correlated errors.";
print "The critical values are no longer valid";
print "-----------------";
endif;
i=1;
dxi=zeros(m,m);
do while i <= m;
@brbeta=0, brv=1@
{dxi[1:i,i],rbetai}=est(maty,matz,n,i,bigt,trm,ones((i+1),1).*.eye(cols(matz)),0,1);
lrtest[i]=suplr(maty,matz,i,n,dxi[1:i,i],rbetai,bigt,0,1);
print "The SupLR test for 0 versus" i "break" lrtest[i,1];
i=i+1;
endo;
print "-----------------";
i=1;
do while i <= m;
cvs[i,.]=cvlr(i,(n+1)*n/2,trm);
i=i+1;
endo;
i=1;
do while i<=4;
print "The critical values at the" clevel[1,i] "%level are (for 0 against 1 to " m " breaks)";
print cvs[1:m,i]';
i=i+1;
endo;
print "-------------------------------------------";
print "b) The WDmax test for up to" m "breaks";
@The WD test@
WD=zeros(1,4);
i=1;
do while i<=m;
wd=(maxc(WD|(lrtest[i]*(cvs[1,.]./cvs[i,.]))))';
i=i+1;
endo;
cvs=cvwd((n+1)*n/2,trm);
i=1;
do while i<=4;
print "-----------------";
print "The WDmax test at the" clevel[1,i] "%level is" Wd[i];
print "(The critical value is" cvs[i] ")";
i=i+1;
endo;
print "-------------------------------------------";
print "c) The output from the seq(l+1|l) test";
i=1;
do while i<=m-1;
seq=seqlr(maty,matz,n,dxi[1:i,i],i+1,bigt,h,brbeta,brv);
cvs=supseq(i,(n+1)*n/2,trm);
print "-----------------";
print "The Seq(" i+1 "|" i ") test is : " seq;
print "(The critical values are" cvs ")";
i=i+1;
endo;
print "-------------------------------------------";
print "-------------------------------------------------------------------------";
print "Output from the estimation procedures allowing" m "breaks";
{rbeta,rvv}=r_estim(maty,matz,bigt,n,m,dxi[.,m],ones((i+1),1).*.eye(cols(matz)),brbeta,brv);
bound2=interval2(maty,matz,m,n,bigt,dxi[.,m],rbeta,rvv,hetq,vauto,brv,brbeta,prewhit);
print "--------------------------------------------";
print " The estimated breaks are:" dxi[.,m]';
print "--------------------------------------------";
print "Confidence intervals for the break dates are:";
i=1;
do while i <= m;
print "The 95% C.I. for the" i "th break is: " bound2[i,1] bound2[i,2];
print "The 90% C.I. for the" i "th break is: " bound2[i,3] bound2[i,4];
i=i+1;
endo;
print "--------------------------- -----------------";
print " The estimated coefficients are:" ;
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rbeta[(jk-1)*(rows(rbeta)/(m+1))+1:jk*(rows(rbeta)/(m+1))]';
jk=jk+1;
endo;
print "---------------------------------------------";
print "The estimated covariance matrices of the errors are";
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rvv[(jk-1)*n+1:jk*n,.];
jk=jk+1;
endo;
print "---------------------------------------------";
endif;
if brv==1 and brbeta==1;
lrtest=zeros(m,1);
print "-------------------------------------------";
print "a) SupLR test for pure structural change in regression coefficietns";
print "(beta) and covariance matrix (sigma)";
print "Number changing parameters" sumq+(n+1)*n/2;
print "-----------------";
if vauto==1;
print "NOTE: The code does not allow to test for structural change in the";
print "covariance matrix with serially correlated errors.";
print "The critical values are no longer valid";
print "-----------------";
endif;
dxi=zeros(m,m);
i=1;
do while i <= m;
@brbeta=1, brv=1@
{dxi[1:i,i],rbetai}=est(maty,matz,n,i,bigt,trm,eye(cols(matz)*(i+1)),1,1);
lrtest[i]=suplr(maty,matz,i,n,dxi[1:i,i],rbetai,bigt,1,1);
print "The SupLR test for 0 versus" i "break" lrtest[i,1];
i=i+1;
endo;
print "-----------------";
i=1;
do while i <= m;
cvs[i,.]=cvlr(i,sumq+(n+1)*n/2,trm);
i=i+1;
endo;
i=1;
do while i<=4;
print "The critical values at the" clevel[1,i] "%level are (for 0 against 1 to " m " breaks)";
print cvs[1:m,i]';
i=i+1;
endo;
@The WD test@
print "---------------------------------------------";
print "b) The WDmax test for up to" m "breaks";
WD=zeros(1,4);
i=1;
do while i<=m;
wd=(maxc(WD|(lrtest[i]*(cvs[1,.]./cvs[i,.]))))';
i=i+1;
endo;
cvs=cvwd(sumq+(n+1)*n/2,trm);
i=1;
do while i<=4;
print "-----------------";
print "The WDmax test at the" clevel[1,i] "%level is" Wd[i];
print "(The critical value is" cvs[i] ")";
i=i+1;
endo;
print "---------------------------------------------";
print "c) The output from the seq(l+1|l) test";
i=1;
do while i<=m-1;
seq=seqlr(maty,matz,n,dxi[1:i,i],i+1,bigt,h,brbeta,brv);
cvs=supseq(i,sumq+(n+1)*n/2,trm);
print "-----------------";
print "The Seq(" i+1 "|" i ") test is : " seq;
print "(The critical values are" cvs ")";
i=i+1;
endo;
print "-------------------------------------------";
print "-------------------------------------------------------------------------";
print "Output from the estimation procedures allowing" m "breaks";
{rbeta,rvv}=r_estim(maty,matz,bigt,n,m,dxi[.,m],eye(cols(matz)*(i+1)),brbeta,brv);
bound2=interval2(maty,matz,m,n,bigt,dxi[.,m],rbeta,rvv,hetq,vauto,brv,brbeta,prewhit);
print "--------------------------------------------";
print " The estimated breaks are:" dxi[.,m]';
print "--------------------------------------------";
print "Confidence intervals for the break dates are:";
i=1;
do while i <= m;
print "The 95% C.I. for the" i "th break is: " bound2[i,1] bound2[i,2];
print "The 90% C.I. for the" i "th break is: " bound2[i,3] bound2[i,4];
i=i+1;
endo;
print "--------------------------- -----------------";
print " The estimated coefficients are:" ;
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rbeta[(jk-1)*(rows(rbeta)/(m+1))+1:jk*(rows(rbeta)/(m+1))]';
jk=jk+1;
endo;
print "---------------------------------------------";
print "The estimated covariance matrices of the errors are";
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rvv[(jk-1)*n+1:jk*n,.];
jk=jk+1;
endo;
print "---------------------------------------------";
endif;
print " ";
print "***********************************************************************";
print "* OUTPUT FROM THE RESTRICTED MODEL: TESTING AND ESTIMATION *";
print "***********************************************************************";
print "Output from the testing procedures:";
print "------------------------------------------------------------------------";
{dx,rbeta}=est(maty,matz,n,m,bigt,trm,R,brbeta,brv);
{rbeta,rvv}=r_estim(maty,matz,bigt,n,m,dx,R,brbeta,brv);
if brv==1 and brbeta==0;
if vauto==0;
rlrtest=suplr(maty,matz,m,n,dx,rbeta,bigt,brbeta,brv);
print "The total number coefficients allowed to change:" n*(n+1)/2;
print "The value of the SUPF test:" rlrtest;
print "-----------------";
print "The critical values are:";
print cvlr(i,n*(n+1)/2,trm);
elseif vauto==1;
print "NOTE: The code does not allow to test for structural change in the";
print "covariance matrix of the errors with serially correlated innovations.";
print "No result is reported";
print "-----------------";
endif;
endif;
if brv==0 and brbeta==1;
cq=(n+1)*n/2*(brv==1)+rows(R)/(m+1)-(rows(R)-cols(r))/(m);
print "The total number coefficients allowed to change:" cq;
if vauto==0;
@no autocorrelation, the SUPLR test is used@
rlrtest=suplr(maty,matz,m,n,dx,rbeta,bigt,brbeta,brv);
print "The value of the SUPLR test:" rlrtest;
print "-----------------";
elseif vauto==1;
@use supf test with restrictions@
rftest=suprft(maty,matz,m,n,dx,rbeta,R,bigt,brbeta,brv,prewhit);
print "The value of the SUPF test:" rftest;
print "-----------------";
endif;
print "The critical values are:";
print cvlr(i,cq,trm);
endif;
if brv==1 and brbeta==1;
if vauto==0;
cq=(n+1)*n/2*(brv==1)+rows(R)/(m+1)-(rows(R)-cols(r))/(m);
print "The total number coefficients allowed to change:" cq;
@no autocorrelation, the SUPLR test is used@
rlrtest=suplr(maty,matz,m,n,dx,rbeta,bigt,brbeta,brv);
print "The value of the SUPLR test:" rlrtest;
print "-----------------";
print "The critical values are:";
print cvlr(i,cq,trm);
elseif vauto==1;
print "NOTE: The code does not allow to test for structural change in the";
print "covariance matrix of the errors with serially correlated innovations.";
print "Instead the test for change in regression coefficients is reported";
print "-----------------";
cq=rows(R)/(m+1)-(rows(R)-cols(r))/(m);
print "the total number of regression coefficients allowed to change:" cq;
@use supf test with restrictions@
rftest=suprft(maty,matz,m,n,dx,rbeta,R,bigt,brbeta,brv,prewhit);
print "The value of the SUPF test:" rftest;
print "-----------------";
print "The critical values are:";
print cvlr(i,cq,trm);
endif;
endif;
print "------------------------------------------------------------------------";
print "Output from the estimation procedures allowing" m "breaks";
print " (The restricted model) ";
print "--------------------------------------------";
print " The estimated breaks are:" dx';
print "--------------------------------------------";
bound2=interval2(maty,matz,m,n,bigt,dx,rbeta,rvv,hetq,vauto,brv,brbeta,prewhit);
print "Confidence intervals for the break dates are:";
i=1;
do while i <= m;
print "The 95% C.I. for the" i "th break is: " bound2[i,1] bound2[i,2];
print "The 90% C.I. for the" i "th break is: " bound2[i,3] bound2[i,4];
i=i+1;
endo;
print "--------------------------- -----------------";
print " The estimated coefficients are:" ;
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rbeta[(jk-1)*(rows(rbeta)/(m+1))+1:jk*(rows(rbeta)/(m+1))]';
jk=jk+1;
endo;
print "---------------------------------------------";
print "The estimated covariance matrices of the errors are";
jk=1;
do while jk<=m+1;
print " for the" jk "th regime: " rvv[(jk-1)*n+1:jk*n,.];
jk=jk+1;
endo;
print "---------------------------------------------";
print "***********************************************************************";
print "* THE END! *";
print "***********************************************************************";
endp;
@*************************************************************************@
proc(2)=transf(y,z);
@
***Procedure that generates the matrix of regressors (matz) and that of
the dependent variable (maty)
***Input: y, dependent variable as a matrix, with the jth row containing
observations at time j
z: regressors as a matrix, with the jth row containing
observations at time j
***Output: maty, dependent variable as a row vector, with the (i-1)*n+jth
row containing ith observation for equation j
matz: regressors as a matrix
***Global: _S, matrix of constants selecting variables appearing in each
regression.
@
local maty,matz,i,_ID;
_ID=eye(cols(y));
maty=y[1,.]';
matz=_ID.*.z[1,.]*_S;
i=2;
do while i<=rows(y);
maty=maty|y[i,.]';
matz=matz|_ID.*.z[i,.]*_S;
i=i+1;
endo;
retp(maty,matz);
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 minus likelhood value (MML).
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);
@ for a given ending date, the index is the inserting break point ,
the cell contains the minus maximized likelhood function (MML)@
ini=(start-1)*bigt-(start-2)*(start-1)/2+1;
@ the starting date of the data@
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];
@ the first term: the MML corresponding to data from date start to j
The second term: the ML corresponding to data from (j+1) and lasts
for last -j
@
j=j+1;
endo;
ssrmin=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(ssrmin,dx);
endp;
@*************************************************************************@
/* Find the optimal partition allowing m breaks, unrestricted model*/
proc(3)=dating_MLE(maty,matz,n,h,m,bigt);
@ Input:
maty: vector dependent variables, with the (i-1)*n+j th element
corresponding to the ith obseration from the jth equation
included in the regression
matz: matrix of independent variables, arranged in the same manner
as maty. The number of columns corresponds to the number of
basic parameters included in the regression
h: mimunum number of observations in each segment
m: number of breaks.
bigt: number of obervations in the sample
Output:
global: The maximized likelhood function under optimal partitons.
datevec: the estimated break dates.
bigvec:the big vector containing the MML for all possible segments
@
local datevec,optdat,optmle,dvec,i,mlemax,datx,j1,ib,jlast;
local global,vecmle,jb,xx,bigvec,zi, yi, k, temp,b,beta;
datevec=zeros(m,m);
@up-tiangular matrix contains the estimated break dates for break
numbers from one to m
@
optdat=zeros(bigt,m);
@row index corresponds to the ending dates, column index corresponds
to the number of breaks permitted before the ending date the cell
contains the optimal final break date
@
optmle=zeros(bigt,m);
@same as above, the cell contains the MML corresponding to that break
sturcture
@
dvec=zeros(bigt,1);
@the index is the date after which we inserting the break point.The
cell contains the corresponding MML
@
global=zeros(m,1);
@Global MML when i breaks are permitted
@
bigvec=zeros(bigt*(bigt+1)/2,1);
@the vector that contains the MML corresponding to all the possible segments.
The MML corresponding to the segment starting at J and lasting for span
corresponds to the index T(j-1)-(j-1)(j-2)/2+span
@
i=1;
do while i <= bigt-h+1;
{vecmle}=mlef(i,maty,matz,h,n,bigt);
bigvec[(i-1)*bigt+i-(i-1)*i/2:i*bigt-(i-1)*i/2,1]=vecmle[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;
{mlemax,datx}=parti(1,h,bigt-h,bigt,bigvec,bigt);
datevec[1,1]=datx;
global[1,1]=mlemax;
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 MML are stored in a vector optssr.
@
j1=2*h;
@First loop. Looking for the optimal one break partitions for break
dates between h and T-h. j1 is the last date of the segments.
@
do while j1 <= bigt;
{mlemax,datx}=parti(1,h,j1-h,j1,bigvec,bigt);
optmle[j1,1]=mlemax; @ again no typo@
optdat[j1,1]=datx;
j1=j1+1;
endo;
global[1,1]=optmle[bigt,1]; @ no typo@
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] = optmle[jb,ib-1]+bigvec[(jb+1)*bigt-jb*(jb+1)/2,1];
jb=jb+1;
endo;
optmle[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] = optmle[jb,ib-1]+bigvec[jb*bigt-jb*(jb-1)/2+jlast-jb,1];
jb=jb+1;
endo;
optmle[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]=optmle[bigt,ib];
ib=ib+1;
endo;
endif; @closing the if for the case m >1@
retp(-global,datevec,bigvec);
endp;
@************************************************************************@
/*Computing MLE based on a recursive formula*/
proc(1)=mlef(start,maty,matz,h,n,last);
local vecmle,v,r,i,j,b,f,g,retcode,t,x1, y1,icstar,gstar,ihstar,bigvar,bstar,
k,tempy,tempz,ystar,zstar,vstar,itr,vvar,bbstar,res,umat,ibigv,segx,segy;
@
***Function: This procedure computes MLE from a data set that
starts at date "start" and ends at date "last".
***return: It returns a vector of MLE of length last-start+1
(stored for convenience in a vector of length T).
***method: iterated FGLS is used to obtain the estimates, the recursive formula
in Tobing and McGilchrist (1992) is applied to reduce
the computational cost.
***Inputs:
start: starting entry of the sample used.
last: ending date of the last segment considered
maty: vector of dependent variables, with the n*(j-1)+i th element
corresponding to the jth observation from ith equation.
matz: matrix of all regressors included in the regression, arranged
in the same way as maty
h: minimal length of a segment
n: number of equations included in the system.
***Note: that the first h-1 elements have no meaning, usually is
initialized as zero
***Note: the maximum number of iterations allowed is 1000
@
vecmle=zeros(last,1);
@ initialize the beta and sigma using the first h observations,
iterated until convergence
@
i=start;
j=start+h-1;
segx=matz[(i-1)*n+1:j*n,.];
segy=maty[(i-1)*n+1:j*n,.];
b=invpd(segx'segx)*segx'segy;
res=segy-segx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/h;
vstar=vvar+1;
itr=1;
do while maxc(abs(vecr(vvar-vstar)))>1e-6 and itr<1000;
vstar=vvar;
ibigv=eye(rows(segy)/n).*.invpd(vstar);
b=invpd(segx'ibigv*segx)*segx'ibigv*segy;
res=segy-segx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/h;
itr=itr+1;
if itr==1000;
print " warning: The number of iterations has reached the limit, FGLS does not converge";
print " consider to increase the minimum length of a segment";
endif;
endo;
bstar=b;
ibigv=eye(h).*.invpd(vvar);
ihstar=invpd(segx'*ibigv*segx);
bstar=b;
ystar=segy;
zstar=segx;
vstar=vvar;
@ start recursion@
do while j<=last;
gstar=ihstar*matz[(j-1)*n+1:j*n,.]';
icstar=invpd(vstar+matz[(j-1)*n+1:j*n,.]*ihstar*matz[(j-1)*n+1:j*n,.]');
b=bstar+gstar*icstar*(maty[(j-1)*n+1:j*n,.]-matz[(j-1)*n+1:j*n,.]*bstar);
tempz=zstar|matz[(j-1)*n+1:j*n,.];
tempy=ystar|maty[(j-1)*n+1:j*n,.];
res=tempy-tempz*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/(j-i+1);
itr=1;
vstar=vvar+10;
do while maxc(abs(vecr(vvar-vstar)))>1e-6 and itr<1000;
vstar=vvar;
@ note that gstar does not need to be updated@
icstar=invpd(vstar+matz[(j-1)*n+1:j*n,.]*ihstar*matz[(j-1)*n+1:j*n,.]');
b=bstar+gstar*icstar*(maty[(j-1)*n+1:j*n,.]-matz[(j-1)*n+1:j*n,.]*bstar);
res=tempy-tempz*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/(j-i+1);
itr=itr+1;
if itr==1000;
print " warning: The number of iterations has reached the limit, FGLS does not converge";
print " consider to increase the minimum length of a segment";
endif;
endo;
vecmle[j,1]= ((j-i+1)*n/2)*(ln(2*pi)+1)+(j-i+1)/2*ln(det(vvar));
@ note that minus the likelihood function is recorded@
bstar=b;
zstar=tempz;
ystar=tempy;
vstar=vvar;
ihstar=ihstar-gstar*icstar*gstar';
j=j+1;
endo;
retp(vecmle);
endp;
@***************************************************************************@
/* Procedure that computes confidence intervals for the break dates based on
the "shrinking shifts asymptotic framework." */
/*
proc(1)= interval(maty,matx,m,n,br,beta,vv);
local a,bound,i,j,A1,A2,Q1,Q2,Pi1,Pi2,vvar1,vvar2,cvf,Omi;
local phi1s,phi2s,eta,omega1,hetq,Delt1,Delt2,Gam1,Gam2;
@now constructing the quantities needed for the limit distribution@
bound=zeros(m,4);
br=br|(rows(maty)/n);
j=1;
do while j<=m;
vvar1=vv[(j-1)*n+1:j*n,.];
vvar2=vv[(j)*n+1:(j+1)*n,.];
A1=sqrm(vvar1)*invpd(vvar2)*(vvar2-vvar1)*invpd(sqrm(vvar1));
A2=sqrm(vvar2)*invpd(vvar1)*(vvar2-vvar1)*invpd(sqrm(vvar2));
Q1=0;
Q2=0;
Pi1=0;
Pi2=0;
if j==1;
i=1;
else;
i=br[j-1]+1;
endif;
do while i<=br[j];
Q1=Q1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
Pi1=Pi1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*vvar1*invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
if j==1;
Q1=Q1/br[j];
Pi1=Pi1/br[j];
else;
Q1=Q1/(br[j]-br[j-1]);
Pi1=Pi1/(br[j]-br[j-1]);
endif;
i=br[j]+1;
do while i<=br[j+1];
Q2=Q2+matx[(i-1)*n+1:i*n,.]'invpd(vvar1)*matx[(i-1)*n+1:i*n,.];
Pi2=Pi2+matx[(i-1)*n+1:(i)*n,.]'invpd(vvar1)*vvar2*invpd(vvar1)*matx[(i-1)*n+1:(i)*n,.];
i=i+1;
endo;
Q2=Q2/(br[j+1]-br[j]);
Pi2=Pi2/(br[j+1]-br[j]);
Delt1=(1/2)*sumc(diag(A1*A1))+(beta[j+1,.]-beta[j,.])*Q1*(beta[j+1,.]-beta[j,.])';
Delt2=(1/2)*sumc(diag(A2*A2))+(beta[j+1,.]-beta[j,.])*Q2*(beta[j+1,.]-beta[j,.])';
Gam1=sqrt((1/2)*vec(A1)'vec(A1)+(beta[j+1,.]-beta[j,.])*Pi1*(beta[j+1,.]-beta[j,.])');
Gam2=sqrt((1/2)*vec(A2)'vec(A2)+(beta[j+1,.]-beta[j,.])*Pi2*(beta[j+1,.]-beta[j,.])');
eta=Delt2/Delt1;
phi1s=Gam1*Gam1/Delt1;
phi2s= Gam2*Gam2/Delt2;
cvf=cvg(eta,phi1s,phi2s);
a=(Delt1/Gam1)^2;
bound[j,1]=br[j,1]-cvf[4,1]/a;
bound[j,2]=br[j,1]-cvf[1,1]/a;
bound[j,3]=br[j,1]-cvf[3,1]/a;
bound[j,4]=br[j,1]-cvf[2,1]/a;
j=j+1;
endo;
bound[.,1]=int(bound[.,1]);
bound[.,2]=int(bound[.,2])+1;
bound[.,3]=int(bound[.,3]);
bound[.,4]=int(bound[.,4])+1;
retp(bound);
endp;
*/
@*****************************************************************@
/*obtain the quantiles of the limiting distribution*/
proc(1)= funcg(x,bet,alph,b,deld,gam);
local aa,xb,g;
if x <= 0;
xb=bet*sqrt(abs(x));
if abs(xb) <= 30 ;
g=-sqrt(-x/(2*pi))*exp(x/8)
-(bet/alph)*exp(-alph*x)*cdfn(-bet*sqrt(abs(x)))
+((2*bet*bet/alph)-2-x/2)*cdfn(-sqrt(abs(x))/2);
else;
aa=ln(bet/alph)-alph*x-xb^2/2-ln(sqrt(2*pi))-ln(xb);
g=-sqrt(-x/(2*pi))*exp(x/8)-exp(aa)*cdfn(-sqrt(abs(x))/2)
+((2*bet*bet/alph)-2-x/2)*cdfn(-sqrt(abs(x))/2);
endif;
else;
xb=deld*sqrt(x);
@print xb bet b deld gam*x -deld*sqrt(x) x;@
if abs(xb) <= 30;
g=1+(b/sqrt(2*pi))*sqrt(x)*exp(-b*b*x/8)
+(b*deld/gam)*exp(gam*x)*cdfn(-deld*sqrt(x))
+(2-b*b*x/2-2*deld*deld/gam)*cdfn(-b*sqrt(x)/2);
else;
aa=ln((b*deld/gam))+gam*x-xb^2/2-ln(sqrt(2*pi))-ln(xb);
g=1+(b/sqrt(2*pi))*sqrt(x)*exp(-b*b*x/8)
+exp(aa)
+(2-b*b*x/2-2*deld*deld/gam)*cdfn(-b*sqrt(x)/2);
endif;
endif;
retp(g);
endp;
@****************************************************************@
proc(1)= cvg(eta,phi1s,phi2s);
@ procedure to get the boundaries for the confidence interval@
@ the default is 90% and 95% level@
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=maxc(200|(100/gam));
lwb=-200.0;
crit=999999.0;
cct=1;
do while abs(crit) >= .00001;
cct=cct+1;
if cct > 1000;
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)/20;
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(1)=sqrm(M);
@compute the square root matrix of a symmetric positive definte matrix@
local sq,s,r,va,ve,dd;
{ va,ve } = eigv(M);
dd=diagrv(eye(rows(m)),va);
sq=ve*sqrt(dd)*ve';
retp(sq);
endp;
@************************************************************@
proc(2)=estim(maty,matz,n,m,br);
@procedure that estimates the coefficients under the
optimal partitions@
local vv,beta,b,j,i,k,res,segx,segy,vvar,umat,itr,vstar,ibigv,rep;
beta=zeros(m+1,cols(matz));
vv=zeros((m+1)*n,n);
br=0|br|(rows(maty)/n);
k=1;
do while k<=m+1;
i=br[k]+1;
j=br[k+1];
segx=matz[(i-1)*n+1:j*n,.];
segy=maty[(i-1)*n+1:j*n,.];
b=invpd(segx'segx)*segx'segy;
res=segy-segx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/(j-i+1);
vstar=vvar+1;
itr=1;
do while maxc(abs(vecr(vvar-vstar)))>1e-6 and itr<200;
vstar=vvar;
ibigv=eye(rows(segy)/n).*.invpd(vstar);
b=invpd(segx'ibigv*segx)*segx'ibigv*segy;
res=segy-segx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/(j-i+1);
itr=itr+1;
if itr==200;
print " warning: the FGLS fail, the result does not converge";
endif;
endo;
beta[k,.]=b';
vv[(k-1)*n+1:k*n,.]=vvar;
k=k+1;
endo;
retp(beta,vv);
endp;
@***********************************************************************@
proc(1)=suplr(maty,matx,m,n,br,rbeta,bigt,brbeta,brv);
@Procedure that constructs the supf test against m breaks.
The test is the LR test under the optimal partitions
computed by dating_MLE.
Control variables:
if brv==1, against breaks in variance
if brbeta==1, against breaks in coefficients
if brv==1 and brbeta==1, allowing changes in both.
@
local i,j,k,rlr,ulr,beta,vv,seg,ibigv,tempx,tempy,res,umat,bbs,vvar,b,bstar,vstar;
local itr,pmatx;
pmatx=pzbar(matx,m,br,bigt);
res=maty-pmatx*rbeta;
res=reshape(res,rows(res)/n,n);
seg=0|br|(rows(maty)/n);
if brv==1 and brbeta==1;
ulr=0;
k=1;
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
vvar=res[i:j,.]'res[i:j,.]/(j-i+1);
ulr=ulr-((j-i+1)*n/2)*(ln(2*pi)+1)-(j-i+1)/2*ln(det(vvar));
k=k+1;
endo;
endif;
if brv==1 and brbeta/=1;
@now, change in the var only, iterative method is used for the
optimization, with the unrestricted estimate as the prior input.
Then the beta is updated and then var. The procedure is iterated
until convergence. Non-linear optimization algorithm could be used
in the place of the iterative procedure
@
ulr=0;
k=1;
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
vvar=res[i:j,.]'res[i:j,.]/(j-i+1);
ulr=ulr-((j-i+1)*n/2)*(ln(2*pi)+1)-(j-i+1)/2*ln(det(vvar));
k=k+1;
endo;
endif;
if brbeta==1 and brv/=1;
ulr=0;
vvar=res[1:bigt,.]'res[1:bigt,.]/(bigt);
ulr=-(bigT*n/2)*(ln(2*pi)+1)-(bigT)/2*ln(det(vvar));
endif;
@now estimate the restricted model@
b=invpd(matx'matx)*matx'maty;
res=maty-matx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/bigT;
vstar=vvar+1;
itr=1;
do while maxc(abs(vecr(vvar-vstar)))>1e-6 and itr<1000;
vstar=vvar;
ibigv=eye(rows(maty)/n).*.invpd(vstar);
b=invpd(matx'ibigv*matx)*matx'ibigv*maty;
res=maty-matx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/bigT;
itr=itr+1;
if itr==1000;
print " warning: the FGLS fail, the result does not converge";
endif;
endo;
rlr=-(bigT*n/2)*(ln(2*pi)+1)-(bigT)/2*ln(det(vvar));
retp(2*(ulr-rlr));
endp;
@***********************************************************************@
/*Procedrue that estimates the model allowing cross regime restrictions*/
/*For the current version of the code, only linear restriction is allowed*/
proc(2)=r_estim(maty,matz,bigt,n,m,br,R,brbeta,brv);
@
Input:
maty: dependent variable
matz: independent variable
n: number of equations
m: number of breaks
bigt: sample size
br: break dates
R: the restriction, taking the form beta=R*teta, where teta is the vector
of basic parameters of the restricted model
brv: brv=1 if the variance-covariance matrix is allowed to change
brv=0 otherwise
brbeta: brbeta=1 if the coefficients are allowed to change
brbeta=0 otherwise
Output:
nbeta: the estimate of coefficients imposing restrictions
Note: the model includes partial structural change and block-partial structural
change models
The estimation of based on iterative feasible GLS
@
local i,j,k,rlr,ulr,beta,vv,seg,ibigv,tempx,tempy,res,umat,bbs,vvar,b,bstar,vstar;
local itr,pmatz,nbeta,nvv;
seg=0|br|(rows(maty)/n);
@now estimate an unrestricted model@
{beta,vv}=estim(maty,matz,n,m,br);
if brv==1 and brbeta==1;
pmatz=pzbar(matz,m,br,bigt);
pmatz=pmatz*R;
k=1;
ibigv=eye(bigT*n);
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
ibigv[(i-1)*n+1:j*n,(i-1)*n+1:j*n]=eye(j-i+1).*.invpd(vv[(k-1)*n+1:k*n,.]);
k=k+1;
endo;
b=invpd(pmatz'ibigv*pmatz)*pmatz'ibigv*maty;
bstar=b+10;
itr=1;
do while maxc(abs(bstar-b))>=1e-6 and itr<=1000;
bstar=b;
k=1;
@update the covariance matrix@
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
tempx=pmatz[(i-1)*n+1:j*n,.];
tempy=maty[(i-1)*n+1:j*n,.];
res=tempy-tempx*bstar;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/(j-i+1);
vv[(k-1)*n+1:k*n,.]=vvar;
k=k+1;
endo;
k=1;
@create the block diagonal matrix of the covariances@
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
ibigv[(i-1)*n+1:j*n,(i-1)*n+1:j*n]=eye(j-i+1).*.invpd(vv[(k-1)*n+1:k*n,.]);
k=k+1;
endo;
b=invpd(pmatz'ibigv*pmatz)*pmatz'ibigv*maty;
itr=itr+1;
if itr==1000;
print " The iteration has reached the upper limit";
endif;
endo;
nbeta=R*b;
nvv=vv;
endif;
if brv==1 and brbeta/=1;
@now, change in the covariance matrix only, iterative method is used for the
optimization, with the unrestricted estimate as the prior input.
The procedure is iterated until convergence. Non-linear optimization algorithm
could be used in place of the iterative procedure
@
ibigv=eye(bigT*n);
k=1;
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
ibigv[(i-1)*n+1:j*n,(i-1)*n+1:j*n]=eye(j-i+1).*.invpd(vv[(k-1)*n+1:k*n,.]);
k=k+1;
endo;
b=invpd(matz'ibigv*matz)*matz'ibigv*maty;
bstar=b+10;
itr=1;
do while maxc(abs(bstar-b))>=1e-6 and itr<=200;
bstar=b;
@update the variance-covariance matrix@
k=1;
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
tempx=matz[(i-1)*n+1:j*n,.];
tempy=maty[(i-1)*n+1:j*n,.];
res=tempy-tempx*bstar;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/(j-i+1);
vv[(k-1)*n+1:k*n,.]=vvar;
k=k+1;
endo;
k=1;
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
ibigv[(i-1)*n+1:j*n,(i-1)*n+1:j*n]=eye(j-i+1).*.invpd(vv[(k-1)*n+1:k*n,.]);
k=k+1;
endo;
b=invpd(matz'ibigv*matz)*matz'ibigv*maty;
itr=itr+1;
if itr==1000;
print "The iteration has reached the upper limit";
endif;
endo;
nbeta=b;
nvv=vv;
k=2;
do while k<=m+1;
nbeta=nbeta|b;
k=k+1;
endo;
endif;
@ the next allows changes in coefficients only@
if brbeta==1 and brv/=1;
k=1;
vvar=0;
do while k<=m+1;
i=seg[k]+1;
j=seg[k+1];
tempx=matz[(i-1)*n+1:j*n,.];
tempy=maty[(i-1)*n+1:j*n,.];
res=tempy-tempx*beta[k,.]';
umat=reshape(res,rows(res)/n,n);
vvar=vvar+umat'umat;
k=k+1;
endo;
vvar=vvar/bigT;
ibigv=eye(Bigt).*.invpd(vvar);
pmatz=pzbar(matz,m,br,bigt);
pmatz=pmatz*R;
b=invpd(pmatz'ibigv*pmatz)*pmatz'ibigv*maty;
@now iterate till convergence@
bstar=b+10;
itr=1;
do while maxc(abs(vec(b-bstar)))>=1e-6 and itr<=1000;
bstar=b;
k=1;
res=maty-pmatz*bstar;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/bigt;
ibigv=eye(bigt).*.invpd(vvar);
b=invpd(pmatz'ibigv*pmatz)*pmatz'ibigv*maty;
itr=itr+1;
if itr==1000;
print " The iteration has reached the upper limit";
endif;
endo;
nbeta=R*B;
nvv=vvar;
k=2;
do while k<=m+1;
nvv=nvv|vvar;
k=k+1;
endo;
endif;
retp(nbeta,nvv);
endp;
@*****************************************************************************@
/* procedure to construct the diagonal partition of matz with m break at dates b */
proc pzbar(matz,m,bb,bigt);
local nt,q1,matzb,i,n;
nt=rows(matz);
q1=cols(matz);
n=nt/bigt;
if m==0;
print " m=0, no break is allowed, procedure terminates";
retp(0);
else;
matzb=zeros(nt,(m+1)*q1);
matzb[1:(bb[1,1]*n),1:q1]=matz[1:(bb[1,1]*n),.];
i=2;
do while i <= m;
matzb[bb[i-1,1]*n+1:bb[i,1]*n,(i-1)*q1+1:i*q1]=matz[bb[i-1,1]*n+1:bb[i,1]*n,.];
i=i+1;
endo;
matzb[bb[m,1]*n+1:nt,m*q1+1:(m+1)*q1]=matz[bb[m,1]*n+1:nt,.];
retp(matzb);
endif;
endp;
@************************************************************************@
proc(1)=residuals(maty,matz,nbeta,q,n,m);
local i, bigvec2,bigt ;
@ input:
y: dependent variable.
z: regressors
b: estimated coefficients.
q: total number of regression coefficients in each regime.
n: number of equations
m: number of breaks
output:
bigvec2: residuals for all the possible segments
@
bigt=rows(maty)/n;
bigvec2=zeros(bigt*(m+1),n);
i=1;
do while i<=m+1;
bigvec2[(i-1)*bigt+1:i*bigt,.]=reshape((maty-matz*nbeta[(i-1)*q+1:i*q]),bigt,n);
i=i+1;
endo;
retp(bigvec2);
endp;
@************************************************************************@
/* the procedure is called only when number of breaks is one */
@
obtain an optimal one break partitions for a segment that starts at
date start and ends at date last,under the given coefficients.
@
proc(2)=parti2(start,b1,b2,last,bigvec2,bigt,n);
local k,dvec,j,m1,m2,optmle,dx,ini,jj;
@ Input:
start: begining of the segment considered
b1: first possible break date
b2: last possible break date
last: end of segment considered
Output: optmle: minus the value of the log-likelihood function
under the optimal partitons
dx: the break dates
Note: bigvec must be a matrix of size 2*bigt by n
@
dvec=zeros(bigt,1);
j=b1;
do while j<=b2;
m1=bigvec2[start:j,.]'bigvec2[start:j,.]/(j-start+1);
m2=bigvec2[1*bigt+j+1:1*bigt+last,.]'bigvec2[1*bigt+j+1:1*bigt+last,.]/(last-j+1);
dvec[j,1]=((last-start+1)*n/2)*(ln(2*pi)+1)+((j-start+1)/2)*ln(det(m1))+((last-j+1)/2)*ln(det(m2));
j=j+1;
endo;
optmle=minc(dvec[b1:b2,1]);
dx=(b1-1)+minindc(dvec[b1:b2,1]);
retp(optmle,dx);
endp;
@************************************************************************@
/* This procedure is called when multiple breaks are present*/
proc(3)=dating_M2(bigvec2,h,m,n,bigt);
@ Input:
h: mimunum number of observations in each segment
m: number of breaks.
bigt: number of obervations in the sample
Output:
global: The maximized likelhood function under optimal partitons.
datevec: the estimated break dates.
@
local datevec,optdat,optmle,dvec,i,mlemax,datx,j1,ib,jlast;
local global,vecmle,jb,xx,bigvec,zi, yi, k, temp,b,beta,pmle,q;
datevec=zeros(m,m);
@up-tiangular matrix contains the estimated break dates for break
numbers from one to m
@
optdat=zeros(bigt,m);
@row index corresponds to the ending dates, column index corresponds
to the number of breaks permitted before the ending date the cell
contains the optimal final break date
@
optmle=zeros(bigt,m);
@same as above, the cell contains the MML corresponding to that break
sturcture
@
dvec=zeros(bigt,1);
@the index is the date after which we inserting the break point.The
cell contains the corresponding MML
@
global=zeros(m,1);
@Global MML when i breaks are permitted
@
bigvec=zeros(bigt*(bigt+1)/2,1);
@the vector that contains the MML corresponding to all the possible segments.
The MML corresponding to the segment starting at J and lasting for span
corresponds to the index T(j-1)-(j-1)(j-2)/2+span
@
if m == 1;
{mlemax,datx}=parti2(1,h,bigt-h,bigt,bigvec2,bigt,n);
datevec[1,1]=datx;
global[1,1]=mlemax;
else;
j1=2*h;
@First loop. Looking for the optimal one break partitions for break
dates between h and T-h. j1 is the last date of the segments.
@
do while j1 <= bigt;
{mlemax,datx}=parti2(1,h,j1-h,j1,bigvec2[1:2*bigt,.],bigt,n);
optmle[j1,1]=mlemax; @ again no typo@
optdat[j1,1]=datx;
j1=j1+1;
endo;
global[1,1]=optmle[bigt,1]; @ no typo@
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;
pmle=bigvec2[bigt*m+jb+1:bigt*(m+1),.]'bigvec2[bigt*m+jb+1:bigt*(m+1),.]/(bigt-jb);
dvec[jb,1] = optmle[jb,ib-1]+((bigt-jb+1)*n/2)*(ln(2*pi)+1)+((bigt-jb+1)/2)*ln(det(pmle));
jb=jb+1;
endo;
optmle[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;
pmle=bigvec2[bigt*ib+jb+1:bigt*ib+jlast,.]'bigvec2[bigt*ib+jb+1:bigt*ib+jlast,.]/(jlast-jb);
dvec[jb,1] = optmle[jb,ib-1]+((jlast-jb+1)*n/2)*(ln(2*pi)+1)+((jlast-jb+1)/2)*ln(det(pmle));
jb=jb+1;
endo;
optmle[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]=optmle[bigt,ib];
ib=ib+1;
endo;
endif; @closing the if for the case m >1@
retp(-global,datevec,bigvec);
endp;
@********************************************************************@
/* This procedure estimates the breaks and the coefficients
using the restrictions, based on the iterative method*/
proc(2)=est(maty,matz,n,m,bigt,trm,R,brbeta,brv);
@ Input:
y: dependent variable
z: independent variable
q: total number of regression coefficients in a single regime
bigt: sample size
Output:
dx: the break dates
delta: the coefficients
@
local h,global,datevec,bigvec,br,tbar,i,t1,t2,zi,delta,bigvec2,tstar,diff,dx,zbar;
local rbeta, rvv,q;
q=cols(matz);
h=round(trm*bigt);
tstar=bigt;
{global,datevec,bigvec}=dating_mle(maty,matz,n,h,m,bigt);
@find the optimal partiton of the sample without restrictions
on the coefficients, serves as the starting value of the
iterative procedure@
br=datevec[.,m]; @break dates@
{rbeta,rvv}=r_estim(maty,matz,bigt,n,m,datevec[.,m],R,brbeta,brv);
@now find optimal partitions under the estimated coefficients@
diff=0;
do while diff/=-1;
bigvec2=residuals(maty,matz,rbeta,q,n,m);
{global,datevec,bigvec}=dating_M2(bigvec2,h,m,n,bigt);
if datevec[.,m]==br;
diff=-1;
dx=br;
else;
br=datevec[.,m];
{rbeta,rvv}=r_estim(maty,matz,bigt,n,m,br,R,brbeta,brv);
diff=diff+1;
endif;
endo;
retp(dx,rbeta);
endp;
@********************************************************************@
/* Procedure that computes confidence intervals for the break dates based on
the "shrinking shifts asymptotic framework. */
proc(1)= interval2(maty,matx,m,n,bigt,br,beta,vv,hetq,vauto,brv,brbeta,prewhit);
local a,bound,i,j,A1,A2,Q1,Q2,Pi1,Pi2,vvar1,vvar2,cvf,Omiga1,omiga2,diagx;
local phi1s,phi2s,eta,Delt1,Delt2,Gam1,Gam2,res,k,tempmat,to,diffb;
@now constructing the quantities needed for the limit distribution@
bound=zeros(m,4);
diagx=pzbar(matx,m,br,bigt);
br=br|(rows(maty)/n);
res=maty-diagx*beta;
res=reshape(res,rows(res)/n,n);
@get standarized residuals@
j=1;
do while j<=m+1;
if j==1;
res[1:br[j],.]=res[1:br[j],.]*invpd(sqrm(vv[1:j*n,.]));
else;
res[br[j-1]+1:br[j],.]=res[br[j-1]+1:br[j],.]*invpd(sqrm(vv[(j-1)*n+1:j*n,.]));
endif;
j=j+1;
endo;
@xy(seqa(1,1,rows(res)),res);@
j=1;
do while j<=m;
@construct confidence inteval for each break@
vvar1=vv[(j-1)*n+1:j*n,.];
vvar2=vv[(j)*n+1:(j+1)*n,.];
A1=sqrm(vvar1)*invpd(vvar2)*(vvar2-vvar1)*invpd(sqrm(vvar1));
A2=sqrm(vvar2)*invpd(vvar1)*(vvar2-vvar1)*invpd(sqrm(vvar2));
Q1=0;
Q2=0;
Pi1=0;
Pi2=0;
if j==1;
i=1;
else;
i=br[j-1]+1;
endif;
if hetq==0;
if (brv==1) and (vauto==1); @ change in variance and/or in the regression coefficients,
allow serial correlation
and change in the distribution of the regressors@
@construct estimate of pi1@
tempmat=zeros(br[j]-i+1,cols(matx));
k=1;
do while k<=br[j]-i+1;
tempmat[k,.]=(matx[(i-1+k-1)*n+1:(i-1+k)*n,.]'invpd(vvar2)*sqrm(vvar1)*(res[i-1+k,.]'))';
k=k+1;
endo;
pi1=correct(tempmat,prewhit); @with scaling@
@construct estimate of omiga@
tempmat=zeros(br[j]-i+1,n^2);
k=1;
do while k<=br[j]-i+1;
tempmat[k,.]=vec(res[i-1+k,.]'res[i-1+k,.]-eye(n))';
k=k+1;
endo;
omiga1=correct(tempmat,prewhit); @with scaling@
do while i<=br[j];
Q1=Q1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
if j==1;
Q1=Q1/br[j];
else;
Q1=Q1/(br[j]-br[j-1]);
endif;
i=br[j]+1;
@construct estimate of pi2@
tempmat=zeros(br[j+1]-i+1,cols(matx));
k=1;
do while k<=br[j+1]-i+1;
tempmat[k,.]=(matx[(i-1+k-1)*n+1:(i-1+k)*n,.]'invpd(vvar1)*sqrm(vvar2)*(res[i-1+k,.]'))';
k=k+1;
endo;
pi2=correct(tempmat,prewhit); @with scaling@
tempmat=zeros(br[j+1]-i+1,n^2);
k=1;
do while k<=br[j+1]-i+1;
tempmat[k,.]=vec(res[i-1+k,.]'res[i-1+k,.]-eye(n))';
k=k+1;
endo;
omiga2=correct(tempmat,prewhit); @with scaling@
do while i<=br[j+1];
Q2=Q2+matx[(i-1)*n+1:i*n,.]'invpd(vvar1)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
Q2=Q2/(br[j+1]-br[j]);
diffb=beta[j*cols(matx)+1:(j+1)*cols(matx)]-beta[(j-1)*cols(matx)+1:(j)*cols(matx)];
Delt1=(1/2)*sumc(diag(A1*A1))+diffb'*Q1*diffb;
Delt2=(1/2)*sumc(diag(A2*A2))+diffb'*Q2*diffb;
Gam1=sqrt((1/4)*vec(A1)'Omiga1*vec(A1)+diffb'*Pi1*diffb);
Gam2=sqrt((1/4)*vec(A2)'Omiga2*vec(A2)+diffb'*Pi2*diffb);
endif;
if (brv==1) and (vauto==0); @change in the variance and/or in the regression coefficients@
omiga1=0;
do while i<=br[j];
omiga1=omiga1+vec(res[i,.]'res[i,.]-eye(n))*vec(res[i,.]'res[i,.]-eye(n))';
Q1=Q1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
Pi1=Pi1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*vvar1*invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
if j==1;
omiga1=omiga1/br[j];
Q1=Q1/br[j];
Pi1=Pi1/br[j];
/*print omiga1 "mmiga";
print Q1 "q1";
print Pi1 "pi1";
*/
else;
omiga1=omiga1/(br[j]-br[j-1]);
Q1=Q1/(br[j]-br[j-1]);
Pi1=Pi1/(br[j]-br[j-1]);
endif;
i=br[j]+1;
omiga2=0;
do while i<=br[j+1];
omiga2=omiga2+vec(res[i,.]'res[i,.]-eye(n))*vec(res[i,.]'res[i,.]-eye(n))';
Pi2=Pi2+matx[(i-1)*n+1:(i)*n,.]'invpd(vvar1)*vvar2*invpd(vvar1)*matx[(i-1)*n+1:(i)*n,.];
Q2=Q2+matx[(i-1)*n+1:i*n,.]'invpd(vvar1)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
omiga2=omiga2/(br[j+1]-br[j]);
Q2=Q2/(br[j+1]-br[j]);
Pi2=Pi2/(br[j+1]-br[j]);
/* print omiga2 "mmiga";
print Q2 "q1";
print Pi2 "pi1";
*/
diffb=beta[j*cols(matx)+1:(j+1)*cols(matx)]-beta[(j-1)*cols(matx)+1:(j)*cols(matx)];
Delt1=(1/2)*sumc(diag(A1*A1))+diffb'*Q1*diffb;
Delt2=(1/2)*sumc(diag(A2*A2))+diffb'*Q2*diffb;
Gam1=sqrt((1/4)*vec(A1)'Omiga1*vec(A1)+diffb'*Pi1*diffb);
Gam2=sqrt((1/4)*vec(A2)'Omiga2*vec(A2)+diffb'*Pi2*diffb);
endif;
if (brv==0) and (brbeta==1) and (vauto==1); @change in the coefficients only@
@construct estimate of pi1@
tempmat=zeros(bigt,cols(matx));
k=1;
do while k<=bigt;
tempmat[k,.]=(matx[(k-1)*n+1:(k)*n,.]'invpd(vvar2)*sqrm(vvar1)*(res[k,.]'))';
k=k+1;
endo;
pi1=correct(tempmat,prewhit); @with scaling@
k=1;
do while k<=bigt;
Q1=Q1+matx[(k-1)*n+1:k*n,.]'invpd(vvar2)*matx[(k-1)*n+1:k*n,.];
k=k+1;
endo;
Q1=Q1/bigt;
pi2=pi1;
Q2=Q1;
diffb=beta[j*cols(matx)+1:(j+1)*cols(matx)]-beta[(j-1)*cols(matx)+1:(j)*cols(matx)];
Delt1=diffb'*Q1*diffb;
Delt2=diffb'*Q2*diffb;
Gam1=sqrt(diffb'*Pi1*diffb);
Gam2=sqrt(diffb'*Pi2*diffb);
endif;
if (brv==0) and (brbeta==1) and (vauto==0); @change in the coefficients only@
k=1;
do while k<=bigt;
Q1=Q1+matx[(k-1)*n+1:k*n,.]'invpd(vvar2)*matx[(k-1)*n+1:k*n,.];
Pi1=Pi1+matx[(k-1)*n+1:k*n,.]'invpd(vvar2)*vvar1*invpd(vvar2)*matx[(k-1)*n+1:k*n,.];
k=k+1;
endo;
Q1=Q1/bigt;
Pi1=Pi1/bigt;
pi2=pi1;
Q2=Q1;
diffb=beta[j*cols(matx)+1:(j+1)*cols(matx)]-beta[(j-1)*cols(matx)+1:(j)*cols(matx)];
Delt1=diffb'*Q1*diffb;
Delt2=diffb'*Q2*diffb;
Gam1=sqrt(diffb'*Pi1*diffb);
Gam2=sqrt(diffb'*Pi2*diffb);
endif;
elseif hetq==1; @ allow the distribution of the regressors to be different, we only need to seperately treat vauto=1/0@
@ for a partial break model, the corresponding coefficients are zero@
if vauto==1;
tempmat=zeros(br[j]-i+1,cols(matx));
k=1;
do while k<=br[j]-i+1;
tempmat[k,.]=(matx[(i-1+k-1)*n+1:(i-1+k)*n,.]'invpd(vvar2)*sqrm(vvar1)*(res[i-1+k,.]'))';
k=k+1;
endo;
pi1=correct(tempmat,prewhit); @with scaling@
@construct estimate of omiga@
tempmat=zeros(br[j]-i+1,n^2);
k=1;
do while k<=br[j]-i+1;
tempmat[k,.]=vec(res[i-1+k,.]'res[i-1+k,.]-eye(n))';
k=k+1;
endo;
omiga1=correct(tempmat,prewhit); @with scaling@
do while i<=br[j];
Q1=Q1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
if j==1;
Q1=Q1/br[j];
else;
Q1=Q1/(br[j]-br[j-1]);
endif;
i=br[j]+1;
tempmat=zeros(br[j+1]-i+1,cols(matx));
k=1;
do while k<=br[j+1]-i+1;
tempmat[k,.]=(matx[(i-1+k-1)*n+1:(i-1+k)*n,.]'invpd(vvar1)*sqrm(vvar2)*(res[i-1+k,.]'))';
k=k+1;
endo;
pi2=correct(tempmat,prewhit); @with scaling@
tempmat=zeros(br[j+1]-i+1,n^2);
k=1;
do while k<=br[j+1]-i+1;
tempmat[k,.]=vec(res[i-1+k,.]'res[i-1+k,.]-eye(n))';
k=k+1;
endo;
omiga2=correct(tempmat,prewhit); @with scaling@
do while i<=br[j+1];
Q2=Q2+matx[(i-1)*n+1:i*n,.]'invpd(vvar1)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
Q2=Q2/(br[j+1]-br[j]);
diffb=beta[j*cols(matx)+1:(j+1)*cols(matx)]-beta[(j-1)*cols(matx)+1:(j)*cols(matx)];
Delt1=(1/2)*sumc(diag(A1*A1))+diffb'*Q1*diffb;
Delt2=(1/2)*sumc(diag(A2*A2))+diffb'*Q2*diffb;
Gam1=sqrt((1/4)*vec(A1)'Omiga1*vec(A1)+diffb'*Pi1*diffb);
Gam2=sqrt((1/4)*vec(A2)'Omiga2*vec(A2)+diffb'*Pi2*diffb);
endif;
if vauto==0; @change in the variance and/or in the regression coefficients@
omiga1=0;
to=0;
do while i<=br[j];
omiga1=omiga1+vec(res[i,.]'res[i,.]-eye(n))*vec(res[i,.]'res[i,.]-eye(n))';
Q1=Q1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
Pi1=Pi1+matx[(i-1)*n+1:i*n,.]'invpd(vvar2)*vvar1*invpd(vvar2)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
if j==1;
omiga1=omiga1/br[j];
Q1=Q1/br[j];
Pi1=Pi1/br[j];
else;
omiga1=omiga1/(br[j]-br[j-1]);
Q1=Q1/(br[j]-br[j-1]);
Pi1=Pi1/(br[j]-br[j-1]);
endif;
i=br[j]+1;
omiga2=0;
do while i<=br[j+1];
omiga2=omiga2+vec(res[i,.]'res[i,.]-eye(n))*vec(res[i,.]'res[i,.]-eye(n))';
Pi2=Pi2+matx[(i-1)*n+1:(i)*n,.]'invpd(vvar1)*vvar2*invpd(vvar1)*matx[(i-1)*n+1:(i)*n,.];
Q2=Q2+matx[(i-1)*n+1:i*n,.]'invpd(vvar1)*matx[(i-1)*n+1:i*n,.];
i=i+1;
endo;
omiga2=omiga2/(br[j+1]-br[j]);
Q2=Q2/(br[j+1]-br[j]);
Pi2=Pi2/(br[j+1]-br[j]);
diffb=beta[j*cols(matx)+1:(j+1)*cols(matx)]-beta[(j-1)*cols(matx)+1:(j)*cols(matx)];
Delt1=(1/2)*sumc(diag(A1*A1))+diffb'*Q1*diffb;
Delt2=(1/2)*sumc(diag(A2*A2))+diffb'*Q2*diffb;
Gam1=sqrt((1/4)*vec(A1)'Omiga1*vec(A1)+diffb'*Pi1*diffb);
Gam2=sqrt((1/4)*vec(A2)'Omiga2*vec(A2)+diffb'*Pi2*diffb);
endif;
endif;
eta=Delt2/Delt1;
phi1s=Gam1*Gam1/Delt1;
phi2s= Gam2*Gam2/Delt2;
cvf=cvg(eta,phi1s,phi2s);
a=(Delt1/Gam1)^2;
@print " a" a;@
bound[j,1]=br[j,1]-cvf[4,1]/a;
bound[j,2]=br[j,1]-cvf[1,1]/a;
bound[j,3]=br[j,1]-cvf[3,1]/a;
bound[j,4]=br[j,1]-cvf[2,1]/a;
j=j+1;
endo;
bound[.,1]=int(bound[.,1]);
bound[.,2]=int(bound[.,2])+1;
bound[.,3]=int(bound[.,3]);
bound[.,4]=int(bound[.,4])+1;
retp(bound);
endp;
@********************************************************************@
/* procedure that obtains the robust standard errors*/
proc correct(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(res);
nt=rows(res);
b=zeros(d,1);
bmat=zeros(d,d);
vstar=zeros(nt-1,d);
vmat=res;
@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(1)=cvlr(k,q,trm);
local cv90, cv95,cv975,cv99;
cv90=(7.551+1.718*q-0.041*q*q-0.610*k-15.846*trm+0.025*(q/trm))*exp(0.338/k-0.014/(trm*k));
cv95=(8.238+1.756*q-0.043*q*q-0.659*k-15.436*trm+0.025*(q/trm))*exp(0.389/k-0.013/(trm*k));
cv975=(8.968+1.788*q-0.045*q*q-0.715*k-15.255*trm+0.025*(q/trm))*exp(0.426/k-0.013/(trm*k));
cv99=(9.879+1.771*q-0.042*q*q-0.777*k-14.551*trm+0.025*(q/trm))*exp(0.471/k-0.012/(trm*k));
retp(k*(cv90~cv95~cv975~cv99));
endp;
proc(1)=cvud(q,trm);
local cv90, cv95,cv975,cv99;
cv90=(6.917+2.930*q-9.275*trm)*exp(-0.028*q-0.406*trm);
cv95=(8.228+3.095*q-9.644*trm)*exp(-0.029*q-0.291*trm);
cv975=(9.436+3.304*q-9.301*trm)*exp(-0.030*q-0.259*trm);
cv99=(11.211+3.366*q-7.279*trm)*exp(-0.027*q-0.268*trm);
retp(cv90~cv95~cv975~cv99);
endp;
proc(1)=cvwd(q,trm);
local cv90, cv95,cv975,cv99;
cv90=(7.316+3.128*q-8.624*trm)*exp(-0.029*q-0.412*trm);
cv95=(9.039+3.318*q-9.969*trm)*exp(-0.030*q-0.327*trm);
cv975=(10.703+3.465*q-11.119*trm)*exp(-0.031*q-0.250*trm);
cv99=(13.189+3.346*q-11.870*trm)*exp(-0.026*q-0.173*trm);
retp(cv90~cv95~cv975~cv99);
endp;
proc(1)=supseq(l,q,trm);
local cv90, cv95,cv975,cv99;
cv90=(8.397+3.702*q-0.209*q*q+0.317*(l+1)-3.736/(l+1)-11.596*trm)*exp(-0.027*q+0.005*q*q);
cv95=(9.879+4.086*q-0.232*q*q+0.322*(l+1)-3.687/(l+1)-11.931*trm)*exp(-0.039*q+0.006*q*q);
cv975=(11.424+4.435*q-0.261*q*q+0.300*(l+1)-3.700/(l+1)-12.292*trm)*exp(-0.047*q+0.006*q*q);
cv99=(13.073+4.954*q-0.317*q*q+0.285*(l+1)-3.419/(l+1)-12.452*trm)*exp(-0.058*q+0.008*q*q);
retp(cv90~cv95~cv975~cv99);
endp;
@********************************************************************@
/*sequential likelihood ratio*/
proc(1)=seqlr(maty,matz,n,dt,nseg,bigt,h,brbeta,brv);
local length,news,dv,is,ssr,ds,ssr0,ftestv,maxf,newd,in,starti,endi;
local dxi,rbetai,lrtest,segy,segz;
@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);
lrtest=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;
starti=dv[is,1]+1;
endi=dv[is+1];
segy=maty[n*(starti-1)+1:n*endi];
segz=matz[n*(starti-1)+1:n*endi,.];
if brbeta==1 and brv==0;
{dxi,rbetai}=est(segy,segz,n,1,endi-starti+1,trm,eye(cols(matz)*2),1,0);
lrtest[is]=suplr(segy,segz,1,n,dxi,rbetai,endi-starti+1,1,0);
endif;
if brbeta==0 and brv==1;
{dxi,rbetai}=est(segy,segz,n,1,endi-starti+1,trm,ones(2,1).*.eye(cols(matz)),0,1);
lrtest[is]=suplr(segy,segz,1,n,dxi,rbetai,endi-starti+1,0,1);
endif;
if brbeta==1 and brv==1;
{dxi,rbetai}=est(segy,segz,n,1,endi-starti+1,trm,eye(cols(matz)*2),1,1);
lrtest[is]=suplr(segy,segz,1,n,dxi,rbetai,endi-starti+1,1,1);
endif;
else;
in=in+1;
lrtest[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(lrtest[1:nseg,1]);
retp(maxf);
endp;
@********************************************************************@
@********************************************************************@
/*sequential f test*/
proc(1)=seqf(maty,matz,n,dt,nseg,bigt,h,brbeta,brv,prewhit);
local length,news,dv,is,ssr,ds,ssr0,ftestv,maxf,newd,in,starti,endi;
local dxi,rbetai,lrtest,segy,segz;
@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);
lrtest=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;
starti=dv[is,1]+1;
endi=dv[is+1];
segy=maty[n*(starti-1)+1:n*endi];
segz=matz[n*(starti-1)+1:n*endi,.];
{dxi,rbetai}=est(segy,segz,n,1,endi-starti+1,trm,eye(cols(matz)*2),1,0);
lrtest[is]=supft(segy,segz,1,n,dxi,rbetai,endi-starti+1,1,0,prewhit);
else;
in=in+1;
lrtest[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(lrtest[1:nseg,1]);
retp(maxf);
endp;
@********************************************************************@
/*supF test allowing m breaks*/
proc(1)=supft(maty,matx,m,n,br,rbeta,bigt,brbeta,brv,prewhit);
local pmatx,res,xbstar,ubstar,ostar,sigma,umat,bigsig,tempmat,k,vmat,vbeta;
local rmat,rsub,j,ftest;
pmatx=pzbar(matx,m,br,bigt);
res=maty-pmatx*rbeta;
umat=reshape(res,rows(res)/n,n);
sigma=umat'umat/bigT;
bigsig=invpd(sqrm(eye(bigt).*.sigma));
xbstar=bigsig*pmatx;
ubstar=bigsig*res;
ubstar=reshape(ubstar,rows(ubstar)/n,n);
tempmat=zeros(bigt,cols(xbstar));
k=1;
do while k<=bigt;
tempmat[k,.]=(xbstar[(k-1)*n+1:(k)*n,.]'(ubstar[k,.])')';
k=k+1;
endo;
vmat=correct(tempmat,prewhit);
vbeta=invpd(xbstar'xbstar/bigt)*(vmat/bigt)*invpd(xbstar'xbstar/bigt);
rsub=zeros(m,m+1);
j=1;
do while j <= m;
rsub[j,j]=-1;
rsub[j,j+1]=1;
j=j+1;
endo;
rmat=rsub.*.eye(cols(matx));
ftest=rbeta'rmat'inv(rmat*vbeta*rmat')*rmat*rbeta;
ftest=(bigt-(m+1)*cols(matx))*ftest/(bigt);
retp(ftest);
endp;
@********************************************************************@
/*restricted supF test allowing m breaks*/
proc(1)=suprft(maty,matx,m,n,br,rbeta,R,bigt,brbeta,brv,prewhit);
local pmatx,res,xbstar,ubstar,ostar,sigma,umat,bigsig,tempmat,k,vmat,vbeta,vdelta;
local rmat,rsub,j,ftest;
pmatx=pzbar(matx,m,br,bigt);
res=maty-pmatx*rbeta;
umat=reshape(res,rows(res)/n,n);
sigma=umat'umat/bigT;
bigsig=invpd(sqrm(eye(bigt).*.sigma));
pmatx=pmatx*R;
xbstar=bigsig*pmatx;
ubstar=bigsig*res;
ubstar=reshape(ubstar,rows(ubstar)/n,n);
tempmat=zeros(bigt,cols(xbstar));
k=1;
do while k<=bigt;
tempmat[k,.]=(xbstar[(k-1)*n+1:(k)*n,.]'(ubstar[k,.])')';
k=k+1;
endo;
vmat=correct(tempmat,prewhit);
vdelta=invpd(xbstar'xbstar/bigt)*(vmat/bigt)*invpd(xbstar'xbstar/bigt);
vbeta=R*vdelta*R';
rsub=zeros(m,m+1);
j=1;
do while j <= m;
rsub[j,j]=-1;
rsub[j,j+1]=1;
j=j+1;
endo;
rmat=rsub.*.eye(cols(matx));
ftest=rbeta'rmat'pinv(rmat*vbeta*rmat')*rmat*rbeta;
ftest=(bigt-(m+1)*cols(matx))*ftest/(bigt);
retp(ftest);
endp;
@********************************************************************@
/* confidence interval for the ordered breaks*/
proc(1)=interval3(matz,maty,n,start,last,dx1,dx2,delt1,delt2,res,segT,bigM,rep);
local ibigv,QQ,p1j,p2j,p12j,p1jp1,p2jp1,p12jp1,step,bigN,bigW,interv,grids,gridN;
local W,k,sv1,sv2,V1,V2,Bv1,Bv2,sbv1,sbv2,z1,z2,z3,z4,z,Qj,Qjp1,obd,simu,count,sigma;
local bound;
step=1000;
obd=zeros(rep,3);
sigma=res'res/rows(res); @covariance matrix@
if hetq==1;
ibigv=eye(dx1-start+1).*.invpd(sigma);
Qj=matz[(start-1)*n+1:dx1*n,.]'ibigv*matz[(start-1)*n+1:dx1*n,.]/(dx1-start+1);
ibigv=eye(last-dx2).*.invpd(sigma);
Qjp1=matz[dx2*n+1:last*n,.]'ibigv*matz[dx2*n+1:last*n,.]/(last-dx2);
p1j=delt1'Qj*delt1;
P2j=delt2'Qj*delt2;
P12j=delt1'Qj*delt2;
p1jp1=delt1'Qjp1*delt1;
P2jp1=delt2'Qjp1*delt2;
P12jp1=delt1'Qjp1*delt2;
elseif hetq==0;
ibigv=eye(last-start+1).*.invpd(sigma);
Qj=matz[(start-1)*n+1:last*n,.]'ibigv*matz[(start-1)*n+1:last*n,.]/(last-start);
P1j=delt1'Qj*delt1;
P2j=delt2'Qj*delt2;
p12j=delt1'Qj*delt2;
endif;
@print ";;" p1j p2j p12j p1jp1 p2jp1 p12jp1;@
bigN=step*bigM;
interv=0.02;
grids=1/interv;
gridN=grids*bigM;
simu=1;
do while simu<=rep;
bigW=rndn(bigN,4);
bigW=cumsumc(bigW);
bigW=bigW/sqrt(step);
W=zeros(gridN,4);
k=1;
do while k<=gridN;
W[k,.]=bigW[step*interv*(k-1)+1,.];
k=k+1;
endo;
sv1=seqa(-bigM,1/grids,2*grids*bigM+1);
sv2=sv1;
V1=sv1.*.ones(rows(sv1),1);
V2=ones(rows(sv2),1).*.sv2;
sbv1=rev(W[.,1])|zeros(1,1)|(W[.,4]);
sbv2=rev(W[.,2])|zeros(1,1)|(W[.,3]); @two sided Brownian motion@
Bv1=sbv1.*.ones(rows(sbv1),1);
Bv2=ones(rows(sbv2),1).*.sbv2;
if hetq==1;
Z1=(V1.<0).*(V2.<=0).*(v1.<=v2).*((Bv1+v1/2)+sqrt(p2j/p1j)*Bv2+v2/2*(p2j+2*P12j)/P1j);
Z2=(V1.<0).*(V2.>0).*(Bv1+V1/2+sqrt(p2jp1/p1j)*Bv2-V2/2*(p2jp1/p1j));
Z3=(v1.>=0).*(v2.>0).*(v1.<=v2).*(sqrt(p1jp1/p1j)*Bv1+sqrt(P2jp1/P1j)*Bv2-(P2jp1/P1j)*v2/2-((P1jp1+2*P12jp1)/P1j)*v1/2);
Z4=(v1.>v2)*(-1);
elseif hetq==0;
Z1=(V1.<0).*(V2.<=0).*(v1.<=v2).*((Bv1+v1/2)+sqrt(p2j/p1j)*Bv2+v2/2*(p2j+2*P12j)/P1j);
Z2=(V1.<0).*(V2.>0).*(Bv1+V1/2+sqrt(p2j/p1j)*Bv2-V2/2*(p2j/p1j));
Z3=(v1.>=0).*(v2.>0).*(v1.<=v2).*(Bv1+sqrt(P2j/P1j)*Bv2-(P2j/P1j)*v2/2-((P1j+2*P12j)/P1j)*v1/2);
Z4=(v1.>v2)*(-1);
endif;
Z=Z1+Z2+Z3+Z4;
Z=v1~v2~Z;
z=sortc(z,3);
obd[simu,.]=z[rows(z),.];
simu=simu+1;
endo;
obd[.,1]=sortc(obd[.,1],1);
obd[.,2]=sortc(obd[.,2],1)@-(dx2-dx1)*p1j@;
count=sumc(abs(abs(obd[.,1])-bigM).<1e-4)+sumc(abs(abs(obd[.,2])-bigM)<1e-4);
if count>=rep*0.05;
print "wrong bound";
endif;
bound=zeros(4,2);
bound[1,.]=int((dx1~dx2)-obd[int(rep*0.025),1:2]/p1j)+ones(1,2);
bound[2,.]=int((dx1~dx2)-obd[int(rep*0.05),1:2]/p1j)+ones(1,2);
bound[3,.]=int((dx1~dx2)-obd[int(rep*0.95),1:2]/p1j);
bound[4,.]=int((dx1~dx2)-obd[int(rep*0.975),1:2]/p1j);
retp(bound);
endp;
@********************************************************************@
/* Testing for ordered breaks, also estimating the break dates */
proc(5)=brorder(maty,matz,n,q1,q2,h,T,trun);
local bigvec,ii,jj,bigT,b,e,fit,vvar,vstar,umat,res,segx,segy,vecmle,itr,ibigv,zz,z1,res1,res2,TT;
local vecmle2,H0like,H1like,lr,brdate,beta1,beta2;
TT=T*T;
vecmle=ones(TT,3+2*(q1+q2))*(-1e12);
vecmle2=ones(TT,3+2*(q1+q2))*(-1e12);
if cols(matz)/=q1+q2;
print "wrong input of columns of regressors";
endif;
zz=matz~matz;
z1=zz;
ii=h;
segx=zz;
do while ii<=T-h;
ii=ii+1;
jj=0;
do while (ii+jj)<=T-h;
z1[ii*n+1:rows(zz),cols(matz)+1:cols(matz)+q1]=zz[ii*n+1:rows(zz),cols(matz)+1:cols(matz)+q1];
z1[1:ii*n,cols(matz)+1:cols(matz)+q1]=0.*zz[1:ii*n,cols(matz)+1:cols(matz)+q1];
z1[(ii+jj)*n+1:rows(zz),cols(matz)+q1+1:2*cols(matz)]=zz[(ii+jj)*n+1:rows(zz),cols(matz)+q1+1:2*cols(matz)];
z1[1:(ii+jj)*n,cols(matz)+q1+1:2*cols(matz)]=0.*zz[1:(ii+jj)*n,cols(matz)+q1+1:2*cols(matz)];
segx=z1;
segy=maty;
if jj==0;
b=inv(segx'segx)*segx'segy;
else;
ibigv=eye(rows(segy)/n).*.invpd(vvar);
b=invpd(segx'ibigv*segx)*segx'ibigv*segy;
endif;
res=segy-segx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/rows(umat);
vstar=vvar+1;
itr=1;
do while maxc(abs(vecr(vvar-vstar)))>1e-8 and itr<1000;
vstar=vvar;
ibigv=eye(rows(segy)/n).*.invpd(vvar);
b=invpd(segx'ibigv*segx)*segx'ibigv*segy;
res=segy-segx*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/rows(umat);
itr=itr+1;
if itr==1000;
print " warning: the FGLS fail, the result does not converge";
endif;
endo;
/*now caculate the likelihood function for the given partition*/
vecmle[(ii-1)*T+jj,1]=ii;
vecmle[(ii-1)*T+jj,2]=jj+ii;
vecmle[(ii-1)*T+jj,3]=-T/2*ln(det(vvar));
vecmle[(ii-1)*T+jj,3+1:3+2*(q1+q2)]=b';
jj=jj+1;
endo;
endo;
vecmle2=vecmle;
vecmle2[.,3]=((vecmle2[.,2]-vecmle2[.,1]).<=trun).*vecmle2[.,3]+((vecmle2[.,2]-vecmle2[.,1]).>trun).*(-1e12);
vecmle=sortc(vecmle,3);
vecmle2=sortc(vecmle2,3);
H1like=vecmle2[rows(vecmle2),3];
@estimate a model with no change@
b=invpd(matz'matz)*matz'maty;
res=maty-matz*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/rows(umat);
vstar=vvar+1;
itr=1;
do while maxc(abs(vecr(vvar-vstar)))>1e-8 and itr<1000;
vstar=vvar;
ibigv=eye(rows(maty)/n).*.invpd(vvar);
b=invpd(matz'ibigv*matz)*matz'ibigv*maty;
res=maty-matz*b;
umat=reshape(res,rows(res)/n,n);
vvar=umat'umat/rows(umat);
itr=itr+1;
if itr==1000;
print " warning: the FGLS fail, the result does not converge";
endif;
endo;
H0like=-T/2*ln(det(vvar));
lr=2*(H1like-H0like);
brdate=vecmle[rows(vecmle),1:2];
beta1=vecmle[rows(vecmle),4:3+(q1+q2)];
beta2=vecmle[rows(vecmle),4+(q1+q2):3+2*(q1+q2)]+vecmle[rows(vecmle),4:3+(q1+q2)];
retp(brdate,beta1',beta2',umat,lr);
endp;
proc(0)=mainop(y,z,n,q1,q2,trmm,bigT,hetq,trun,getcon,bigM,rep);
local lrtest,cvs,clevel,brdate,beta1,beta2,res,delt1,delt2,bound,h;
local maty,matz,i;
h=int(trmm*bigt);
print "------------------------------------------------------------------------";
print "The basic specifications for testing and estimation:";
print "-------------------------------------------";
print "(1) Trimming=" trmm;
print "(2) T=" bigt;
print "(3) The break in the first equation occurs first";
if hetq==1;
print "(4) the distribution of the regressors is allowed to change";
else;
print "(4) the distributions of the regressors is NOT allowed to change";
endif;
@transform the data@
{maty,matz}=transf(y,z);
@testing and estimation@
{brdate,beta1,beta2,res,lrtest}=brorder(maty,matz,n,q1,q2,h,bigT,trun);
delt1=(beta2[1:q1]-beta1[1:q1])|zeros(q2,1);
delt2=zeros(q1,1)|(beta2[q1+1:q1+q2]-beta1[q1+1:q1+q2]);
if getcon==1;
bound=interval3(matz,maty,n,1,T,brdate[1],brdate[2],delt1,delt2,res,bigT,bigM,rep);
endif;
print " ";
print "Output from the testing procedures:";
print "------------------------------------------------------------------------";
clevel=10~5~2.5~1; @significance level@
cvs=cvlr(1,q1+q2,h/bigT);
print "The SupLR test for ordered breaks is" lrtest;
print "-----------------";
i=1;
do while i<=4;
print "The critical values at the" clevel[1,i] "% level are";
print cvs[i];
i=i+1;
endo;
print "------------------------------------------------------------------------";
print "Output from the estimation procedure";
print "--------------------------------------------";
print " The estimated breaks are:" brdate;
if getcon==1;
print "--------------------------------------------";
print "Confidence intervals for the break dates are:";
print "The 95% C.I. for the first break is: " bound[4,1] bound[1,1];
print "The 95% C.I. for the second break is: " bound[4,2] bound[1,2];
print "The 90% C.I. for the first break is: " bound[3,1] bound[2,1];
print "The 90% C.I. for the second break is: " bound[3,2] bound[2,2];
endif;
print "-------------------------------------------";
print "The estimated coefficients are:";
print "For the first regime:" beta1';
print "For the second regime:" beta2';
print "-------------------------------------------";
endp;6 Answers
0
The matrix compatibility issue you are facing is arising because of the kronecker product being computed on line 559 of your original code
matz=_ID.*.z[1,.]*_S;
The kronecker product between _ID, a 3x3 matrix, and the first row of the Z matrix, Z[1,.], results in a (3*1) x (3*2) matrix. This 3 x 6 matrix is not compatible for matrix multiplication with _S, a 2 x 2 matrix.
As I am not exactly sure the purpose of this line without further understanding of the model, it is difficult to suggest a suitable solution for your specific model. However, one possible solution would be to replace line 559 with
matz=_ID.*.(z[1,.]*_S);
In this case, the matrix product between Z[1,.], a 1x2 matrix, and _S, a 2x2 matrix, occurs prior to the kronecker multiplication. You will also have to change line 563
matz=matz|_ID.*.z[i,.]*_S;
to
matz=matz|_ID.*.(z[i,.]*_S);.
These revisions should solve the matrix incompatibility.
0
hi, I have figured out where does my error come from, that I should set up_S=eye(6) not 2x2; thanks for helping.
also the original code is from Estimating and testing structural changes in multivariate regressions. (Econometrica, 2007) (developed by Zhongjun Qu). Revised January 2007.
0
when I run the procedure with my new data, it cames with error G0121 : Matrix not positive definite.
is there anything wrong with my data?
| 0.016401343 | 0.004698341 | 0.010611354 |
| 0.010907711 | 0.004499951 | 0.00620937 |
| 0.000593074 | 0.004605865 | -0.003906877 |
| 0.034496269 | 0.003998587 | 0.029890405 |
| 0.019209886 | 0.003422499 | 0.015211299 |
| 0.007585408 | 0.003273966 | 0.004162909 |
| -0.0085012 | 0.003533595 | -0.011775165 |
| 0.030979637 | 0.003212406 | 0.027446042 |
| 0.010152538 | 0.00268135 | 0.006940132 |
| 0.001686743 | 0.002887214 | -0.000994607 |
| 0.020203772 | 0.002439332 | 0.017316559 |
| -0.007452367 | 0.002851368 | -0.009891699 |
| 0.016140456 | 0.002641396 | 0.013289088 |
| 0.012645047 | 0.002504387 | 0.01000365 |
| 0.011499131 | 0.002498013 | 0.008994743 |
| -0.007078074 | 0.002742743 | -0.009576087 |
| 0.021471088 | 0.002491184 | 0.018728345 |
| 0.011765185 | 0.002406283 | 0.009274 |
| -0.00871663 | 0.002539652 | -0.011122912 |
| 0.001888882 | 0.002968714 | -0.00065077 |
| 0.013268911 | 0.002454432 | 0.010300196 |
| 0.013505226 | 0.002236705 | 0.011050794 |
| -0.016217885 | 0.002482635 | -0.01845459 |
| -0.025173032 | 0.003654861 | -0.027655667 |
| 0.027732054 | 0.003427403 | 0.024077193 |
| 0.020482706 | 0.003224572 | 0.017055304 |
| -0.002171272 | 0.003808907 | -0.005395844 |
| 0.0080987 | 0.004143698 | 0.004289793 |
| 0.019021195 | 0.004203931 | 0.014877496 |
| 0.025473933 | 0.003919736 | 0.021270002 |
| -0.002633047 | 0.003887835 | -0.006552783 |
| 0.016860722 | 0.004384871 | 0.012972888 |
| 0.017429634 | 0.003302067 | 0.013044764 |
| 0.00334764 | 0.003561497 | 4.55738E-05 |
| 0.023342596 | 0.003201791 | 0.019781099 |
| 0.013152459 | 0.002785796 | 0.009950669 |
| 0.005416554 | 0.002630605 | 0.002630758 |
| 0.013828938 | 0.003300886 | 0.011198333 |
| 0.021305664 | 0.002905654 | 0.018004778 |
| -0.003925277 | 0.002999645 | -0.006830931 |
| 0.025001833 | 0.002540244 | 0.022002188 |
| 0.007432064 | 0.002280439 | 0.00489182 |
| 0.00929789 | 0.002113675 | 0.007017451 |
| -0.000237278 | 0.002142468 | -0.002350954 |
| 0.025530879 | 0.001968987 | 0.023388411 |
| 0.016445446 | 0.00167556 | 0.014476459 |
| 0.022909153 | 0.001393941 | 0.021233594 |
| 0.018845306 | 0.001203836 | 0.017451365 |
| 0.01953282 | 0.001079305 | 0.018328983 |
| -0.008374468 | 0.000924551 | -0.009453773 |
| -0.008151046 | 0.001110964 | -0.009075597 |
| -0.017965501 | 0.001343473 | -0.019076465 |
| 0.024327993 | 0.001487583 | 0.02298452 |
| 0.010451443 | 0.001609234 | 0.00896386 |
| 0.006186574 | 0.001513426 | 0.00457734 |
| 0.012686271 | 0.001562207 | 0.011172845 |
| 0.006261605 | 0.001580993 | 0.004699398 |
| -0.037836006 | 0.001962901 | -0.039416999 |
| 0.024621228 | 0.002052984 | 0.022658327 |
| 0.014991891 | 0.001765415 | 0.012938907 |
| 0.000224058 | 0.001722681 | -0.001541356 |
| 0.012650642 | 0.002065529 | 0.010927962 |
| 0.022416691 | 0.001777976 | 0.020351162 |
0
Your first and third column are very similar -- which can contribute to your issue of "Matrix Not Positive Definitive". For more information on diagnosing these errors, please see our blog "Diagnosing a singular matrix."
Your Answer
6 Answers
0
The matrix compatibility issue you are facing is arising because of the kronecker product being computed on line 559 of your original code
matz=_ID.*.z[1,.]*_S;
The kronecker product between _ID, a 3x3 matrix, and the first row of the Z matrix, Z[1,.], results in a (3*1) x (3*2) matrix. This 3 x 6 matrix is not compatible for matrix multiplication with _S, a 2 x 2 matrix.
As I am not exactly sure the purpose of this line without further understanding of the model, it is difficult to suggest a suitable solution for your specific model. However, one possible solution would be to replace line 559 with
matz=_ID.*.(z[1,.]*_S);
In this case, the matrix product between Z[1,.], a 1x2 matrix, and _S, a 2x2 matrix, occurs prior to the kronecker multiplication. You will also have to change line 563
matz=matz|_ID.*.z[i,.]*_S;
to
matz=matz|_ID.*.(z[i,.]*_S);.
These revisions should solve the matrix incompatibility.
0
hi, I have figured out where does my error come from, that I should set up_S=eye(6) not 2x2; thanks for helping.
also the original code is from Estimating and testing structural changes in multivariate regressions. (Econometrica, 2007) (developed by Zhongjun Qu). Revised January 2007.
0
when I run the procedure with my new data, it cames with error G0121 : Matrix not positive definite.
is there anything wrong with my data?
| 0.016401343 | 0.004698341 | 0.010611354 |
| 0.010907711 | 0.004499951 | 0.00620937 |
| 0.000593074 | 0.004605865 | -0.003906877 |
| 0.034496269 | 0.003998587 | 0.029890405 |
| 0.019209886 | 0.003422499 | 0.015211299 |
| 0.007585408 | 0.003273966 | 0.004162909 |
| -0.0085012 | 0.003533595 | -0.011775165 |
| 0.030979637 | 0.003212406 | 0.027446042 |
| 0.010152538 | 0.00268135 | 0.006940132 |
| 0.001686743 | 0.002887214 | -0.000994607 |
| 0.020203772 | 0.002439332 | 0.017316559 |
| -0.007452367 | 0.002851368 | -0.009891699 |
| 0.016140456 | 0.002641396 | 0.013289088 |
| 0.012645047 | 0.002504387 | 0.01000365 |
| 0.011499131 | 0.002498013 | 0.008994743 |
| -0.007078074 | 0.002742743 | -0.009576087 |
| 0.021471088 | 0.002491184 | 0.018728345 |
| 0.011765185 | 0.002406283 | 0.009274 |
| -0.00871663 | 0.002539652 | -0.011122912 |
| 0.001888882 | 0.002968714 | -0.00065077 |
| 0.013268911 | 0.002454432 | 0.010300196 |
| 0.013505226 | 0.002236705 | 0.011050794 |
| -0.016217885 | 0.002482635 | -0.01845459 |
| -0.025173032 | 0.003654861 | -0.027655667 |
| 0.027732054 | 0.003427403 | 0.024077193 |
| 0.020482706 | 0.003224572 | 0.017055304 |
| -0.002171272 | 0.003808907 | -0.005395844 |
| 0.0080987 | 0.004143698 | 0.004289793 |
| 0.019021195 | 0.004203931 | 0.014877496 |
| 0.025473933 | 0.003919736 | 0.021270002 |
| -0.002633047 | 0.003887835 | -0.006552783 |
| 0.016860722 | 0.004384871 | 0.012972888 |
| 0.017429634 | 0.003302067 | 0.013044764 |
| 0.00334764 | 0.003561497 | 4.55738E-05 |
| 0.023342596 | 0.003201791 | 0.019781099 |
| 0.013152459 | 0.002785796 | 0.009950669 |
| 0.005416554 | 0.002630605 | 0.002630758 |
| 0.013828938 | 0.003300886 | 0.011198333 |
| 0.021305664 | 0.002905654 | 0.018004778 |
| -0.003925277 | 0.002999645 | -0.006830931 |
| 0.025001833 | 0.002540244 | 0.022002188 |
| 0.007432064 | 0.002280439 | 0.00489182 |
| 0.00929789 | 0.002113675 | 0.007017451 |
| -0.000237278 | 0.002142468 | -0.002350954 |
| 0.025530879 | 0.001968987 | 0.023388411 |
| 0.016445446 | 0.00167556 | 0.014476459 |
| 0.022909153 | 0.001393941 | 0.021233594 |
| 0.018845306 | 0.001203836 | 0.017451365 |
| 0.01953282 | 0.001079305 | 0.018328983 |
| -0.008374468 | 0.000924551 | -0.009453773 |
| -0.008151046 | 0.001110964 | -0.009075597 |
| -0.017965501 | 0.001343473 | -0.019076465 |
| 0.024327993 | 0.001487583 | 0.02298452 |
| 0.010451443 | 0.001609234 | 0.00896386 |
| 0.006186574 | 0.001513426 | 0.00457734 |
| 0.012686271 | 0.001562207 | 0.011172845 |
| 0.006261605 | 0.001580993 | 0.004699398 |
| -0.037836006 | 0.001962901 | -0.039416999 |
| 0.024621228 | 0.002052984 | 0.022658327 |
| 0.014991891 | 0.001765415 | 0.012938907 |
| 0.000224058 | 0.001722681 | -0.001541356 |
| 0.012650642 | 0.002065529 | 0.010927962 |
| 0.022416691 | 0.001777976 | 0.020351162 |
0
Your first and third column are very similar -- which can contribute to your issue of "Matrix Not Positive Definitive". For more information on diagnosing these errors, please see our blog "Diagnosing a singular matrix."