Chapter 20. Nonlinear least squares
180
Table 20.1: Nonlinearregression: the NIST tests
Analytical derivatives Numerical derivatives
Failures in 54 tests
4
5
Average iterations
32
127
Mean of min. correct figures,
8.120
6.980
parameters
Worst of min. correct figures,
4
3
parameters
Mean of min. correct figures,
8.000
5.673
standarderrors
Worst of min. correct figures,
5
2
standarderrors
Percent correct to at least 6 figures,
96.5
91.9
parameters
Percent correct to at least 6 figures,
97.7
77.3
standarderrors
Convert pdf to ppt online without email - application software cloud: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
Convert pdf to ppt online without email - application software cloud: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
21.1 Generic ML estimation with gretl
Maximum likelihood estimation is a cornerstone of modern inferential procedures. Gretl provides
away to implement this method for a wide range of estimation problems, by use of the mle com-
mand. We give here a few examples.
To give a foundation for the examples that follow, we start from a brief reminder on the basics
of ML estimation. Given a sample of size T, it is possible to define the density function
1
for the
whole sample, namely the joint distribution ofall the observations fY;,where Y
y
1
;:::;y
T
 
.
Its shape is determined by a k-vector of unknown parameters , which we assume is contained in
aset Ò, and which can be used to evaluate the probability of observing a sample with any given
characteristics.
After observing the data, the values Y are given, and this function can be evaluated for any legiti-
mate value of . In this case, we prefer to call it the likelihood function; the needfor another name
stems from the fact that this function works as a density when we use the y
t
sas arguments and 
as parameters, whereas in this context  is taken as the function’s argument, and the data Y only
have the role ofdetermining its shape.
In standardcases,this function hasa unique maximum. The location ofthe maximum isunaffected
if we consider the logarithm of the likelihood (or log-likelihood for short): this function will be
denoted as
‘  logfY;
The log-likelihood functions that gretl can handle are those where ‘ can be written as
‘ 
XT
t1
t

which is true in most cases of interest. The functions ‘
t
 are called the log-likelihood contribu-
tions.
Moreover, the location of the maximum is obviously determinedby the data Y. This means that the
value
ˆ
Y Argmax
2Ò
‘
(21.1)
is some function of the observed data (a statistic), which has the property, under mild conditions,
of being a consistent, asymptotically normal and asymptotically efficient estimator of .
Sometimes it is possible to write down explicitly the function
ˆ
Y; in general, it need not be so. In
these circumstances, the maximum can be found by means of numerical techniques. These often
rely on the fact that the log-likelihood is a smooth function of , and therefore on the maximum
its partial derivatives shouldall be 0. The gradient vector, or score vector, is a function that enjoys
many interesting statistical properties in its own right; it will be denoted here as g. It is a
1
Wearesupposing herethatourdataarearealizationofcontinuousrandom variables. Fordiscreterandomvariables,
everythingcontinuestoapplybyreferringtotheprobability functioninsteadofthedensity. Inbothcases,thedistribution
may be conditionalonsome exogenousvariables.
181
Chapter 21. Maximum likelihood estimation
182
k-vector with typical element
g
i
 
@‘
@
i
XT
t1
@‘
t

@
i
Gradient-basedmethods can be shortly illustrated as follows:
1. pick a point 
0
2Ò;
2. evaluate g
0
;
3. if g
0
is “small”, stop. Otherwise, compute a direction vector dg
0
;
4. evaluate 
1

