%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "Using OLS to estimate and test for structural changes in models with
% endogenous regressors" by Pierre Perron & Yohei Yamamoto
% Nov 30, 2012

% This code must be run with function files of Bai and Perron (1998) procedure 
% in this folder or provided by the authors website.
% http://people.bu.edu/perron/code.html

% The data is stored in py-data.csv (25KB) and was originally downloaded 
% from Andre Kurmann's web site (Kurmann, 2007) 

% This program consists of the following sections
% 0. Data 
% 0-1. parameter settings
% 0-2. data loading
% 0-3. data transformation
% 1. Structural change tests in the first stage
% 2. Estimation and tests  
% 2-1. first stage regression
% 2-2. second stage estimation
% 2-3. HHB procedure
% 2-4. subsample coefficients and moments
% 2-5. fixed regressors bootstrap (supressed initially)
% 2-6. CI for break dates

% To replicate Table 4, set numbr=0 and check the results in 2-2
% To replicate Table 5, set numbr=1,2,3 and checkthe results in 1
% To replicate Table 6, set numbr=1,2 and check the results in 2-4,2-6
% Please identify the name of the variable for each result by yourself.

clear all;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 0. Data
%%% 0-1. parameter settings
% parameters for Bai and Perron (1998) procedure
robust   = 1;
prewhit  = 0;
eps      = .1;
siglev   = [10;5;2.5;1];
signif   = 2;
hetdat   = 1;
hetvar   = 1;
fixb     = 0;

% other parameters
zlag     = 1;
order    = 2;
daty1    = 1947;     % First Year of Data Set 
datq1    = 1;        % First Quarter of Data Set 
daty2    = 1998;     % Last Year of Data Set 
datq2    = 2;        % Last Quarter of Data Set
esty1    = 1960;
estq1    = 1;
esty2    = 1997;
estq2    = 4;
estq1    = estq1-(zlag+1);
first    = (esty1-daty1-1)*4+estq1+(5-datq1); % first Obs for est. 
last     = (esty2-daty1-1)*4+estq2+(5-datq1); % last Obs for est.

numbr    = 1;      % number of reduced form breaks 

%%% 0-2. data loading
data     = csvread('py-data.csv');				

gdp      = data(first:last,1);  %nominal total GDP, 47:1 - 98:2
gdpq     = data(first:last,2);  %real total GDP, 47:1 - 98:2 
gbpu     = data(first:last,3);  %nominal GDP nonfarm business, 1947:1 - 1998:2
gbpuq    = data(first:last,4);  %real GDP nonfarm business, 1947:1 - 1998:2
gcq      = data(first:last,5);  %(real) personal consumption expenditures, total, 1947:1 -1997:2
gcnq     = data(first:last,6);  %real personal consumption expenditures, nondurables, 1947:1 -1997:2
gcsq     = data(first:last,7);  %real personal consumption expenditures, services, 1947:1 -1997:2
gifq     = data(first:last,8);  %real fixed investment, total, 1947:1 -1997:2
girq     = data(first:last,9);  %(real) fixed investment (app. private domestic), residential, 1947:1 -1997:2
ggeq     = data(first:last,10); %real government consumption exp & gross investm., 1947:1 -1997:2
lbmnu    = data(first:last,11); %hours of all persons, nonfarm business 47:1 - 98:2
lbcpu    = data(first:last,12); %nominal wage rate, nonfarm business, 47:1 -98:2
lblcpu   = data(first:last,13); %nominal unit labor cost, nonfarm business, 47:1 - 98:2
fyff     = data(first:last,14); %fed funds rate (effective), % per annum, 54:3 - 98:2
fygm3    = data(first:last,15); %U.S. Tbill 3-months rate, % per annum, 47:1 - 98:2
fygt10   = data(first:last,16); %U.S. Tbill 10-year rate, % per annum, 53:2 - 98:2
fm2      = data(first:last,17); %FM2, nominal money stock, 59:1 - 98:2
lbgdpu   = data(first:last,18); %implicit GDP deflator, nonfarm business, 47:1 - 98:2
gdpd     = data(first:last,19); %implicit GDP deflator, total, 47:1 - 97:2
psccom   = data(first:last,20); %spot market price index, 48:1 - 98:2

