(* :Title: EmpiricalLikelihood *)

(* :Author: Mark Fisher *)

(* :Context: EmpricalLikelihood` *)

(* :Mathematica Version: 5.0 *)

(* :History:
		Version 0.0 October 2002
		Version 0.1 May 2003
			Rewrote the syntax for EL.
			Streamlined the iteration loop (now in ELMainEngine).
			Added SandwichEstimator (with same syntax as EL).
		Version 0.2 October 2003
			Reorganized. Added Euclidiean likelihood and ELScore,
			among others.
		Version 0.3 January 2004
			Completely reorganized.
		Version 0.4 March 2004
			Changes in what ELMakeObject does.
		Version 0.5 April 2005
			Various tweaks.
*)

(* :Keywords: empirical likelihood, nonparametric confidence regions *)

(* :Summary: Compute empirical likelihood from data *)

(* :References:
	Art B. Owen, Empirical Likelihood, Chapman & Hall, 2001.
	Ron C. Mittelhammer, George G. Judge, and Douglas J. Miller,
		Econometric Foundations, Cambridge University Press, 2000
		[Chapter 12: Empirical Likelihood Estimation and Inference].
	*)

(* :Discussion:

The main function in the package is EL, which computes the value of
the "profile likelihood ratio function" at a specified point in the
"parameter" space. (Actually, with the default setting Deviance ->
True, EL returns the "deviance" (i.e., -2 * Log[el], where el is the
value of the profile likelihood ratio function).

The first step is to compute an ELMomentObject via

	mobj = ELMakeMomentObject[moments, vars, parms],

where moments is a list of moment conditions, vars is a list of
variables in the moment conditions, parms is a list of parameters in
the moment conditions. ELMakeMomentObject computes the function

	fun = Compile @@ {Join[parms, vars], moments}

(An option to ELMakeMomentObject determines whether the internal
function is compiled or not.) ELMakeMomentObject returns an
ELMomentObject that contains the parameters, the dimensions of the
moment matrix, the BlockingParameters, the optional weights, and the
function. These "internals" are accessible via the selectors
ELParameters, ELMomentDimensions, BlockingParameters, ELWeights, and
ELMomentObjectFunction.

Given a matrix of observations (data), a vector of numerical values (parvals),

	mobj[data, parmvals]

evaluates the function over the data. It computes

	vals = fun @@@ (Join[parmvals, #]& /@ data)

If required (as determined by options to ELMakeMomentObject),
it also computes

	vals = (Plus @@@ Partition[vals, M, L])/M

where M and L are the blocking parameters and

	vals = w * # & /@ vals

where w is a vector of weights.

To compute the deviance, use EL[mobj, data, parmvals]. To compute the
profile likelihood, use EL[mobj, data, parmvals, Deviance -> False]. Other
likelihoods can be specified using the option Likelihood. Valid
likelihoods are EL (empricial likelihood, the default), ET
(exponential tilting), CU (continuously updated GMM), and ETEL
(exponentially tilted empirical likelihood). Other options to EL
include GradientTolerance and UsePseudoInverse.

The function SandwichEstimator[mobj, gobj, data, parmvals] computes
the variance matrix evaluated at the parameter values parmvals,
where

gobj = ELMakeGradientObject[moments, vars, parmsyms]

is an ELMomentObject that contains the gradient tensor.
ELScore[mobj, gobj, data, parmvals] has the same syntax.


To draw from the likelihood as a Bayesian posterior, use
ELMetropolis[mobj, start, step, n], where step can be either a
scalar or a vector of the same dimension as start. ELMetropolis
takes the option RandomFunction; the default is RandomFunction ->
(RandomArray[NormalDistribution[0, 1], #]&).

The number of moments may be greater than, less than, or equal to the
number of parameters. With fewer estimating functions than parameters,
there will be underidentification, while with more estimating
functions than parameters there will be overidentification.

Given Likelihood -> EL, the code uses Owen's pseudo log function,
which equals the log for values greater than or equal to 1/n (where n
is the number of observations) and is quadratic for values less than n
(including negative values). Consequently, the code occasionally
returns negative weights (and complex values for the deviance) at
intermediate iterations.

The package also contains a number of functions that produce
estimating functions (functions of moments): ELRegressionMoments,
ELCorrelationMoments, and ELQuantileMoment. (See the usage
statements.)

*)

(* :Examples:
Needs["Statistics`NormalDistribution`"]
seedRandom[2]
data = RandomArray[NormalDistribution[0, 1], {30, 1}];
mean = Mean[data[[All, 1]]];
var = Variance[data[[All,1]]];

mobj1 = ELMakeMomentObject[{x - m}, {x}, {m}];
Plot[EL[mobj1, data, {m}, Deviance -> False],
	{m, -1, 1}, PlotRange -> All, GridLines -> {{mean}, {}}];
int1 = NIntegrate[EL[mobj1, data, {m}, Deviance -> False], {m, -1, 1}];
Plot[EL[mobj1, data, {m}, Deviance -> False]/int1,
	{m, -1, 1}, PlotRange -> All, GridLines -> {{mean}, {}}];

mobj2 = ELMakeMomentObject[{x - m, (x - m)^2 - v}, {x}, {m, v}];
ContourPlot[EL[mobj2, data, {m, v}, Deviance -> False,
	UsePseudoInverse -> True],
	{m, -1, 1}, {v, .25, 2}, PlotRange -> All, ContourShading -> False,
	PlotPoints -> 60, Epilog -> {PointSize[.02], RGBColor[0, 0, 0],
		Point[{mean, var}], RGBColor[1, 0, 0], Point[{0, 1}]}]
ContourPlot[EL[mobj2, data, {m, v}, Deviance -> True,
	UsePseudoInverse -> True],
	{m, -1, 1}, {v, .25, 2}, PlotRange -> All, ContourShading -> False,
	PlotPoints -> 60, Epilog -> {PointSize[.02], RGBColor[0, 0, 0],
		Point[{mean, var}], RGBColor[1, 0, 0], Point[{0, 1}]}]
*)

BeginPackage["EmpiricalLikelihood`",
	{"Statistics`DescriptiveStatistics`",
	"Statistics`NormalDistribution`",
	"Utilities`FilterOptions`",
	"MetropolisMCMC`"}]