0
dg
0
;
5. substitute 
0
with 
1
;
6. restart from 2.
Many algorithms of this kind exist; they basically differ from one another in the way they compute
the direction vector dg
0
, to ensure that ‘
1
> ‘
0
(so that we eventually end up on the
maximum).
The default method gretl uses to maximize the log-likelihood is a gradient-basedalgorithm known
as the BFGS (Broyden, Fletcher, Goldfarb and Shanno) method. This technique is used in most
econometric and statistical packages, as it is well-established and remarkably powerful. Clearly,
in order to make this technique operational, it must be possible to compute the vector g for
any value of . In some cases this vector can be written explicitly as a function of Y. If this is
not possible or too difficult the gradient may be evaluated numerically. The alternative Newton-
Raphson algorithm is also available, which is more effective under some circumstances but is also
more fragile; see section21.8 and chapter31 for details.
The choice of the starting value, 
0
,is crucial in some contexts and inconsequential in others. In
general, however, it is advisable to start the algorithm from “sensible” values whenever possible. If
aconsistent estimator is available, this is usually a safe and efficient choice: this ensures that in
large samples the starting point will be likely close to
ˆ
 and convergence can be achieved in few
iterations.
The maximum number of iterations allowed for the BFGS procedure, and the relative tolerance
for assessing convergence, can be adjusted using the set command: the relevant variables are
bfgs_maxiter (default value 500) and bfgs_toler (default value, the machine precision to the
power 3/4).
Covariance matrix and standard errors
By default the covariance matrix of the parameter estimates is based on the Outer Product of the
Gradient. That is,
d
Var
OPG
ˆ
 
G
0
ˆ
G
ˆ

1
(21.2)
where G
ˆ
 is the T k matrixof contributions to the gradient. Two other options are available. If
the --hessian flag is given, the covariance matrix is computed from a numerical approximation to
the Hessian at convergence. If the --robust option is selected, the quasi-ML “sandwich” estimator
is used:
d
Var
QML
ˆ
  H
ˆ

1
G
0
ˆ
G
ˆ
H
ˆ

1
where H denotes the numerical approximation to the Hessian.
Chapter 21. Maximum likelihood estimation
183
21.2 Gamma estimation
Supposewe havea sample ofT independent andidentically distributedobservationsfrom a Gamma
distribution. The density function for each observation x
t
is
fx
t

p
—p
x
p 1
t
exp x
t
(21.3)
The log-likelihood for the entire sample can be written as the logarithm of the joint density of all
the observations. Since these are independent and identical, the joint density is the product of the
individual densities, and hence its log is
‘;p 
XT
t1
log
"
p
—p
x
p 1
t
exp x
t
#
XT
t1
t
(21.4)
where
t
p logx
t
  p  logx
t
x
t
and   is the log of the gamma function. In order to estimate the parameters  and p via ML, we
need to maximize (21.4) with respect to them. The corresponding gretl code snippet is
scalar alpha = 1
scalar p = 1
mle logl =
p*ln(alpha * x) - lngamma(p) - ln(x) - alpha * x
params alpha p
end mle
The first two statements
alpha = 1
p = 1
are necessary to ensure that the variables alpha and p exist before the computation of logl is
attempted. Inside the mle block these variables (which could be either scalars, vectors or a com-
bination of the two — see below for an example) are identified as the parameters that should be
adjusted to maximize the likelihood via the params keyword. Their values will be changed by the
execution of the mle command; upon successful completion, they will be replaced by the ML esti-
mates. The starting value is 1 for both; this is arbitrary and does not matter much in this example
(more on this later).
The above code can be made more readable, and marginally more efficient, by defining a variable
to hold x
t
.This command can be embedded in the mle block as follows:
mle logl =
p*ln(ax) - lngamma(p) - ln(x) - ax
series ax = alpha*x
params alpha p
end mle
The variable ax is not added to the params list, of course, since it is just an auxiliary variable to
facilitate the calculations. You can insert as many such auxiliary lines as you require before the
params line, with the restriction that they must contain either (a) commands to generate series,
scalars or matrices or (b) print commands (which may be used to aid in debugging).
In a simple example like this, the choice of the starting values is almost inconsequential; the algo-
rithm is likely to converge no matter what the starting values are. However, consistent method-of-
moments estimators of p and can be simply recoveredfrom the sample mean m and variance V:
since it can be shown that
Ex
t
 p=
