/* EL_TEST.PRG 

Empirical Likelihood Estimation and Testing and Confidence Interval Program

written by:

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 program calculates empirical likelihood estimates and tests of nonlinear hypotheses in GMM-type moment models.
It also constructs profile confidence intervals.

The empirical likelihood is optimized using the CONSTRAINED OPTIMIZATION package.  
The CO and PGRACH libraries must be active.  To do so, type "library co,pgraph" at the Gauss prompt. 

The user needs to provide a procedure "_moments" which computes the moments,
and provide a vector of starting values "_beta0".

To test a hypothesis, the user needs to provide a procedure "h_test" which is a (nonlinear) function
of the parameter b, which is a zero vector under the null hypothesis.

To construct a confidence interval for a parameter of interest, the use needs to provide a procedure "h_theta" which
computes this parameter from the entire parameter vector.

Format of Procedures:

m=_moments(b);
Input: 	b	kx1 parameter vector
Output:	m	nxm matrix of "moment conditions".  Rows (columns) correspond to observations (moment conditions).
Example:		proc _moments(b); retp(x.*(y-z*b)); endp;

h=h_test(b);
Input: 	b	kx1 parameter vector
Output:	h	rx1 vector of constraints.  Should be zero under H0
Example:		proc h_test(b); retp(b[1]+b[2]-2); endp;

theta=h_theta(b);
Input: 	b	kx1 parameter vector
Output:	theta	1x1 parameter
Example:		proc h_theta(b); retp(b[1]+b[2]); endp;

_beta0= kx1 	initial value for parameter vector


The user may also select the following :   **************/

_testh0=1;		@ Set to 1 to test the hypothesis in "h_test", else set to 0 to not do test @
_output = "el_test.out";	@ Name of output file @
_bnames="";	@ Optional.  kx1 vector of names for parameters
		     If left blank, the program will substitute "Beta1" etc @
_graph=0;		@ Set to 1 to compute and display profile criterion, else set to 0 @
_cregion=1;	@ Set to 1 to compute confidence region for "theta" @
_level=.90;		@ Level of confidence region, if computed @
_grid=21;		@ Number of gridpoints for graph of profile criterion, if computed @



@  Define procedures and initial conditions @

n=500;
x=ones(n,1)~rndn(n,3);
z1=x[.,2]+rndn(n,1)*3;
z2=x[.,3]+x[.,4]+rndn(n,1)*3;
z=z1~z2;
y=z*ones(2,1)+rndn(n,1)*10;

proc _moments(b);
local e;
  e=x.*(y-z*b);
retp(e);
endp;

@  Starting value for beta  @
_beta0=((z'x)*invpd(moment(x,0))*(x'y))/((z'x)*invpd(moment(x,0))*(x'z));

@ Test of Non-Linear Restrictions @
proc h_test(b);
local h;
  if _testh0==1;
     h=b[1]+b[2]-2;
  else;
    h=miss(0,0);
  endif;
retp(h);
endp;

@ Parameter of Interest: Theta @
proc h_theta(b);
local h;
  if _graph==1 or _cregion==1;
     h=b[1]+b[2];
  else;
    h=miss(0,0);
  endif;
retp(h);
endp;




/********************************************************************
Program Starts Here
********************************************************************/

output file =^_output reset;
output off;

e=_moments(_beta0);
n=rows(e);m=cols(e);k=rows(_beta0);
if (_bnames.=="") or (rows(_bnames) .ne k);
_names=0 $+ "Beta" $+ ftocv(seqa(1,1,k),1,0);
else;
_names=_bnames;
endif;
_names2=0 $+ "Lamba" $+ ftocv(seqa(1,1,m),1,0);
_names2=_names|_names2;
lambda=zeros(m,1);
_alpha0=_beta0|lambda;
coset;
trap 4;
_co_ParNames=_names2;
_co_MaxIters=500; 
_co_IneqProc=&_pp;
_co_GradProc=&_grad_el;
string _co_Options={stepbt newton forward none};

