Chapter 31. Numerical methods
280
relevant data. On successful completion it returns an m  n matrix holding the Jacobian. For
example,
matrix Jac = fdjac(theta, SumOC(&theta, &X))
where we assume that SumOC is a user-definedfunction with the following structure:
function matrix SumOC (matrix *theta, matrix *X)
matrix V = ...
# do some computation
return V
end function
This may come in handy in several cases: for example,if you use BFGSmax to estimate a model, you
may wish to calculate a numerical approximation to the relevant Jacobian to construct a covariance
matrix for your estimates.
Another example is the delta method: if you have a consistent estimator of a vector of parameters
ˆ
, and a consistent estimate of its covariance matrix Ö, you may need to compute estimates for a
nonlinearcontinuous transformation   g. In this case, a standardresult in asymptotic theory
is that
8
<
:
ˆ
p
! 
p
T
ˆ
 
d
! N0;Ö
9
=
;
)
8
<
:
ˆ
 g
ˆ

p
!    g
p
T
ˆ
d
! N0;JÖJ0
9
=
;
where T is the sample size and J is the Jacobian
@gx
@x
x
.
Script31.6 exemplifies such a case: the example is taken fromGreene(2003), section 9.3.1. The
slight differencesbetween the results reportedin the original source and what gretl returns are due
to the fact that the Jacobian is computed numerically, rather than analytically as in the book.
On the subject of numerical versusanalytical derivatives, one may wonder what difference it makes
to use one method or another. Simply put, the answer is: analytical derivatives may be painful to
derive and to translate into code, but in most cases they are much faster than using fdjac; as a
consequence, if you need to use derivatives as part of an algorithm that requires iteration (such
as numerical optimization, or a Monte Carlo experiment), you’ll definitely want to use analytical
derivatives.
Analytical derivativesare also,in most cases, more precise than numerical ones, but this advantage
may or may not be negligible in practice depending on the practical details: the two fundamental
aspects to take in consideration are nonlinearity and machine precision.
As an example,consider the derivative of a highly nonlinear function such as the matrixinverse. In
order to keep the example simple, let’s focus on 22 matrices and define the function
function matrix vecinv(matrix x)
A = mshape(x,2,2)
return vec(inv(A))’
end function
which, given vecA, returns vecA 1. As is well known (see for exampleMagnusandNeudecker
(1988)),
@vecA
1
@vecA
 A
