Chapter 21. Maximum likelihood estimation
190
The mle block in the above script can therefore be modified as follows:
mle logl = y*ln(P) + (1-y)*ln(1-P)
series ndx = b0 + b1*x1 + b2*x2 + b3*x3
series P = cnorm(ndx)
series m = dnorm(ndx)*(y/P - (1-y)/(1-P))
deriv b0 = m
deriv b1 = m*x1
deriv b2 = m*x2
deriv b3 = m*x3
end mle --verbose
Note that the params statement has been replaced by a series of deriv statements; these have the
double function of identifying the parameters over which to optimize and providing an analytical
expression for their respective score elements.
21.6 Debugging ML scripts
We have discussed above the main sorts of statements that are permitted within an mle block,
namely
 auxiliary commands to generate helper variables;
 deriv statements to specify the gradient with respect to each of the parameters; and
 a params statement to identify the parameters in case analytical derivatives are not given.
For the purpose of debugging ML estimators one additional sort of statement is allowed: you can
print the value of a relevant variable at each step of the iteration. This facility is more restricted
then the regular print command. The command word print should be followed by the name of
just one variable (a scalar, series or matrix).
In the last exampleabovea key variablenamedm wasgenerated, formingthe basisfor the analytical
derivatives. To track the progress of this variable one could add a print statement within the ML
block, as in
series m = dnorm(ndx)*(y/P - (1-y)/(1-P))
print m
21.7 Using functions
The mle command allows you to estimate models that gretl does not provide natively: in some
cases, it may be a good idea to wrap up the mle block in a user-defined function (see Chapter13),
so as to extend gretl’s capabilities in a modular and flexible way.
As an example, we will take a simple case of a model that gretl does not yet provide natively:
the zero-inflated Poisson model, or ZIP for short.
4
In this model, we assume that we observe a
mixed population: for some individuals, the variable y
t
is (conditionally on a vector of exogenous
covariates x
t
)distributed as a Poisson random variate; for some others, y
t
is identically 0. The
trouble is, we don’t know which category a given individual belongs to.
For instance, suppose we have a sample of women, and the variable y
t
represents the number of
children that woman t has. There may be a certain proportion, , of women for whom y
t
0 with
certainty (maybe out of a personal choice, or due to physical impossibility). But there may be other
4
Theactual ZIP modelis in facta bitmoregeneralthan theone presented here. Thespecializedversiondiscussed in
this section was chosen for thesake ofsimplicity. For futher details, seeGreene(2003).
Converting pdf to ppt online - control application platform:C# Create PDF from PowerPoint Library to convert pptx, ppt to PDF in C#.net, ASP.NET MVC, WinForms, WPF
Online C# Tutorial for Creating PDF from Microsoft PowerPoint Presentation
www.rasteredge.com
Converting pdf to ppt online - control application platform:VB.NET Create PDF from PowerPoint Library to convert pptx, ppt to PDF in vb.net, ASP.NET MVC, WinForms, WPF
VB.NET Tutorial for Export PDF file from Microsoft Office PowerPoint
www.rasteredge.com
Chapter 21. Maximum likelihood estimation
191
women for whom y
t
0 just as a matter of chance — they haven’t happened to have any children
at the time of observation.
In formulae:
Py
t
kjx
t
 
d
t
1 
"
e
t
y
t
t
y
t
!
#
t
expx
t