ELMomentObject::usage = "ELMomentObject is an object that contains
the observations expressed in terms of the blocked moment conditions
as a function of the parameters."

ELMakeMomentObject::usage = "ELMakeMomentObject[moments, vars, parmsyms]
returns an ELMomentObject."

ELMomentObjectData::usage = "ELMomentObjectData[mobj] returns the data
contained in the ELMomentObject mobj."

ELMomentObjectFunction::usage = "ELMomentObjectFunction[mobj] returns
the function contained in the ELMomentObject mobj."

ComputeGradient::usage = "ComputeGradient is an option for ELMakeMomentObject,
which specifies whether to compute the gradient of the moments with respect to
the parameters. The default setting is ComputeGradient -> False."

ELParameters::usage = "ELParameters is an ELMomentObject selector.
ELParameters[mobj] returns the list of parameters in the ELMomentObject mobj."

ELMomentDimensions::usage = "ELMomentDimensions is an ELMomentObject selector.
ELMomentDimensions[mobj] returns the dimensions of the moment matrix in the
ELMomentObject mobj."

ELMakeGradientObject::usage = "ELMakeGradientObject[moments,
vars, parmsyms] returns an ELMomentObject that contains the gradient
tensor. It is equivalent to ELMakeMomentObject[moments, vars,
parmsyms, ComputeGradient -> True]."

ELMakeBothObjects::usage = "ELMakeBothObjects[moments, vars,
parmsyms] returns the pair of moment and gradient objects."

ELDevianceFromWeights::usage "ELDevianceFromWeights[type][weights] returns
the deviance for the specified type, where 0 indicates empirical likelihood,
1 indicates exponential tilting, and 2 indicates continuously updated GMM."

BlockingParameters::usage = "BlockingParameters is an option for
ELMakeMomentObject which specifies the block size M and the offset L
(between the start of adjacent blocks) for averaging over (i.e., blocking)
dependendent observations. The usage is BlockingParameters -> {M, L}, where
M and L are a positive integers. The default setting is BlockingParameters
-> {1, 1}."

BlockingFactor::usage = "BlockingFactor[T, {M, L}] returns T/(M *
Floor[(T - (M - L))/L]), where T is the number of observations, M is
the block size, and L is the offset between blocks. The chi-square
value is automatically multiplied by the blocking factor."

EL::usage = "EL[mobj, data, parmvals] returns -2*Log[el] where el is the
empirical likelihood function as determined by the estimating
functions embodied in the ELMomentObject mobj at the specified
parameter values parmvals. (By setting the option Deviance -> False,
EL will return el.) Computing the empirical likelihood involves
solving a minimization problem iteratively, the stopping of which is
controlled by the option GradientTolerance. The default setting is
GradientTolerance -> 10.^(-12). By setting the option Likelihood ->
ET, EL computes the weights using exponential tilting (a.k.a. maximum
entropy) rather than their empirical likelihood. Other valid settings
are CU for continuously updated GMM and ETEL for exponentially tilted
empirical likelihood. The default setting is Likelihood -> EL. The
other option for EL is UsePseudoInverse."

GradientTolerance::usage = "GradientTolerance is an option for EL that
specifies the norm of gradient for convergence. The default setting is
GradientTolerance -> 10.^(-12)."

UsePseudoInverse::usage = "UsePseudoInverse is an option for EL that
specifies whether to use PseudoInverse instead of Inverse. The default
setting is UsePseudoInverse -> False. If the setting is True, the
option Tolerance can be passed to PseudoInverse."

Deviance::usage "Deviance is an option for EL that specifies whether
to return the deviance (-2 log profile likelihood) or the profile
likelihood. The default is Deviance -> True."

Likelihood::usage = "Likelihood is an option for EL that specifies
which nonparametric likelihood to use. The default is Likelihood ->
EL. Other valid likelihoods are ET, CU, and ETEL."

ET::usage = "ET is a setting for the option Likelihood, which
specifies the use of exponential tilting (maximum entropy)."

CU::usage = "CU is a setting for the option Likelihood, which
specifies the use of continuous updating GMM (Euclidean likelihood)."

ETEL::usage = "ETEL is a setting for the option Likelihood, which
specifies the use of exponentially-tilted empirical likelihood."

GMM::usage = "GMM is a setting for the option Likelihood, which
specifies the use of GMM."

ELWeights::usage = "ELWeights is used when ReportAll -> True is set to
indicate the weights. ELWeights[mobj, data, parmvals] returns the weights
computed from the given ELMomentObject at the specified parameter values."

