ivregress liml, does not return the likelihood e(ll), and I am not able to do likelihood ratio testing. Why is that?

[Joro Kolev]

When I fit Limited Information Maximum Likelihood regression model, say

ivregress liml price (mpg =head trunk weight)
ereturn list

there is no e(ll), the maximised likelihood is missing in what ivregress liml returns.

Then likelihood ratio testing becomes impossible:

. qui ivregress liml price (mpg =head trunk weight)

. est sto full

. qui ivregress liml price

. lrtest full .
full does not contain scalar e(ll)
r(498);


This is puzzling. The estimator is called Limited Information Maximum Likelihood, and yet it does not return a value for the maximused Likelihood...

( I understand that there are at least two ways how to solve the problem. One of them is a closed form solution to an eigen value problem. But still we are finding a solution to a maximum likelihood problem and the maximised likelihood is an important quantity that should be recoverable.)

My question is how can I recover the maximised value of the likelihood? How can I do likelihood ratio testing after ivregress liml?

[Malte Hessenius]

It is definitely not very intuitive that no e(ll) is stored.

What about ivreg2 as an alternative?

ssc install ivreg2
ssc install ranktest // necessary for the following command
ivreg2 depvar X1 X2 .... (endogenousvar = Z1 Z2), liml
ereturn list // yields e(ll)

[Joro Kolev]

Thank you Max,

ivreg2 does provide e(ll), but I do not think that this likelihood is the correct system likelihood. I spent a couple of hours studying the documentation of ivreg2, I did not find anything on how they calculate e(ll). I know that the ivreg2 e(ll) is not what I need because Pagan in a paper in Economics Letters shows that Iterated Seemingly Unrelated regression is equivalent to Limited information maximum likelihood, here it is

. ivreg2 price (mpg =head trunk weight), liml noheader
------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |  -342.7629   67.19733    -5.10   0.000    -474.4673   -211.0586
       _cons |   13465.18   1464.046     9.20   0.000      10595.7    16334.66
------------------------------------------------------------------------------
. dis e(ll)
-688.45705

. sureg (price mpg)  (mpg head trunk weight), isure nolog noheader

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
price        |
         mpg |  -342.7628   52.21744    -6.56   0.000    -445.1071   -240.4185
       _cons |   13465.18   1154.152    11.67   0.000     11203.08    15727.27
-------------+----------------------------------------------------------------
mpg          |
    headroom |   .3575597   .5797288     0.62   0.537    -.7786879    1.493807
       trunk |  -.0944815   .1355905    -0.70   0.486    -.3602339     .171271
      weight |  -.0058398   .0006548    -8.92   0.000    -.0071231   -.0045565
       _cons |   39.15977   1.700472    23.03   0.000     35.82691    42.49263
------------------------------------------------------------------------------

. dis e(ll)
-877.96867

As we can see the parameter estimates in the structural equation are (almost) the same. I guess the almost comes from numerical precision, they have to be the same according to Pagan.

But the log likelihoods are different, in fact very different -688.45705 vs -877.96867 ! (This difference spills over to the Likelihood ratio test, not shown here, but the likelihood ratio test the two produce are very different, in this example LR chi2(1) = 14.51 (ivreg2) vs LR chi2(1) = 25.24 (sureg, isure) for Ho: mpg=0. )

To preempt why am I not just using sureg when I know how to get the e(ll) via sureg: Sureg is an iterative estimator and takes ages when I do simulations.

So probably I should modify my question and say: How can I get the correct likelihood from a one step closed form (not iterative estimator). So that I can speed up my simulations to manageable time.

[David Roodman]

I'm just adding to this old thread because I hit the same issue. It looks to me like the e(ll) returned by ivreg2 uses the same formula as for OLS, and should not be used. It is computed as if you'd gotten the same fit through OLS.

In principle, one can compute the proper ll without iteration. But you need to go back to the formula for the LIML likelihood. This contains the determinant of the covariance matrix of the errors across the two equations (or more equations if you have multiple instrumented variables). Unfortunately it looks to me like neither -ivregress liml- nor -ivreg ..., liml- supplies this matrix.

