function result = QLR_hac(y,x,z,w,heter,hac)
%This file does Andrews' (1993) test for structural breaks
%Barbara Rossi, Duke University, 2003

%OUTPUT is: [QLR test stat; pvQLR (its pvalue)]; 

%INPUT is:
%The constant is included automatically (so it cannot be included in p1)
%y is (T*1) dependent variable
%x is (T*p) independent variable whose param are time-varying
%z is (T*q) independent variable whose param are not time-varying
%w is (T*k) vector of instruments
%if no subsets then set z=0;
%heter=1 if robust to cond heteroskedasticity; heter=0 if cond homoskedasticity is imposed

load pvqlrsb.txt;   
[T,restr]=size(x);  
t2v=round(T*0.15):1:round(T*0.85); 
Andrews=zeros(length(t2v),1);  
for t2=t2v;
   Andrews(t2-round(T*0.15)+1,1)=chowgmmhac(y,x,z,w,t2,heter,hac);
end;
SupLR=max(Andrews);
pvSupLR=pvcalc(SupLR,pvqlrsb,restr);

bd = find(Andrews == SupLR);
brdate = bd + round(T*0.15);
result=[SupLR; pvSupLR; brdate];

%-------------------------------------------
function b = gmmbeta(y,z,x,heter)
T=length(y); Sxx=x'*x/T; Sxz=x'*z/T; Sxy=x'*y/T; [a,k]=size(x); [a,l]=size(z);
if k<l; disp('order condition violated'); 
else
W1=inv(Sxx);
b1=(inv(Sxz'*W1*Sxz))*(Sxz'*W1*Sxy);
e=y-z*b1;
S=x'*(diag(e.^2))*x./T;
W2=inv(S);
b2=(inv(Sxz'*W2*Sxz))*(Sxz'*W2*Sxy); b=b2;
end;

%-------------------------------------------
function res = gmmres(y,z,x,heter)
T=length(y); Sxx=x'*x/T; Sxz=x'*z/T; Sxy=x'*y/T; 
[a,k]=size(x); [a,l]=size(z);
if k<l; disp('order condition violated'); 
else
W1=inv(Sxx);
b1=(inv(Sxz'*W1*Sxz))*(Sxz'*W1*Sxy);
e=y-z*b1;
S=zeros(k,k); for i=1:T; S=S+(x(i,:)'*x(i,:)).*(e(i)^2); end; S=S./T;
W2=inv(S);
b2=(inv(Sxz'*W2*Sxz))*(Sxz'*W2*Sxy);
avarb=inv(Sxz'*W2*Sxz);
if heter==1; b=b2; elseif heter==0; b=b1; end;
end;
res=y-z*b;
%-------------------------------------------
function wald = chowgmmhac(y,x,z,w,t,heter,hac)

if z==0;
[T,p]=size(x); pi=t/T; [T,k]=size(w);
y1=y(1:t,:); x1=x(1:t,:); w1=w(1:t,:); 
y2=y(t+1:T,:); x2=x(t+1:T,:); w2=w(t+1:T,:); 
X=[x1,zeros(t,p);zeros(T-t,p),x2]; W=[w1,zeros(t,k);zeros(T-t,k),w2];
b=gmmbeta(y,X,W,heter); b1=b(1:p); b2=b(p+1:2*p);
e=gmmres(y,X,W,heter); 
M1=w1'*x1./t; 
M2=w2'*x2./(T-t); 
ew1=mmult(e(1:t),w1); 
ew2=mmult(e(t+1:T),w2); 
if heter==1;
S1=w1'*(diag(e(1:t).^2))*w1/t; 
S2=w2'*(diag(e(t+1:T).^2))*w2/(T-t); 
    if hac==1; S1=nw(ew1,t^(1/3)); S2=nw(ew2,(T-t)^(1/3)); end;
elseif heter==0; s2=(sum(e(1:T).^2))/(T); 
    if hac==1;   s2=nw(e(1:T),T^(1/3)); end; 
    S1=s2*w1'*w1/t; S2=s2*w2'*w2/(T-t); 
end;
V1=inv(M1'*(inv(S1))*M1); V2=inv(M2'*(inv(S2))*M2);
wald=T*(b1-b2)'*(inv(V1/pi+V2/(1-pi)))*(b1-b2);
else;
[T,p]=size(x); [T,q]=size(z); pi=t/T; [T,k]=size(w);
y1=y(1:t,:); x1=x(1:t,:); z1=z(1:t,:); w1=w(1:t,:);
y2=y(t+1:T,:); x2=x(t+1:T,:); z2=z(t+1:T,:); w2=w(t+1:T,:);
X=[x1,zeros(t,p),z1;zeros(T-t,p),x2,z2]; W=[w1,zeros(t,k);zeros(T-t,k),w2];
b=gmmbeta(y,X,W,heter); b1=b(1:p); b2=b(p+1:2*p);
e=gmmres(y,X,W,heter);
M1=w1'*x1./t; 
M2=w2'*x2./(T-t);
ew1=mmult(e(1:t),w1); 
ew2=mmult(e(t+1:T),w2); 
if heter==1; S1=w1'*(diag(e(1:t).^2))*w1/t; S2=w2'*(diag(e(t+1:T).^2))*w2/(T-t);
    if hac==1; S1=nw(ew1,t^(1/3)); S2=nw(ew2,(T-t)^(1/3)); end;
elseif heter==0; s2=(sum(e(1:T).^2))/(T);
    if hac==1;   s2=nw(e(1:T),T^(1/3)); end; 
S1=s2*w1'*w1/t;S2=s2*w2'*w2/(T-t);
end;
V1=inv(M1'*(inv(S1))*M1); V2=inv(M2'*(inv(S2))*M2);
wald=T*(b1-b2)'*(inv(V1/pi+V2/(1-pi)))*(b1-b2);
end;

%---------------------------------------------------
function pvalue=pvcalc(osser,tavola,kk)
numero=tavola(:,[1,kk+1+1]);
uno=0; 
if osser<=tavola(1,kk+2) pv=1; uno=uno+1; end;
if osser>=tavola(34,kk+2) pv=0; uno=uno+1; end;
if uno==0;
   rigal=find(numero(:,2)<=osser);
   riga=[max(rigal);max(rigal)+1]; sel=tavola(riga,[1,kk+1+1]);
   pv=sel(2,1)+(sel(2,2)-osser)*(sel(1,1)-sel(2,1))/(sel(2,2)-sel(1,2));
   end;
pvalue=pv;

%---------------------------------------------------------
function result = nw(y,qn)
[T,k]=size(y); ybar=ones(T,1)*((sum(y))/T);
dy=y-ybar;
G0=dy'*dy/T;
for j=1:qn-1;
   gamma=(dy(j+1:T,:)'*dy(1:T-j,:))./(T-1);
   G0=G0+(gamma+gamma').*(1-abs(j/qn));
end;
result=G0;