ELGradient::usage = "ELGradient[mobj, gobj, data, parmals] returns the
gradient. Does not work properly with Likelihood -> ETEL. Useful in FindMinimum."

ELOptimize::usage = "ELOptimize[mobj, gobj, data, start] returns the optimal
value. Not intended for ETEL."

SandwichEstimator::usage = "SandwichEstimator[mobj, gobj, data, parmvals]
returns the sandwich estimator of the parameters (given the moment
matrix in mobj and the gradient tensor in gobj) evaluated at the list
of parmvals. SandwichEstimator takes the option ELWeights -> weights;
the default is ELWeights -> Automatic, which sets the weights to
Table[1/n, {n}] where n is the number of observations.
SandwichEstimator also takes the option UsePseudoInverse."

ELRegressionMoments::usage = "ELRegressionMoments[model, variables,
parameters] returns a moment function for a (possibly nonlinear)
regression based on the given residual. The last variable listed is
assumed to be the dependent variable. Here is an example of a
nonlinear regression model: ELRegressionMoments[a(1-Exp[-b x]), {x,
y}, {a, b}]."

ELCorrelationMoments::usage "ELCorrelationMoments[] returns a moments
function for the means (m1 and m2), standard deviations (s1 and s2),
and correlation (rho) of two variables."

ELQuantileMoment::usage = "ELQuantileMoment[q, x, v] returns {{UnitStep[v -
x] - q}, {x}, {t}}, where q \[Element] (0,1) is numeric and x and v are
symbolic. Fixing q (as in this case), the estimating equation {UnitStep[v -
x] - q} produces the likelihood that value v is the q quantile. If we
reverse the roles of q and v (making v numeric and leaving q symbolic), the
estimating equation produces the likelihood that quantile of the specified
value v is q."

WeightedKernelDensityEstimate::usage =
"WeightedKernelDensityEstimate[x, h, data, weights] returns the
weighted kernel density estimate at x given bandwidth h and the data,
which are a list of pairs of data points with their weights.
Unweighted data (a vector) are weighted equally."

WKDEFunction::usage = "WKDEFunction[x, h, data, weights] is
the compiled version that WeightedKernelDensityEstimate calls."

WeightedRandomIntegerFunction::usage =
"WeightedRandomIntegerFunction[weights] returns a function that
returns a random integer in the range {1, Length[weights]}
with probability proportional to the weights, and which
can be used for resampling data."

WeightedResampling::usage = "WeightedResampling[data, weights, n]
returns a sample of n drawn from the data using probabilities
proportional to the weights."

AssembleData::usage = "AssembleData[{data1, data2,...}] merges a
number of (univariate) data sets into a single data set with dummy
variables that distinguish the individual data sets. The resulting
merged data set is useful for testing the equivalence of means, for
example."

ELScore::usage = "ELScore[mobj, gobj, data, parmvals]
returns the generalized empirical likelihood score statistic."

ELEMStatistic::usage = "ELEMStatistic[mobj, gobj, data, parmvals]
returns the generalized empirical likelihood EM statistic."

ELScoreBootstrap::usage = "ELScoresootStrap[mobj, mparmvals, n]."

ELBootstrap::usage = "ELBootstrap[mobj, parmvals, n]."

ELBias::usage = "ELBias[data, efuns, vars, parms] returns the second-order
bias as given by Newey and Smith.\n*** This function is currently out of service. ***"

ELLambda::usage = "ELLambda[vals, usePseudo, gradtol] returns the
vector of Lagrange multipliers. ELLambda take two optional arguments:
the Boolean that indicates the inverse function and the gradient
tolerance. The default Boolean is False, which indicates to use
Inverse, and the default tolerace is 10^(-12). For example:
ELLambda[vals, True, 10^(-6)]."

ELWeightsFromLambda::usage = "ELWeightsFromLambda[vals, lambda]."

ELWeightsFromLambdaEntropy::usage = "ELWeightsFromLambdaEntropy[vals,
lambda]."

ELWeightsFromLambdaEuclidean::usage =
"ELWeightsFromLambdaEuclidean[vals, lambda]."

ELWeightsEntropy::usage = "ELWeightsEntropy[vals, usePseudo, gradtol]."

ELWeightsEuclidean::usage = "ELWeightsEuclidean[vals, usePseudo, gradtol]."

ELWeightsEntropyCompiled::usage =
"ELWeightsEntropyCompiled[vals, usePseudo, gradtol]."

ELWeightsEuclideanCompiled::usage =
"ELWeightsEuclideanCompiled[vals, usePseudo, gradtol]."

ELValue::usage = "ELValue[vals, usePseudo, gradtol] returns the EL value."

ELValueEntropy::usage = "ELValueEntropy[vals, usePseudo, gradtol] returns the EL
entropy value."

ELValueEuclidean::usage = "ELValueEuclidean[vals, usePseudo, gradtol] returns the
EL Euclidean value. (gradtol plays no role.)"

ELValueCompiled::usage = "ELValueCompiled[vals, usePseudo, gradtol] returns the EL value."

ELValueEntropyCompiled::usage = "ELValueEntropyCompiled[vals, usePseudo, gradtol]
returns the EL value."

ELLambdaCompiled::usage = "ELLambdaCompiled[vals, usePseudo,
tol] is a compiled version of ELLambda, where usePseudo is either
True or False. All three arguments must be given."