Vx
t
 p=
2
Chapter 21. Maximum likelihood estimation
184
it follows that the following estimators
¯
m=V
¯
p 
m
¯
are consistent, andtherefore suitable to be used asstartingpoint for the algorithm. The gretl script
code then becomes
scalar m = mean(x)
scalar alpha = m/var(x)
scalar p = m*alpha
mle logl =
p*ln(ax) - lngamma(p) - ln(x) - ax
series ax = alpha*x
params alpha p
end mle
Another thing to note is that sometimes parameters are constrained within certain boundaries: in
this case, for example, both  and p must be positive numbers. Gretl does not check for this: it
is the user’s responsibility to ensure that the function is always evaluated at an admissible point
in the parameter space during the iterative search for the maximum. An effective technique is to
define a variable for checking that the parameters are admissible and setting the log-likelihood as
undefined if the check fails. An example, which uses the conditional assignment operator, follows:
scalar m = mean(x)
scalar alpha = m/var(x)
scalar p = m*alpha
mle logl = check ? p*ln(ax) - lngamma(p) - ln(x) - ax : NA
series ax = alpha*x
scalar check = (alpha>0) && (p>0)
params alpha p
end mle
21.3 Stochastic frontier cost function
When modeling a cost function, it is sometimes worthwhile to incorporate explicitly into the sta-
tistical model the notion that firms may be inefficient, so that the observed cost deviates from the
theoretical figure not only because of unobserved heterogeneity between firms, but also because
two firms could be operating at a different efficiency level, despite being identical under all other
respects. In this case we may write
C
i
C
i
u
i
v
i
where C
i
issome variable cost indicator, C
i
is its“theoretical” value, u
i
isa zero-mean disturbance
term and v
i
is the inefficiency term, which is supposed to be nonnegative by its very nature.
Alinear specification for C
i
is often chosen. For example, the Cobb–Douglas cost function arises
when C
i
is a linear function of the logarithms of the input prices andthe output quantities.
The stochastic frontier model is a linear model of the form y
i
x
i
"
i
in which the error term
"
i
is the sum of u
i
and v
i
. A common postulate is that u
i
N0;2
u
and v
i
N0;2
v
. If
independence between u
i
and v
i
is also assumed, then it is possible to show that the density
function of "
i
has the form:
f"
i

s
2
Ø
"
i
1
"
i
(21.5)
where Ø and are, respectively,the distribution anddensity function ofthe standardnormal,

q
2
u

2
v
and  
u
v
.
Chapter 21. Maximum likelihood estimation
185
As a consequence, the log-likelihood for one observation takes the form (apart form an irrelevant
constant)
t
logØ
"
i
"
log 
"
2
i
22
#
Therefore, a Cobb–Douglas cost function with stochastic frontier is the model described by the
following equations:
logC
i
logC
i
"
i
logC
i
c
mX
j1
j
logy
ij
Xn
j1
j
logp
ij
"
i
u
i
v
i
u
i
N0;
2
u
v
i
N0;
2
v
In most cases, one wants to ensure that the homogeneity of the cost function with respect to
the prices holds by construction. Since this requirement is equivalent to
P
n
j1
j
1, the above
equation for C
i
can be rewritten as
logC
i
logp
in
c 
mX
j1
j
logy
ij
Xn
j2
j
logp
ij
logp
in
"
i
(21.6)
The above equation could be estimated by OLS, but it would suffer from two drawbacks: first,
the OLS estimator for the intercept c is inconsistent because the disturbance term has a non-zero
expected value; second, the OLS estimators for the other parameters are consistent, but inefficient
in view of the non-normality of "
i
. Both issues can be addressed by estimating (21.6) by maximum
likelihood. Nevertheless, OLS estimation is a quick and convenient way to provide starting values
for the MLE algorithm.
Example21.1 shows how to implement the model described so far. The banks91 file contains part
of the data used inLucchetti,PapiandZazzaro (2001).
The script in example21.1 is relatively easy to modify to show how one can use vectors (that is,
1-dimensional matrices) for storingthe parameters to optimize: example21.2 holds essentially the
same script in which the parametersof the cost function are stored together in a vector. Ofcourse,
this makes also possible to use variable lists and other refinements which make the code more
compact andreadable.
21.4 GARCH models
GARCH models are handled by gretl via a native function. However, it is instructive to see how they
can be estimated through the mle command.
2
The following equations provide the simplest example of a GARCH(1,1) model:
y
t
"
t
"
t
u
t