%%% 0-3.  data transformation 
T        = size(gdp,1);
inf      = log(gdpd(2:T,1))-log(gdpd(1:T-1,:));
s        = log(lblcpu)-log(gdpd);
s        = s(2:T);
dw       = log(lbcpu(2:T))-log(lbcpu(1:T-1))...
             -(log(lbgdpu(2:T))-log(lbgdpu(1:T-1)));
y        = log(gdpq(2:T));
spread   = fygt10(2:T)-fyff(2:T);
dcp      = log(psccom(2:T))-log(psccom(1:T-1))...
             -(log(lbgdpu(2:T))-log(lbgdpu(1:T-1)));

time     = (1:size(y,1))';
[pol,dm] = polyfit(time,y,order); 
ytrend   = polyval(pol,time); 
ygap     = y-ytrend;

Y        = inf(zlag+1:T-2,:);
Ylag     = inf(zlag  :T-3,:);
Yfut     = inf(zlag+2:T-1,:);
X        = [Yfut];
% here we can change the mc measure
S        = ygap(zlag+1:T-2,:);
% Note: use s for labor income share
%       use ygap for GDP gap

%  set the instruments
z        = [inf s ygap spread dw dcp];

z        = z(1:T-2,:);
[T,n]    = size(z);
Z        = [];
for i=1:zlag;
   Z     = [Z z(zlag+1-i:T-i,:)];
end
K        = size(Z,2);
Z        = [ones(T-zlag,1) Z];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1. Structural change test in the first stage
T        = size(X,1);
epsz     = 0.1;
hz       = round(T*epsz);
mz       = 4;
pz       = 0;
qz       = size(Z,2);
if pz==0
   [globz,bdz,bigvecz] = dating(X,Z,hz,mz,qz,T);   
else
   [globz,bdz,bigvecz] = nldat(X,Z,[],hz,mz,pz,qz,T,fixb,epsz,20,0,0);
end

RFsf     = zeros(mz,1);
RFseq    = zeros(mz-1,2);
for i=1:mz
  RFsf(i,1)   = pftest(X,Z,i,qz,T,bdz,prewhit,robust,[],pz,hetdat,hetvar);
end
for i=1:mz-1           
 [seqrf,ndat] = spflp1(bigvecz,bdz(1:i,i),i+1,X,Z,hz,qz,prewhit,robust,[],pz,hetdat,hetvar);
  RFseq(i,:)  = [seqrf ndat];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2. Estimation and tests