ELValueEuclideanCompiled::usage = "ELValueEuclideanCompiled[vals, usePseudo, gradtol]
returns the EL Euclidean value."

ELLambdaEntropy::usage = "ELLambdaEntropy[vals, usePseudo, tol]."

ELLambdaEntropyCompiled::usage = "ELLambdaEntropyCompiled[vals, usePseudo, tol]."

ELLambdaEuclidean::usage = "ELLambdaEuclidean[vals, invfun, tol]."

ELLambdaEuclideanCompiled::usage = "ELLambdaEuclideanCompiled[vals, usePseudo, tol]."

ELParameterFunction::usage = "ELParameterFunction[data, moments, vars,
parms] returns a pure function of the parameters (compiled function
with the default option setting Comiled -> True)."

ELMetropolis::usage = "ELMetropolis[mobj, data, start, lambda, n]
makes n draws from the empirical likelihood by computing a Metropolis
Markov Chain Monte Carlo. The average step size is controled by
lambda. Options can be passed to EL."

NumberOfBlockedObservations::usage = "NumberOfBlockedObservations[T, {M, L}]
returns the number of blocked observations given the number of raw observations
T and the blocking parameters M and L, which is Floor[(T - (M - L))/L]."

ParameterDataJoin::usage = "ParameterDataJoin[parmvals, data] joins the
vector of parameter values with each row of the data matrix."

ELBlock::usage = "ELBlock[gvals, {M, L}] blocks the moment matrix gvals.
ELBlock[gvals, {M, L}, 2] blocks the gradient tensor gvals."

WeightValues::usage = "WeightValues[gvals, weights] weights
the matrix of gvals."

Begin["`Private`"]

BlockingFactor[T_, {M_, L_}] := T/(M * Floor[(T - (M - L))/L])

NumberOfBlockedObservations[T_, {M_, L_}] := Floor[(T - (M - L))/L]

ELMakeMomentObject::badblock = "BlockingParameters are not positive integers."

Options[ELMakeMomentObject] =
	{BlockingParameters -> {1, 1}, Compiled -> True, ELWeights -> Automatic,
	ComputeGradient -> False}

