%%%%%%%%%%%
% IVTAR.M %
%%%%%%%%%%%
% for questions, contact

% Bruce E. Hansen
% Department of Economics
% Social Science Building
% University of Wisconsin
% Madison, WI 53706-1393
% bhansen@ssc.wisc.edu
% http://www.ssc.wisc.edu/~bhansen/


% This is a MATLAB procedure.  It computes estimates and confidence
% intervals for threshold models with endogenous regressors.

% It computes the estimator described in
% "Instrumental Variable Estimation of a Threshold Model"
% written by Mehmet Caner and Bruce E. Hansen

% Be sure to have the Gauss graphics library active.

% If you run this program, it is currently set up to generate some simulated data 
% and estimate the model on this data.  
% For your own application, you can substitute your own data.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Example using Simulated Data %
function ivtar
n=100;
kx=4;
sig=[1 .5;.5 1];
e=normrnd(0,1,n,2)*chol(sig);
q=normrnd(0,1,n,1);
x=normrnd(0,1,n,kx);
z2=[ones(n,1),q];
z1=x*ones(kx,1).*3+e(:,2);
zz=(z1+z2*ones(2,1))*.1;
y=zz.*(q<0)-zz.*(q>0)+e(:,1);

qhat=gmm_thresh(y,z1,z2,x,q,.9,.8,1,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GMM_THRESH

% Computes threshold model with endogenous variables

% Format

% qhat = gmm_thresh(y,z1,z2,x,q,_conf,_conf1,_conf2,_reduced);

% Inputs
% y	nx1	dependent variable
% z1	nxk1	endogenous rhs variables
% z2	nxk2	exogenous rhs variables
% x	nxm	instrumental variables not included in z2, m>=k1
% q	nx1	threshold variable
% conf_	1x1	level for confidence interval for threshold, e.g. conf_=.9
% conf1_	1x1	confidence level for threshold interval used as input to slope intervals, e.g. conf2_=.8 
% conf2_	1x1	interval dummy - to determine method for threshold interval used as input to slope intervals
%		set to 0 for uncorrected (homoskedastic) interval
%		set to 1 for heteroskedastic correction by quadratic
%		set to 2 for heteroskedastic correction by nonparametric kernel
% reduced_	1x1	reduced form dummy
%		set to 0 if reduced form is linear (no threshold)
%		set to 1 if reduced form is a linear threshold model

% Output
% qhat	1x1	threshold estimate

% Most of the output is written to the screen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function qhat=gmm_thresh(y,z1,z2,x,q,conf_,conf1_,conf2_,reduced_);
xx=[x,z2];
if reduced_==0
    z1hat=xx*regress(z1,xx);
else
    [z1hat,rhohat]=joint_thresh(z1,xx,q);
end;
zhat=[z1hat,z2];
[yhat,qhat,qcf0,qcf1,qcf2] = joint_thresh_ci(y,zhat,q,conf_,1);

z=[z1,z2];
da=(q<=qhat);
db=1-da;
ind=0;
for i=1:length(da)
    if da(i)==1
        if ind==0
            ya=y(i,:);
            xa=xx(i,:);
            za=z(i,:);
            ind=1;
        else
            ya=[ya;y(i,:)];
            xa=[xa;xx(i,:)];
            za=[za;z(i,:)];
        end;
    end;
end;
ind=0;
for i=1:length(db)
    if db(i)==1
        if ind==0
            yb=y(i,:);
            xb=xx(i,:);
            zb=z(i,:);
            ind=1;
        else
            yb=[yb;y(i,:)];
            xb=[xb;xx(i,:)];
            zb=[zb;z(i,:)];
        end;
    end;
end;

[beta1,se1,jstat1] = gmm_linear(ya,za,xa);
[beta2,se2,jstat2] = gmm_linear(yb,zb,xb);


beta1l=beta1-se1*1.96;
beta1u=beta1+se1*1.96;
beta2l=beta2-se2*1.96;
beta2u=beta2+se2*1.96;

[yhat,qhat,qcf0i,qcf1i,qcf2i] = joint_thresh_ci(y,zhat,q,conf1_,0);
if conf2_==0
    qcf=qcf0i;
elseif conf2_==1
    qcf=qcf1i;
elseif conf2_==2
    qcf=qcf2i;
end;
qq=unique(q);
qqcf2=(qq<=qcf(2));
ind=0;
for i=1:length(qqcf2)
    if qqcf2(i)==1
        if ind==0
            temp=qq(i);
            ind=1;
        else
            temp=[temp;qq(i)];
        end;
    end;
end;
qq=temp;
qqcf1=(qq>=qcf(1)); 
ind=0;
for i=1:length(qqcf1)
    if qqcf1(i)==1
        if ind==0
            temp=qq(i);
            ind=1;
        else
            temp=[temp;qq(i)];
        end;
    end;
end;
qq=temp;
clear qqcf1;
clear qqcf2;
ind=0;
i=1;
while i<=length(qq(:,1));
    qi=qq(i);
    dai=(q<=qi);
    dbi=1-dai;
    for ii=1:length(dai)
        if dai(ii)==1
            if ind==0
                ya=y(ii,:);
                xa=xx(ii,:);
                za=z(ii,:);
                ind=1;
            else
                ya=[ya;y(ii,:)];
                xa=[xa;xx(ii,:)];
                za=[za;z(ii,:)];
            end;
        end;
    end;
    ind=0;
    for ii=1:length(dbi)
        if dbi(ii)==1
            if ind==0
                yb=y(ii,:);
                xb=xx(ii,:);
                zb=z(ii,:);
                ind=1;
            else
                yb=[yb;y(ii,:)];
                xb=[xb;xx(ii,:)];
                zb=[zb;z(ii,:)];
            end;
        end;
    end;
    [beta1i,se1i,jstat1i] = gmm_linear(ya,za,xa);
    [beta2i,se2i,jstat2i] = gmm_linear(yb,zb,xb);
    beta1l=min([(beta1i-se1i*1.96),beta1l]')';
    beta1u=max([(beta1i+se1i*1.96),beta1u]')';
    beta2l=min([(beta2i-se2i*1.96),beta2l]')';
    beta2u=max([(beta2i+se2i*1.96),beta2u]')';
    i=i+1;
end;
disp(' ');
disp(' ');
fprintf('Threshold Estimate:                       %f\n',qhat);
fprintf('Confidence Interval - Uncorrected:        %f  %f\n',qcf0(1),qcf0(2));
fprintf('Confidence Interval - Het Corrected Quad: %f  %f\n',qcf1(1),qcf1(2));
fprintf('Confidence Interval - Het Corrected NP:   %f  %f\n',qcf2(1),qcf2(2));
disp(' ');
disp(' ');
fprintf('Regime 1 : Threshold variable less than %f\n',qhat);
fprintf('Number of observations:                   %f\n',sum(da));
disp(' ');
disp('         Estimates          S.E.          Lower');
for i=1:length(se2)
    fprintf('%f    %f     %f    %f\n',beta1(i),se1(i),beta1l(i),beta1u(i));
end;
disp(' ');
fprintf('Regime 2 : Threshold variable greater than %f\n',qhat);
fprintf('Number of observations:                   %f\n',sum(db));
disp(' ');
disp('         Estimates          S.E.          Lower');
for i=1:length(se2)
    fprintf('%f    %f     %f    %f\n',beta2(i),se2(i),beta2l(i),beta2u(i));
end;
disp(' ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JOINT_THRESH

% Estimates the threshold in a multivariate threshold model

% Format
% [yhat,qhat] = joint_thresh(y,x,q);

% Inputs
% y	nxm	dependent variable(s)
% x	nxk	independent variables (should include constant)
% q	nx1	threshold variable

% Outputs
% yhat	nxm	Predicted values for y
% qhat	1x1	Threshold Estimate 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yhat,qhat]=joint_thresh(y,x,q);
n=length(y(:,1));
k=length(y(1,:));
e=y-x*regress(y,x);
s0=det(e'*e);
n1=round(.05*n)+k;
n2=round(.95*n)-k;
qs=sortrows(q,1);
qs=qs(n1:n2);
qs=unique(qs);
qn=length(qs(:,1));
sn=zeros(qn,1);
r=1;
while r<=qn
    d=(q<=qs(r));
    xx=x.*(d*ones(1,length(x(1,:))));
    xx=xx-x*regress(xx,x);
    ex=e-xx*regress(e,xx);
    sn(r)=det(ex'*ex);
    r=r+1;
end;
[temp,r]=min(sn);
smin=sn(r);
qhat=qs(r);
d=(q<=qhat);
x1=x.*(d*ones(1,length(x(1,:))));
x2=x.*((1-d)*ones(1,length(x(1,:))));
beta1=regress(y,x1);
beta2=regress(y,x2);
yhat=x1*beta1+x2*beta2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JOINT_THRESH_CI

% Estimates the threshold in a multivariate threshold model 
% and computes asymptotic confidence intervals

% Format
% [yhat,qhat,qcf_0,qcf_h1,qcf_h2] = joint_thresh_ci(y,x,q,_conf,_graph);

% Inputs
% y	nxm	dependent variable(s)
% x	nxk	independent variables (should include constant)
% q	nx1	threshold variable
% conf_	1x1	level for confidence interval, e.g. conf_=.9
% graph	1x1	Set _graph=1 to view graph of likelihood ratio
%		Set graph_=0 to not display graph

% Outputs
% yhat	nxm	Predicted values for y
% qhat	1x1	Threshold Estimate 
% qcf_0	1x2	Confidence interval for threshold, no heteroskedasticity correction
% qcf_h1	1x2	Confidence interval for threshold, heteroskedasticity correction 
% 		using quadratic variance estimate 
% 		(only if m=1)
% qcf_h2	1x2	Confidence interval for threshold, heteroskedasticity correction 
% 		using variance estimated by Epanechnikov kernel with automatic bandwidth
% 		 (only if m=1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [yhat,qhat,qcf_h0,qcf_h1,qcf_h2]=joint_thresh_ci(y,x,q,conf_,graph_);
n=length(y(:,1));
k=length(x(1,:));
e=y-x*regress(y,x);
s0=det(e'*e);
n1=round(.05*n)+k;
n2=round(.95*n)-k;
qs=sortrows(q,1);
qs=qs(n1:n2);
qs=unique(qs);
qn=length(qs(:,1));
sn=zeros(qn,1);
r=1;
while r<=qn
    d=(q<=qs(r));
    xx=x.*(d*ones(1,length(x(1,:))));
    xx=xx-x*regress(xx,x);
    ex=e-xx*regress(e,xx);
    sn(r)=det(ex'*ex);
    r=r+1;
end;
[temp,r]=min(sn);
smin=sn(r);
qhat=qs(r);
d=(q<=qhat);
x1=x.*(d*ones(1,length(x(1,:))));
x2=x.*((1-d)*ones(1,length(x(1,:))));
beta1=regress(y,x1);
beta2=regress(y,x2);
yhat=x1*beta1+x2*beta2;
e=y-yhat;
lr=n*(sn/smin-1);
sig2=smin/n;

if length(y(1,:))>1
    eta1=1;
    eta2=1;
else
    r1=(x*(beta1-beta2)).^2;
    r2=r1.*(e.^2);
    qx=[q.^0,q.^1,q.^2];
    qh=[qhat.^0,qhat.^1,qhat.^2];
    m1=(r1'/qx')';
    m2=(r2'/qx')';
    g1=qh*m1;
    g2=qh*m2;
    eta1=((g2'/g1')/sig2')';
    sigq=sqrt(mean((q-mean(q)').^2)');
    hband=2.344*sigq/(n^(.2));
    u=(qhat-q)/hband;
    u2=u.^2;
    f=mean((1-u2).*(u2<=1))'*(.75/hband);
    df=-mean(-u.*(u2<=1))'*(1.5/(hband^2));
    eps=r1-qx*m1;
    sige=(eps'*eps)/(n-3);
    hband=sige/(4*f*((m1(3)+(m1(2)+2*m1(3)*qhat)*df/f)^2));
    u2=((qhat-q)/hband).^2;
    kh=((1-u2)*.75/hband).*(u2<=1);
    g1=mean(kh.*r1)';
    g2=mean(kh.*r2)';
    eta2=((g2'/g1')/sig2')';
end;
c1=-2*log(1-sqrt(conf_));
lr0=(lr>=c1);
lr1=(lr>=(c1*eta1));
lr2=(lr>=(c1*eta2));
if max(lr0)==1
    [temp,indtemp1]=min(lr0);
    revlr0=lr0(1);
    for ii=2:length(lr0)
        revlr0=[lr0(ii);revlr0];
    end;
    [temp,indtemp2]=min(revlr0);
    qcf_h0=[qs(indtemp1),qs(qn+1-indtemp2)];
else
    qcf_h0=[qs(1),qs(qn)];
end;
if max(lr1)==1
    [temp,indtemp1]=min(lr0);
    revlr1=lr1(1);
    for ii=2:length(lr1)
        revlr1=[lr1(ii);revlr1];
    end;
    [temp,indtemp2]=min(revlr1);
    qcf_h1=[qs(indtemp1),qs(qn+1-(indtemp2))];
else
    qcf_h1=[qs(1),qs(qn)];
end;
if max(lr2)==1
    [temp,indtemp1]=min(lr2);
    revlr2=lr2(1);
    for ii=2:length(lr2)
        revlr2=[lr2(ii);revlr2];
    end;
    [temp,indtemp2]=min(revlr2);
    qcf_h2=[qs(indtemp1),qs(qn+1-indtemp2)];
else
    qcf_h2=[qs(1),qs(qn)];
end;

if graph_==1;
    figure;
    clr = ones(qn,1)*c1;
    plot(qs,lr,qs,clr,qs,clr*eta1,qs,clr*eta2);
    grid on;
    title('Confidence Interval Construction for Threshold');
    xlabel('Threshold Variable');
    ylabel('Likelihood Ratio Sequence in Gama');
    legend('LR_n(Gama)','90% Critical','Hetero Corrected-1','Hetero Corrected-2');
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GMM_LINEAR

% Computes the GMM estimator of a linear model

% Format
% [beta,se,jstat] = gmm_linear(y,z,x);

% Inputs
% y	nx1	dependent variable
% z	nxk	rhs variables
% x	nxl	instruments variables (should include constant and exogenous parts of z), l>=k

% Outputs
% beta	kx1	Regression slope estimates
% se	kx1	standard errors
% jstat	1x1	J Statistic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [beta,se,jstat]=gmm_linear(y,z,x);
pihat=regress(z,x);
xz=x'*z;
xy=x'*y;
beta=((pihat'*xy)'/(pihat'*xz)')';
e=y-z*beta;
xe=x.*(e*ones(1,length(x(1,:))));
g=inv(xe'*xe);
v=inv(xz'*g*xz);
beta=v*(xz'*g*xy);
se=sqrt(diag(v));
m=x'*(y-z*beta);
jstat=m'*g*m;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REGRESS

% Computes a linear regression. Uses generalized inverse if X'X is singular

% Format
% beta = regress(y,x);

% Inputs
% y	nxm	dependent variable(s)
% x	nxk	independent variables (should include constant)

% Output
% beta	kxm	Regression slope estimates

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function beta=regress(y,x);
warning off;
mw=' ';
beta=(y'/x')';
[mw,idw] = lastwarn;
lastwarn(' ');
warning on;
if (1-(mw==' '))
    beta=pinv(x'*x)*(x'*y);
end;