%%% 2-1. first stage regression
if numbr==0
  Xh     = Z*inv(Z'*Z)*Z'*X;  
elseif numbr==1
  bdz1 = bdz(1,1);
  Z1   = Z(1:bdz1,:);
  Z2   = Z(bdz1+1:end,:);
  X1   = X(1:bdz1,:);
  X2   = X(bdz1+1:end,:);
  Xh1  = Z1*inv(Z1'*Z1)*Z1'*X1; 
  Xh2  = Z2*inv(Z2'*Z2)*Z2'*X2;
  Xh   = [Xh1;Xh2];
elseif numbr==2
  bdz1 = bdz(1,2);
  bdz2 = bdz(2,2);
  Z1   = Z(1:bdz1,:);
  Z2   = Z(bdz1+1:bdz2,:);
  Z3   = Z(bdz2+1:end,:);
  X1   = X(1:bdz1,:);
  X2   = X(bdz1+1:bdz2,:);
  X3   = X(bdz2+1:end,:);
  Xh1  = Z1*inv(Z1'*Z1)*Z1'*X1; 
  Xh2  = Z2*inv(Z2'*Z2)*Z2'*X2;
  Xh3  = Z3*inv(Z3'*Z3)*Z3'*X3;
  Xh   = [Xh1;Xh2;Xh3];
elseif numbr==3
  bdz1 = bdz(1,3);
  bdz2 = bdz(2,3);
  bdz3 = bdz(3,3);
  Z1   = Z(1:bdz1,:);
  Z2   = Z(bdz1+1:bdz2,:);
  Z3   = Z(bdz2+1:bdz3,:);
  Z4   = Z(bdz3+1:end,:);  
  X1   = X(1:bdz1,:);
  X2   = X(bdz1+1:bdz2,:);
  X3   = X(bdz2+1:bdz3,:);
  X4   = X(bdz3+1:end,:);
  Xh1  = Z1*inv(Z1'*Z1)*Z1'*X1; 
  Xh2  = Z2*inv(Z2'*Z2)*Z2'*X2;
  Xh3  = Z3*inv(Z3'*Z3)*Z3'*X3;
  Xh4  = Z4*inv(Z4'*Z4)*Z4'*X4;
  Xh   = [Xh1;Xh2;Xh3;Xh4];  
end        

%%% 2-2. second stage estimation
C        = ones(T,1);
x        = [ C Ylag S X];
xh       = [ C Ylag S Xh];
w        = [ ];
wh       = [ ];
% NOTE: here we set variables with break in (x)
%                             without break in (w)
W        = [w x];
Wh       = [w xh]; 
k        = size(W,2);

%%% 2-2-1. full sample coefficient estimation 
betao    = inv(W'*W)*W'*Y;
betat    = inv(Wh'*Wh)*Wh'*Y;
reso     = Y-W*betao;
rest     = Y-Wh*betat;
esuo     = reso'*reso/T;
esut     = rest'*rest/T;
vo       = esuo*inv(W'*W);
vt       = esut*inv(Wh'*Wh);
stdo     = sqrt(diag(vo));
stdt     = sqrt(diag(vt));
parao    = [betao stdo];
parat    = [betat stdt];

haco     = correct(W, reso,0);
hact     = correct(Wh,rest,0);
vhaco    = inv(W'*W)*haco*T*inv(W'*W);
vhact    = inv(Wh'*Wh)*hact*T*inv(Wh'*Wh);
stdhaco  = sqrt(diag(vhaco));
stdhact  = sqrt(diag(vhact));

%%% 2-2-2. structural change tests
T        = size(Y,1);
h        = round(T*eps);
m        = 5;
p        = size(w,2);
q        = size(x,2);
if p==0
   [Oglob,Odatevec,Obigvec] = dating(Y,x, h,m,q,T);
   [Tglob,Tdatevec,Tbigvec] = dating(Y,xh,h,m,q,T);
else
   [Oglob,Odatevec,Obigvec] = nldat(Y,x, w, h,m,p,q,T,fixb,eps,20,0,0);
   [Tglob,Tdatevec,Tbigvec] = nldat(Y,xh,wh,h,m,p,q,T,fixb,eps,20,0,0);
end
OLSsf = pftest(Y,x, 1,q,T,Odatevec,prewhit,robust,w,p,hetdat,hetvar);
TLSsf = pftest(Y,xh,1,q,T,Tdatevec,prewhit,robust,wh,p,hetdat,hetvar);

Oseqf = zeros(m-1,2);
Tseqf = zeros(m-1,2);
for i=1:m-1           
 [Osupfl,Ondat] = spflp1(Obigvec,Odatevec(1:i,i),i+1,Y,x,h,q,prewhit,robust,[],p,hetdat,hetvar);
  Oseqf(i,:)    = [Osupfl Ondat];
 [Tsupfl,Tndat] = spflp1(Tbigvec,Tdatevec(1:i,i),i+1,Y,xh,h,q,prewhit,robust,[],p,hetdat,hetvar);
  Tseqf(i,:)    = [Tsupfl Tndat]; 
end

%%% 2-3. HHB procedure
if numbr==1
 Y1  = Y(1:bdz1,:);
 Y2  = Y(bdz1+1:end,:);
 xh1 = [ C(1:bdz1,:) Ylag(1:bdz1,:) S(1:bdz1,:) Xh1];   
 xh2 = [ C(bdz1+1:end,:) Ylag(bdz1+1:end,:) S(bdz1+1:end,:) Xh2];
 T1  = size(Y1,1);
 T2  = size(Y2,1);
 h1  = round(T1*eps);
 h2  = round(T2*eps);
 m1  = 1;
 m2  = 1;
 %p   = size(w,2);
 p   = 0;
 q   = size(xh2,2);
 if p==0
   [~,HHBbd1,~] = dating(Y1,xh1,h,m1,q,T1);  
   [~,HHBbd2,~] = dating(Y2,xh2,h,m2,q,T2);
 else
   [~,HHBbd1,~] = nldat(Y1,xh1,Wh1,h1,m1,p,q,T1,fixb,eps,20,0,0); 
   [~,HHBbd2,~] = nldat(Y2,xh2,Wh2,h2,m2,p,q,T2,fixb,eps,20,0,0);
 end
 HHBsf1 = pftest(Y1,xh1,1,q,T1,HHBbd1,prewhit,robust,[],p,hetdat,0);
 HHBsf2 = pftest(Y2,xh2,1,q,T2,HHBbd2,prewhit,robust,[],p,hetdat,0);

 xh1bar = pzbar(xh1,1,HHBbd1);
 CIHHB1 = interval(Y1,xh1,xh1bar,HHBbd1,q,m1,robust,prewhit,0,0,[],0);
 xh2bar = pzbar(xh2,1,HHBbd2);
 CIHHB2 = interval(Y2,xh2,xh2bar,HHBbd2,q,m2,robust,prewhit,0,0,[],0);
elseif numbr==2
 Y1  = Y(1:bdz1,:);
 Y2  = Y(bdz1+1:bdz2,:);
 Y3  = Y(bdz2+1:end,:);
 xh1 = [ C(1:bdz1,:) Ylag(1:bdz1,:) S(1:bdz1,:) Xh1];
 xh2 = [ C(bdz1+1:bdz2,:) Ylag(bdz1+1:bdz2,:) S(bdz1+1:bdz2,:) Xh2];
 xh3 = [ C(bdz2+1:end,:) Ylag(bdz2+1:end,:) S(bdz2+1:end,:) Xh3];
 T1  = size(Y1,1);
 T2  = size(Y2,1);
 T3  = size(Y3,1);
 h1  = round(T1*eps);
 h2  = round(T2*eps);
 h3  = round(T3*eps);
 m1  = 1;
 m2  = 1;
 m3  = 1;
 %p        = size(w,2);
 p   = 0;
 q        = size(xh1,2);
 if p==0
   [~,HHBbd1,~] = dating(Y1,xh1,h1,m1,q,T1);
   [~,HHBbd2,~] = dating(Y2,xh2,h2,m2,q,T2);
   [~,HHBbd3,~] = dating(Y3,xh3,h3,m3,q,T3);   
 else
   [~,HHBbd1,~] = nldat(Y1,xh1,Wh1,h1,m1,p,q,T,fixb,eps,20,0,0);
   [~,HHBbd2,~] = nldat(Y2,xh2,Wh2,h2,m2,p,q,T,fixb,eps,20,0,0);
   [~,HHBbd3,~] = nldat(Y3,xh3,Wh3,h3,m3,p,q,T,fixb,eps,20,0,0);
 end
 HHBsf1  = pftest(Y1,xh1,1,q,T1,HHBbd1,prewhit,robust,[],p,hetdat,0);   
 HHBsf2  = pftest(Y2,xh2,1,q,T2,HHBbd2,prewhit,robust,[],p,hetdat,0);   
 HHBsf3  = pftest(Y3,xh3,1,q,T3,HHBbd3,prewhit,robust,[],p,hetdat,0);    

 xh1bar = pzbar(xh1,1,HHBbd1(:,1));
 CIHHB1 = interval(Y1,xh1,xh1bar,HHBbd1(:,1),q,1,robust,prewhit,0,0,[],0);
 xh2bar = pzbar(xh2,1,HHBbd2(:,1));
 CIHHB2 = interval(Y2,xh2,xh2bar,HHBbd2(:,1),q,1,robust,prewhit,0,0,[],0);
 xh3bar = pzbar(xh3,1,HHBbd3(:,1));
 CIHHB3 = interval(Y3,xh3,xh3bar,HHBbd3(:,1),q,1,robust,prewhit,0,0,[],0);

elseif numbr==3
 Y1  = Y(1:bdz1,:);
 Y2  = Y(bdz1+1:bdz2,:);
 Y3  = Y(bdz2+1:bdz3,:);
 Y4  = Y(bdz3+1:end,:);
 xh1 = [ C(1:bdz1,:) Ylag(1:bdz1,:) S(1:bdz1,:) Xh1];
 xh2 = [ C(bdz1+1:bdz2,:) Ylag(bdz1+1:bdz2,:) S(bdz1+1:bdz2,:) Xh2];
 xh3 = [ C(bdz2+1:bdz3,:) Ylag(bdz2+1:bdz3,:) S(bdz2+1:bdz3,:) Xh3];
 xh4 = [ C(bdz3+1:end,:) Ylag(bdz3+1:end,:) S(bdz3+1:end,:) Xh4];

 T1  = size(Y1,1);
 T2  = size(Y2,1);
 T3  = size(Y3,1);
 T4  = size(Y4,1);
 h1  = round(T1*eps);
 h2  = round(T2*eps);
 h3  = round(T3*eps);
 h4  = round(T4*eps);
 m1  = 1;
 m2  = 1;
 m3  = 1;
 m4  = 1;
 %p        = size(w,2);
 p   = 0;
 q   = size(xh1,2);
 if p==0
   [~,HHBbd1,~] = dating(Y1,xh1,h1,m1,q,T1);
   [~,HHBbd2,~] = dating(Y2,xh2,h2,m2,q,T2);
   [~,HHBbd3,~] = dating(Y3,xh3,h3,m3,q,T3);   
   [~,HHBbd4,~] = dating(Y4,xh4,h4,m4,q,T4);   
 else
   [~,HHBbd1,~] = nldat(Y1,xh1,Wh1,h1,m1,p,q,T,fixb,eps,20,0,0);
   [~,HHBbd2,~] = nldat(Y2,xh2,Wh2,h2,m2,p,q,T,fixb,eps,20,0,0);
   [~,HHBbd3,~] = nldat(Y3,xh3,Wh3,h3,m3,p,q,T,fixb,eps,20,0,0);
   [~,HHBbd4,~] = nldat(Y4,xh4,Wh4,h4,m4,p,q,T,fixb,eps,20,0,0);
 end
 HHBsf1  = pftest(Y1,xh1,1,q,T1,HHBbd1,prewhit,robust,[],p,hetdat,0);   
 HHBsf2  = pftest(Y2,xh2,1,q,T2,HHBbd2,prewhit,robust,[],p,hetdat,0);   
 HHBsf3  = pftest(Y3,xh3,1,q,T3,HHBbd3,prewhit,robust,[],p,hetdat,0);  
 HHBsf4  = pftest(Y4,xh4,1,q,T4,HHBbd4,prewhit,robust,[],p,hetdat,0);  

 xh1bar = pzbar(xh1,1,HHBbd1(:,1));
 CIHHB1 = interval(Y1,xh1,xh1bar,HHBbd1(:,1),q,1,robust,prewhit,0,0,[],0);
 xh2bar = pzbar(xh2,1,HHBbd2(:,1));
 CIHHB2 = interval(Y2,xh2,xh2bar,HHBbd2(:,1),q,1,robust,prewhit,0,0,[],0);
 xh3bar = pzbar(xh3,1,HHBbd3(:,1));
 CIHHB3 = interval(Y3,xh3,xh3bar,HHBbd3(:,1),q,1,robust,prewhit,0,0,[],0);
 xh4bar = pzbar(xh4,1,HHBbd4(:,1));
 CIHHB4 = interval(Y4,xh4,xh4bar,HHBbd4(:,1),q,1,robust,prewhit,0,0,[],0);
end

%%% 2-4. subsample coefficients and moments
 odate   =125;
 xho     = [xh(1:odate,:) zeros(odate,k);zeros(T-odate,k) xh(odate+1:T,:)];
 xo      = [x(1:odate,:) zeros(odate,k);zeros(T-odate,k) x(odate+1:T,:)];
 
 Who     = [w xho];
 Wo      = [w xo];
 if p==0
   Who1  = [xh(1:odate,:)];
   Who2  = [xh(odate+1:end,:)];
   Wo1   = [x(1:odate,:)];
   Wo2   = [x(odate+1:end,:)]; 
 else
   Who1  = [w(1:odate,:) xh(1:odate,:)];
   Who2  = [w(odate+1:end,:) xh(odate+1:end,:)];
   Wo1   = [w(1:odate,:) x(1:odate,:)];
   Wo2   = [w(odate+1:end,:) x(odate+1:end,:)]; 
 end 
 
 beta     = inv(Who'*Who)*Who'*Y;
 beta_ols = inv(Wo'*Wo)*Wo'*Y;
 res      = Y-Who*beta;
 res1     = res(1:odate,:);
 res2     = res(odate+1:end,:);
 esu1     = res1'*res1/(odate-k);
 esu2     = res2'*res2/(T-odate-k);
 v1       = esu1*inv(Who1'*Who1);
 v2       = esu2*inv(Who2'*Who2);
 std1     = sqrt(diag(v1));
 std2     = sqrt(diag(v2));
 hac1     = correct(Who1,res1,0);
 hac2     = correct(Who2,res2,0);
 vhac1    = inv(Who1'*Who1)*hac1*(odate-k)*inv(Who1'*Who1);
 vhac2    = inv(Who2'*Who2)*hac2*(T-odate-k)*inv(Who2'*Who2); 
 stdhaco1 = sqrt(diag(vhac1));
 stdhaco2 = sqrt(diag(vhac2));

% moments in each regime
Qo1    = x(1:odate,:)'*x(1:odate,:)/odate;
Qo2    = x(odate+1:T,:)'*x(odate+1:T,:)/(T-odate);
Q1     = Qo1(end,end);
Q2     = Qo2(end,end);
reses  = Y-Who*beta;
reses1 = reses(1:odate,:);
reses2 = reses(odate+1:end,:);

% add eta
reses1 = reses1+(X(1:odate,:)-Xh(1:odate,:))*beta(4,1);
reses2 = reses2+(X(odate+1:T,:)-Xh(odate+1:T,:))*beta(8,1);

Phi1   = X(1:odate,1)'*reses1/odate;
Phi2   = X(odate+1:end,1)'*reses2/(T-odate);
bias   = Phi1/Q1-Phi2/Q2;

%%% 2-5. Fixed regressors bootstrap
%  numbt=1000;
%  simBTR = zeros(numbt,1);
%  simBTR21=zeros(numbt,1);
%  simBTR32=zeros(numbt,1);
%  simBTO = zeros(numbt,1);
%  simBTT = zeros(numbt,1);
%  m=1;
%   for simb=1:numbt
%      disp(simb)
%      Yb     = randn(T,1);
%      [globRb,datevecRb,bigvecRb] = dating(Yb,Z,hz,3,qz,T);     
%      [globOb,datevecOb,bigvecOb] = dating(Yb,x,h,m,size(x,2),T);
%      [globTb,datevecTb,bigvecTb] = dating(Yb,xh,h,m,size(xh,2),T);
%      sfRb    = pftest(Yb,Z,1,qz,T,datevecRb,prewhit,robust,[],p,hetdat,hetvar);     
%      sfOb    = pftest(Yb,x,1,size(x,2),T,datevecOb,prewhit,robust,[],p,hetdat,hetvar);
%      sfTb    = pftest(Yb,xh,1,size(xh,2),T,datevecTb,prewhit,robust,[],p,hetdat,hetvar);
%      simBTR(simb,1)=sfRb;
%      simBTO(simb,1)=sfOb;
%      simBTT(simb,1)=sfTb;    
%             
%      [sf21,Ondat] = spflp1(bigvecRb,datevecRb(1:1,1),1+1,X,Z,h,qz,prewhit,robust,[],p,hetdat,hetvar);
%      [sf32,Tndat] = spflp1(bigvecRb,datevecRb(1:2,2),2+1,X,Z,h,qz,prewhit,robust,[],p,hetdat,hetvar);
%      simBTR21(simb,1)=sf21;
%      simBTR32(simb,1)=sf32;
%  
%   end
%   cvBTR1  = quantile(simBTR,0.99);
%   cvBTR5  = quantile(simBTR,0.95);
%   cvBTR10 = quantile(simBTR,0.9);
%   cvBTO1  = quantile(simBTO,0.99);
%   cvBTO5  = quantile(simBTO,0.95);
%   cvBTO10 = quantile(simBTO,0.9);
%   cvBTT1  = quantile(simBTT,0.99);
%   cvBTT5  = quantile(simBTT,0.95);
%   cvBTT10 = quantile(simBTT,0.9);

%%% 2-6. CI for break date
[T,k]  = size(x);
% we have to manually put these!!!!!
odate  = 125;
ostart = 84;
oend   = 151;
tdate  = 125;
tstart = 84;
tend   = 151;
xb     = [x(ostart:odate,:) zeros(odate-ostart+1,k);zeros(oend-odate,k) x(odate+1:oend,:)];
xhb    = [xh(tstart:tdate,:) zeros(tdate-tstart+1,k);zeros(tend-tdate,k) xh(tdate+1:tend,:)];
if p==0
Wb     = xb;
Whb    = xhb;
W      = x;
Wh     = xh;
Yb     = Y(ostart:oend,:);
Yhb    = Y(tstart:tend,:);
else
Wb     = xb-w*inv(w'*w)*w'*xb;
Whb    = xhb-wh*inv(wh'*wh)*wh'*xhb;
W      = x-w*inv(w'*w)*w'*x;
Wh     = xh-wh*inv(wh'*wh)*wh'*xh;
Yb     = Y-w*inv(w'*w)*w'*Y;
Yhb    = Y-w*inv(wh'*wh)*wh'*Y;
end

gammao = inv(Wb'*Wb)*Wb'*Yb;
gammat = inv(Whb'*Whb)*Whb'*Yhb;
RR     = [eye(k)*(-1) eye(k)];
deltao = RR*gammao;
deltat = RR*gammat;
Qo1    = x(ostart:odate,:)'*x(ostart:odate,:)/(odate-ostart+1);
Qo2    = x(odate+1:oend,:)'*x(odate+1:oend,:)/(oend-odate);
Qt1    = xh(tstart:tdate,:)'*xh(tstart:tdate,:)/(tdate-tstart+1);
Qt2    = xh(tdate+1:tend,:)'*xh(tdate+1:tend,:)/(tend-tdate);
Uo     = reso(ostart:oend,:)'*reso(ostart:oend,:)/(oend-ostart+1-k);
Ut     = rest(tstart:tend,:)'*rest(tstart:tend,:)/(tend-tstart+1-k);
Uo1    = reso(ostart:odate,:)'*reso(ostart:odate,:)/(odate-ostart+1);
Uo2    = reso(odate+1:oend,:)'*reso(odate+1:oend,:)/(oend-odate);
Ut1    = rest(tstart:tdate,:)'*rest(tstart:tdate,:)/(tdate-tstart+1);
Ut2    = rest(tdate+1:tend,:)'*rest(tdate+1:tend,:)/(tend-tdate);

xio    = (deltao'*Qo2*deltao)/(deltao'*Qo1*deltao);
xit    = (deltat'*Qt2*deltat)/(deltat'*Qt1*deltat);
phio   = Uo2/Uo1*xio;
phit   = Ut2/Ut1*xit;

Lo     = deltao'*Qo1*deltao/Uo;
Lt     = deltat'*Qt1*deltat/Ut;

CIOy   = [odate-floor(11/Lo)-1 odate+floor(11/Lo)+1];
CITy   = [tdate-floor(11/Lt)-1 tdate+floor(11/Lt)+1];