t
u
t
N0;1
h
t
!"
2
t 1
h
t 1
:
Since the variance of y
t
depends on past values, writing down the log-likelihood function is not
simply a matter of summing the log densities for individual observations. As is common in time
2
Thegigaddon,whichhandlesothervariantsofconditionally heteroskedasticmodels, usesmleasitsinternalengine.
Chapter 21. Maximum likelihood estimation
186
Example 21.1: Estimation of stochastic frontiercost function (with scalarparameters)
open banks91.gdt
# transformations
series cost = ln(VC)
series q1 = ln(Q1)
series q2 = ln(Q2)
series p1 = ln(P1)
series p2 = ln(P2)
series p3 = ln(P3)
# Cobb-Douglas cost function with homogeneity restrictions
# (for initialization)
genr rcost = cost - p1
genr rp2 = p2 - p1
genr rp3 = p3 - p1
ols rcost const q1 q2 rp2 rp3
# Cobb-Douglas cost function with homogeneity restrictions
# and inefficiency
scalar b0 = $coeff(const)
scalar b1 = $coeff(q1)
scalar b2 = $coeff(q2)
scalar b3 = $coeff(rp2)
scalar b4 = $coeff(rp3)
scalar su = 0.1
scalar sv = 0.1
mle logl = ln(cnorm(e*lambda/ss)) - (ln(ss) + 0.5*(e/ss)^2)
scalar ss = sqrt(su^2 + sv^2)
scalar lambda = su/sv
series e = rcost - b0*const - b1*q1 - b2*q2 - b3*rp2 - b4*rp3
params b0 b1 b2 b3 b4 su sv
end mle
Chapter 21. Maximum likelihood estimation
187
Example 21.2: Estimationof stochastic frontier cost function (withmatrix parameters)
open banks91.gdt
# transformations
series cost = ln(VC)
series q1 = ln(Q1)
series q2 = ln(Q2)
series p1 = ln(P1)
series p2 = ln(P2)
series p3 = ln(P3)
# Cobb-Douglas cost function with homogeneity restrictions
# (for initialization)
genr rcost = cost - p1
genr rp2 = p2 - p1
genr rp3 = p3 - p1
list X = const q1 q2 rp2 rp3
ols rcost X
X = const q1 q2 rp2 rp3
# Cobb-Douglas cost function with homogeneity restrictions
# and inefficiency
matrix b = $coeff
scalar su = 0.1
scalar sv = 0.1
mle logl = ln(cnorm(e*lambda/ss)) - (ln(ss) + 0.5*(e/ss)^2)
scalar ss = sqrt(su^2 + sv^2)
scalar lambda = su/sv
series e = rcost - lincomb(X, b)
params b su sv
end mle
Chapter 21. Maximum likelihood estimation
188
series models, y
t
cannot be considered independent of the other observations in our sample, and
consequently the density function for the whole sample (the joint density for all observations) is
not just the product of the marginal densities.
Maximum likelihood estimation, in these cases, is achieved by considering conditional densities, so
what we maximize is a conditional likelihood function. Ifwe define the information set at time t as
F
t
y
t
;y
t 1
;:::
 