One other option is my -cmp- command. It also fits the model iteratively, but in a quick test I found it twice as fast as -sureg ..., isure- for what that's worth. If using cmp, you may need to make sure that all equations have the same number of observations if you want to exactly match ivreg2. cmp will add observations to an equation if they are complete for that equation but not others. Check the e(N1), e(N2), ... return values from cmp to get equation-specific sample sizes.

Example:

cmp (price  = mpg) (mpg = head trunk weight) if price<., ind(1 1) nolr

Note that if you have exogenous controls, they need to be listed in both equations in the cmp and sureg command lines.

[Joro Kolev]

David's -cmp-, Stata's -sem-, and Stata's -sureg, isure- use different optimisation algorithms. It depends on the application which of them would calculate fastest, when they all can calculate the same thing as in Seemingly Unrelated (Triangular) Regressions.

In simulations I tend to find what David says, that -cmp- and -sem- are faster than -sureg, isure-.

Note also that both -cmp- and -sem- use Mata, so on first run when they compile their Mata objects they are substantially slower, and then on second run and on they speed up a lot.

Here are the timings in the example above, and I did verify that the estimated parameters are the same:

. timer clear

. timeit 1: qui cmp (price  = mpg) (mpg = head trunk weight) if price<., ind(1 1) nolr

. timeit 2: qui sem (price  <- mpg) (mpg <- head trunk weight), cov(e.price*e.mpg)

. timeit 3: qui sureg (price  = mpg) (mpg = head trunk weight), isure nolog

. timer list
   1:      0.05 /        1 =       0.0530
   2:      0.03 /        1 =       0.0250
   3:      0.07 /        1 =       0.0730

[Joro Kolev]

As for my solution to calculate the Likelihood Ratio test without iteration, it goes to what David said about the formula for the Likelihood Ratio Test in a Seemingly Unrelated Regressions model. The formula for the test statistic is:

LR = N*[log(det(Sr)) - log(det(Su))], where N is the number of observations, Su is the unrestricted covariance matrix of the errors, and Sr is the restricted covariance matrix of the errors.

So the following should work:

. clear

. sysuse auto
(1978 automobile data)

. qui ivregress liml price (mpg =head trunk weight)

. predict double eprice, resid

. qui reg mpg head trunk weight

. predict double empg, resid

. qui correlate eprice empg, cov

. mat Su = r(C)

. qui reg price

. predict double epriceR, resid

. qui correlate epriceR empg, cov

. mat Sr = r(C)

. matlist 74*(log(det(Sr)) - log(det(Su)))

             |        c1
-------------+----------
          r1 |  24.87883

So according to the non-iterated procedure, the LR test statistic for the Ho: b_mpg = 0, is LR = 24.87883 .

Now with the -sureg, isure- we get:

. qui sureg (price  = mpg) (mpg = head trunk weight), isure nolog

. est sto full

. qui sureg (price ) (mpg = head trunk weight), isure nolog

. lrtest full .

Likelihood-ratio test
Assumption: . nested within full

 LR chi2(1) =  25.24
Prob > chi2 = 0.0000

So we are pretty close with the LR statistic, 25.24 vs 24.87883.

I think the small difference comes from ambiguity as to what the estimation method should be under Ho. The first procedure does OLS equation by equation, and the second procedure does -sureg, isure- under Ho.

[David Roodman]

Here's a much better solution. It goes back to the math of LIML. As I mentioned, to compute the log likelihood, you need the covariance matrix of the residuals from all equations. Generating the residuals from the second-stage/structural equation is easy. To get the residuals from the first-stage/reduced-form equations, you need to run those regressions manually while including the second-stage residuals as an extra control. See generally appendix B of the Roodman et al. (2019) boottest article in the Stata Journal and Davidson and MacKinnon (1993), chapter 18.

Here's a little program that's fairly general. It runs after an -ivregress liml- estimate. It could be extended to support ivreg2 , liml, but I think ivreg2 is slower. It handles observations weights and cases of multiple instrumented variables, as well as overidentified models. I haven't tested it thoroughly.