d
t
(
1 for y
t
0
0 for y
t
>0
Writing a mle block for this model is not difficult:
mle ll = logprob
series xb = exp(b0 + b1 * x)
series d = (y=0)
series poiprob = exp(-xb) * xb^y / gamma(y+1)
series logprob = (alpha>0) && (alpha<1) ? \
log(alpha*d + (1-alpha)*poiprob) : NA
params alpha b0 b1
end mle -v
However, the code above has to be modified each time we change our specification by, say, adding
an explanatory variable. Using functions, we can simplify this task considerably and eventually be
able to write something easy like
list X = const x
zip(y, X)
Example 21.3: Zero-inflated Poisson Model — user-level function
/*
user-level function: estimate the model and print out
the results
*/
function void zip(series y, list X)
matrix coef_stde = zip_estimate(y, X)
printf "\nZero-inflated Poisson model:\n"
string parnames = "alpha,"
string parnames += varname(X)
modprint coef_stde parnames
end function
Let’ssee how this can be done. First we needto define a function called zip() that will take two ar-
guments: a dependent variable y and a list of explanatory variables X. An example of such function
can be seen in script21.3. By inspecting the function code, you can see that the actual estimation
does not happen here: rather, the zip() function merely uses the built-in modprint command to
print out the results coming from another user-written function, namely zip_estimate().
The function zip_estimate() is not meant to be executed directly; it just contains the number-
crunching part of the job, whose results are then picked up by the end function zip(). In turn,
zip_estimate() calls other user-written functions to perform other tasks. The whole set of “in-
ternal” functions is shown in the panel21.4.
control application platform:Online Convert PowerPoint to PDF file. Best free online export
Your PDF file is converted to look just the same as it does in your office software. Creating a PDF from PPTX/PPT has never been so easy! Easy converting!
www.rasteredge.com
control application platform:How to C#: Convert PDF, Excel, PPT to Word
How to C#: Convert PDF, Excel, PPT to Word. Online C# Tutorial for Converting PDF, MS-Excel, MS-PPT to Word. PDF, MS-Excel, MS-PPT to Word Conversion Overview.
www.rasteredge.com
Chapter 21. Maximum likelihood estimation
192
Example 21.4: Zero-inflatedPoisson Model — internal functions
/* compute log probabilities for the plain Poisson model */
function series ln_poi_prob(series y, list X, matrix beta)
series xb = lincomb(X, beta)
return -exp(xb) + y*xb - lngamma(y+1)
end function
/* compute log probabilities for the zero-inflated Poisson model */
function series ln_zip_prob(series y, list X, matrix beta, scalar p0)
# check if the probability is in [0,1]; otherwise, return NA
if (p0>1) || (p0<0)
series ret = NA
else
series ret = ln_poi_prob(y, X, beta) + ln(1-p0)
series ret = (y=0) ? ln(p0 + exp(ret)) : ret
endif
return ret
end function
/* do the actual estimation (silently) */
function matrix zip_estimate(series y, list X)
# initialize alpha to a "sensible" value: half the frequency
# of zeros in the sample
scalar alpha = mean(y=0)/2
# initialize the coeffs (we assume the first explanatory
# variable is the constant here)
matrix coef = zeros(nelem(X), 1)
coef[1] = mean(y) / (1-alpha)
# do the actual ML estimation
mle ll = ln_zip_prob(y, X, coef, alpha)
params alpha coef
end mle --hessian --quiet
return $coeff ~ $stderr
end function
control application platform:VB.NET PowerPoint: Convert & Render PPT into PDF Document
This VB.NET PowerPoint to PDF conversion tutorial will illustrate our effective PPT to PDF converting control SDK from following aspects.
www.rasteredge.com
control application platform:VB.NET PowerPoint: Complete PowerPoint Document Conversion in VB.
image or document formats, such as PDF, BMP, TIFF that can be converted from PPT document, please corresponding VB.NET guide for converting PowerPoint document
www.rasteredge.com
Chapter 21. Maximum likelihood estimation
193
All the functions shown in21.3 and21.4 can be stored in a separate inp file and executed once, at
the beginning of our job, by means ofthe include command. Assumingthe name ofthis script file
is zip_est.inp, the following is an example script which (a) includes the script file, (b) generates a
simulated dataset, and (c) performs the estimation of a ZIP model on the artificial data.
set echo off
set messages off
# include the user-written functions
include zip_est.inp
# generate the artificial data
nulldata 1000
set seed 732237
scalar truep = 0.2
scalar b0 = 0.2
scalar b1 = 0.5
series x = normal()
series y = (uniform()<truep) ? 0 : genpois(exp(b0 + b1*x))
list X = const x
# estimate the zero-inflated Poisson model
zip(y, X)
The results are as follows:
Zero-inflated Poisson model:
coefficient
std. error
z-stat
p-value
-------------------------------------------------------
alpha
0.203069
0.0238035
8.531
1.45e-17 ***
const
0.257014
0.0417129
6.161
7.21e-10 ***
x
0.466657
0.0321235
14.53
8.17e-48 ***
Afurther step may then be creating a function package for accessing your new zip() function via
gretl’s graphical interface. For details on how to do this, see section13.5.
21.8 Advanced use of mle: functions, analytical derivatives, algorithm choice
All the techniques decribed in the previous sections may be combined, and mle can be used for
solving non-standard estimation problems (provided, of course, that one chooses maximum likeli-
hood as the preferred inference method).
The strategy that, as of this writing, has proven most successful in designing scripts for this pur-
pose is:
 Modularize your code as much as possible.
 Use analytical derivatives whenever possible.
 Choose your optimization method wisely.
In the rest of this section, we will expand on the probit example of section21.5 to give the reader
an idea of what a “heavy-duty” application ofmle lookslike. Most ofthe code fragmentscome from
mle-advanced.inp, which is one of the sample scripts supplied with the standard installation of
gretl (see under File > Script files > Practice File).
control application platform:VB.NET PowerPoint: Customize PPT Document Rendering Options in VB.
to render and convert PPT slide to various formats, including PDF, BMP, TIFF, SVG, PNG, JPEG, GIF and JBIG2. In the process of converting PPT slide to any of
www.rasteredge.com
control application platform:VB.NET PowerPoint: Process & Manipulate PPT (.pptx) Slide(s)
control add-on can do PPT creating, loading controls, PDF document, image to pdf files and for capturing, viewing, processing, converting, compressing and
www.rasteredge.com
Chapter 21. Maximum likelihood estimation
194
BFGSwith and without analytical derivatives
The example in section21.5 can be made more general by using matrices and user-written func-
tions. Consider the following code fragment:
list X = const x1 x2 x3
matrix b = zeros(nelem(X),1)
mle logl = y*ln(P) + (1-y)*ln(1-P)
series ndx = lincomb(X, b)
series P = cnorm(ndx)
params b
end mle
In this context, the fact that the model we are estimating has four explanatory variables is totally
incidental: the code is written in such a way that we could change the content of the list X without
having to make any other modification. This was made possible by:
1. gathering the parameters to estimate into a single vector b rather than using separate scalars;
2. using the nelem() function to initialize b, so that its dimension is kept track ofautomatically;
3. using the lincomb() function to compute the index function.
Aparallel enhancement could be achieved in the case of analytically computed derivatives: since
bis now a vector, mle expects the argument to the deriv keyword to be a matrix, in which each
column is the partial derivative to the corresponding element of b. It is useful to re-write the score
for the i-th observation as
@‘
i
@
m
i
x
0
i
(21.8)
where m
i
is the “signed Mills’ ratio”, that is
m
i
y
i
x
0
i

Øx
0
i

1  y
i
x
0
i

1 Øx
0
i

;
which was computedin section21.5 via
series P = cnorm(ndx)
series m = dnorm(ndx)*(y/P - (1-y)/(1-P))
Here, we will code it in a somewhat terser way as
series m = y ? invmills(-ndx) : -invmills(ndx)
and make use of the conditional assignment operator and of the specialized function invmills()
for efficiency. Building the score matrix is now easily achieved via
mle logl = y*ln(P) + (1-y)*ln(1-P)
series ndx = lincomb(X, b)
series P = cnorm(ndx)
series m = y ? invmills(-ndx) : -invmills(ndx)
matrix mX = {X}
deriv b = mX .* {m}
end mle
in which the {} operator was used to turn series and lists into matrices (see chapter15). However,
proceeding in this way for more complex models than probit may imply inserting into the mle
control application platform:VB.NET PowerPoint: Convert PowerPoint to BMP Image with VB PPT
in VB class for rendering and converting PowerPoint presentations converters, such as VB.NET PDF Converter, Excel to the corresponding guide on C# PPT to BMP
www.rasteredge.com
control application platform:C# TIFF: Learn to Convert MS Word, Excel, and PPT to TIFF Image
doc.ConvertToDocument(DocumentType.TIFF, @"output.tif"); C# Demo for Converting PowerPoint to TIFF. Add references (Extra); Load your PPT (.pptx) document.
www.rasteredge.com
Chapter 21. Maximum likelihood estimation
195
block a long series of instructions; the example above merely happens to be short because the
score matrix for the probit model is very easy to write in matrix form.
Abetter solution is writing a user-level function to compute the score andusing that inside the mle
block, as in
function matrix score(matrix b, series y, list X)
series ndx = lincomb(X, b)
series m = y ? invmills(-ndx) : -invmills(ndx)
return {m} .* {X}
end function
[...]
mle logl = y*ln(P) + (1-y)*ln(1-P)
series ndx = lincomb(X, b)
series P = cnorm(ndx)
deriv b = score(b, y, X)
end mle
In this way, no matter how complex the computation of the score is, the mle block remains nicely
compact.
Newton’s method and the analytical Hessian
Since version 1.9.7, gretl offers the user the option of using Newton’s method for maximizing the
log-likelihood. In terms of the notation used in section21.1, the direction for updating the inital
parameter vector 
0
is given by
dg
0
   H
0
1
g
0
;
(21.9)
where H is the Hessian of the total loglikelihood computed at  and 0 < < 1 is a scalar called
the step length.
The above expression makes a few points clear:
1. At each step, it must be possible to compute not only the score g, but also its derivative
H;
2. the matrix H should be nonsingular;
3. it is assumed that for some positive value of , ‘
1
> ‘
0
; in other words, that going in
the direction dg
0
 leads upwards for some step length.
The strength of Newton’s method lies in the fact that, if the loglikelihood is globally concave,
then (21.9) enjoys certain optimality properties and the number of iterations required to reach the
maximum is often much smaller than it would be with other methods, such as BFGS. However, it
may have some disadvantages: for a start, the Hessian H may be difficult or very expensive to
compute; moreover, the loglikelihood may not be globally concave, so for some values of , the
matrix H is not negative definite or perhaps even singular. Those cases are handled by gretl’s
implementation of Newton’s algorithm by means of several heuristic techniques
5
,but a number of
adverse consequences may occur, which range from longer computation time for optimization to
non-convergence of the algorithm.
As a consequence, using Newton’s method is advisable only when the computation of the Hessian
is not too CPU-intensive and the nature of the estimator is such that it is known in advance that
5
The gist to it is that, if H is not negative definite, itis substituted by k dgH 1 k H, where k is a suitable
scalar; however, ifyou’re interested in the precise details, you’ll be much better off looking at the source code: the file
you’llwanttolookatislib/src/gretl_bfgs.c.
control application platform:How to C#: Convert Word, Excel and PPT to PDF
How to C#: Convert Word, Excel and PPT to PDF. Online C# Tutorial for Converting MS Office Word, Excel and PowerPoint to PDF. MS Office
www.rasteredge.com
Chapter 21. Maximum likelihood estimation
196
the loglikelihoodis globally concave. The probit models satisfies both requisites, so we will expand
the preceding example to illustrate howto use Newton’s method in gretl.
Afirst example may be given simply by issuing the command
set optimizer newton
before the mle block.
6
This will instruct gretl to use Newton’s method insteadofBFGS. Ifthe deriv
keyword is used, gretl will differentiate the score function numerically; otherwise, if the score
has to be computed itself numerically, gretl will calculate H by differentiating the loglikelihood
numerically twice. The latter solution, though, is generally to be avoided, as may be extremely
time-consuming and may yield imprecise results.
Amuch better option is to calculate the Hessian analytically andhave gretl use its true value rather
than a numerical approximation. In most cases, this is both much faster and numerically stable,
but of course comes at the price of having to differentiate the loglikelihood twice to respect with
the parameters and translate the resulting expressions into efficient hansl code.
Luckily, both tasks are relatively easy in the probit case: the matrix ofsecond derivatives of ‘
i
may
be written as
@
2
i
@@0
 m
i
m
i
x
0
i
x
i
x
0
i
so the total Hessian is
Xn
i1
@
2
i
@@0
 X
0
2
6
6
6
6
6
4
w
1
w
2
.
.
.
w
n
3
7
7
7
7
7
5
X
(21.10)
where w
i
m
i
m
i
x
0
i
.It can be shown that w
i
>0,so the Hessian is guaranteedto be negative
definite in all sensible cases and the conditions are ideal for applying Newton’s method.
Ahansl translation of equation (21.10) may look like
function void Hess(matrix *H, matrix b, series y, list X)
/* computes the negative Hessian for a Probit model */
series ndx = lincomb(X, b)
series m = y ? invmills(-ndx) : -invmills(ndx)
series w = m*(m+ndx)
matrix mX = {X}
H = (mX .* {w})’mX
end function
There are two characteristics worth noting of the function above. For a start, it doesn’t return
anything: the result of the computation is simply stored in the matrix pointed at by the first
argument ofthe function. Second, the result is not the Hessian proper, but rather its negative. This
function becomes usable from within an mle block by the keyword hessian. The syntax is
mle ...
...
hessian funcname(&mat_addr, ...)
end mle
In other words, the hessian keyword must be followed by the call to a function whose first argu-
ment is a matrix pointer which is supposed to be filled with the negative of the Hessian at .
Another feature worth noting is that gretl does not perform any numerical check on whether the
function computes the Hessian correctly or not. On the one hand, this means that you can trick
6
To gobackto BFGS, youuse set optimizer bfgs.
Chapter 21. Maximum likelihood estimation
197
mle into using alternatives to the Hessian and thereby implement other optimization methods. For
example, if you substitutein equation21.9 the Hessian Hwith the negative ofthe OPGmatrix  G
0
G,
as defined in (21.2), you get the so-called BHHH optimization method (see Berndtet al.(1974)).
Again, the sample file mle-advanced.inp provides an example. On the other hand, you may want
to perform a check ofyour analytically-computed H matrix versus a numerical approximation.
If you have a function that computes the score, this is relatively simple to do by using the fdjac
function, briefly described in section31.4, which computes a numerical approximation to a deriva-
tive. In practice, you need a function computing g as a row vector and then use fdjac to
differentiate it numerically with respect to . The result can then be compared to your analytically-
computed Hessian. The code fragment below shows an example of how this can be done in the
probit case:
function matrix totalscore(matrix *b, series y, list X)
/* computes the total score */
return sumc(score(b, y, X))
end function
function void check(matrix b, series y, list X)
/* compares the analytical Hessian to its numerical
approximation obtained via fdjac */
matrix aH
Hess(&aH, b, y, X) # stores the analytical Hessian into aH
matrix nH = fdjac(b, "totalscore(&b, y, X)")
nH = 0.5*(nH + nH’) # force symmetry
printf "Numerical Hessian\n%16.6f\n", nH
printf "Analytical Hessian (negative)\n%16.6f\n", aH
printf "Check (should be zero)\n%16.6f\n", nH + aH
end function
Chapter 22
GMM estimation
22.1 Introduction and terminology
The Generalized Method of Moments (GMM) is a very powerful and general estimation method,
which encompasses practically all the parametric estimation techniques used in econometrics. It
was introduced in Hansen (1982) and Hansen and Singleton (1982); an excellent and thorough
treatment is given in chapter 17 ofDavidsonandMacKinnon(1993).
The basic principle on which GMM is built is rather straightforward. Suppose we wish to estimate
ascalar parameter  based on a sample x
1
;x
2
;:::;x
T
. Let 
0
indicate the “true” value of . Theo-
retical considerations (either of statistical or economic nature) may suggest that a relationship like
the following holds:
E
x
t
g
0 a   
0
;
(22.1)
with g a continuous and invertible function. That is to say, there exists a function of the data
and the parameter, with the property that it has expectation zero ifand only ifit is evaluated at the
true parametervalue. For example, economic models with rational expectationsleadto expressions
like (22.1) quite naturally.
If the sampling model for the x
t
sis such that some version of the Law of Large Numbers holds,
then
¯
X
1
T
XT
t1
x
t
p
! g
0
;
hence, since g is invertible, the statistic
ˆ
g
1
¯
X
p
! 
0
;
so
ˆ
is a consistent estimator of . A different way to obtain the same outcome is to choose, as an
estimator of , the value that minimizes the objective function
F 
2
4
1
T
XT
t1
x
t
g
3
5
2
¯
X g
2
;
(22.2)
the minimum is trivially reached at
ˆ
 g 1
¯
X, since the expression in square brackets equals 0.
The above reasoning can be generalized as follows: suppose  is an n-vector and we have m
relations like
Ef
i
x
t
; 0 for i  1:::m;
(22.3)
where E is a conditional expectation on a set of p variables z
t
, called the instruments. In the
above simple example, m  1 and fx
t
;  x
t
g, and the only instrument used is z
t
1.
Then, it must also be true that
E
h
f
i
x
t
; z
j;t
i
E
h
f
i;j;t

i
0 for i 1:::m
and j  1:::p;
(22.4)
equation (22.4) is known as an orthogonality condition, or moment condition. The GMM estimator is
defined as the minimum of the quadratic form
F;W 
¯
f
0
W
¯
f;
(22.5)
198
Chapter 22. GMM estimation
199
where
¯
fis a 1 m p vector holding the average of the orthogonality conditions and W is some
symmetric, positive definite matrix, known as the weights matrix. A necessary condition for the
minimum to exist is the order condition n  m p.
The statistic
ˆ
 Argmin
F;W
(22.6)
is a consistent estimator of  whatever the choice of W. However, to achieve maximum asymp-
totic efficiency W must be proportional to the inverse of the long-run covariance matrix of the
orthogonality conditions; if W is not known, a consistent estimator will suffice.
These considerations lead to the following empirical strategy:
1. Choose a positive definite W andcompute the one-step GMM estimator
ˆ
1
.Customary choices
for W are I
mp
or I
m
Z
0
Z
1
.
2. Use
ˆ
1
to estimate Vf
i;j;t
 and use its inverse as the weights matrix. The resulting esti-
mator
ˆ
2
is called the two-step estimator.
3. Re-estimate Vf
i;j;t
 by means of
ˆ
2
and obtain
ˆ
3
; iterate until convergence. Asymp-
totically, these extra steps are unnecessary, since the two-step estimator is consistent and
efficient; however, the iterated estimator often hasbetter small-sample properties andshould
be independent of the choice of W made at step 1.
In the special case when the number of parameters n is equal to the total number of orthogonality
conditions m  p, the GMM estimator
ˆ
is the same for any choice of the weights matrix W, so the
first step is sufficient; in this case, the objective function is 0 at the minimum.
If, on the contrary, n < m  p, the second step (or successive iterations) is needed to achieve
efficiency, and the estimator so obtained can be very different, in finite samples, from the one-
step estimator. Moreover, the value of the objective function at the minimum, suitably scaled by
the number of observations, yields Hansen’s J statistic; this statistic can be interpreted as a test
statistic that has a 
2
distribution with m p n degrees of freedom under the null hypothesis of
correct specification. See Davidson and MacKinnon (1993, section 17.6) for details.
In the following sections we will show how these ideas are implemented in gretl through some
examples.
22.2 GMM as Method of Moments
This section draws from a kind contribution by Alecos Papadopoulos, whom we thank.
Avery simple illustration of GMM can be given by dropping the “G”, via an example of the time-
honored statistical technique known as“method ofmoments”; let’s see howto estimate the param-
eters of a gamma distribution, which we also used asan example for ML estimation in section21.2.
Assume that we have an i.i.d. sample of size T from a gamma distribution. The gamma density can
be parameterized in terms of the two parameters p (shape) and  (scale), both real and positive.
1
In order to estimate them by the method of moments, we need two moment conditions so that we
have two equations and two unknowns (in the GMM jargon, this amounts to exact identification).
The two relations we need are
Ex
i
p 
Vx
i
p 
2
1
In section21.2 we used a slightly different, perhaps more common, parametrization, employing   1=. We are
switching totheshape/scale parametrization here for thesake ofconvenience.
Documents you may be interested
Documents you may be interested