/*
** Author: Langche Zeng, (C) 1999
** langche@everest.fas.harvard.edu
**
** Hlogit: Heteroscedastic Logit Model.
**
** This procedure estimates the heteroscedastic logit model developed
** in Langche Zeng, "A Heteroscedastic GEV Discrete Choice Model,"
** forthcoming, \em{Socialogical Methods and Research}. 
** The heteroscedasticity is across individuals only, and
** the parameterization of \theta is the exponential form. Standard
** logit or "weighted logit" can also be estimated using this
procedure.
** The paper was completed with the support of NSF grant SBR-9753126.
** The code is provided here for public non-commercial use and standard
** disclaimers apply.
**
** USAGE
**
** library maxlik;
** {b,vc,llik,ret} = hlogit(y,x,w);  
**
** INPUTS to hlogit
**
** y = nx1 dependent variable vector
** x = nxk explanatory variable matrix (with constant in column 1)
** w = nxm weighting variable vector
**
** GLOBALS:
**
** _hlogit_out = 0 turn off unnecessary printing
**               1 print intermediate results of estimation
**
** _hlogit_w   =0, no weighting, standard logit
**             =1, if weighting llik instead of util. function
**             =2, if weighting util function (hlogit)
** 
** OUTPUTS
**
** b    = kx1 (if logit or weighted logit) or (k+m)x1 (if hlogit) 
**          vector of coefficients
** vc   = kxk (if logit or weighted logit) or (k+m)x(k+m) (if hlogit)
**          variance matrix
** llik = scalar log-likelihood at the maximum
** ret  = maxlik return code
*/

declare matrix _hlogit_out=0;
declare matrix _hlogit_w=0;
declare matrix k=0;                     @internal global, cols(x)@
declare matrix m=0;                     @internal global, cols(w)@

proc (4)=hlogit(y,x,w);
  local b,vc,llik,ret,g,stval;

  if _hlogit_out==0;
    _max_diagnostic=miss(0,0);          @ turn off unnecessary printing
@   
    __output=0;                 
  else;
    _max_diagnostic=3; 
    __output=1;
  endif;
  
  k=cols(x);
  m=cols(w);
  
  if _hlogit_w==1 and m>1;
    "Only one variable can be used to weight the llik function.";
    end;
  endif;

  if _hlogit_w/=2;                      @ logit, or weighted llik@
    stval=zeros(k,1);
  else;
    stval=zeros(k+m,1);
  endif;
  
  _max_covpar=3;                        @ white's se @
  {b,llik,g,vc,ret}=maxprt(maxlik(y~x~w,0,&ll,stval)); 
  
  retp(b,vc,llik,ret);
endp;

/* llik proc for hlogit */
proc ll(b,dta);
  local x,y,w,xb,wc,llik;
  y=dta[.,1];
  x=dta[.,2:(k+1)];
  xb=x*b[1:k];
  w=dta[.,(k+2):cols(dta)];
  if _hlogit_w==0;
    w=1;
    wc=0;
  elseif _hlogit_w==1;
    wc=0;
  else;
    wc=w*b[(k+1):rows(b)];
    w=1;
  endif;
  llik=w.*y.*xb.*exp(wc)-w.*ln(1+exp(xb.*exp(wc))); 
  retp(llik);
endp;  

