help mata ghk()
-------------------------------------------------------------------------------
Title
[M-5] ghk() -- Geweke-Hajivassiliou-Keane (GHK) multivariate normal
simulator
Syntax
S = ghk_init(real scalar npts)
(varies) ghk_init_method(S [, string scalar method])
(varies) ghk_init_start(S [, real scalar start])
(varies) ghk_init_pivot(S [, real scalar pivot])
(varies) ghk_init_antithetics(S [, real scalar anti])
real scalar ghk_query_npts(S)
real scalar ghk(S, real vector x, V)
real scalar ghk(S, real vector x, V, real rowvector dfdx, dfdv)
where S, if declared, should be declared
transmorphic S
and where method, optionally specified in ghk_init_method(), is
method description
------------------------------------------------------------
"halton" Halton sequences
"hammersley" Hammersley's variation of the Halton set
"random" pseudorandom uniforms
------------------------------------------------------------
Description
S = ghk_init(npts) constructs the transmorphic object S. Calls to
ghk_init_method(S, method), ghk_init_start(S, start), ghk_init_pivot(S,
pivot), and ghk_init_antithetics(S, anti) prior to calling ghk(S, ...)
allow you to modify the simulation algorithm through the object S.
ghk(S, x, V) returns a real scalar containing the simulated value of the
multivariate normal (MVN) distribution with variance-covariance V at the
point x. First, code S = ghk_init(npts) and then use ghk(S, ...) to
obtain the simulated value based on npts simulation points.
ghk(S, x, V, dfdx, dfdv) does the same thing but also returns the
first-order derivatives of the simulated probability with respect to x in
dfdx and the simulated probability derivatives with respect to vech(V) in
dfdv. See vech() in [M-5] vec() for details of the half-vectorized
operator.
The ghk_query_npts(S) function returns the number of simulation points,
the same value given in the construction of the transmorphic object S.
Remarks
Halton and Hammersley point sets are composed of deterministic sequences
on [0,1] and, for sets of dimension less than 10, generally have better
coverage than that of the uniform pseudorandom sequences.
Antithetic draws effectively double the number of points and reduce the
variability of the simulated probability. For draw u, the antithetic
draw is 1 - u. To use antithetic draws, call ghk_init_antithetic(S, 1)
prior to executing ghk().
By default, ghk() will pivot the wider intervals of integration (and
associated rows/columns of the covariance matrix) to the interior of the
multivariate integration. This improves the accuracy of the quadrature
estimate. When ghk() is used in a likelihood evaluator for [R] ml or
[M-5] optimize(), discontinuities may result in the computation of
numerical second-order derivatives using finite differencing (for the
Newton-Raphson optimize technique) when few simulation points are used,
resulting in a nonpositive-definite Hessian. To turn off the interval
pivoting, call ghk_init_pivot(S, 0) prior to executing ghk().
If you are using ghk() in a likelihood evaluator, be sure to use the same
sequence with each call to the likelihood evaluator. For a uniform
pseudorandom sequence, ghk_init_method("random"), set the seed of the
uniform random-variate generator -- see rseed() in [M-5] rseed() -- to
the same value with each call to the likelihood evaluator.
If you are using the Halton or Hammersley point sets, you will want to
keep the sequences going with each call to ghk() within one likelihood
evaluation. This can be done in one expression executed after each call
to ghk(): ghk_init_start(S, ghk_init_start(S)+ghk_query_npts(S)). With
each call to the likelihood evaluator, you will need to reset the
starting index to 1. This last point assumes that the transmorphic
object S is not recreated on each call to the likelihood evaluator.
Unlike [M-5] ghkfast_init(), the transmorphic object S created by
ghk_init() is inexpensive to create, so it is possible to re-create it
with each call to your likelihood evaluator instead of storing it as
external global and reusing the object with each likelihood evaluation.
Alternatively, the initialization function for optimize(),
optimize_init_arguments(), can be used.
Conformability
All initialization functions have 1 x 1 inputs and have 1 x 1 or void
outputs except
ghk_init(npts)
input:
npts: 1 x 1
output:
S: transmorphic
ghk_query_npts(S)
input:
S: transmorphic
output:
result: 1 x 1
ghk(S, x, V):
input:
S: transmorphic
x: 1 x m or m x 1
V: m x m (symmetric, positive definite)
output:
result: 1 x 1
ghk(S, x, V, dfdx, dfdv):
input:
S: transmorphic
x: 1 x m or m x 1
V: m x m (symmetric, positive definite)
output:
result: 1 x 1
dfdx: 1 x m
dfdv: 1 x m(m+1)/2
Diagnostics
The maximum dimension, m, is 20.
V must be symmetric and positive definite. ghk() will return a missing
value when V is not positive definite. When ghk() is used in an ml (or
optimize()) likelihood evaluator, return a missing likelihood to ml and
let ml take the appropriate action (i.e., step halving).
Source code
ghk.mata
Also see
Manual: [M-5] ghk()
Help: [M-5] ghkfast(), [M-5] halton(); [M-4] statistical
