Find (or calculate) log-likelihood value, AIC, and BIC for SUR model (for each equation) with systemfit

Question

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)?

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:

library(systemfit)
library(LaplacesDemon)

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)