*
* @ELFCalc( options ) start end
* # list of series with moment conditions
*
* Computes the empirical likelihood for a set of moment conditions,
* that is, it finds the maximum of sum log w(i) subject to
*
*   (a) sum w(i)=1
*   (b) sum w(i)f(i)=0
*
* where f are the moment conditions input as a set of series
*
* This operates by minimizing the Lagrangean with respect to
* the Lagrange multipliers on the (b) conditions. The Lagrange
* multiplier on the adding up constraint is pegged at 1, with
* the w's corrected for this by scaling by their sum over the
* sample. The w's are concentrated out by solving their FONC's
* (w(i)**-1 - 1+lambda*f(i))=0, i.e. w(i)=1.0/(1+lambda*f(i)).
*
* Options:
*   ITERS=maximum number of iterations [40]
*   SUBITERS=maximum number of subiterations [10]
*   CVCRIT=convergence criterion on Lagrange multipliers [1.e-6]
*   TRACE/[NOTRACE]
*   W=series of weights
*
* Returns %FUNCVAL as the empirical log likelihood
*
* Revision Schedule
*  09/05 Written by Tom Doan, Estima
*
procedure ELFCalc start end
*
option integer iters	  40
option integer subiters	  10
option real    cvcrit     1.e-6
option switch  trace      0
option series  *w
*
local vect[integer] moments
local equation temp
local vector lambda lambdax g d
local integer time
local series wx
local integer n it subits
local real f0 scale func lmvalue sumw
local integer startl endl
*
* Get the list of moment conditions. Figure out the
* range based upon them.
*
enter(varying) moments
compute n=%rows(moments)
inquire(reglist) startl<<start endl<<end
# moments
compute %nobs=endl-startl+1
*
* This is used to make it easy to pull out the observations
* for a given time period.
*
equation temp
# moments
*
* Initialize the lambdas to zero. This gives rise to
* equal weights.
*
dim lambda(n)
compute lambda=%const(0.0)
set wx startl endl = 1.0
*
* Compute the concentrated Lagrangean at lambda=0.
*
compute f0=-%nobs*log(%nobs)+%nobs
do it=1,iters
   *
   * Compute the cross product and mean of the weighted moment
   * conditions
   *
   mcov(mean=g) startl endl wx
   # moments
   *
   * Compute the direction vector. Because g is a mean while %cmom
   * is a raw cross product, we need to multiply by %nobs to correct
   * for the scaling of g.
   *
   compute d=%nobs*inv(%cmom)*g
   compute scale=1.0
   do subits=1,subiters
      compute lambdax=lambda+scale*d
      *
      * Compute w's at the test set of lambdas
      *
      set wx startl endl = 1.0/(1+%dot(lambdax,%eqnxvector(temp,t)))
      sstats startl endl log(wx)>>func wx>>sumw
      *
      * Compute the concentrated Lagrangean
      *
      compute lmvalue=func-%nobs*log(sumw)+sumw
      *
      * Take the first valid set of values which improves (reduces)
      * the function value.
      *
      if %valid(lmvalue).and.lmvalue<f0
         break
      compute scale=scale*.5
   end do subits
   if %testdiff(lambda,lambdax)<cvcrit
      break
   compute lambda=lambdax,f0=lmvalue
   set wx startl endl = wx/sumw
end do iters
compute %funcval=func-%nobs*log(sumw)
if %defined(w)
   set w startl endl = wx
end


