(* :Title: GMM *)

(* :Author: Mark Fisher *)

(* :Summary:
This package provides functions to estimate GMM
(Generalized Method of Moments).
*)

(* :Notes to myself:
	Need to deal with
	ELMomentObject.
	Estimating functions.
*)

(* :Discussion:

This package remains in a somewhat unfinished state.

It does not establish its own Context to encapsulate its functions.

You should think of the GMM code here as an adjunct to my package
EmpiricalLikelihood, which creates the tools to deal with moment-based
estimating equations.

In particular, in order to use the GMM (the main engine here) one must use
ELMakeMomentObject and ELMakeGradientObject (or ELMakeBothObjects) to
produce the ELMomentObjects that contain the appropriate moment and
gradient information that GMM requires as input. These functions are
defined in EmpiricalLikelihood. (That package further requires my package
MetropolisMCMC.) You should read the usage notes for the functions
(and perhaps the discussion at the top of the package itself).

Here is a simple example that illustrates the basic usage.
In this example, the list of moment conditions is

{x - m, (x - m)^2 - v}

where {x} is list of symbols that represent columns in the data matrix, and
{m, v} is the list of symbols that represent the parameters to be
estimated.

Clear[m, v, x]
{mobj, gobj} = ELMakeBothObjects[{x - m, (x - m)^2 - v}, {x}, {m, v}];
data = Table[Random[], {30}, {1}]; (* a matrix with one column *)
GMM[mobj, gobj, data, {.5, .1}]

*)

(* :Mathematica Version: 5.0 or greater *)

Needs["EmpiricalLikelihood`"]
Needs["Utilities`FilterOptions`"]