program define get_ll, rclass
  _assert e(cmd)=="ivregress" & e(estimator)=="liml", msg(only works after -ivregress liml-) rc(198)
  tempvar V e0
  local es `e0'
  predict double `e0' if e(sample), resid
  mvreg `e(endog)' = `e(exog)' `e0' [aweight `e(wexp)']
  forvalues i=1/`e(k_eq)' {
    tempname e`i'
    local es `es' `e`i''
    predict double `e`i'' if e(sample), resid eq(#`i')
  }
  mat accum `V' = `es' [`e(wtype)'`e(wexp)'], nocons
  return scalar ll = -e(N) / 2 * ((1+e(k_eq)) * (ln(2*_pi / e(N)) + 1) + ln(det(`V')))
end

Example:

ivregress liml price gear_ratio (mpg length = headroom trunk weight) [aw=mpg]
get_ll
scalar ll = r(ll)

cmp (price = gear_ratio mpg length) (mpg = headroom trunk weight gear_ratio) (length = headroom trunk weight gear_ratio) [aw=mpg] if e(sample), nolr ind(1 1 1) qui tech(bhhh)

di ll, e(ll)

At the end I get a match between the result of get_ll and cmp:

. di ll, e(ll)
-1121.3231 -1121.3231

[David Roodman]

Oh, here's a better version that will also run after -ivregress 2sls- when it is exactly identified. For when the number of instruments = the number of instrumented variables, 2SLS and LIML are the same.

program define get_ll, rclass
  _assert e(cmd)=="ivregress" & (e(estimator)=="liml" | `:word count `e(exog)''==`:word count `e(endog)' `e(exogr)''), ///
      msg(only works after -ivregress liml- or exactly identified -ivregress 2sls-) rc(198)
  tempvar V e0
  local es `e0'
  predict double `e0' if e(sample), resid
  mvreg `e(endog)' = `e(exog)' `e0' [aweight `e(wexp)']
  forvalues i=1/`e(k_eq)' {
    tempname e`i'
    local es `es' `e`i''
    predict double `e`i'' if e(sample), resid eq(#`i')
  }
  mat accum `V' = `es' [`e(wtype)'`e(wexp)'], nocons
  return scalar ll = -e(N) / 2 * ((1+e(k_eq)) * (ln(2*_pi / e(N)) + 1) + ln(det(`V')))
end

[David Roodman]

Even better. :-)
This version works after ivreg2, restores the user's estimation results at exit, and also does the math slightly more elegantly.

As an alternative to the procedure I described above, adding the second-stage residuals as controls in the first stage, one can just regress all the endogenous variables, including the dependent variable, on all the exogenous variables; get the residuals; form their covariance matrix; get its determinant; and then multiply by the value of the kappa parameter returned by ivregress (or "lambda" from ivreg2). That's the same as the determinant described ealier. Then one plugs that into the formula for the log likelihood.

program define get_ll2, rclass
  tempname V hold kappa
  if e(cmd)=="ivregress" {
    local mvregeq `e(depvar)' `e(endog)' = `e(exog)'
    scalar `kappa' = e(kappa)
    local k_endog = e(k_endog)
  }
  else {
    local mvregeq `e(depvar)' `e(instd)' = `e(insts)'
    scalar `kappa' = e(lambda)
    local k_endog = e(endog_ct)
  }
  _assert e(cmd)=="ivregress" & (`"`e(estimator)'"'=="liml" | `:word count `e(exog)' '==`:word count `e(endog)' `e(exogr)'') | ///
          e(cmd)=="ivreg2"    & (`"`e(model)'"'    =="liml" | e(endog_ct)==e(exexog_ct)) , ///
      msg(only works after -ivregress liml-, -ivreg2 ..., liml- or exactly identified 2SLS) rc(198)
  local wexp `"`e(wexp)'"'
  local weight = cond(`"`wexp'"'=="", "", "aweight")
  _est hold `hold', restore
  mvreg `mvregeq' [`weight' `wexp']
  forvalues i=1/`e(k_eq)' {
    tempvar e`i'
    local es `es' `e`i''
    predict double `e`i'' if e(sample), resid eq(#`i')
  }
  _est unhold `hold'
  mat accum `V' = `es' [`e(wtype)'`e(wexp)'], nocons
  return scalar ll = -e(N) / 2 * ((1+`k_endog') * (ln(2*_pi / e(N)) + 1) + ln(det(`V') * `kappa'))
end

[David Roodman]

I'm getting weird errors from the Stata forum, and edits to posts are getting lost. In the penultimate line above, "`e(k_endog)'`e(endog_ct)'" should be "`k_endog'".