%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This code replicates the empirical example in 
%%% Yamamoto and Perron (2013) "Estimating and Testing Multiple Structural 
%%% Changes in Linear Models Using Band Spectral Regressions 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
clear all;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 0) Parameter values
% They are common to structural change tests (Bai and Perron (1998) code)
eps1=.15; 
hetdat=1; 
hetvar=0;
robust=1;
prewhit=0;
fixb=0;
m=5;

% Band spectral regressions
wl=pi*(1/16);  % lower band limit
wh=pi*(1/2);   % higher band limit
% The truncation rule is set to log(T) but you can change it in trun_rule.m

% The followings are specific to this example
data_transform=0;   % 0:level, 1 first-differences
method=2;  % 1:full, 2:truncation, 3:business-cycle
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1) Dataset
% Data set for hours-productivity
data0=csvread('data_hour.csv');
  % 1st column ln(productivity)
  % 2nd column ln(hours)
  prodata=data0(2:end,1)-data0(1:end-1,1);
  data=[prodata data0(2:end,2)]; 
  
% y: LHS variable  
yy=data(1:end,2);    
if data_transform==0 % linear detrend
   yl=yy(2:end,1);
   c=[1:1:size(yl,1)]';    
   gamma=inv(c'*c)*c'*yl;
   y0=yl-c*gamma;
elseif data_transform==1 % first difference
   yd=yy(2:end,1)-yy(1:end-1,1); 
   y0=yd;
end
y=y0(5:end,1);

% x: variables with breaking coefficients (with number p)
% z: variables with constant coefficients (with number k)
x0=data(2:end,1);     
x=[x0(4:end-1,1) x0(3:end-2,1) x0(2:end-3,1) x0(1:end-4,1)];
z=[];
bigt=size(y,1);
y=y-mean(y); 
x=x-ones(bigt,1)*mean(x); 

p=size(x,2);
k=size(z,2);
h=ceil(eps1*bigt);

disp('***********************************************************')  
disp('The options chosen are:')
disp('---------------------------------------------------------')     
if method==1
  disp('full spectral band')
elseif method==2
  disp('truncation')
elseif method==3
  disp(['spectral band from  ' num2str(wl/pi) '*pi to  ' num2str(wh/pi) '*pi'])  
end
disp('---------------------------------------------------------')     
disp(['eps1 = ' num2str(eps1)])
disp(['hetdat = ' num2str(hetdat)])
disp(['hetvar = ' num2str(hetvar)])
disp(['robust = ' num2str(robust)])
disp(['hetvar = ' num2str(hetvar)])
disp(['The maximum number of breaks is: ' num2str(m)])
disp('***********************************************************')     
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2) Break date estimation
% The following calls the procedure to obtain the break dates and the
% associated sum of squared residuals for all numbers of breaks between 1
% and m.
siglev=[10,5,2.5,1]';
if method==1 
   [glb,datevec,bigvec]=dating(y,x,h,m,p,bigt);
elseif method==2
   [glb,datevec,bigvec]=dating_trun(y,x,h,m,p,bigt);
elseif method==3
   [glb,datevec,bigvec]=dating_bsr(y,x,h,m,p,bigt,wl,wh);
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 3) Number of breaks
% The following calls the procedure to estimate various tests for the
% number of breaks. It also returns the UDmax and WDmax test.
disp('************************************')
disp('Output from the testing procedures')
disp('************************************')
disp('a) supF tests against a fixed number of breaks')
disp('-----------------------------------------')

%%% a) supF tests against a fixed number of breaks
ftest=zeros(m,1); 
wftest=zeros(m,1);
 
j=2; cv=getcv1(j,eps1);       
for i=1:m
   if method==1 
   ftest(i,1)=pftest(y,x,i,p,bigt,datevec,prewhit,robust,z,k,hetdat,hetvar); 
   wftest(i,1)=cv(p,1)*ftest(i,1)/cv(p,i);
   elseif method==2
   ftest(i,1)=pftest_trun(y,x,i,p,bigt,datevec,robust,z,k,hetdat,hetvar); 
   wftest(i,1)=cv(p,1)*ftest(i,1)/cv(p,i);
   elseif method==3
   ftest(i,1)=pftest_bsr(y,x,i,p,bigt,datevec,robust,z,k,hetdat,hetvar,wl,wh); 
   wftest(i,1)=cv(p,1)*ftest(i,1)/cv(p,i);    
   end    
   disp(['The supF test for 0 versus ' num2str(i) ' breaks (scaled by q) is: ' num2str(ftest(i,1))]) 
   
end
disp('--------------------------------------------------------') 
for j=1:4
 cv=getcv1(j,eps1);    
 disp(['The critical values at the ' num2str(siglev(j,1)) '% level are (for k=1 to ' num2str(m) '):' num2str(cv(p,1:m))])
end
    
%%% b) Dmax tests against an unknown number of breaks
UDmax=max(ftest);
WDmax=max(wftest); 
disp('***************************************************************')
disp('b) Dmax tests against an unknown number of breaks')
disp('---------------------------------------------------------')
disp(['The UDmax test is: ' num2str(max(ftest))])
for j=1:4
  cvm=getdmax(j,eps1);
  disp(['(the critical value at the ' num2str(siglev(j,1)) '% level is: ' num2str(cvm(p,1)) ')'])
end
disp('---------------------------------------------------------')
disp(['The WDmax test at the ' num2str(siglev(j,1)) '% level is: ' num2str(max(wftest))])
 
%%% 3) SupF(l+1|l) tests
disp('**************************************************************')
disp('c) supF(l+1|l) tests using global optimizers under the null')
disp('----------------------------------------------------------------')

supflval=zeros(m,2);
ndat0=datevec(1,1);
disp(['The break date for SupF is: ' num2str(round((ndat0-1)/4)+1949) ':Q' num2str(mod(ndat0+1,4)+1)])
for i=1:m-2           
  if method==1 
   [supfl,ndat]=spflp1(bigvec,datevec(1:i,i),i+1,y,x,h,p,prewhit,robust,z,k,hetdat,hetvar);
  elseif method==2
   [supfl,ndat]=spflp1_trun(bigvec,datevec(1:i,i),i+1,y,x,h,p,robust,z,k,hetdat,hetvar);  
  elseif method==3
   [supfl,ndat]=spflp1_bsr(bigvec,datevec(1:i,i),i+1,y,x,h,p,robust,z,k,hetdat,hetvar,wl,wh);   
  end
  supflval(i,:)=[supfl ndat];
  disp(['The supF(' num2str(i+1) '|' num2str(i) ') test is ' num2str(supfl)])
  disp(['It corresponds to a new break at: ' num2str(round((ndat-1)/4)+1949) ':Q' num2str(mod(ndat+1,4)+1)])
end
 
for j=1:4
  cv=getcv2(j,eps1);
  disp(['The critical values of supF(l+1|l) at the ' num2str(siglev(j,1)) '% level are (for i=1 to ' num2str(m) ' are: ' num2str(cv(p,1:m))])
end
disp('**************************************************************')
 