ELMakeMomentObject[moments_List, vars_List, parmsyms_List,
		opts___?OptionQ] :=
	Module[{compiled, gradient, w, M, L, moms, fun},
	compiled = TrueQ[Compiled /. {opts} /. Options[ELMakeMomentObject]];
	gradient = TrueQ[ComputeGradient /. {opts} /. Options[ELMakeMomentObject]];
	w = ELWeights /. {opts} /. Options[ELMakeMomentObject];
	{M, L} = BlockingParameters /. {opts} /. Options[ELMakeMomentObject];
	Switch[{M, L},
		{_Integer?Positive, _Integer?Positive},
			Null, (* OK *)
		{_, _},
			Message[ELMakeMomentObject::badblock]; Return[$Failed]
		];
	moms = If[gradient, D[moments, #]& /@ parmsyms, moments];
	fun = If[compiled, Compile, Function] @@ {Join[parmsyms, vars], moms};
	ELMomentObject[parmsyms, Dimensions[moms], {M, L}, w, fun]
	]



(* ELMomentObject formating *)
Format[mobj_ELMomentObject] :=
	StringForm["ELMomentObject[`1`, `2`, `3` <>]",
		mobj[[1]], mobj[[2]], mobj[[3]]
		]

(* ELMomentObject selectors *)
ELParameters[mobj_ELMomentObject] := mobj[[1]]

ELMomentDimensions[mobj_ELMomentObject] := mobj[[2]]

BlockingParameters[mobj_ELMomentObject] := mobj[[3]]

ELWeights[mobj_ELMomentObject] := mobj[[4]]

ELMomentObjectFunction[mobj_ELMomentObject] := mobj[[-1]]

BlockingFactor[mobj_ELMomentObject] :=
	BlockingFactor[ELMomentDimensions[mobj][[1]], BlockingParameters[mobj]]


(* ELMomentObject evaluation rules *)
ELMomentObject[_, _, {1, 1}, Automatic, fun_][data_List, parmvals_List] :=
	Transpose[fun @@@ ParameterDataJoin[parmvals, data]]

ELMomentObject[_, dims_, {M_, L_}, w_, fun_][data_List, parmvals_List] :=
	Module[{vals},
	vals = fun @@@ ParameterDataJoin[parmvals, data];
	If[M !=1 || L != 1,
		vals = ELBlock[vals, {M, L}, Length[dims]];
		];
	If[w =!= Automatic,
		vals = WeightValues[vals, w];
		];
	Transpose[vals]
	]

(* helpers for ELMomentObject[{args}] *)
ParameterDataJoin = Compile[{{parmvals, _Real, 1}, {data, _Real, 2}},
	Join[parmvals, #] & /@ data]

ELBlock[vals_, {M_, L_}, dim_:1] :=
	If[dim == 1, ELblock1[vals, M, L], ELblock2[vals, M, L]]

ELblock1 = Compile[{{vals, _Real, 2}, {M, _Integer, 0}, {L, _Integer, 0}},
    Plus @@@ Partition[vals, M, L]/M
    ]

ELblock2 = Compile[{{vals, _Real, 3}, {M, _Integer, 0}, {L, _Integer, 0}},
    Plus @@@ Partition[vals, M, L]/M
    ]

WeightValues = Compile[{{vals, _Real, 2}, {w, _Real, 1}},
	Transpose[w * # & /@ Transpose[vals]]
	]

ELMakeGradientObject[moments_List, vars_List, parmsyms_List,
		opts___?OptionQ] :=
	ELMakeMomentObject[moments, vars, parmsyms,
		ComputeGradient -> True, opts]

ELMakeBothObjects[moments_List, vars_List, parmsyms_List,
		opts___?OptionQ] :=
	Through[{ELMakeMomentObject, ELMakeGradientObject}[
		moments, vars, parmsyms, opts]]

(* private helper function *)
Likelihood2Integers[likelihood_] :=
	Switch[likelihood,
		EL,   {0, 0},
		ET,   {1, 1},
		CU,   {2, 2},
		ETEL, {1, 0},
		GMM,  {2, 3},
		_,    Message[EL::badlike]; {0, 0} (* just in case *)
		]

EL::badlike = "Bad Likelihood option. Likelihood being set to EL."

Options[EL] = {
	Deviance -> True,
	Likelihood -> EL,
	UsePseudoInverse -> False,
	GradientTolerance -> 10.^(-12)
	}

EL[mobj_ELMomentObject, data_?(MatrixQ[#, NumericQ]&),
		args_?(VectorQ[#, NumericQ]&), opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, usedev, wtype, liketype, weights, dev},
	{tol, pseudo, likelihood, usedev} = {GradientTolerance, UsePseudoInverse,
		Likelihood, Deviance} /. {opts} /. Options[EL];
	{wtype, liketype} = Likelihood2Integers[likelihood];
	weights = ELWeights[wtype][mobj[data, args], TrueQ[pseudo], tol];
	dev = BlockingFactor[Length[data], BlockingParameters[mobj]] *
		ELDevianceFromWeights[liketype][weights];
	If[TrueQ[usedev],
		dev,
		If[dev < 1416.6, Exp[-dev/2], 0.]
		]
	]

ELWeights[mobj_ELMomentObject, data_?(MatrixQ[#, NumericQ]&),
		args_?(VectorQ[#, NumericQ]&), opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, wtype, liketype},
	{tol, pseudo, likelihood} = {GradientTolerance, UsePseudoInverse,
		Likelihood} /. {opts} /. Options[EL];
	{wtype, liketype} = Likelihood2Integers[likelihood];
	ELWeights[wtype][mobj[data, args], TrueQ[pseudo], tol]
	]

(* not for ETEL *)

ELGradient[mobj_ELMomentObject, gobj_ELMomentObject,
		data_?(MatrixQ[#, NumericQ] &), parmvals_?(VectorQ[#, NumericQ]&),
		opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, type, vals, lam},
	{tol, pseudo, likelihood} = {GradientTolerance, UsePseudoInverse,
		Likelihood} /. {opts} /. Options[EL];
	type = Last @ Likelihood2Integers[likelihood];
	vals = mobj[data, parmvals];
	lam = ELLambda[type][vals, pseudo, tol];
	2 * Length[data] *
		(gobj[data, parmvals].lam).ELWeightsFromLambda[type][vals, lam]
	]

(* mainly for EL and CU, also for ET, but not for ETEL *)
ELOptimize[mobj_ELMomentObject, gobj_ELMomentObject,
		data_?(MatrixQ[#, NumericQ] &), start_?(VectorQ[#, NumericQ]&),
		opts___?OptionQ] :=
	Module[{fmopts, elopts, step, parms, fm},
	fmopts = FilterOptions[FindMinimum, opts];
	elopts = FilterOptions[EL, opts];
	Switch[Likelihood /. {elopts} /. Options[EL],
		ET, step = .05,
		_,  step = 1
		];
	parms = ELParameters[mobj];
	fm = FindMinimum @@ {EL[mobj, data, parms, elopts],
		Transpose[{parms, start}],
		fmopts,
		Gradient -> ELGradient[mobj, gobj, data, parms, elopts],
		Method -> {"QuasiNewton",
			"StepMemory" -> Automatic,
			"StepControl" -> {"LineSearch",
				"Method" -> Automatic,
				"MaxRelativeStepSize" -> step,
				"CurvatureFactor" -> .5,
				"DecreaseFactor" -> .01
					}
				}
		};
	{fm[[1]], parms /. fm[[2]]}
	]

ELBootstrap[mobj_ELMomentObject, data_?(MatrixQ[#, NumericQ]&),
		args_?(VectorQ[#, NumericQ]&), n_Integer?Positive, opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, wtype, liketype,
		vals, nobs, weights, rfun, bf, wfun, draws, bsvals, bsweights},
	{tol, pseudo, likelihood} = {GradientTolerance, UsePseudoInverse,
		Likelihood} /. {opts} /. Options[EL];
	{wtype, liketype} = Likelihood2Integers[likelihood];
	vals = mobj[data, args];
	weights = ELWeights[wtype][vals, pseudo, tol];
	rfun = WeightedRandomIntegerFunction[weights];
	bf = BlockingFactor[mobj];
	vals = Transpose[vals];
	nobs = Length[vals];
	pseudo = TrueQ[pseudo];
	wfun = ELWeights[wtype];
	dfun = ELDevianceFromWeights[liketype];
	Table[
		draws = Table[rfun[], {nobs}];
		bsvals = Transpose[ vals[[ draws ]] ];
		bsweights = wfun[bsvals, pseudo, tol];
		bf * dfun[bsweights],
		{n}]
	]

SandwichEstimator[mobj_ELMomentObject, gobj_ELMomentObject,
		data_?(MatrixQ[#, NumericQ]&), parmvals_?(VectorQ[#, NumericQ]&),
		opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, invfun, wtype, liketype, vals,
		n, weights, G, Omega},
	{tol, pseudo, likelihood} = {GradientTolerance, UsePseudoInverse,
		Likelihood} /. {opts} /. Options[EL];
	invfun = If[TrueQ[pseudo], PseudoInverse, Inverse];
	{wtype, liketype} = Likelihood2Integers[likelihood];
	vals = mobj[data, parmvals];
	n = Length[data];
	weights = If[liketype == 3, (* GMM *)
			Table[1/n, {n}],
			ELWeights[wtype][vals, TrueQ[pseudo], tol]
			];
	G = Transpose[gobj[data, parmvals], {1, 3, 2}].weights;
	Omega = If[wtype == 2, (* CU and GMM *)
			vals.Transpose[vals]/n,
			vals.(weights * Transpose[vals])
		];
	invfun[G.invfun[Omega].Transpose[G]]/n
	]

ELScore[mobj_ELMomentObject, gobj_ELMomentObject, data_?(MatrixQ[#, NumericQ]&),
		parmvals_?(VectorQ[#, NumericQ]&), opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, invfun, wtype, liketype, vals,
		weights, G, Omega, score},
	{tol, pseudo, likelihood} = {GradientTolerance, UsePseudoInverse,
		Likelihood} /. {opts} /. Options[EL];
	invfun = If[TrueQ[pseudo], PseudoInverse, Inverse];
	{wtype, liketype} = Likelihood2Integers[likelihood];
	vals = mobj[data, parmvals];
	lambda = ELLambda[wtype][vals, TrueQ[pseudo], tol];
	weights =  ELWeightsFromLambda[wtype][vals, lambda];
	G = Transpose[gobj[data, parmvals], {1, 3, 2}].weights;
	Omega = vals.(weights * Transpose[vals]);
	score = G.lambda;
	Dimensions[vals][[2]] * score.invfun[G.invfun[Omega].Transpose[G]].score
	]

ELEMStatistic[mobj_ELMomentObject, gobj_ELMomentObject, data_?(MatrixQ[#, NumericQ]&),
		parmvals_?(VectorQ[#, NumericQ]&), opts___?OptionQ] :=
	Module[{tol, pseudo, likelihood, invfun, wtype, liketype, vals,
		weights, G, Omega, score},
	{tol, pseudo, likelihood} = {GradientTolerance, UsePseudoInverse,
		Likelihood} /. {opts} /. Options[EL];
	invfun = If[TrueQ[pseudo], PseudoInverse, Inverse];
	{wtype, liketype} = Likelihood2Integers[likelihood];
	vals = mobj[data, parmvals];
	gbar = (Tr /@ vals)/Length[data];
	lambda = ELLambda[wtype][vals, TrueQ[pseudo], tol];
	weights =  ELWeightsFromLambda[wtype][vals, lambda];
	G = Transpose[gobj[data, parmvals], {1, 3, 2}].weights;
	Omega = vals.(weights * Transpose[vals]);
	score = G.gbar;
	Dimensions[vals][[2]] * score.invfun[G.invfun[Omega].Transpose[G]].score
	]


ELMetropolis[mobj_ELMomentObject, data_, start_List, step_, n_Integer?Positive,
		opts___?OptionQ] :=
	MetropolisMCMC[EL[mobj, data, #, Deviance -> False, opts]&, start, step, n, opts]

(************* beginning of various engines ***********)

(* for easy reference *)
(* 0 <=> EL, 1 <=> ET, 2 <=> CU *)

(*** compiled versions ***)

(* used extensively *)
ELWeights[0] = ELWeightsCompiled
ELWeights[1] = ELWeightsEntropyCompiled
ELWeights[2] = ELWeightsEuclideanCompiled

ELWeightsCompiled =
	Compile[{{vals, _Real, 2}, {usePseudo, True | False}, {tol, _Real, 0}},
	Module[{m = 1, n = 1, lambda = {1.}, d1 = {1.}, d2 = {1.}, grad = {1.}},
	{m, n} = Dimensions[vals];
	lambda = Table[0., {m}];
	If[usePseudo,
		(* then use PseudoInverse *)
		While[
			{d1, d2} = Transpose[
				If[# >= 1/n, {1/#, #^(-2)}, {n(2 - n #), n^2}] & /@
					(1 + lambda.vals)];
			grad = vals.d1;
			grad.grad > tol,
			(*update lambda*)
			lambda += PseudoInverse[vals.(d2 * Transpose[vals])].grad
			],
		(* else use Inverse *)
		While[
			{d1, d2} = Transpose[
				If[# >= 1/n, {1/#, #^(-2)}, {n(2 - n #), n^2}] & /@
					(1 + lambda.vals)];
			grad = vals.d1;
			grad.grad > tol,
			(*update lambda*)
			lambda += LinearSolve[vals.(d2 * Transpose[vals]), grad]
			]
		];
	(1/n)/(1 + lambda.vals)
	]]

ELWeightsEntropyCompiled =
	Compile[{{vals, _Real, 2}, {usePseudo, True | False}, {tol, _Real, 0}},
	Module[{
		lambda = Table[0., {Length[vals]}],
		n = Dimensions[vals][[2]],
		d = {1.}, grad = {1.}},
	If[usePseudo,
		(* then use PseudoInverse *)
		While[
			(* d = Exp[-lambda.vals]; *)
			d = If[# < 708.3, Exp[-#], 0]& /@ (lambda.vals);
			grad = vals.d;
			grad.grad > tol,
			(* update lambda *)
			lambda += PseudoInverse[vals.(d * Transpose[vals])].grad
			],
		(* else use LinearSolve *)
		While[
			(* d = Exp[-lambda.vals]; *)
			d = If[# < 708.3, Exp[-#], 0]& /@ (lambda.vals);
			grad = vals.d;
			grad.grad > tol,
			(* update lambda *)
			lambda += LinearSolve[vals.(d * Transpose[vals]), grad]
			]
		];
	d/(Plus @@ d)
	]]

ELWeightsEuclideanCompiled =
	Compile[{{vals, _Real, 2}, {usePseudo, True | False}, {tol, _Real}},
	Module[{n = Dimensions[vals][[2]], mbarhat = {1.}, devs = {{1.}}},
	mbarhat = vals.Table[1., {n}]/n;
	devs = vals - mbarhat;
	lambda = If[usePseudo,
		(* then use PseudoInverse *)
		n * PseudoInverse[devs.Transpose[devs]].mbarhat,
		(* else use Inverse *)
		n * LinearSolve[devs.Transpose[devs],mbarhat]
		];
	(1 - lambda.devs)/n
	]]

(* used extensively *)
ELDevianceFromWeights[0] = ELDevianceFromWeightsCompiled
ELDevianceFromWeights[1] = ELDevianceFromWeightsEntropyCompiled
ELDevianceFromWeights[2] = ELDevianceFromWeightsEuclideanCompiled

(*
ELDevianceFromWeightsCompiled =
	Compile[{{w, _Real, 1}},
	-2 * Tr[ Log[Length[w] * w] ]
	]

ELDevianceFromWeightsEntropyCompiled =
	Compile[{{w, _Real, 1}},
	2 * Tr[(# * Log[#])& [ Length[w] * w] ]
	]
*)

ELDevianceFromWeightsCompiled =
	With[{m = .1 * $MaxMachineNumber},
		Compile[{{w, _Real, 1}},
			If[Min[w] > 0.,
				-2 * Tr[ Log[Length[w] * w] ],
				m
			]
		]
	]

ELDevianceFromWeightsEntropyCompiled =
	Compile[{{w, _Real, 1}},
		2 * Tr[
			If[# == 0., 0, #*Log[#]]& /@ (Length[w]*w)
		]
	]

ELDevianceFromWeightsEuclideanCompiled =
	Compile[{{w, _Real, 1}},
	(#.#)& [(Length[w] * w) - 1]
	]

(* used occasionally *)
ELLambda[0] = ELLambdaCompiled
ELLambda[1] = ELLambdaEntropyCompiled
ELLambda[2] = ELLambdaEuclideanCompiled

ELLambdaCompiled =
	Compile[{{vals, _Real, 2}, {usePseudo, True | False}, {tol, _Real, 0}},
	Module[{m = 1, n = 1, lambda = {1.}, d1 = {1.}, d2 = {1.}, grad = {1.}},
	{m, n} = Dimensions[vals];
	lambda = Table[0., {m}];
	If[usePseudo,
		(* then use PseudoInverse *)
		While[
			{d1, d2} = Transpose[
				If[# >= 1/n, {1/#, #^(-2)}, {n(2 - n #), n^2}] & /@
					(1 + lambda.vals)];
			grad = vals.d1;
			grad.grad > tol,
			(*update lambda*)
			lambda += PseudoInverse[vals.(d2 * Transpose[vals])].grad
			],
		(* else use Inverse *)
		While[
			{d1, d2} = Transpose[
				If[# >= 1/n, {1/#, #^(-2)}, {n(2 - n #), n^2}] & /@
					(1 + lambda.vals)];
			grad = vals.d1;
			grad.grad > tol,
			(*update lambda*)
			lambda += LinearSolve[vals.(d2 * Transpose[vals]), grad]
			]
		];
	lambda]
	]

ELLambdaEntropyCompiled =
	Compile[{{vals, _Real, 2}, {usePseudo, True | False}, {tol, _Real, 0}},
	Module[{lambda = Table[0., {Length[vals]}], d = {1.}, grad = {1.}},
	If[usePseudo,
		(* then use PseudoInverse *)
		While[
			(* d = Exp[-lambda.vals]; *)
			d = If[# < 708.3, Exp[-#], 0]& /@ (lambda.vals);
			grad = vals.d;
			grad.grad > tol,
			(* update lambda *)
			lambda += PseudoInverse[vals.(d * Transpose[vals])].grad
			],
		(* else use Inverse *)
		While[
			(* d = Exp[-lambda.vals]; *)
			d = If[# < 708.3, Exp[-#], 0]& /@ (lambda.vals);
			grad = vals.d;
			grad.grad > tol,
			(* update lambda *)
			lambda += LinearSolve[vals.(d * Transpose[vals]), grad]
			]
		];
	-lambda]
	]

ELLambdaEuclideanCompiled =
	Compile[{{vals, _Real, 2}, {usePseudo, True | False}, {tol, _Real}},
	Module[{n = Dimensions[vals][[2]], mbarhat = {1.}, devs = {{1.}}},
	mbarhat = vals.Table[1., {n}]/n;
	devs = vals - mbarhat;
	If[usePseudo,
		(* then use PseudoInverse *)
		n * PseudoInverse[devs.Transpose[devs]].mbarhat,
		(* else use Inverse *)
		n * Inverse[devs.Transpose[devs]].mbarhat
		]
	]]

(* only used by ELScore; not compiled *)
ELWeightsFromLambda[0] = ELWeightsFromLambda
ELWeightsFromLambda[1] = ELWeightsFromLambdaEntropy
ELWeightsFromLambda[2] = ELWeightsFromLambdaEuclidean

ELWeightsFromLambda[vals_, lambda_] :=
	(1/Dimensions[vals][[2]])/(1 + lambda.vals)

ELWeightsFromLambdaEntropy[vals_, lambda_] := #/Tr[#]& @ Exp[lambda.vals]

ELWeightsFromLambdaEuclidean[vals_, lambda_] :=
	With[{n = Dimensions[vals][[2]]},
	(1 - lambda.(vals - vals.Table[1.,{n}]/n))/n
	]


(***************************** other stuff ********************************)

(* weighted resampling *)
WeightedResampling[data_, weights_, n_] :=
	With[{r = WeightedRandomIntegerFunction[weights]},
	data[[#]]& /@ Table[r[], {n}, {Length[data]}]
	]

WeightedRandomIntegerFunction[weights_] :=
	Evaluate[
		Ceiling[Interpolation[
		Transpose[{
			FoldList[Plus, 0, weights]/Tr[weights],
			Range[0, Length[weights]]}],
		InterpolationOrder -> 1][x]]
		]& /. x :> Random[]


(* WeightedKernelDensityEstimate *)

(* the following expression is modified to avoid non-machine numbers:
		(1/(E^((data - x)^2/(2*h^2))*Sqrt[2*Pi])).weights/h
*)
WKDEFunction = Compile[
	{{x, _Real, 0}, {h, _Real, 0}, {data, _Real, 1}, {weights, _Real, 1}},
		(If[# < 708.3, Exp[-#]/Sqrt[2*Pi]/h, 0]& /@ ((data - x)^2/(2 * h^2))) . weights
    ]


WeightedKernelDensityEstimate[x_?NumericQ, h_?NumericQ,
		data_?(VectorQ[#, NumericQ]&), weights_?(VectorQ[#, NumericQ]&)] :=
	WKDEFunction[x, h, data, weights]

WeightedKernelDensityEstimate[x_?NumericQ, h_?NumericQ,
		data_?(VectorQ[#, NumericQ]&)] :=
	With[{n = Length[data]},
	WeightedKernelDensityEstimate[x, h, data, Table[1/n, {n}]]
	]

WeightedKernelDensityEstimate[x_?(VectorQ[#, NumericQ]&), h_?NumericQ,
		data_?(MatrixQ[#, NumericQ]&), weights_?(VectorQ[#, NumericQ]&)] :=
	Times @@ MapThread[WeightedKernelDensityEstimate[#1, h, #2, weights]&,
		{x, Transpose[data]}]

WeightedKernelDensityEstimate[x_?(VectorQ[#, NumericQ]&),
		h_?(VectorQ[#, NumericQ]&), data_?(MatrixQ[#, NumericQ]&),
		weights_?(VectorQ[#, NumericQ]&)] :=
	Times @@ MapThread[WeightedKernelDensityEstimate[##, weights]&,
		{x, h, Transpose[data]}]

WeightedKernelDensityEstimate[x_?(VectorQ[#, NumericQ]&), h_?NumericQ,
		data_?(MatrixQ[#, NumericQ]&)] :=
	With[{n = Length[data]},
	WeightedKernelDensityEstimate[x, h, data, Table[1/n, {n}]]
	]

WeightedKernelDensityEstimate[x_?(VectorQ[#, NumericQ]&), h_?(VectorQ[#, NumericQ]&),
		data_?(MatrixQ[#, NumericQ]&)] :=
	With[{n = Length[data]},
	Times @@ MapThread[WeightedKernelDensityEstimate[##, Table[1/n, {n}]]&,
		{x, Transpose[data], h}]
	]

(***** functions that assemble various things *****)

AssembleData[data:{__?VectorQ}] := AssembleData[Partition[#, 1]& /@ data]

AssembleData[data:{__?MatrixQ}] :=
	Flatten[MapIndexed[Append[ #1, #2[[1]] ]&, data, {2}], 1]

(* make moment function *)
ELRegressionMoments[model_, variables_List, parameters_List] :=
	{(Last[variables] - model) * (D[model, #]& /@ parameters), variables, parameters}

ELCorrelationMoments[] :=
	{{Global`x1 - Global`m1, Global`x2 - Global`m2,
	(Global`x1 - Global`m1)^2 - Global`s1^2,
	(Global`x2 - Global`m2)^2 - Global`s2^2,
	(Global`x1 - Global`m1)(Global`x2 - Global`m2) -
	Global`s1 Global`s2 Global`rho},
	{Global`x1, Global`x2},
	{Global`m1, Global`m2, Global`s1, Global`s2, Global`rho}}

ELQuantileMoment[q_, x_, v_] := {{UnitStep[v - x] - q}, {x}, {t}}


End[]
EndPackage[]