\chapter{Maximum likelihood estimation}
\label{chap:mle}

\section{Generic ML estimation with gretl}
\label{sec:mle-intro}

Maximum likelihood estimation is a cornerstone of modern inferential
procedures. Gretl provides a way to implement this method for a wide
range of estimation problems, by use of the \texttt{mle} command. 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\footnote{We are
  supposing here that our data are a realization of continuous random
  variables. For discrete random variables, everything continues to
  apply by referring to the probability function instead of the
  density. In both cases, the distribution may be conditional on some
  exogenous variables.} for the whole sample, namely the joint
distribution of all the observations $f(\mathbf{Y} ; \theta)$, where
$\mathbf{Y} = \left\{ y_1, \ldots, y_T \right\}$.  Its shape is
determined by a $k$-vector of unknown parameters $\theta$, which we
assume is contained in a set $\Theta$, and which can be used to
evaluate the probability of observing a sample with any given
characteristics.

After observing the data, the values $\mathbf{Y}$ are given, and this
function can be evaluated for any legitimate value of $\theta$. In
this case, we prefer to call it the \emph{likelihood} function; the
need for another name stems from the fact that this function works as
a density when we use the $y_t$s as arguments and $\theta$ as
parameters, whereas in this context $\theta$ is taken as the
function's argument, and the data $\mathbf{Y}$ only have the role of
determining its shape.

In standard cases, this function has a unique maximum.  The location
of the maximum is unaffected if we consider the logarithm of the
likelihood (or log-likelihood for short): this function will be
denoted as
\[
  \LogLik(\theta) = \log  f(\mathbf{Y}; \theta)
\] 
The log-likelihood functions that gretl can handle are those where
$\LogLik(\theta)$ can be written as
\[
  \LogLik(\theta) = \sum_{t=1}^T \ell_t(\theta)
\] 
which is true in most cases of interest. The functions $\ell_t(\theta)$
are called the log-likelihood contributions.

Moreover, the location of the maximum is obviously determined by the
data $\mathbf{Y}$. This means that the value
\begin{equation}
  \label{eq:maxlik}
  \hat{\theta}(\mathbf{Y}) = \argmax_{\theta \in \Theta} \LogLik(\theta)
\end{equation}
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 $\theta$.

