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)