Programming and Numerical Analysis: User Written Programs and Estimators

LIMDEP and NLOGIT’s components are designed to work together in a fully integrated package. The estimation programs, matrix and scalar calculators, variable transformation commands and programming tools can all be used together to design estimation and hypothesis testing routines. The following example shows how one might use a number of these features to fit a constrained model then test the constraint using a likelihood ratio and a Lagrange multiplier test.

Example: Poisson Model with a Fixed Value Restriction

This program estimates a Poisson regression model with one of the coefficients fixed at 1.0 and tests the restriction in three ways (the standard trinity of tests, Wald, LM and LR). This is an iterative procedure that uses Newton’s method. (There is actually a single line, ‘hardwired’ command that does the same thing.)

?===================================================================
? Set up the matrix procedure.  We strip off the last element of 
? b to obtain the vector beta. delta is the update vector, 
? initialized at 0, so the first iteration uses the starting values.
?===================================================================
NAMELIST	; x 	= one,a,c,d,e,c67,c72,c77,p77 $
POISSON	; Lhs 	= num ; Rhs = x,logmth $
MATRIX	; beta	= B(1:9) ; delta = [9|0] $   
CALC	; conv	= 1 ; loglu = logl $
?===================================================================
? This is the iterative procedure to fit the constrained estimator.
? 1.  Update beta.
? 2.  Compute expected values and residuals, imposing constraint
?     on slope on log(months)=1.
? 3.  Iterate using Newton's method
?===================================================================
PROCEDURE
MATRIX	; beta	= beta - delta 	$ update coefficient vector
CREATE	; index	= x'beta + logmth 
  	; ey   	= Exp(index) 
  	; uy 	= ey - num 
  	; logpy	= num*index - ey - lgm(num+1) $
MATRIX	; g   	= x'uy 	? first derivatives vector
	; h   	=   	? negative of 2nd derivatives
	; delta	= h * g 	$ update vector
CALC	; List 	; conv = g'delta $ convergence measure.
ENDPROCEDURE
?===================================================================
? Execute procedure until convergence, then display final results. 
?===================================================================
EXECUTE	; until conv < .00001$ 
MATRIX	; Stat(beta,h,X) $
?===================================================================
? Likelihood ratio test then Wald test then LM test
?===================================================================
CALC	; List	; loglr = sum(logpy) ; lrtest = 2*loglu - loglr)
	; critical = ctb(.95,1) $
WALD	; Start	= beta ; Var = h 
  	; Labels	= 8_b, bmth ; fcn = bmth $
POISSON	; Lhs	= num ; Rhs = x,logmth ; Start = beta, 1 
  	; Maxit	= 0 $

The program components operate in concert to allow you to create new estimators and models.  You can use the random number generators with other parts of the program to construct MCMC estimators and Gibbs samplers.

Example: MCMC Estimation of the Binary Probit Model

This program uses data augmentation to estimate the posterior means of the coefficients in a probit model.  We write the program in generic form so it could be modified for any data set.

? This is the only application specific part of the program:
?
NAMELIST	; x = the independent variables $
CREATE	; y = the dependent variable $
?
? From this point forward, the routine is completely generic
?
CALC  	; rep 	= 2000 ; burnin = 500 $
CALC  	; k 	= col(x)$
MATRIX	; xx 	= x'x ; xxi =  $
PROBIT	; Lhs 	= y ; Rhs = x $ (Classical estimator)
CALC	; iter	= 0 $
MATRIX	; beta	= init(k,1,0) ; bbar = beta ; bv = init(k,k,0)$
PROC 	= Gibbs $
?
? Gibbs sampler
?
DO FOR	; Simulate ; r = 1,rep $
CALC	; iter = iter+1$
?
? Data augmentation step:  E[yi*|beta,y,X]
?
CREATE	; mui = x'beta ; f = rnu(0,1) 
      	; if(y = 1) ysg = mui + inp(1-(1-f)*phi( mui)); 
      (else)  ysg = mui + inp(     f *phi(-mui)) $
?
? Gibbs step: E[beta|yi*,y,X]
?
MATRIX	; mb = xxi*x'ysg ; beta = rndm(mb,xxi) $
MATRIX	; if(iter > burnin) ; bbar = bbar+beta ; bv = bv+beta*beta'$
ENDDO	; Simulate $
ENDPROC	$
EXECUTE	; Proc = Gibbs $
CALC	; ri = rep - burnin ; ri = 1/ri $
MATRIX	; bbar = ri*bbar  ; bv = ri*bv-bbar*bbar' $
MATRIX	; Stat(bbar,bv,x) ; Stat(b,varb,x) $