Sometimes it is possible to write down explicitly the function
$\hat{\theta}(\mathbf{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 $\theta$, and therefore on the maximum its partial
derivatives should all be 0.  The \textsl{gradient vector}, or
\textsl{score vector}, is a function that enjoys many interesting
statistical properties in its own right; it will be denoted here as
$\mathbf{g}(\theta)$.  It is a $k$-vector with typical element
\[
g_i(\theta) = \frac{\partial\LogLik(\theta)}{\partial\theta_i} 
  = \sum_{t=1}^T \frac{\partial\ell_t(\theta)}{\partial\theta_i}
\]

Gradient-based methods can be briefly illustrated as follows:

\begin{enumerate}
\item pick a point $\theta_0 \in \Theta$;
\item evaluate $\mathbf{g}(\theta_0)$;
\item if $\mathbf{g}(\theta_0)$ is ``small'', stop. Otherwise, compute
  a direction vector $d(\mathbf{g}(\theta_0))$;
\item evaluate $\theta_1 = \theta_0 + d(\mathbf{g}(\theta_0))$;
\item substitute $\theta_0$ with $\theta_1$;
\item restart from 2.
\end{enumerate}

Many algorithms of this kind exist; they basically differ from one
another in the way they compute the direction vector
$d(\mathbf{g}(\theta_0))$, to ensure that $\LogLik(\theta_1) >
\LogLik(\theta_0)$ (so that we eventually end up on the maximum).

The default method gretl uses to maximize the log-likelihood is a
gradient-based algorithm known as the \textbf{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
$\mathbf{g}(\theta)$ for any value of $\theta$. In some cases this
vector can be written explicitly as a function of $\mathbf{Y}$. If
this is not possible or too difficult the gradient may be evaluated
numerically. The alternative \textbf{Newton-Raphson} algorithm is also
available. This method is more effective under some circumstances but
is also more fragile; see section \ref{sec:mle-adv} and chapter
\ref{chap:numerical} for details.\footnote{Note that some of the
  statements made below (for example, regarding estimation of the
  covariance matrix) have to be modified when Newton's method is
  used.}

The choice of the starting value, $\theta_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 a consistent 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 $\hat{\theta}$ 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 \cmd{set} command: the relevant variables are
\verb+bfgs_maxiter+ (default value 500) and \verb+bfgs_toler+ (default
value, the machine precision to the power 3/4).

\section{Syntax}
\label{sec:ml-syntax}

ML estimation in gretl is supported by the \cmd{mle} command
block. This consists of an initial line holding the keyword
\texttt{mle} plus an equation for the loglikelihood; one or more
statements within the block (details below); and a trailer line to
close the block: \texttt{end mle}. Option flags may be appended to the
trailer line.

Listing~\ref{ex:ols-mle} gives a simple but complete example which serves
to illustrate the equivalence of MLE and OLS in the context of the
normal linear model.

\begin{script}[htbp]
  \scriptinfo{ols-mle}{OLS and MLE}
\begin{scode}
open data9-7
list X = const INCOME PRICE
ols QNC X
matrix b = $coeff
scalar s2 = $sigma^2
scalar l2pi = log(2*$pi)
scalar n = $nobs
mle lt = -0.5*l2pi -0.5*log(s2) - 1/(2*s2) * uhat^2
   series uhat = QNC - lincomb(X, b)
   s2 = sum(uhat^2)/n
   params b
end mle
\end{scode}
\end{script}

\subsection{Initial line of block}

If possible the given expression should evaluate to a series or vector
(contribution to the loglikelihood per observation). Failing that, it
must evaluate to a scalar (the total loglikelihood).  The identifier
on the left-hand side (\texttt{lt} in Listing~\ref{ex:ols-mle}) is up to
the user. If the variable in question is defined prior to the
\texttt{mle} block it can be referenced after ML estimation; otherwise
it is treated as a temporary variable and is destroyed after
estimation.

\subsection{Lines within the block}

These may take three forms:
\begin{enumerate}
\item ``Helper'' statements that calculate auxiliary quantities (in
  the example, \texttt{uhat} and \texttt{s2}). Such statements will be
  evaluated before the loglikelihood and then re-evaluated on each
  iteration.
\item Keyword plus parameter, as in ``\texttt{params b}'', which tells
  \texttt{mle} that the parameter to be adjusted to maximize the
  loglikelihood is the vector \texttt{b}. This sort of statement can
  also be used to specify analytical derivatives of the loglikelihood
  with respect to the parameters; see section~\ref{sec:anal-der} for
  discussion and examples.
\item Statements employing \texttt{print} or \texttt{printf} to track
  the progress of calculation, which can be useful for debugging.
\end{enumerate}

\subsection{Final line}

In the example above this merely terminates the block, but if one
wanted standard errors calculated via a numerical approximation to the
Hessian (for instance) one could substitute
\begin{code}
end mle --hessian
\end{code}
For a full listing of applicable options see the \cmd{mle} entry in
the \GCR.

\section{Covariance matrix and standard errors}
\label{sec:ml-vcv}

By default the covariance matrix of the parameter estimates is
based on the Outer Product of the Gradient (OPG).  That is,
\begin{equation}
  \label{eq:OPGmat}
  \widehat{\mbox{Var}}_{\mbox{\scriptsize OPG}}(\hat{\theta}) =
  \left(G'(\hat{\theta}) G(\hat{\theta}) \right)^{-1}
\end{equation}
where $G(\hat{\theta})$ is the $T \times k$ matrix of contributions to
the gradient.  Other options are available.  If the \option{hessian}
flag is given, the covariance matrix is computed from a numerical
approximation to the Hessian at convergence.  If the \option{robust}
option is given the quasi-ML ``sandwich'' estimator is used:
\[
\widehat{\mbox{Var}}_{\mbox{\scriptsize QML}}(\hat{\theta}) = H(\hat{\theta})^{-1}
  G'(\hat{\theta}) G(\hat{\theta}) H(\hat{\theta})^{-1}
\]
where $H$ denotes the numerical approximation to the Hessian. A
refinement here is that if the \texttt{hac} parameter is appended to
the \option{robust} option, as in
\begin{code}
end mle --robust=hac
\end{code}
the sandwich estimator is augmented in the manner of
\cite{newey-west87} to allow for serial correlation in the
gradient. (Note that this only makes sense for time-series data.) In
that case the details of the HAC estimator can be controlled via the
\cmd{set} command, as described in chapter~\ref{chap:robust_vcv}.

Cluster-robust estimation is also available: in order to activate it,
use the \option{cluster=}\emph{\texttt{clustvar}}, where \cmd{clustvar}
should be a discrete series. See section \ref{sec:vcv-cluster} for
more details.

Note, however, that if the log-likelihood function supplied by the
user just returns a scalar value---as opposed to a series or vector
holding per-observation contributions---then the OPG method is not
applicable and so the covariance matrix must be estimated via a
numerical approximation to the Hessian.

\section{Gamma estimation}
\label{sec:ml-gamma}

Suppose we have a sample of $T$ independent and identically
distributed observations from a Gamma distribution with shape
parameter $k$ and scale parameter $\theta$. The density
function for each observation $x_t$ is
\begin{equation}
  \label{eq:gammadens}
  f(x_t) = \frac{1}{\Gamma(k) \theta^k} \cdot \,x_t^{k-1} \,e^{-x_t/\theta}
\end{equation}
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
\begin{equation}
  \label{eq:gammaloglik}
  \LogLik(k, \theta) = \sum_{t=1}^T \log \left[ \frac{x_t^{k-1}
      e^{-x_t/\theta}}{\Gamma(k) \theta^k} \right] =
      \sum_{t=1}^T \ell_t
\end{equation}
where 
\[
  \ell_t = k \, \log (x_t/\theta) - \gamma(k) - \log x_t - x_t/\theta
\]
and $\gamma(\cdot)$ is the log of the gamma function.  In order to
estimate the parameters $k$ and $\theta$ via ML, we need to maximize
(\ref{eq:gammaloglik}) with respect to them. Here's a simple snippet
of gretl code to do the job.

\begin{code}
scalar k = 1
scalar theta = 1

mle logl =  k*ln(x/theta) - lngamma(k) - ln(x) - x/theta 
  params k theta
end mle 
\end{code}

The first two statements above are necessary to ensure that the
variables \texttt{k} and \texttt{theta} exist before the computation
of \texttt{logl} is attempted. Inside the \texttt{mle} block these
variables (which could be either scalars, vectors or a combination of
the two---see below for examples) are identified via the
\texttt{params} keyword as the parameters that should be adjusted to
maximize the likelihood.  Their values will be changed by the
execution of the \texttt{mle} command; upon successful completion,
they will be replaced by the ML estimates. We set the starting value
to 1 for both; this is arbitrary, but does not matter much in this
example (more on this later).

The code above can be made more readable, and marginally more
efficient, by defining a series to hold $x_t/\theta$, which we'll call
\texttt{y}. This command can be embedded in the \texttt{mle} block as
follows:
\begin{code}
mle logl = k*ln(y) - lngamma(k) - ln(x) - y
  series y = x/theta
  params k theta
end mle 
\end{code}
You can insert as many such auxiliary lines as you require before the
\texttt{params} line, with the restriction that they must contain
either (a) commands to generate series, scalars or matrices or (b)
\texttt{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 algorithm is likely to converge regardless
of the initialization. However, consistent method-of-moments
estimators of $k$ and $\theta$ can easily be recovered from the sample
mean $m$ and variance $V$. It can be shown that
\[
  E(x_t) = k\theta \qquad  V(x_t) = k\theta^2
\]
hence the following estimators 
\begin{eqnarray*}
  \bar{\theta} & = & V/m \\
  \bar{k} & = & m/\theta
\end{eqnarray*}
are consistent, and therefore suitable to be used as a starting point
for the algorithm.  The original initializers for $k$ and $\theta$
could then be replaced by
\begin{code}
scalar m = mean(x)
scalar theta = var(x)/m
scalar k = m/theta
\end{code}

Another thing to note is that sometimes parameters are constrained
within certain boundaries. In this case, for example, both Gamma
parameters must be positive numbers. Gretl can't check for this
automatically; it's 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 scalar variable that checks the validity of
the parameters, and set the log-likelihood as undefined (\texttt{NA})
if the check fails. To implement this, the mle block above could be
modified as follows:
\begin{code}
mle logl = check ? k*ln(y) - lngamma(k) - ln(x) - y : NA
  series y = x/theta
  scalar check = k > 0 && theta > 0
  params k theta
end mle 
\end{code}

For reference, Listing~\ref{ex:gamma-mle} presents a complete script that
generates artificial Gamma data and obtains ML estimates of the
parameters, employing the various refinements described above.

\begin{script}[htbp]
  \scriptinfo{gamma-mle}{ML estimation of Gamma parameters}
\begin{scode}
# create an empty data set with 200 observations
nulldata 200

# fix a random seed for replicability
set seed 1707138404

# generate a Gamma random variable x with shape k = 3 and scale theta = 2
series x = randgen(G, 3, 2)

# initialize estimates via sample moments
m = mean(x)
theta = var(x) / m
k = m / theta

mle logl = check ? k*ln(y) - lngamma(k) - ln(x) - y : NA
   series y = x/theta
   check = k > 0 && theta > 0
   params k theta
end mle
\end{scode}
  The \texttt{mle} output is:
\begin{scodebot}
Model 1: ML, using observations 1-200
logl = check ? k*ln(y) - lngamma(k) - ln(x) - y : NA
Standard errors based on Outer Products matrix

             estimate   std. error     z     p-value 
  ---------------------------------------------------
  k          3.28159     0.309764    10.59   3.18e-26 ***
  theta      1.86066     0.181566    10.25   1.21e-24 ***

Log-likelihood      -504.8970   Akaike criterion     1013.794
Schwarz criterion    1020.391   Hannan-Quinn         1016.463
\end{scodebot}
\end{script}

\section{Stochastic frontier cost function}
\label{sec:frontier}

\emph{%
Note: this section has the sole purpose of illustrating the \cmd{mle}
command. For the estimation of stochastic frontier cost or production
functions, you may want to use the \package{frontier} function package.}

When modeling a cost function, it is sometimes worthwhile to
incorporate explicitly into the statistical 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 in all other
respects. In this case we may write
\[
  C_i = C^*_i + u_i + v_i
\]
where $C_i$ is some variable cost indicator, $C_i^*$ is its
``theoretical'' value, $u_i$ is a zero-mean disturbance term and $v_i$
is the inefficiency term, which is supposed to be nonnegative by its
very nature. A linear 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 and the output
quantities.

The \emph{stochastic frontier} model is a linear model of the form
$y_i = x_i \beta + \varepsilon_i$ in which the error term
$\varepsilon_i$ is the sum of $u_i$ and $v_i$.  

A common postulate is that $u_i \sim N(0,\sigma_u^2)$ and
$v_i \sim \left|N(0,\sigma_v^2)\right|$. If independence between $u_i$
and $v_i$ is also assumed, then it is possible to show that the
density function of $\varepsilon_i$ has the form:
\begin{equation}
  \label{eq:frontdens}
  f(\varepsilon_i) = 
   \sqrt{\frac{2}{\pi}} 
   \Phi\left(\frac{\lambda \varepsilon_i}{\sigma}\right)
   \frac{1}{\sigma} \phi\left(\frac{\varepsilon_i}{\sigma}\right)
\end{equation}
where $\Phi(\cdot)$ and $\phi(\cdot)$ are, respectively, the distribution and density
function of the standard normal, $\sigma =
\sqrt{\sigma^2_u + \sigma^2_v}$ and $\lambda = \frac{\sigma_u}{\sigma_v}$.

As a consequence, the log-likelihood for one observation takes the
form (apart form an irrelevant constant)
\[
  \ell_t = 
  \log\Phi\left(\frac{\lambda \varepsilon_i}{\sigma}\right) -
  \left[ \log(\sigma) + \frac{\varepsilon_i^2}{2 \sigma^2} \right]
\]
Therefore, a Cobb--Douglas cost function with stochastic frontier is the
model described by the following equations: 
\begin{eqnarray*}
  \log C_i & = & \log C^*_i + \varepsilon_i \\
  \log C^*_i & = & c + \sum_{j=1}^m \beta_j \log y_{ij} + \sum_{j=1}^n \alpha_j \log p_{ij} \\
  \varepsilon_i & = & u_i + v_i \\
  u_i & \sim & N(0,\sigma_u^2) \\
  v_i & \sim & \left|N(0,\sigma_v^2)\right| 
\end{eqnarray*}

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 $\sum_{j=1}^n \alpha_j = 1$, the above
equation for $C^*_i$ can be rewritten as

\begin{equation}
  \label{eq:CobbDouglasFrontier}
  \log C_i - \log p_{in}  = c + \sum_{j=1}^m \beta_j \log y_{ij} +
  \sum_{j=2}^n \alpha_j (\log p_{ij} - \log p_{in})  + \varepsilon_i
\end{equation}

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
$\varepsilon_i$. Both issues can be addressed by estimating
(\ref{eq:CobbDouglasFrontier}) by maximum likelihood. Nevertheless,
OLS estimation is a quick and convenient way to provide starting
values for the MLE algorithm.

Listing \ref{ex:frontier1} shows how to implement the model
described so far. The \texttt{banks91} file contains part of the data
used in \citet*{lucchetti01}.

\begin{script}[htbp]
  \scriptinfo{frontier1}{Estimation of stochastic frontier cost
    function (with scalar parameters)}
\begin{scode}
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)
series rcost = cost - p1
series rp2 = p2 - p1
series 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
\end{scode}
\end{script}

The script in example \ref{ex:frontier1} is relatively easy to
modify to show how one can use vectors (that is, 1-dimensional
matrices) for storing the parameters to optimize: example
\ref{ex:frontier2} holds essentially the same script in which
the parameters of the cost function are stored together in a
vector. Of course, this makes also possible to use variable lists and
other refinements which make the code more compact and readable.

\begin{script}[htbp]
  \scriptinfo{frontier2}{Estimation of stochastic frontier cost
    function (with matrix parameters)}
\begin{scode}
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)
series rcost = cost - p1
series rp2 = p2 - p1
series 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
\end{scode}
\end{script}

\section{GARCH models}
\label{sec:garch}

GARCH models are handled by gretl via a native function. However, it is
instructive to see how they can be estimated through the \texttt{mle}
command.\footnote{The \package{gig} addon, which handles other variants
of conditionally heteroskedastic models, uses \texttt{mle} as its
internal engine.}

The following equations provide the simplest example of a GARCH(1,1)
model:
\begin{eqnarray*}
  y_t & = & \mu + \varepsilon_t \\
  \varepsilon_t & = & u_t \cdot \sigma_t \\
  u_t & \sim & N(0,1) \\
  h_t & = & \omega + \alpha \varepsilon^2_{t-1} + \beta h_{t-1}.
\end{eqnarray*}
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 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 \emph{conditional} densities, so what we maximize is a
conditional likelihood function. If we define the information set at
time $t$ as
\[
  F_t = \left\{ y_t, y_{t-1}, \ldots \right\} ,
\]
then the density of $y_t$ conditional on $F_{t-1}$ is normal:
\[
  y_t | F_{t-1} \sim N\left[ \mu, h_{t} \right].
\]

By means of the properties of conditional distributions, the joint
density can be factorized as follows
\[
  f(y_t, y_{t-1}, \ldots) = \left[ \prod_{t=1}^T f(y_t |F_{t-1})
  \right] \cdot f(y_0)
\]
If we treat $y_0$ as fixed, then the term $f(y_0)$ does not depend on
the unknown parameters, and therefore the conditional log-likelihood
can then be written as the sum of the individual contributions as
\begin{equation}
  \label{eq:garchloglik}
  \LogLik(\mu,\omega,\alpha,\beta) = \sum_{t=1}^T \ell_t
\end{equation}
where 
\[
  \ell_t = \log \left[ \frac{1}{\sqrt{h_t}} \phi\left( \frac{y_t - \mu}{\sqrt{h_t}}
    \right) \right] = 
    - \frac{1}{2} \left[ \log(h_t) + \frac{(y_t - \mu)^2}{h_t} \right]
\]

The following script shows a simple application of this technique,
which uses the data file \texttt{djclose}; it is one of the example
dataset supplied with gretl and contains daily data from the Dow Jones
stock index.

\begin{code}
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
\end{code}

\section{Analytical derivatives}
\label{sec:anal-der}

Computation of the score vector is essential for the working of the
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
\[
  \pder{\LogLik(\theta_1, \ldots, \theta_n)}{\theta_i} \simeq 
  \frac{\LogLik(\theta_1, \ldots, \theta_i + h, \ldots, \theta_n) -
    \LogLik(\theta_1, \ldots, \theta_i - h, \ldots, \theta_n)}{2h}
\]
where $h$ is a small number. 

In many situations, this is rather efficient and accurate. A better
approximation to the true derivative may be obtained by forcing
\cmd{mle} to use a technique known as \emph{Richardson Extrapolation},
which gives extremely precise results, but is considerably more
CPU-intensive. This feature may be turned on by using the \cmd{set}
command as in
\begin{code}
  set bfgs_richardson on
\end{code}

However, one might want to avoid the approximation and specify an
exact function for the derivatives. As an example, consider the
following script:
%
\begin{code}
nulldata 1000

series x1 = normal()
series x2 = normal()
series x3 = normal()

series ystar = x1 + x2 + x3 + normal()
series 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
\end{code}

Here, 1000 data points are artificially generated for an ordinary
probit model:\footnote{Again, gretl does provide a native
  \texttt{probit} command (see section \ref{sec:logit-probit}), but a
  probit model makes for a nice example here.} $y_t$ is a binary
variable, which takes the value 1 if $y_t^* = \beta_1 x_{1t} + \beta_2
x_{2t} + \beta_3 x_{3t} + \varepsilon_t > 0$ and 0 otherwise.
Therefore, $y_t = 1$ with probability $\Phi(\beta_1 x_{1t} + \beta_2
x_{2t} + \beta_3 x_{3t}) = \pi_t$.  The probability function for one
observation can be written as
\[
  P(y_t) = \pi_t^{y_t} ( 1 -\pi_t )^{1-y_t}
\]
Since the observations are independent and identically distributed,
the log-likelihood is simply the sum of the individual
contributions. Hence
\[
  \LogLik = \sum_{t=1}^T y_t \log(\pi_t) + (1 - y_t) \log(1 - \pi_t)
\]
The \option{verbose} switch at the end of the \texttt{end mle}
statement produces a detailed account of the iterations done by the
BFGS algorithm.

In this case, numerical differentiation works rather well;
nevertheless, computation of the analytical score is straightforward,
since the derivative $\pder{\LogLik}{\beta_i}$ can be written as
\[
  \pder{\LogLik}{\beta_i} = \pder{\LogLik}{\pi_t} \cdot \pder{\pi_t}{\beta_i}
\]
via the chain rule, and it is easy to see that
\begin{eqnarray*}
  \pder{\LogLik}{\pi_t} & = & \frac{y_t}{\pi_t} - \frac{1 - y_t}{1 -
    \pi_t} \\
  \pder{\pi_t}{\beta_i} & = & \phi(\beta_1 x_{1t} + \beta_2 x_{2t} +
  \beta_3 x_{3t}) \cdot x_{it}
\end{eqnarray*}

The \texttt{mle} block in the above script can therefore be modified
as follows:
%
\begin{code}
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
\end{code}

Note that the \texttt{params} statement has been replaced by a series
of \texttt{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.

\section{Debugging ML scripts}
\label{sec:mle-debug}

We have discussed above the main sorts of statements that are
permitted within an \texttt{mle} block, namely 
%
\begin{itemize}
\item auxiliary commands to generate helper variables;
\item \texttt{deriv} statements to specify the gradient with respect
  to each of the parameters; and
\item a \texttt{params} statement to identify the parameters in case
  analytical derivatives are not given.
\end{itemize}

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 \texttt{print} command.  The command word \texttt{print}
should be followed by the name of just one variable (a scalar, series
or matrix).

In the last example above a key variable named \texttt{m} was
generated, forming the basis for the analytical derivatives.  To track
the progress of this variable one could add a print statement within
the ML block, as in
%
\begin{code}
series m = dnorm(ndx)*(y/P - (1-y)/(1-P))
print m
\end{code}

\section{Using functions}
\label{sec:mle-func}

The \texttt{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 \texttt{mle} block in a user-defined function (see
Chapter \ref{chap:functions}), 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.\footnote{The actual ZIP model is in fact a bit more general
than the one presented here. The specialized version discussed in this
section was chosen for the sake of simplicity. For futher details, see
\cite{greene03}.} 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, $\alpha$, 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 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:
\begin{eqnarray*}
  P(y_t = k | x_t) & = & \alpha d_t + (1 - \alpha) 
  \left[e^{-\mu_t} \frac{\mu_t^{y_t}}{y_t!}\right] \\
    \mu_t & = & \exp(x_t \beta) \\
    d_t & = & 
    \left\{ 
      \begin{array}{ll} 
        1 & \mathrm{for} \quad y_t = 0 \\ 
        0 & \mathrm{for} \quad y_t > 0 
      \end{array}
    \right. 
\end{eqnarray*}

Writing a \texttt{mle} block for this model is not difficult:
\begin{code}
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
\end{code}

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
\begin{code}
list X = const x
zip(y, X)
\end{code}

\begin{script}[htbp]
  \scriptinfo{poisson-mle}{Zero-inflated Poisson Model -- user-level function}
\begin{scode}
/*
  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
\end{scode}
\end{script}

Let's see how this can be done.  First we need to define a function
called \texttt{zip()} that will take two arguments: a dependent
variable \texttt{y} and a list of explanatory variables \texttt{X}. An
example of such function can be seen in script~\ref{ex:poisson-mle}. By
inspecting the function code, you can see that the actual estimation
does not happen here: rather, the \texttt{zip()} function merely
uses the built-in \cmd{modprint} command to print out the results
coming from another user-written function, namely
\dtk{zip_estimate()}.

\begin{script}[htbp]
  \scriptinfo{poisson-functions}{Zero-inflated Poisson Model --- internal functions}
\begin{scode}
/* 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
\end{scode}
\end{script}

The function \dtk{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 \texttt{zip()}. In
turn, \dtk{zip_estimate()} calls other user-written functions to
perform other tasks. The whole set of ``internal'' functions is shown
in the panel \ref{ex:poisson-functions}.

All the functions shown in \ref{ex:poisson-mle} and \ref{ex:poisson-functions} can
be stored in a separate \texttt{inp} file and executed once, at the
beginning of our job, by means of the \texttt{include}
command.  Assuming the name of this script file is
\dtk{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.

\begin{code}
set verbose 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 : randgen(p, exp(b0 + b1*x))
list X = const x

# estimate the zero-inflated Poisson model
zip(y, X)
\end{code}

The results are as follows:

\begin{code}
Zero-inflated Poisson model:

             coefficient   std. error   z-stat   p-value 
  -------------------------------------------------------
  alpha       0.209738     0.0261746     8.013   1.12e-15 ***
  const       0.167847     0.0449693     3.732   0.0002   ***
  x           0.452390     0.0340836    13.27    3.32e-40 ***
\end{code}

A further step may then be creating a function package for accessing
your new \texttt{zip()} function via gretl's graphical interface. For
details on how to do this, see section \ref{sec:func-packages}.

\section{Advanced use of \texttt{mle}: functions, analytical
  derivatives, algorithm choice}
\label{sec:mle-adv}

All the techniques decribed in the previous sections may be combined,
and \cmd{mle} can be used for solving non-standard estimation problems
(provided, of course, that one chooses maximum likelihood as the
preferred inference method).

The strategy that, as of this writing, has proven most successful in
designing scripts for this purpose is:
\begin{itemize}
\item Modularize your code as much as possible.
\item Use analytical derivatives whenever possible.
\item Choose your optimization method wisely.
\end{itemize}

In the rest of this section, we will expand on the probit example of
section \ref{sec:anal-der} to give the reader an idea of what a
``heavy-duty'' application of \cmd{mle} looks like. Most of the code
fragments come from \verb|mle-advanced.inp|, which is one of the
sample scripts supplied with the standard installation of gretl
(see under \emph{File $>$ Script files $>$ Practice File}).

\subsection{BFGS with and without analytical derivatives}
\label{sec:mle-adv-bfgs}

The example in section \ref{sec:anal-der} can be made more general by
using matrices and user-written functions. Consider the following code
fragment:
\begin{code}
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
\end{code}

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 \texttt{X}
without having to make any other modification. This was made possible
by:
\begin{enumerate}
\item gathering the parameters to estimate into a single vector $b$
  rather than using separate scalars;
\item using the \texttt{nelem()} function to initialize $b$, so that
  its dimension is kept track of automatically;
\item using the \texttt{lincomb()} function to compute the index
  function.
\end{enumerate}

A parallel enhancement could be achieved in the case of analytically
computed derivatives: since $b$ is now a vector, \cmd{mle} expects the
argument to the \cmd{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
\begin{equation}
  \label{eq:mle-probscore}
  \pder{\LogLik_i}{\beta} = m_i \mathbf{x}_i'
\end{equation}
where $m_i$ is the ``signed Mills' ratio'', that is 
\[
m_i = y_i \frac{\phi(\mathbf{x}_i'\beta)}{\Phi(\mathbf{x}_i'\beta)} - 
(1-y_i) \frac{\phi(\mathbf{x}_i'\beta)}{1 - \Phi(\mathbf{x}_i'\beta)} ,
\]
which was computed in section \ref{sec:anal-der} via
\begin{code}
  series P = cnorm(ndx)
  series m = dnorm(ndx)*(y/P - (1-y)/(1-P))
\end{code}
Here, we will code it in a somewhat terser way as
\begin{code}
  series m = y ? invmills(-ndx) : -invmills(ndx)
\end{code}
and make use of the conditional assignment operator and of the
specialized function \cmd{invmills()} for efficiency.  Building the
score matrix is now easily achieved via
\begin{code}
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
\end{code}
in which the \verb|{}| operator was used to turn series and lists into
matrices (see chapter \ref{chap:matrices}). However, proceeding in
this way for more complex models than probit may imply inserting into
the \cmd{mle} 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.

A better solution is writing a user-level function to compute the
score and using that inside the \cmd{mle} block, as in
\begin{code}
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
\end{code}
In this way, no matter how complex the computation of the score is,
the \cmd{mle} block remains nicely compact.

\subsection{Newton's method and the analytical Hessian}
\label{sec:mle-adv-hessian}

As mentioned above, gretl offers the user the option of using
Newton's method for maximizing the log-likelihood. In terms of the
notation used in section \ref{sec:mle-intro}, the direction for
updating the inital parameter vector $\theta_0$ is given by
\begin{equation}
  \label{eq:mle-newton}
  d\left[\mathbf{g}(\theta_0)\right] = -\lambda
  \mathbf{H}(\theta_0)^{-1}\mathbf{g}(\theta_0) ,
\end{equation}
where $\mathbf{H}(\theta)$ is the Hessian of the total loglikelihood
computed at $\theta$ and $0 < \lambda < 1$ is a scalar called the
\emph{step length}.

The above expression makes a few points clear:
\begin{enumerate}
\item At each step, it must be possible to compute not only the score
  $\mathbf{g}(\theta)$, but also its derivative $\mathbf{H}(\theta)$;
\item the matrix $\mathbf{H}(\theta)$ should be nonsingular;
\item it is assumed that for some positive value of $\lambda$,
  $\LogLik(\theta_1) > \LogLik(\theta_0)$; in other words, that going
  in the direction $d\left[\mathbf{g}(\theta_0)\right]$ leads upwards
  for some step length.
\end{enumerate}

The strength of Newton's method lies in the fact that, if the
loglikelihood is globally concave, then \eqref{eq:mle-newton} 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 $\mathbf{H}(\theta)$ may be difficult or very
expensive to compute; moreover, the loglikelihood may not be globally
concave, so for some values of $\theta$, the matrix
$\mathbf{H}(\theta)$ 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\footnote{The gist to it is that, if $\mathbf{H}$ is not
  negative definite, it is substituted by $k \cdot
  \mathrm{dg}(\mathbf{H}) + (1-k) \cdot \mathbf{H}$, where $k$ is a
  suitable scalar; however, if you're interested in the precise
  details, you'll be much better off looking at the source code: the
  file you'll want to look at is \dtk{lib/src/gretl_bfgs.c}.}, but
a number of adverse consequences may occur, which range from
\emph{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 the
loglikelihood is globally concave. The probit models satisfies both
requisites, so we will expand the preceding example to illustrate how
to use Newton's method in gretl.

A first example may be given simply by issuing the command
\begin{code}
  set optimizer newton
\end{code}
before the \cmd{mle} block.\footnote{To go back to BFGS, you use
  \cmd{set optimizer bfgs}.} This will instruct gretl to use
Newton's method instead of BFGS. If the \cmd{deriv} keyword is used,
gretl will differentiate the score function numerically;
otherwise, if the score has to be computed itself numerically,
gretl will calculate $\mathbf{H}(\theta)$ 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.

A much better option is to calculate the Hessian analytically and have
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 \app{hansl} code.

Luckily, both tasks are relatively easy in the probit case: the matrix
of second derivatives of $\LogLik_i$ may be written as
\[
  \pder{{}^2\LogLik_i}{\beta \partial \beta'} = 
  - m_i \left( m_i + \mathbf{x}_i'\beta \right) \mathbf{x}_i \mathbf{x}_i'
\]
so the total Hessian is 
\begin{equation}
  \label{eq:mle-tothess}
  \sum_{i=1}^n \pder{{}^2\LogLik_i}{\beta \partial \beta'} = 
  - X' \left[
    \begin{array}{cccc}
      w_1 & & & \\
      & w_2 & & \\
      & & \ddots & \\
      & & & w_n
    \end{array}
  \right] X 
\end{equation}
where $w_i = m_i \left( m_i + \mathbf{x}_i'\beta \right)$. It can be
shown that $w_i > 0$, so the Hessian is guaranteed to be negative
definite in all sensible cases and the conditions are ideal for
applying Newton's method.

A \app{hansl} translation of equation \eqref{eq:mle-tothess} may look
like
\begin{code}
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
\end{code}

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 of the
function. Second, the result is not the Hessian proper, but rather its
negative. This function becomes usable from within an \cmd{mle} block
by the keyword \cmd{hessian}. The syntax is
\begin{code}
mle ...
    ...
    hessian funcname(&mat_addr, ...)
end mle
\end{code}
In other words, the \cmd{hessian} keyword must be followed by the call
to a function whose first argument is a matrix pointer which is
supposed to be filled with the \emph{negative} of the Hessian at
$\theta$.

We said above (section~\ref{sec:mle-intro}) that the covariance matrix
of the parameter estimates is by default estimated using the Outer
Product of the Gradient (so long as the log-likelihood function
returns the per-observation contributions). However, if you supply a
function that computes the Hessian then by default it is used in
estimating the covariance matrix. If you wish to impose use of OPG
instead, append the \option{opg} option to the end of the \cmd{mle}
block.

Note that gretl does not perform any numerical check on whether a
user-supplied function computes the Hessian correctly. On the one
hand, this means that you can trick \cmd{mle} into using alternatives
to the Hessian and thereby implement other optimization methods. For
example, if you substitute in equation \ref{eq:mle-newton} the Hessian
$\mathbf{H}$ with the negative of the OPG matrix
$-\mathbf{G}'\mathbf{G}$, as defined in \eqref{eq:OPGmat}, you get the
so-called BHHH optimization method (see \cite{bhhh74}). Again, the
sample file \verb|mle-advanced.inp| provides an example. On the other
hand, you may want to perform a check of your analytically-computed
$\mathbf{H}$ matrix versus a numerical approximation.

If you have a function that computes the score, this is relatively
simple to do by using the \cmd{fdjac} function, briefly described in
section \ref{sec:fdjac}, which computes a numerical approximation to a
derivative. In practice, you need a function computing
$\mathbf{g}(\theta)$ as a row vector and then use \cmd{fdjac} to
differentiate it numerically with respect to $\theta$. 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:
\begin{code}
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
\end{code}

\section{Estimating constrained models}
\label{sec:mle-constr}

In many cases, you may want to perform ML estimation of a model under
some kind of constraint. Mathematically, this amounts to maximizing the
log-likelihood $\LogLik(\theta)$ under some restriction. Assume
that the restriction can be represented as $g(\theta) = 0$, where
the function $g(\cdot)$ is differentiable. On paper, the most
straightforward way to accomplish this task is to set up a Lagrangean
\[
  \mathcal{L}(\theta) = \LogLik(\theta) + \lambda' g(\theta)
\]
and solve the first-order conditions that arise from differentiating
the Lagrangean with respect to $\theta$ and $\lambda$.

If an explicit solution can be found, then all is well; but in many
cases the resulting system of equations cannot be solved explicitly,
so that numerical optimisation in necessary. In such cases the
approach above is not particularly useful; a different strategy is
much more convenient.

The idea is to find an alternative parametrization---a means of
expressing the vector $\theta$ as a (differentiable) function of a
smaller set of parameters $\psi$. In other words, find a function
$h(\cdot)$ such that any admissible value of $\theta$ can be written
as $\theta = h(\psi)$ and $g[h(\psi)] = 0$ for any value of
$\psi$. Then maximization of the log-likelihood is simply a question
of operating on $\LogLik^*(\psi) = \LogLik[h(\psi)]$ using an ordinary
unconstrained numerical optimization routine.

Once the ML estimate $\hat{\psi}$ is available, it is easy to recover
the corresponding constrained vector $\hat{\theta} =
h(\hat{\psi})$. Computing the covariance matrix involves an extra
step, known as the \emph{delta method}: the asymptotic covariance
matrix of $\hat{\theta}$ can be computed as

\begin{equation}
  \label{eq:mle-constvar}
  V\left( \hat{\theta} \right) = J(\hat{\psi})' V\left( \hat{\psi}
  \right) J(\hat{\psi})
\end{equation}

where $J$ is the Jacobian matrix, holding the partial derivatives of
$h(\psi)$. It is recommended that the Jacobian matrix should be
computed analytically whenever possible, but as a fallback strategy,
numerical differentiation (available via the function \cmd{fdjac} ---
see section \ref{sec:fdjac}) is a viable alternative. Note that the
matrix produced by this method will be singular by construction.

The example reported in script \ref{ex:constrained-mle} is perhaps a
little contrived, but useful to elucidate the technique. Suppose we
wish to estimate mean and variance of an iid sample of Gaussian random
variables, under the constraint that
$V(x_i) = \sigma^2 = \exp[E(x_i)] = \exp(\mu)$. Of course the
unconstrained ML estimators $\hat{\mu} = \bar{X}$ and
$\hat{\sigma}^2 = n^{-1} \sum_i (x_i - \bar{X})^2$ are not guaranteed
to satisfy the constraints (in fact, the probability that they do is
0).

The Lagrangean in this case would be
\[
  \mathcal{L}(\theta) =
  K - \frac{n}{2} \log \sigma^2 - \frac{1}{2 \sigma^2} \sum_i (x_i -
  \mu)^2 + \lambda (e^\mu - \sigma^2)
\]
and finding an explicit solution by solving the first-order conditions
is not at all easy. Fortunately, numerical optimization becomes
straightforward by expressing the constrained parameters as
\[
  \theta = \left[ \mu, \sigma^2 \right]' = [\psi, \exp(\psi)]' =
  h(\psi) ;
\]
after maximizing the log-likelihood, the covariance matrix for
$\theta$ can be recovered by computing the Jacobian as
\[
  J(\psi) =
  \left[ \begin{array}{cc}
      \frac{d \mu}{d \psi} &
           \frac{d \sigma^2}{d \psi} 
  \end{array}\right] =          
  \left[ \begin{array}{cc}
      1 & \exp(\psi)
  \end{array}\right]           
\]
and applying formula \eqref{eq:mle-constvar}.

\begin{script}[htbp]
  \scriptinfo{constrained-mle}{Example of ML estimation of a model under constraints}
  \begin{scode}
set verbose off
set seed 7120

function matrix h(matrix psi)
    ret = psi[1] | exp(psi[1])
    return ret
end function

function matrix anJacob(matrix psi)
    # the derivative of h
    return 1 ~ exp(psi[1])
end function

nulldata 1000
# generate artificial data from a N(1, e) distribution
series x = 1 + normal() * exp(0.5)

# show that the unconstrained estimates don't satisfy the restriction

scalar muhat = mean(x)
scalar s2hat = sst(x)/$nobs
printf "unconstrained estimates:  mean = %g, variance = %g\n", muhat, s2hat
printf "check:  vhat - exp(muhat) = %g\n\n", s2hat - exp(muhat)

# now estimate under the constraint exp(mean) = variance

psi = {1}
mle loglik = -0.5*log(2*$pi) - 0.5*log(s2) - 0.5*(x-m)^2/s2
    matrix par = h(psi)
    scalar m = par[1]
    scalar s2 = par[2]
    params psi
end mle

# now map psi to the constrained parametrisation
matrix par = h(psi)

# show that now the constraint holds
printf "check:  vhat - exp(muhat) = %g\n\n", par[2] - exp(par[1])

# take care of the covariance matrix
matrix vpar = qform(anJacob(psi)', $vcv)
# alternatively, one could use the numerical Jacobian, as in
# matrix vpar = qform(fdjac(psi, "h(psi)"), $vcv)

# finally, print out the constrained parameters via "modprint"
matrix cs = par ~ sqrt(diag(vpar))
modprint cs "mean variance"
\end{scode}
\end{script}

Running the example script should produce the following output:

\begin{code}
unconstrained estimates:  mean = 1.00314, variance = 2.8903
check:  vhat - exp(muhat) = 0.163481

Model 1: ML, using observations 1-1000
loglik = -0.5*log(2*$pi) - 0.5*log(s2) - 0.5*(x-m)^2/s2
Standard errors based on Outer Products matrix

             estimate   std. error     z      p-value 
  ----------------------------------------------------
  psi[1]     1.03763    0.0357311    29.04   2.07e-185 ***

Log-likelihood      -1949.972   Akaike criterion     3901.943
Schwarz criterion    3906.851   Hannan-Quinn         3903.808

check:  vhat - exp(muhat) = 0

             coefficient   std. error     z      p-value 
  -------------------------------------------------------
  mean         1.03763     0.0357311    29.04   2.07e-185 ***
  variance     2.82251     0.100851     27.99   2.35e-172 ***
\end{code}
%$

\section{Handling non-convergence gracefully}
\label{sec:mle-nonconv}

If the numerical aspects of the estimation procedure are complex, it
is possible that \cmd{mle} fails to find the maximum within the number
of iterations stipulated via the \verb|bfgs_maxiter| state variable
(which defaults to 500).

In these cases, \cmd{mle} will exit with error and it's up to the user
to handle the situation appropriately. For example, it is possible
that \cmd{mle} is used inside a loop and you don't want the loop to
stop in case convergence is not achieved. The \cmd{catch} command
modifier (see also the \GCR) is an excellent tool for this purpose.

The example provided in listing \ref{ex:catch-mle} illustrates the
usage of \cmd{catch} in an artificially simple context: we use the
\cmd{mle} command for estimating mean and variance of a Gaussian rv
(of course you don't need the \cmd{mle} apparatus for this, but it
makes for a nice example). The gist of the example is using the
\verb|set bfgs_maxiter| command to force \cmd{mle} to abort after a
very small number of iterations, so that you can have an idea on how
to use the \cmd{catch} modifier and the associated \dollar{error}
accessor to handle the situation.

You may want to increase the maximum number if BFGS iterations in the
example to check what happens if the algorithm is allowed to
converge. Note that, upon successful completion of \cmd{mle}, a bundle
named \dollar{model} is available, containing several quantities that
may be of your interest, including the total number of function
evaluations.

\begin{script}[htbp]
  \scriptinfo{catch-mle}{Handling non-convergence via \cmd{catch}}
\begin{scode}
set verbose off
nulldata 200
set seed 8118

# generate simulated data from a N(3,4) variate
series x = normal(3,2)

# set starting values
scalar m = 0
scalar s2 = 1

# set iteration limit to a ridiculously low value
set bfgs_maxiter 10 

# perform ML estimation; note the "catch" modifier
catch mle loglik = -0.5* (log(2*$pi) + log(s2) + e2/s2)
    series e2 = (x - m)^2
    params m s2
end mle --quiet

# grab the error and proceed as needed
err = $error
if err
    printf "Not converged! (m = %g, s2 = %g)\n", m, s2
else
    printf "Converged after %d iterations\n", $model.grcount
    cs = $coeff ~ sqrt(diag($vcv))
    pn = "m s2"
    modprint cs pn
endif
\end{scode}
\end{script}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "gretl-guide"
%%% End: 