GMM::usage = "GMM[mobj, gobj, data, start] computes the GMM estimator of
the parameters in the moment conditions (and gradients) in the
ELMomentObjects given the data and the starting values. GMM returns a list
containing the chi-square value and the estimated parameters. GMM takes the
option MaxFixedPointSteps. For a two-step estimator use MaxFixedPointSteps
-> 2. The fixed-point loop iterates on the pair {parvals, chol}, where
parvals are the parameter values and chol is the Cholesky decomposition of
the weighting matrix. Convergence is controlled by the option SameTest. The
default setting is SameTest -> (Max[Abs[#1[[1]] - #2[[1]]]] < 10^(-6) &),
which ignors chol."

MaxFixedPointSteps::usage = "MaxFixedPointSteps is an option for GMM which
specifies the maximum number of steps for the fixed point. The default
setting is MaxFixedPointSteps -> 100. For a two-step estimator use
MaxFixedPointSteps -> 2."

InitialWeightingMatrix::usage = "InitialWeightingMatrix is an option for
GMM. The default setting is InitialWeightingMatrix -> IdentityMatrix. The
other valid setting is Automatic which computes the \"optimal\" weighting
matrix from the starting values."

NeweyWestLagLength::usage = "NeweyWestLagLength is an option for GMM which
specifies the number of lags to include in the Newey-West correction. The
default setting is NeweyWestLagLength -> 0."

Options[GMM] = {
	InitialWeightingMatrix -> IdentityMatrix,
	MaxFixedPointSteps -> 100,
	SameTest -> (Max[Abs[#1[[1]] - #2[[1]]]] < 10^(-6) &),
	NeweyWestLagLength -> 0,
	Center -> False
	}

GMM::maxsteps = "MaxFixedPointSteps must be a positive integer."

(* returns the (uncentered) chi-square value *)
GMM[mobj_ELMomentObject, gobj_ELMomentObject, data_?(MatrixQ[#, NumericQ]&),
		start_?(VectorQ[#, NumericQ]&), opts___?OptionQ] :=
	Module[{fmopts, sameQ, maxsteps, w, lags, center, parmsyms, mnum, chol,
		fp, parvals, residual},
	fmopts = FilterOptions[FindMinimum, opts];
	{sameQ, maxsteps, w, lags, center} = {SameTest, MaxFixedPointSteps,
		InitialWeightingMatrix, NeweyWestLagLength, Center} /.
		{opts} /. Options[GMM];
	center = TrueQ[center];
	If[maxsteps < 1 || !IntegerQ[maxsteps],
		Message[GMM::maxsteps]; Return[$Failed]];
	parmsyms = ELParameters[mobj];
	mnum = ELMomentDimensions[mobj][[1]];
	Switch[w,
		IdentityMatrix,
			chol = IdentityMatrix[mnum],
		Automatic,
			vals = mobj[data, start];
			If[center, vals = vals - (Mean /@ vals)];
			w = Inverse[ GMMNeweyWest[vals, lags] ];
			chol = CholeskyDecomposition[(w + Transpose[w])/2],
		_,
			chol = CholeskyDecomposition[(w + Transpose[w])/2]
		];
	fp = FixedPointList[GMMLoop[#,
			mobj, gobj, parmsyms, data, lags, center, fmopts] &,
		{start, chol}, maxsteps, SameTest -> sameQ];
	parvals = fp[[-1, 1]];
	chol = fp[[-2, 2]]; (* use previous iteration *)
	residual = GMMResidual[{parvals, chol}, mobj, data];
	{residual.residual, parvals}
	]

(* 	Note: GMMFunction is not needed as long as the LevenbergMarquardt
	method is used and the Residual is computed. *)
GMMLoop[{parvals_, chol_},
		mobj_, gobj_, parmsyms_, data_, lags_, center_, opts___?OptionQ] :=
	Module[{parmsol, vals, w},
	parmsol = parmsyms /. Last @
		FindMinimum[Null, (* GMMFunction[{parmsyms, chol}, mobj, data] *)
			Transpose[{parmsyms, parvals}],
			Method -> {"LevenbergMarquardt",
				"Residual" -> GMMResidual[{parmsyms, chol}, mobj, data],
				"Jacobian" -> GMMJacobian[{parmsyms, chol}, gobj, data],
				"StepControl" -> {
					"TrustRegion", (* required step control *)
					(* these are the defaults; see advanced documentation *)
					"StartingScaledStepSize" -> Automatic,
					"MaxScaledStepSize" -> Infinity,
					"AcceptableStepRatio" -> 1/10000}},
				opts];
	vals = mobj[data, parmsol];
	If[center, vals = vals - (Mean /@ vals)];
	w = Inverse[ GMMNeweyWest[vals, lags] ];
	{parmsol, CholeskyDecomposition[(w + Transpose[w])/2]}
	]

GMMResidual[{parvals_?(VectorQ[#, NumericQ]&), chol_}, mobj_, data_] :=
	chol.(Tr /@ mobj[data, parvals])

GMMJacobian[{parvals_?(VectorQ[#, NumericQ]&), chol_}, gobj_, data_] :=
	chol.Map[Tr, Transpose[gobj[data, parvals], {2, 3, 1}], {2}]

GMMFunction[{parvals_?(VectorQ[#, NumericQ]&), chol_}, mobj_, data_] :=
	With[{residual = chol.(Tr /@ mobj[data, parvals])},
	(1/2) * residual.residual]

GMMNeweyWest[vals_, 0] := vals.Transpose[vals]

GMMNeweyWest[vals_, q_] := vals.Transpose[vals] +
	Sum[With[{t = (Drop[#, j]& /@ vals).Transpose[Drop[#, -j]& /@ vals]},
	(1 - j/(q + 1))(t + Transpose[t])], {j, 1, q}]

(* not a standard estimator *)
GMMOptimize[mobj_ELMomentObject, data_?(MatrixQ[#, NumericQ] &),
		start_?(VectorQ[#, NumericQ] &), opts___?OptionQ] :=
	Module[{parmsyms, fm},
	parmsyms = ELParameters[mobj];
	fm = FindMinimum[Null,
		Transpose[{parmsyms, start}],
		Evaluate[FilterOptions[FindMinimum, opts]],
		Method -> {"LevenbergMarquardt",
			"Residual" -> GMMOptimizeResidual[mobj, data, parmsyms, opts]}];
	{fm[[1]], parmsyms /. fm[[2]]}
	]

GMMOptimizeResidual[mobj_ELMomentObject, data_?(MatrixQ[#, NumericQ] &),
		parvals_?(VectorQ[#, NumericQ] &), opts___?OptionQ] :=
	Module[{center, lags, vals, g, w},
	{center, lags} = {Center, NeweyWestLagLength} /. {opts} /.
		Options[GMMChiSquareValue];
	vals = mobj[data, parvals];
	g = Tr /@ vals;
	If[TrueQ[center], vals = vals - g/Length[data]];
	w = Inverse[ GMMNeweyWest[vals, lags] ];
	Sqrt[2] * CholeskyDecomposition[(w + Transpose[w])/2].g
	]

Options[GMMCovarianceMatrix] = {UsePseudoInverse -> False}

(* can use SandwichEstimator with Likelihood -> GMM *)
GMMCovarianceMatrix[mobj_ELMomentObject, gobj_ELMomentObject,
		data_?(MatrixQ[#, NumericQ]&), parvals_?(VectorQ[#, NumericQ]&),
		opts___?OptionQ] :=
	Module[{pseudo, invfun, vals, G, Omega},
	pseudo = UsePseudoInverse /. {opts} /. Options[GMMCovarianceMatrix];
	invfun = If[TrueQ[pseudo], PseudoInverse, Inverse];
	vals = mobj[data, parvals];
	Omega = vals.Transpose[vals];
	G = Transpose[gobj[data, parvals], {1, 3, 2}].Table[1, {Length[data]}];
	invfun[G.invfun[Omega].Transpose[G]]
	]