1
0
A
1
;
which is rather easy to code in hansl as
function matrix grad(matrix x)
iA = inv(mshape(x,2,2))
return -iA’ ** iA
end function
Adding pdf to powerpoint slide - control software 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
Adding pdf to powerpoint slide - control software 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 31. Numerical methods
281
Example 31.6: Delta Method
function matrix MPC(matrix *param, matrix *Y)
beta = param[2]
gamma = param[3]
y = Y[1]
return beta*gamma*y^(gamma-1)
end function
# William Greene, Econometric Analysis, 5e, Chapter 9
set echo off
set messages off
open greene5_1.gdt
# Use OLS to initialize the parameters
ols realcons 0 realdpi --quiet
genr a = $coeff(0)
genr b = $coeff(realdpi)
genr g = 1.0
# Run NLS with analytical derivatives
nls realcons = a + b * (realdpi^g)
deriv a = 1
deriv b = realdpi^g
deriv g = b * realdpi^g * log(realdpi)
end nls
matrix Y = realdpi[2000:4]
matrix theta = $coeff
matrix V = $vcv
mpc = MPC(&theta, &Y)
matrix Jac = fdjac(theta, MPC(&theta, &Y))
Sigma = qform(Jac, V)
printf "\nmpc = %g, std.err = %g\n", mpc, sqrt(Sigma)
scalar teststat = (mpc-1)/sqrt(Sigma)
printf "\nTest for MPC = 1: %g (p-value = %g)\n", \
teststat, pvalue(n,abs(teststat))
control software platform:VB.NET PowerPoint: Add Image to PowerPoint Document Slide/Page
insert or delete any certain PowerPoint slide without affecting on C#.NET PPT image adding library. powerful & profession imaging controls, PDF document, tiff
www.rasteredge.com
control software platform:VB.NET PowerPoint: Edit PowerPoint Slide; Insert, Add or Delete
To view C# code for adding, inserting or To view more VB.NET PowerPoint slide processing functions powerful & profession imaging controls, PDF document, image
www.rasteredge.com
Chapter 31. Numerical methods
282
Using the fdjac function to obtain the same result is even easier: you just invoke it like
fdjac(a, "vecinv(a)")
In order to see what the difference is, in terms of precision, between analytical and numerical
Jacobians, let’s start from A 
"
2 1
1 1
#
.The following code
a = {2; 1; 1; 1}
ia = vecinv(a)
ag = grad(a)
ng = fdjac(a, "vecinv(a)")
dg = ag - ng
print ag ng dg
gives
ag (4 x 4)
-1
1
1
-1
1
-2
-1
2
1
-1
-2
2
-1
2
2
-4
ng (4 x 4)
-1
1
1
-1
1
-2
-1
2
1
-1
-2
2
-1
2
2
-4
dg (4 x 4)
-3.3530e-08
-3.7251e-08
-3.7251e-08
-3.7255e-08
2.6079e-08
5.2150e-08
3.7251e-08
6.7060e-08
2.6079e-08
3.7251e-08
5.2150e-08
6.7060e-08
-2.2354e-08
-5.9600e-08
-5.9600e-08
-1.4902e-07
in which the analytically-computed derivative and its numerical approximation are essentially the
same. If, however, you set A 
"
1:0001 1
1
1
#
you end up evaluating the function at a point in
which the function itself is considerably more nonlinear since the matrix is much closer to being
singular. As a consequence, the numerical approximation becomes much less satisfactory:
ag (4 x 4)
-1.0000e+08
1.0000e+08
1.0000e+08
-1.0000e+08
1.0000e+08
-1.0001e+08
-1.0000e+08
1.0001e+08
1.0000e+08
-1.0000e+08
-1.0001e+08
1.0001e+08
-1.0000e+08
1.0001e+08
1.0001e+08
-1.0002e+08
ng (4 x 4)
-9.9985e+07
1.0001e+08
1.0001e+08
-9.9985e+07
9.9985e+07
-1.0002e+08
-1.0001e+08
9.9995e+07
9.9985e+07
-1.0001e+08
-1.0002e+08
9.9995e+07
-9.9985e+07
1.0002e+08
1.0002e+08
-1.0001e+08
control software platform:VB.NET PowerPoint: VB Code to Draw and Create Annotation on PPT
PDF, TIFF, MS Word and Excel). Most of end users would like to install and use Microsoft PowerPoint software and create PPT slide annotation through adding a
www.rasteredge.com
control software platform:C# PowerPoint - How to Process PowerPoint
Use the provided easy to call and write APIs programmed in C# class to develop user-defined PowerPoint slide adding and inserting projects.
www.rasteredge.com
Chapter 31. Numerical methods
283
dg (4 x 4)
-14899.
-14901.
-14901.
-14900.
14899.
14903.
14901.
14902.
14899.
14901.
14903.
14902.
-14899.
-14903.
-14903.
-14903.
Moreover, machine precision may have its impact: if you take A  0:00001
"
2 1
1 1
#
,the matrix
itself is not singular at all, but the order of magnitude of its elements is close enough to machine
precision to provoke problems:
ag (4 x 4)
-1.0000e+10
1.0000e+10
1.0000e+10
-1.0000e+10
1.0000e+10
-2.0000e+10
-1.0000e+10
2.0000e+10
1.0000e+10
-1.0000e+10
-2.0000e+10
2.0000e+10
-1.0000e+10
2.0000e+10
2.0000e+10
-4.0000e+10
ng (4 x 4)
-1.0000e+10
1.0000e+10
1.0000e+10
-1.0000e+10
1.0000e+10
-2.0000e+10
-1.0000e+10
2.0000e+10
1.0000e+10
-1.0000e+10
-2.0000e+10
2.0000e+10
-1.0000e+10
2.0000e+10
2.0000e+10
-4.0000e+10
dg (4 x 4)
-488.30
-390.60
-390.60
-195.33
634.79
781.21
390.60
585.98
634.79
488.26
683.55
585.98
-781.27
-976.52
-781.21
-1367.3
control software platform:VB.NET PowerPoint: VB Codes to Create Linear and 2D Barcodes on
PowerPoint PDF 417 barcode library is a mature and This PowerPoint ISSN barcode adding control is compatible ITF-14 barcode on any PowerPoint document slide
www.rasteredge.com
control software platform:VB.NET PowerPoint: Read, Edit and Process PPTX File
SDK into VB.NET class application by adding several compact well, like reading Excel in VB.NET, Reading PDF in VB Independent from Microsoft PowerPoint Product.
www.rasteredge.com
Chapter 32
Discrete and censored dependent variables
This chapter deals with models for dependent variables that are discrete or censored or otherwise
limited (as in event counts or durations, which must be positive) and that therefore call for estima-
tion methods other than the classical linear model. We discuss several estimators(mostly based on
the Maximum Likelihood principle), adding some details andexamples to complement the material
on these methods in the Gretl Command Reference.
32.1 Logit and probit models
It often happens that one wants to specify and estimate a model in which the dependent variable
is not continuous, but discrete. A typical example is a model in which the dependent variable is
the occupational status of an individual (1 = employed, 0 = unemployed). A convenient way of
formalizing this situation is to consider the variable y
i
as a Bernoulli random variable and analyze
its distribution conditional on the explanatory variables x
i
.That is,
y
i
(
1 P
i
0 1  P
i
(32.1)
where P
i
Py
i
1jx
i
is a given function of the explanatory variables x
i
.
In most cases, the function P
i
is a cumulative distribution function F, applied to a linear combi-
nation of the x
i
s. In the probit model, the normal cdf is used, while the logit model employs the
logistic function Ó. Therefore, we have
probit
P
i
Fz
i
 Øz
i
(32.2)
logit
P
i
Fz
i
 Óz
i

1
1e z
i
(32.3)
z
i
Xk
j1
x
ij
j
(32.4)
where z
i
is commonly known as the index function. Note that in this case the coefficients 
j
cannot
be interpreted as the partial derivatives of Ey
i
jx
i
with respect to x
ij
.However, for a given value
of x
i
it is possible to compute the vector of “slopes”, that is
slope
j
¯x 
@Fz
@x
j
z¯z
Gretl automatically computes the slopes, setting each explanatory variable at its sample mean.
Another, equivalent way of thinking about this model is in terms of an unobserved variable y
i
which can be described thus:
y
i
Xk
j1
x
ij
j
"
i
z
i
"
i
(32.5)
We observe y
i
1 whenever y
i
>0 and y
i
0 otherwise. If "
i
is assumed to be normal, then we
have the probit model. The logit model arises if we assume that the density function of "
i
is
"
i

@Ó"
i
@"
i
e
"
i
1e "
i
2
284
control software platform:C# PowerPoint: C# Guide to Add, Insert and Delete PPT Slide(s)
offer this C#.NET PowerPoint slide adding, inserting and guide for each PowerPoint slide processing operation & profession imaging controls, PDF document, tiff
www.rasteredge.com
control software platform:VB.NET PowerPoint: Sort and Reorder PowerPoint Slides by Using VB.
easily VB.NET PPT image adding and inserting clip art or screenshot to PowerPoint document slide at powerful & profession imaging controls, PDF document, image
www.rasteredge.com
Chapter 32. Discrete and censored dependent variables
285
Both the probit and logit model are estimated in gretl via maximum likelihood, where the log-
likelihood can be written as
L 
X
y
i
0
ln1 Fz
i

X
y
i
1
lnFz
i
;
(32.6)
which is always negative, since 0 < F < 1. Since the score equations do not have a closed form
solution, numerical optimization is used. However, in most cases this is totally transparent to the
user, since usually only a few iterations are needed to ensure convergence. The --verbose switch
can be used to track the maximization algorithm.
Example 32.1: Estimation of simple logit and probit models
open greene19_1
logit GRADE const GPA TUCE PSI
probit GRADE const GPA TUCE PSI
As an example, we reproduce the results given in chapter 21 ofGreene(2000), where the effective-
ness of a program for teaching economics is evaluated by the improvements of students’ grades.
Running the code in example32.1 gives the output reported in Table32.1; note that, for the probit
model, a conditional moment test on skewness and kurtosis is printed out automatically as a test
for normality.
In this context, the $uhat accessor function takes a special meaning: it returns generalized resid-
uals as defined inGourieroux,Monfort,RenaultandTrognon (1987), which can be interpreted as
unbiased estimators ofthe latent disturbances "
i
.These are defined as
u
i
8
<
:
y
i
ˆ
P
i
for the logit model
y
i
ˆz
i
Øˆz
i
1 y
i

ˆz
i
1 Øˆz
i
for the probit model
(32.7)
Among other uses, generalized residuals are often used for diagnostic purposes. For example, it is
very easy to set up an omitted variables test equivalent to the familiar LM test in the context of a
linear regression; example32.2 shows how to perform a variable addition test.
Example 32.2: Variableadditiontest in a probit model
open greene19_1
probit GRADE const GPA PSI
series u = $uhat
ols u const GPA PSI TUCE -q
printf "Variable addition test for TUCE:\n"
printf "Rsq * T = %g (p. val. = %g)\n", $trsq, pvalue(X,1,$trsq)
The perfect prediction problem
One curious characteristic of logit and probit models is that (quite paradoxically) estimation is not
feasible if a model fits the data perfectly; this is called the perfect prediction problem. The reason
control software platform:VB.NET PowerPoint: Extract & Collect PPT Slide(s) Using VB Sample
functions, like VB.NET PPT slide adding/removing, PPT read this VB.NET PowerPoint slide processing tutorial & profession imaging controls, PDF document, image
www.rasteredge.com
control software platform:VB.NET PowerPoint: PPTX to SVG Conversion; Render PPT to Vector
into VB.NET project by adding project reference PowerPoint files that end with .pptx file suffix can powerful & profession imaging controls, PDF document, tiff
www.rasteredge.com
Chapter 32. Discrete and censored dependent variables
286
Model 1: Logit estimates using the 32 observations 1-32
Dependent variable: GRADE
VARIABLE
COEFFICIENT
STDERROR
T STAT
SLOPE
(at mean)
const
-13.0213
4.93132
-2.641
GPA
2.82611
1.26294
2.238
0.533859
TUCE
0.0951577
0.141554
0.672
0.0179755
PSI
2.37869
1.06456
2.234
0.449339
Mean of GRADE = 0.344
Number of cases ’correctly predicted’ = 26 (81.2%)
f(beta’x) at mean of independent vars = 0.189
McFadden’s pseudo-R-squared = 0.374038
Log-likelihood = -12.8896
Likelihood ratio test: Chi-square(3) = 15.4042 (p-value 0.001502)
Akaike information criterion (AIC) = 33.7793
Schwarz Bayesian criterion (BIC) = 39.6422
Hannan-Quinn criterion (HQC) = 35.7227
Predicted
0
1
Actual 0
18
3
1
3
8
Model 2: Probit estimates using the 32 observations 1-32
Dependent variable: GRADE
VARIABLE
COEFFICIENT
STDERROR
T STAT
SLOPE
(at mean)
const
-7.45232
2.54247
-2.931
GPA
1.62581
0.693883
2.343
0.533347
TUCE
0.0517288
0.0838903
0.617
0.0169697
PSI
1.42633
0.595038
2.397
0.467908
Mean of GRADE = 0.344
Number of cases ’correctly predicted’ = 26 (81.2%)
f(beta’x) at mean of independent vars = 0.328
McFadden’s pseudo-R-squared = 0.377478
Log-likelihood = -12.8188
Likelihood ratio test: Chi-square(3) = 15.5459 (p-value 0.001405)
Akaike information criterion (AIC) = 33.6376
Schwarz Bayesian criterion (BIC) = 39.5006
Hannan-Quinn criterion (HQC) = 35.581
Predicted
0
1
Actual 0
18
3
1
3
8
Test for normality of residual -
Null hypothesis: error is normally distributed
Test statistic: Chi-square(2) = 3.61059
with p-value = 0.164426
Table 32.1: Examplelogit and probit
Chapter 32. Discrete and censored dependent variables
287
why this problem arisesis easy tosee by considering equation (32.6): if for some vector andscalar
kit’s the case that z
i
<k whenever y
i
0 and z
i
>k whenever y
i
1, the same thing is true
for any multiple of . Hence, L can be made arbitrarily close to 0 simply by choosing enormous
values for . As a consequence, the log-likelihood has no maximum, despite being bounded.
Gretl has a mechanism for preventing the algorithm from iterating endlessly in search of a non-
existent maximum. One sub-case of interest is when the perfect prediction problem arises because
of a single binary explanatory variable. In this case, the offending variable is dropped from the
model and estimation proceeds with the reduced specification. Nevertheless, it may happen that
no single “perfect classifier” exists among the regressors,in which case estimation is simply impos-
sible and the algorithm stops with an error. This behavior is triggered during the iteration process
if
maxz
i
i:y
i
0
<minz
i
i:y
i
1
If this happens, unless your model is trivially mis-specified (like predicting if a country is an oil
exporter on the basis of oil revenues), it is normally a small-sample problem: you probably just
don’t have enough data to estimate your model. You may want to drop some of your explanatory
variables.
This problem is well analyzed inStokes (2004); the results therein are replicated in the example
script murder_rates.inp.
32.2 Ordered response models
These models constitutea simplevariation on ordinary logit/probit models,andare usually applied
when the dependent variable is a discrete and ordered measurement—not simply binary, but on
an ordinal rather than an interval scale. For example, this sort of model may be applied when the
dependent variable is a qualitative assessment such as “Good”, “Average” and “Bad”.
In thegeneral case,consider an orderedresponsevariable,y, that can take on any ofthe J1 values
0;1;2;:::;J. We suppose, as before, that underlying the observed response is a latent variable,
y
X"  z "
Now define “cut points”, 
1
<
2
< <
J
,such that
y 0 if y

1
y 1 if 
1
<y 
2
.
.
.
y J if y >
J
For example, if the response takes on three values there will be two such cut points, 
1
and 
2
.
The probability that individual i exhibits response j, conditional on the characteristics x
i
,is then
given by
Py
i
jjx
i

8
>
>
<
>
>
:
Py

1
jx
i
 F
1
z
i
for j  0
P
j
<y

j1
jx
i
F
j1
z
i
 F
j
z
i
 for 0 <j <J
Py
>
J
jx
i
 1 F
J
z
i
for j  J
(32.8)
The unknown parameters 
j
are estimated jointly with the s via maximum likelihood. The ˆ
j
estimates are reported by gretl as cut1, cut2 and so on. For the probit variant, a conditional
moment test for normality constructed in the spirit ofChesherandIrish (1987) is also included.
Note that the 
j
parameters can be shifted arbitrarily by adding a constant to z
i
,so the model is
under-identified if there is some linear combination of the explanatory variables which is constant.
The most obvious case in which this occurs is when the model contains a constant term; for this
Chapter 32. Discrete and censored dependent variables
288
reason, gretl drops automatically the intercept if present. However, it may happen that the user in-
adventently specifies a list ofregressorsthat may be combinedin such a way to produce a constant
(for example,by using a full set of dummy variables for a discrete factor). If this happens, gretl will
also drop any offending regressors.
In order to apply these models in gretl, the dependent variable must either take on only non-
negative integer values, or be explicitly marked as discrete. (In case the variable has non-integer
values, it will be recoded internally.) Note that gretl does not provide a separate command for
ordered models: the logit and probit commands automatically estimate the ordered version if
the dependent variable is acceptable, but not binary.
Example32.3 reproduces the results presented in section 15.10 ofWooldridge(2002a). The ques-
tion of interest in this analysis is what difference it makes, to the allocation of assets in pension
funds, whether individual plan participants have a choice in the matter. The response variable is
an ordinal measure of the weight of stocks in the pension portfolio. Having reported the results
of estimation ofthe ordered model, Wooldridge illustrates the effect of the choice variable by ref-
erence to an “average” participant. The example script shows how one can compute this effect in
gretl.
After estimating orderedmodels, the $uhat accessor yields generalizedresiduals as in binary mod-
els; additionally, the $yhat accessor function returns ˆz
i
,so it is possible to compute an unbiased
estimator of the latent variable y
i
simply by adding the two together.
32.3 Multinomial logit
When the dependent variable is not binary and does not have a natural ordering, multinomial
models are used. Multinomial logit is supported in gretl via the --multinomial option to the
logit command. Simple models can also be handled via the mle command (see chapter21). We
give here an example of such a model. Let the dependent variable, y
i
, take on integer values
0;1;:::p. The probability that y
i
k is given by
Py
i
kjx
i

expx
i
k
P
p
j0
expx
i
j
For the purpose of identification one of the outcomes must be taken as the “baseline”; it is usually
assumed that 
0
0, in which case
Py
i
kjx
i

expx
i
k
1
P
p
j1
expx
i
j
and
Py
i
0jx
i

1
1
P
p
j1
expx
i
j
:
Example32.4 reproduces Table 15.2 inWooldridge(2002a), based on data on career choice from
Keane and Wolpin(1997).Thedependentvariableistheoccupationalstatusofanindividual(0=in
school;1 = not in school andnot working; 2 =working), andthe explanatory variables are education
and work experience (linear and square) plus a “black” binary variable. The full data set is a panel;
here the analysis is confined to a cross-section for 1987.
32.4 Bivariate probit
The bivariate probit model is simply a two-equation system in which each equation is a probit
model, but the two disturbance terms may not be independent. In formulae,
y
1;i
k
1
X
j1
x
ij
j
"
1;i
y
1;i
1 () y
1;i
>0
(32.9)
Chapter 32. Discrete and censored dependent variables
289
Example 32.3: Ordered probit model
/*
Replicate the results in Wooldridge, Econometric Analysis of Cross
Section and Panel Data, section 15.10, using pension-plan data from
Papke (AER, 1998).
The dependent variable, pctstck (percent stocks), codes the asset
allocation responses of "mostly bonds", "mixed" and "mostly stocks"
as {0, 50, 100}.
The independent variable of interest is "choice", a dummy indicating
whether individuals are able to choose their own asset allocations.
*/
open pension.gdt
# demographic characteristics of participant
list DEMOG = age educ female black married
# dummies coding for income level
list INCOME = finc25 finc35 finc50 finc75 finc100 finc101
# Papke’s OLS approach
ols pctstck const choice DEMOG INCOME wealth89 prftshr
# save the OLS choice coefficient
choice_ols = $coeff(choice)
# estimate ordered probit
probit pctstck choice DEMOG INCOME wealth89 prftshr
k = $ncoeff
matrix b = $coeff[1:k-2]
a1 = $coeff[k-1]
a2 = $coeff[k]
/*
Wooldridge illustrates the ’choice’ effect in the ordered probit
by reference to a single, non-black male aged 60, with 13.5 years
of education, income in the range $50K - $75K and wealth of $200K,
participating in a plan with profit sharing.
*/
matrix X = {60, 13.5, 0, 0, 0, 0, 0, 0, 1, 0, 0, 200, 1}
# with ’choice’ = 0
scalar Xb = (0 ~ X) * b
P0 = cdf(N, a1 - Xb)
P50 = cdf(N, a2 - Xb) - P0
P100 = 1 - cdf(N, a2 - Xb)
E0 = 50 * P50 + 100 * P100
# with ’choice’ = 1
Xb = (1 ~ X) * b
P0 = cdf(N, a1 - Xb)
P50 = cdf(N, a2 - Xb) - P0
P100 = 1 - cdf(N, a2 - Xb)
E1 = 50 * P50 + 100 * P100
printf "\nWith choice, E(y) = %.2f, without E(y) = %.2f\n", E1, E0
printf "Estimated choice effect via ML = %.2f (OLS = %.2f)\n", E1 - E0,
choice_ols
Documents you may be interested
Documents you may be interested