;
then the density of y
t
conditional on F
t 1
is normal:
y
t
jF
t 1
N ;h
t
:
By means of the properties of conditional distributions, the joint density can be factorized as
follows
fy
t
;y
t 1
;::: 
2
4
YT
t1
fy
t
jF
t 1
3
5
fy
0
Ifwe treat y
0
as fixed,then the term fy
0
does not dependon the unknown parameters,andthere-
fore the conditional log-likelihood can then be written as the sum of the individual contributions
as
‘;!;; 
XT
t1
t
(21.7)
where
t
log
"
1
p
h
t
y
t
p
h
t
!#
 
1
2
"
logh
t

y
t
2
h
t
#
The following script shows a simple application of thistechnique,which uses the data file djclose;
it is one of the example dataset supplied with gretl and contains daily data from the Dow Jones
stock index.
open djclose
series y = 100*ldiff(djclose)
scalar mu = 0.0
scalar omega = 1
scalar alpha = 0.4
scalar beta = 0.0
mle ll = -0.5*(log(h) + (e^2)/h)
series e = y - mu
series h = var(y)
series h = omega + alpha*(e(-1))^2 + beta*h(-1)
params mu omega alpha beta
end mle
21.5 Analytical derivatives
Computation ofthe score vectoris essential forthe workingofthe BFGS method. In all the previous
examples, no explicit formula for the computation of the score was given, so the algorithm was fed
numerically evaluated gradients. Numerical computation of the score for the i-th parameter is
performed via a finite approximation of the derivative, namely
@‘
1
;:::;
n
@
i
‘
1
;:::;
i
h;:::;
n
 ‘
1
;:::;
i
h;:::;
n
2h
Chapter 21. Maximum likelihood estimation
189
where h is a small number.
In many situations, this is rather efficient and accurate. A better approximation to the true deriva-
tive may be obtained by forcing mle to use a technique known as Richardson Extrapolation, which
gives extremely precise results, but is considerably more CPU-intensive. Thisfeature may be turned
on by using the set command as in
set bfgs_richardson on
However, one might want to avoid the approximation and specify an exact function for the deriva-
tives. As an example, consider the following script:
nulldata 1000
genr x1 = normal()
genr x2 = normal()
genr x3 = normal()
genr ystar = x1 + x2 + x3 + normal()
genr y = (ystar > 0)
scalar b0 = 0
scalar b1 = 0
scalar b2 = 0
scalar b3 = 0
mle logl = y*ln(P) + (1-y)*ln(1-P)
series ndx = b0 + b1*x1 + b2*x2 + b3*x3
series P = cnorm(ndx)
params b0 b1 b2 b3
end mle --verbose
Here, 1000 data points are artificially generated for an ordinary probit model:
3
y
t
is a binary
variable, which takes the value 1 if y
t

1
x
1t

2
x
2t

3
x
3t
"
t
>0 and 0 otherwise. Therefore,
y
t
1 with probability Ø
1
x
1t

2
x
2t

3
x
3t

t
.The probability function for one observation
can be written as
Py
t

y
t
t
1  
t
1 y
t
Since the observations are independent and identically distributed, the log-likelihood is simply the
sum of the individual contributions. Hence
‘
XT
t1
y
t
log
t
1 y
t
log1  
t
The --verbose switch at the end of the end mle statement produces a detailed account of the
iterations done by the BFGS algorithm.
In thiscase, numerical differentiation worksrather well;nevertheless, computation ofthe analytical
score is straightforward, since the derivative
@‘
@
i
can be written as
@‘
@
i
@‘
@
t
@
t
@
i
via the chain rule, and it is easy to see that
@‘
@
t
y
t
t
1 y
t
1 
t
@
t
@
i

1
x
1t

2
x
2t

3
x
3t
 x
it
3
Again, gretldoesprovidea native probitcommand (seesection32.1), buta probit modelmakesfora nice example
here.
Documents you may be interested
Documents you may be interested