mvProbitLogLik: Log Likelihood Values for Multivariate Probit Models

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/mvProbitLogLik.R

Description

Function mvProbitLogLik calculates log likelihood values of multivariate probit models.

The logLik model returns or calculates log likelihood values of multivariate probit models estimated by mvProbit.

Usage

1
2
3
4
5
6
7
8
mvProbitLogLik( formula, coef, sigma = NULL, data,
   algorithm = "GHK", nGHK = 1000, 
   returnGrad = oneSidedGrad, oneSidedGrad = FALSE, eps = 1e-6, 
   random.seed = 123, ... )

## S3 method for class 'mvProbit'
logLik( object, coef = NULL, data = NULL, 
   algorithm = NULL, nGHK = NULL, random.seed = NULL, ... )

Arguments

formula

a "formula": a symbolic description of the model (currently, all binary outcome variables must have the same regressors).

coef

a numeric vector of the model coefficients; if argument sigma is not specified, this vector must also include the correlation coefficients; the order of elements is explained in the section “details”.

sigma

optional argument for specifying the covariance/correlation matrix of the residuals (must be symmetric and have ones on its diagonal); if this argument is not specified, the correlation coefficients must be specified by argument coef.

data

a data.frame containing the data.

algorithm

algorithm for computing integrals of the multivariate normal distribution, either function GenzBretz(), Miwa(), or TVPACK() (see documentation of pmvnorm) or character string "GHK" (see documentation of ghkvec).

nGHK

numeric value specifying the number of simulation draws of the GHK algorithm for computing integrals of the multivariate normal distribution.

returnGrad

logical. If TRUE, the returned object has an attribute "gradient", which is a matrix and provides the gradients of the log-likelihood function with respect to all parameters (coef and upper triangle of sigma) evaluated at each observation and obtained by (two-sided) numeric finite-difference differentiation.

oneSidedGrad

logical. If TRUE, attribute "gradient" of the returned object is obtained by one-sided numeric finite-difference differentiation.

eps

numeric. The step size for the numeric finite-distance differentiation.

random.seed

an integer used to seed R's random number generator; this is to ensure replicability when computing (cumulative) probabilities of the multivariate normal distribution which is required to calculate the log likelihood values; set.seed( random.seed ) is called each time before a (cumulative) probability of the multivariate normal distribution is computed; defaults to 123.

object

an object of class "mvProbit" (returned by mvProbit.

...

additional arguments are passed to pmvnorm when calculating conditional expectations.

Details

If the logLik method is called with object as the only argument, it returns the log-likelihood value found in the maximum likelihood estimation. If any other argument is not NULL, the logLik method calculates the log-likelihood value by calling mvProbitLogLik. All arguments that are NULL, are taken from argument object.

If the model has n dependent variables (equations) and k explanatory variables in each equation, the order of the model coefficients in argument coef must be as follows: b_{1,1}, ..., b_{1,k}, b_{2,1}, ..., b_{2,k}, ..., b_{n,1}, ..., b_{n,k}, where b_{i,j} is the coefficient of the jth explanatory variable in the ith equation. If argument sigma is not specified, argument coef must additionally include following elements: R_{1,2}, R_{1,3}, R_{1,4}, ..., R_{1,n}, R_{2,3}, R_{2,4}, ..., R_{2,n}, ..., R_{n-1,n}, where R_{i,j} is the correlation coefficient corresponding to the ith and jth equation.

The ‘state’ (or ‘seed’) of R's random number generator is saved at the beginning of the mvProbitLogLik function and restored at the end of this function so that this function does not affect the generation of random numbers outside this function although the random seed is set to argument random.seed and the calculation of the (cumulative) multivariate normal distribution uses random numbers.

Value

mvProbitLogLik returns a vector containing the log likelihood values for each observation.

If argument returnGrad is TRUE, the vector returned by mvProbitLogLik has an attribute "gradient", which is a matrix and provides the gradients of the log-likelihood function with respect to all parameters (coef and upper triangle of sigma) evaluated at each observation and obtained by numeric finite-difference differentiation.

The logLik method returns the total log likelihood value (sum over all observations) with attribute df equal to the number of estimated parameters (model coefficients and correlation coefficients).

Author(s)

Arne Henningsen

References

Greene, W.H. (1996): Marginal Effects in the Bivariate Probit Model, NYU Working Paper No. EC-96-11. Available at http://ssrn.com/abstract=1293106.

See Also

mvProbit, mvProbitMargEff, probit, glm

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
## generate a simulated data set
set.seed( 123 )
# number of observations
nObs <- 10

# generate explanatory variables
xMat <- cbind( 
   const = rep( 1, nObs ),
   x1 = as.numeric( rnorm( nObs ) > 0 ),
   x2 = as.numeric( rnorm( nObs ) > 0 ),
   x3 = rnorm( nObs ),
   x4 = rnorm( nObs ) )

# model coefficients
beta <- cbind( c(  0.8,  1.2, -1.0,  1.4, -0.8 ),
               c( -0.6,  1.0,  0.6, -1.2, -1.6 ),
               c(  0.5, -0.6, -0.7,  1.1,  1.2 ) )

# covariance matrix of error terms
library( miscTools )
sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )

# generate dependent variables
yMatLin <- xMat %*% beta 
yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
colnames( yMat ) <- paste( "y", 1:3, sep = "" )

# log likelihood values
myData <- as.data.frame( cbind( xMat, yMat ) )
logLikVal <- mvProbitLogLik( cbind( y1, y2, y3 ) ~ x1 + x2 + x3 + x4, 
   coef = c( beta ), sigma = sigma, data = myData )
print( logLikVal )

mvProbit documentation built on May 30, 2017, 5:47 a.m.