/**************
IVTAR.PRG

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 GAUSS 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 */
n=100;
kx=4; 
let sig[2,2]=1 .5 .5 1;
e=rndn(n,2)*chol(sig);
q=rndn(n,1);
x=rndn(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.
*******************************************************************/

proc gmm_thresh(y,z1,z2,x,q,_conf,_conf1,_conf2,_reduced);
local xx,z1hat,rhohat,rhocf0,rhocf1,rhocf2,beta1,beta2,zhat,yhat,qhat,
qcf0,qcf1,qcf2,z,da,db,ya,yb,xa,xb,za,zb,se1,se2,jstat1,jstat2,qq,
beta1l,beta1u,beta2l,beta2u,qcf0i,qcf1i,qcf2i,qcf,i,qi,dai,dbi,beta1i,beta2i,se1i,se2i,jstat1i,jstat2i;

xx=x~z2; 
if _reduced==0;
  z1hat=xx*regress(z1,xx);
elseif _reduced==1;
  {z1hat,rhohat} = joint_thresh(z1,xx,q);
endif;
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;
ya=selif(y,da); xa=selif(xx,da); za=selif(z,da);
yb=selif(y,db); xb=selif(xx,db); zb=selif(z,db);
{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; 
endif;
qq=unique(q,1);
qq=selif(qq,(qq.<=qcf[2]));
qq=selif(qq,(qq.>=qcf[1]));

i=1; do while i<=rows(qq);
  qi=qq[i];
  dai=(q.<=qi);
  dbi=1-dai;
  ya=selif(y,dai); xa=selif(xx,dai); za=selif(z,dai);
  yb=selif(y,dbi); xb=selif(xx,dbi); zb=selif(z,dbi);
  {beta1i,se1i,jstat1i} = gmm_linear(ya,za,xa);
  {beta2i,se2i,jstat2i} = gmm_linear(yb,zb,xb);
  beta1l=minc(((beta1i-se1i*1.96)~beta1l)');
  beta1u=maxc(((beta1i+se1i*1.96)~beta1u)');
  beta2l=minc(((beta2i-se2i*1.96)~beta2l)');
  beta2u=maxc(((beta2i+se2i*1.96)~beta2u)');
i=i+1;endo;

"";"";
"Threshold Estimate	" qhat;
"Confidence Interval - Uncorrected             " qcf0;
"Confidence Interval - Het Corrected Quad " qcf1;
"Confidence Interval - Het Corrected NP    " qcf2;
"";"";
"Regime 1 : Threshold variable less than " qhat;
"Number of observations " sumc(da);"";
"      Estimates              S.E.              Lower                 Upper";
beta1~se1~beta1l~beta1u;
"";
"Regime 2 : Threshold variable greater than " qhat;
"Number of observations " sumc(db);"";
"      Estimates              S.E.              Lower                 Upper";
beta2~se2~beta2l~beta2u;
"";
retp(qhat);
endp;

/***************************************
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 

*******************************************************************/

proc (2) = joint_thresh(y,x,q);
local n,k,e,s0,qs,qn,sn,n1,n2,r,d,xx,ex,qhat,x1,x2,beta1,beta2,yhat,smin;

n = rows(y);
k=cols(x);
e=y-x*regress(y,x);
s0=det(moment(e,0));
n1=round(.05*n)+k;
n2=round(.95*n)-k;
qs=sortc(q,1);
qs=qs[n1:n2];
qs = unique(qs,1);
qn = rows(qs);
sn = zeros(qn,1);
r=1; do while r<=qn;
  d=(q.<=qs[r]);
  xx=x.*d;
  xx=xx-x*regress(xx,x);
  ex=e-xx*regress(e,xx);
  sn[r]=det(moment(ex,0));
r=r+1;endo;
r=minindc(sn);
smin=sn[r];
qhat=qs[r];
d=(q.<=qhat);
x1=x.*d;
x2=x.*(1-d);
beta1=regress(y,x1);
beta2=regress(y,x2);
yhat=x1*beta1+x2*beta2;

retp(yhat,qhat);
endp;

/***************************************
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(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)

*******************************************************************/

proc (5) = joint_thresh_ci(y,x,q,_conf,_graph);
local n,k,e,s0,qs,qn,sn,n1,n2,r,d,xx,ex,qhat,x1,x2,beta1,beta2,yhat,
smin,lr,r1,r2,qx,qh,m1,m2,g1,g2,sigq,hband,u,u2,f,df,eps,sige,kh,eta2,eta1,
c1,lr0,lr1,lr2,clr,tit,qcf_0,qcf_h1,qcf_h2,sig2;

n = rows(y);
k=cols(x);
e=y-x*regress(y,x);
s0=det(moment(e,0));
n1=round(.05*n)+k;
n2=round(.95*n)-k;
qs=sortc(q,1);
qs=qs[n1:n2];
qs = unique(qs,1);
qn = rows(qs);
sn = zeros(qn,1);
r=1; do while r<=qn;
  d=(q.<=qs[r]);
  xx=x.*d;
  xx=xx-x*regress(xx,x);
  ex=e-xx*regress(e,xx);
  sn[r]=det(moment(ex,0));
r=r+1;endo;
r=minindc(sn);
smin=sn[r];
qhat=qs[r];
d=(q.<=qhat);
x1=x.*d;
x2=x.*(1-d);
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 cols(y) .> 1;
  eta1=1; eta2=1;
else;
  r1 = (x*(beta1-beta2)).^2;
  r2 = r1.*(e.^2);
  qx = q.^(0~1~2);
  qh = qhat.^(0~1~2);
  m1 = r1/qx;
  m2 = r2/qx;
  g1 = qh*m1;
  g2 = qh*m2;
  eta1=(g2/g1)/sig2;
  sigq = sqrt(meanc((q-meanc(q)).^2));
  hband = 2.344*sigq/(n^(.2));
  u = (qhat-q)/hband;
  u2 = u.^2;
  f = meanc((1-u2).*(u2.<=1))*(.75/hband);
  df = -meanc(-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 = meanc(kh.*r1);
  g2 = meanc(kh.*r2);
  eta2 = (g2/g1)/sig2;
endif;
c1 = -2*ln(1-sqrt(_conf));
lr0 = (lr .>= c1);
lr1 = (lr .>= (c1*eta1));
lr2 = (lr .>= (c1*eta2));
if maxc(lr0)==1;
  qcf_0=qs[minindc(lr0)]~qs[qn+1-minindc(rev(lr0))];
else;
  qcf_0=qs[1]~qs[qn];
endif;
if maxc(lr1)==1;
  qcf_h1=qs[minindc(lr1)]~qs[qn+1-minindc(rev(lr1))];
else;
  qcf_h1=qs[1]~qs[qn];
endif;
if maxc(lr2)==1;
  qcf_h2=qs[minindc(lr2)]~qs[qn+1-minindc(rev(lr2))];
else;
  qcf_h2=qs[1]~qs[qn];
endif;

if _graph==1;
   graphset;
   _pdate = "";
   fonts("complex simgrma");
   title("Confidence Interval Construction for Threshold");
   ylabel("Likelihood Ratio Sequence in \202g\201"); 
   tit = "Threshold Variable";
   xlabel(tit);
   _pxsci = 8;
   _pysci = 8;
   clr = ones(qn,1)*c1;
   _plegstr="" $+ "LR]n[(\202g\201)\000" $+ ftocv(_conf*100,2,0) $+ "% Critical";
   if cols(y) .== 1;
    clr=clr~(clr*eta1)~(clr*eta2);
    _plegstr = _plegstr $+ "\000Hetero Corrected - 1\000Hetero Corrected - 2";
  endif;
   _plegctl = 1;
  let _pltype = 6 4 3 1;
  _plwidth =2;
   xy(qs,lr~clr);
endif;

retp(yhat,qhat,qcf_0,qcf_h1,qcf_h2);
endp;

/***************************************
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

*******************************************************************/
proc (3) = gmm_linear(y,z,x);
local pihat,xz,xy,beta,e,g,v,se,m,jstat;
 pihat=regress(z,x);
 xz=x'z; xy=x'y;
 beta=(pihat'xy)/(pihat'xz);
 e=y-z*beta;
 g=invpd(moment(x.*e,0)); 
 v=invpd(xz'g*xz);
 beta=v*(xz'g*xy);
 se=sqrt(diag(v));
 m=x'(y-z*beta);
 jstat=m'g*m;
retp(beta,se,jstat);
endp;

/***************************************
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

*******************************************************************/
proc regress(y,x);
local beta;
trap 1;
  beta=y/x;
  if scalerr(beta); 
    beta=pinv(moment(x,0))*(x'y);
  endif;
trap 0;
retp(beta);
endp;