I have estimated SUR model with systemfit (R package). With the estimated results, I am trying to get logLik, AIC and BIC for each of SUR model. For the whole logLik, I can just use logLik function. However, how can I get the goodness-of-fit (logLik, AIC, and BIC) for each equation of SUR model?

Here is a sample of SUR estimation with systemfit.

data( "Kmenta" )
eqDemand <- consump ~ price + income
eqSupply <- consump ~ price + farmPrice + trend
system <- list( demand = eqDemand, supply = eqSupply )

fitsur <- systemfit( system, "SUR", data = Kmenta )    
logLik( fitsur )

How can I get the logLik, AIC and BIC for each equation (eqDemand,eqSupply)?


1 Answer 1


Here is how I would proceed, inspired by, e.g., this and this thread. In principle, one could of course consider other likelihoods than the normal, but there seems to be no good reason to do so here.

eq2 <- summary(fitsur$eq[[2]])                              
n <- eq2$df[1] + eq2$df[2]                                  # get sample size
sigma.sq.hat.MLE <- eq2$ssr/n                               # the MLE of the error variance 
logLiK.eq2 <- -n/2*(log(2*pi) + 1 + log(sigma.sq.hat.MLE))  # the log-likelihood
p <- eq2$df[1] + 1                                          # the number of parameters, including the error variance
aic.eq2 <- -2*logLiK.eq2 + 2*p                              # definition of AIC

The generalization to BIC should be obvious, and surely also how to get equation 1.

The reasoning underlying the above calculations is that the likelihood of such a system of equations is matrix normal with, at the MLE, means equal to the fitted values and column-covariance matrix equal to the covariance matrix of the residuals across equations, say, $V$. Equivalently, the likelihood then is multivariate normal, where the observations belonging to equation $i$ are normal with variance $V_{ii}I_n$, thanks to the properties of marginal distributions of multivariate normals.

To illustrate the point about the matrix normal, let us compare that log-likelihood to logLik. I generate a somewhat larger dataset, as the small Kmenta dataset leads to small discrepancies I haven't traced back:


n <- 1e3  # be careful with larger n - dmatrixnorm takes a long time!!
y1 <- rnorm(n)
y2 <- rnorm(n)
x <- rnorm(n)

eq1 <- y1 ~ x
eq2 <- y2 ~ x
system <- list(eq1, eq2)

fitsur <- systemfit(system, "SUR")    
M <- cbind(fitted(fitsur$eq[[1]]), fitted(fitsur$eq[[2]]))
U <- diag(n)
V <- fitsur$residCov

> dmatrixnorm(cbind(y1,y2), M, U, V, log=T) 
[1] -2848.044

> logLik(fitsur)
'log Lik.' -2848.042 (df=7)
  • $\begingroup$ I just checked your code. Sorry for the late feedback.Your answer is the exactly same as I did before. My concern was the correlation between equations. With the function 'logLik' from the 'systemfit' package, we can get the log-likelihood for the whole equation. Do you think logLik(fitsur) = logLiK.eq2 + logLiK.eq1? $\endgroup$ Mar 16, 2019 at 19:34
  • $\begingroup$ @Johnlegend, I spell out my reasoning in mroe detail - I am not so sure about the additive part, though. $\endgroup$ Mar 17, 2019 at 17:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.