_co_EqProc = &_p_mo;
{alpha,l0,g,ret}=co(&_el,_alpha0);
beta=alpha[1:k];
lambda=alpha[k+1:k+m];

e=_moments(beta);
p=1./(n*(1+e*lambda));
g=gradp(&v_moment,beta);
gg=zeros(m,k);ee=zeros(m,m);
j=1; do while j<=n;
  gg=gg+g[(j-1)*m+1:(j-1)*m+m,.].*p[j];
  ee=ee+(e[j,.]'e[j,.]).*p[j];
j=j+1;endo;
v=invpd(gg'invpd(ee)*gg)/n;
se=sqrt(diag(v));
lr=2*(l0-n*ln(n));

format 12,6;
let fmt1[1,3]="*.*lf" 12 8;
let fmt2[1,3]="*.*lf" 12 6;
output on;
"Unconstrained Estimation";
"";"";
pr=printfm("Variable"~"Estimate"~"St Error",0~0~0,fmt1);"";
pr=printfm(_names~beta~se,0~1~1,fmt2);"";"";
"Negative Log Empirical Likelihood";
 l0;
"";
"ELR Test of Over-Identifying Restrictions";
pr=printfm("Test"~"dof"~"P_Value",0~0~0,fmt1);"";
pr=printfm(lr~(m-k)~cdfchic(lr,m-k),1~1~1,fmt2);"";"";
output off;


if _testh0 ==1;

h=h_test(beta);
r=rows(h);
hh=gradp(&h_test,beta)';
seh=sqrt(diag(hh'v*hh));

_co_EqProc = &_p_mo2;
{alpha1,l1,g,ret}=co(&_el,alpha);
beta1=alpha1[1:k];
lambda1=alpha1[k+1:k+m];

e=_moments(beta1);
p=1./(n*(1+e*lambda1));
g=gradp(&v_moment,beta1);
gg=zeros(m,k);ee=zeros(m,m);
j=1; do while j<=n;
  gg=gg+g[(j-1)*m+1:(j-1)*m+m,.].*p[j];
  ee=ee+(e[j,.]'e[j,.]).*p[j];
j=j+1;endo;
v1=invpd(gg'invpd(ee)*gg)/n;
hh=gradp(&h_test,beta1)';
v2=v1-v1*hh*invpd(hh'v1*hh)*(hh'v1);
se=sqrt(diag(v2));
lr=2*(l1-n*ln(n));
lr1=2*(l1-l0);

format 12,6;
let fmt1[1,3]="*.*lf" 12 8;
let fmt2[1,3]="*.*lf" 12 6;
output on;
"Estimated value of h(beta), s.e.'s";
 h~seh;"";"";
"**************************************************";
"";"";
"Constrained Estimation";
"";"";
pr=printfm("Variable"~"Estimate"~"St Error",0~0~0,fmt1);"";
pr=printfm(_names~beta1~se,0~1~1,fmt2);"";"";
"Negative Log Empirical Likelihood";
 l1;
"";
"ELR Test of Over-Identifying Restrictions";
pr=printfm("Test"~"dof"~"P_Value",0~0~0,fmt1);"";
pr=printfm(lr~(m-k)~cdfchic(lr,m-k),1~1~1,fmt2);"";"";
"ELR Test of H0";
pr=printfm("Test"~"dof"~"P_Value",0~0~0,fmt1);"";
pr=printfm(lr1~(m-k+r)~cdfchic(lr1,m-k+r),1~1~1,fmt2);"";"";
output off;

endif;

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

if (_graph==1) or (_cregion==1);
  thetahat=h_theta(beta);
  hh=gradp(&h_theta,beta)';
  setheta=sqrt(diag(hh'v*hh));
  if _graph==1;
    {c1,c2} = graph_el(&_elr,thetahat,setheta,alpha);
  else;
    c1=c_solve(&_elr,thetahat,setheta,alpha,1);
    c2=c_solve(&_elr,thetahat,setheta,alpha,2);
  endif;
  output on;
  "";
  "Estimates and Confidence Interval for Parameter of Interest";
  pr=printfm("Theta"~"s.e."~"C1"~"C2",0~0~0~0,fmt1);"";
  pr=printfm(thetahat~setheta~c1~c2,1~1~1~1,fmt2);"";"";
  output off;
endif;

/********************************************************************
PROCS
********************************************************************/

proc _el(alpha);
local pv,l;
  pv=(_moments(alpha[1:k])*alpha[k+1:k+m]+1)*n;
  if minc(pv) .> 0;
    l=sumc(ln(pv));
  else;
    l=n*ln(n);
  endif;
retp(l);
endp;

proc _el2(beta);
local pv,l,lambda0;
  lambda0=varget("lambda");
  pv=(_moments(beta)*lambda0+1)*n;
  if minc(pv) .> 0;
    l=sumc(ln(pv));
  else;
    l=n*ln(n);
  endif;
retp(l);
endp;

proc _grad_el(alpha);
local beta,g1,g2,g,e,lambda0,lambda_ok;
  beta=alpha[1:k];
  lambda0=alpha[k+1:k+m];
  lambda_ok=varput(lambda0,"lambda");
  g1=gradp(&_el2,beta)';
  e=_moments(beta);
  g2=sumc(e./(e*lambda0+1));
  g=g1|g2;
retp(g);
endp;

proc _p_mo(alpha);
local e,p,g,pv;
  e=_moments(alpha[1:k]);
  pv=(e*alpha[k+1:k+m]+1)*n;
  p=1./pv;
  g=e'p;
retp(g);
endp;

proc _p_mo2(alpha);
local g;
  g=_p_mo(alpha);
  if _testh0==1;
    g=g|h_test(alpha[1:k]);
  endif;
retp(g);
endp;

proc _p_mo3(alpha);
local g,theta0;
  theta0=varget("theta");
  g=(_p_mo(alpha))|(h_theta(alpha[1:k])-theta0);
retp(g);
endp;

proc v_moment(b);
local v;
  v=vecr(_moments(b));
  retp(v);
endp;

proc _pp(alpha);
local pv;
  pv=(_moments(alpha[1:k])*alpha[k+1:k+m]+1)*n-1;
retp(pv);
endp;

proc (3) = _elr(theta0,alpha0,a0);
local theta_ok,l0,l,lr,a,g,ret;
  _co_EqProc = &_p_mo;
  l0=_el(alpha0);
  theta_ok=varput(theta0,"theta");
  _co_EqProc = &_p_mo3;
  {a,l,g,ret}=co(&_el,a0);
  lr=2*(l-l0);
retp(lr,a,ret);
endp;

proc (2) = graph_el(&like,thetahat,setheta,alpha0);
local like:proc,cr,a0,l0,lr,j,theta0,theta_ok,l,g,ret,theta1,theta2,aa,lrs,d,d1,d2,c1,c2,
c1a,c1b,c2a,c2b,lev,tit1,tit2,z,a1;
  cr=cdfchii(_level,1);
  a0=alpha0;
  lr=0;j=0;
  do while lr < (cr+1); 
    j=j+1;
    theta0=thetahat-setheta*j;
    {lr,a0,ret}=like(theta0,alpha0,a0);
     theta0~lr;
  endo;
  theta1=theta0;
  a0=alpha0;
  lr=0;j=0;
  do while lr < (cr+1); 
    j=j+1;
    theta0=thetahat+setheta*j;
    {lr,a0,ret}=like(theta0,alpha0,a0);
     theta0~lr;
  endo;
  theta2=theta0;
  aa=seqa(theta1,(theta2-theta1)/(_grid-1),_grid);
  lrs=zeros(_grid,1);
  for i (_grid,1,-1);
    theta0=aa[i];
    {lr,a0,ret}=like(theta0,alpha0,a0);
     lrs[i]=lr;
     theta0~lr;
  endfor;
  z=ones(_grid,1)*cr;
  d=(lrs.<cr);
  d1=maxindc(d); 
  d2=_grid+1-maxindc(rev(d));
  c1a=aa[d1-1];
  c1b=aa[d1];
  c1=((lrs[d1-1]-cr)*c1b+(cr-lrs[d1])*c1a)/(lrs[d1-1]-lrs[d1]);
  c2a=aa[d2];
  c2b=aa[d2+1];
  c2=((lrs[d2+1]-cr)*c2a+(cr-lrs[d2])*c2b)/(lrs[d2+1]-lrs[d2]);
  graphset;
  fonts("complex");
  _plwidth=2;
  _pcolor=15;
  lev=ftocv(_level*100,1,0);
  tit1="" $+ lev $+ "% Empirical Likelihood Confidence Interval for Theta";
  tit2="" $+ "Profile EL Criterion and " $+ lev $+ "% Critical Value";
  title(tit1);
  xlabel("Theta");
  ylabel(tit2);
  _plwidth=2;
  _pdate="";
  _psym=
  (c1~cr~2~3~15~1~0)|
  (c2~cr~2~3~15~1~0);
  xy(aa,lrs~z);
retp(c1,c2);
endp;


proc c_solve(&like,thetahat,setheta,alpha0,endpoint);
local like:proc,gg,g_l,g_r,theta_l,theta_r,theta,a0,cr,sn,dtheta,nogo,j,lr,ret,a1,try;
  gg=0; g_l=0;  theta_l=thetahat; g_r=0;  theta_r=thetahat; theta=thetahat;a0=alpha0;
  cr=cdfchii(_level,1);
  if endpoint==1; sn=-1; elseif endpoint==2; sn=1;endif;
  dtheta=sqrt(cr)*setheta*.6*sn;  
  nogo=1;
  for j (1,100,1);
    theta=theta+dtheta;
    {lr,a0,ret}=like(theta,alpha0,a0);
    if ret>2;
       theta=theta-dtheta;
       dtheta=dtheta/2;
       a0=alpha0;
    else;
      gg=lr; 
      theta_l=theta_r;
      g_l=g_r;
      theta_r=theta;
      g_r=gg;
    endif;
    try = (abs(cr-gg) .< .01) +(gg>cr);
    if try==1;
      nogo=0; break;
    endif;
  endfor;
  if nogo==1; goto _end; endif;
  try= (abs(cr-g_r) .< .01)+(abs(theta_l-theta_r) .< .000001);
  if try==1; goto _end; endif;
  a1=a0; nogo=1;
  for j (1,20,1);
    theta=(theta_l+theta_r)/2;
   {lr,a1,ret}=like(theta,alpha0,a1);
    if ret>2;
     {lr,a1,ret}=like(theta,alpha0,a0);
    endif;
    if lr>cr;
      theta_r=theta; g_r=lr;
    else;
      theta_l=theta; g_l=lr;
    endif;
    try= (abs(cr-lr) .< .01)+(abs(theta_l-theta_r) .< .000001);
    if try==1; nogo=0; break; endif;
    theta=(theta_l*(g_r-cr)+theta_r*(cr-g_l))/(g_r-g_l);
   {lr,a1,ret}=like(theta,alpha0,a1);
    if ret>2;
     {lr,a1,ret}=like(theta,alpha0,a0);
    endif;
    if lr>cr;
      theta_r=theta; g_r=lr;
    else;
      theta_l=theta; g_l=lr;
    endif;
    try= (abs(cr-lr) .< .01)+(abs(theta_l-theta_r) .< .000001);
    if try==1; nogo=0; break; endif;
  endfor;
  _end:
retp(theta);
endp;