2
$\begingroup$

I am following up on the question "Why does systemfit yield identical results for OLS and WLS?". It deals with estimating a system of linear equations where each of them has one unique regressor and one common regressor. The estimation is done via OLS and WLS using the systemfit package in R.

It turns out that the point estimates are equal, $\hat\beta_\text{OLS}=\hat\beta_\text{WLS}$ and the estimated coefficient covariance matrices are identical, too, $\widehat{\text{Cov}}_\text{OLS}(\hat\beta) = (X'\hat\Omega_\text{OLS}^{-1}X)^{-1} =(X'\hat\Omega_\text{WLS}^{-1}X)^{-1} = \widehat{\text{Cov}}_\text{WLS}(\hat\beta)$ because the estimated error covariance matrices are identical, $\hat\Omega_\text{OLS}=\hat\Omega_\text{WLS}$ (unless equal error variances across equations are imposed by control = systemfit.control(singleEqSigma = FALSE), but this is not the default option). I can also replicate the findings manually:

library(systemfit)

# Generate and prepare the data
n <-  8  # sample size
m <-  2  # length of the "second part" of the sample
N <-  3  # number of equations
set.seed(321); x <- matrix(rnorm(n*N),ncol=N); colnames(x) <- paste0("x",1:N) # generate regressors
dummy <- c( rep(0,n-m), rep(1,m) ) # generate a common regressor
x <- cbind(x,dummy)                # include the common regressor with the rest of the regressors
set.seed(123); y <- matrix(rnorm(n*N),ncol=N); colnames(y) <- paste0("y",1:N) # a placeholder for dependent variables
for(i in 1:N){
 y[,i] <- i + sqrt(i)*x[,i] - i*dummy + y[,i]*15*sqrt(i)   
 # y[,i] is a linear function of x[,i] and dummy, 
 # plus an error term with equation-specific variance - just what WLS is made for
}
data1 <- as.data.frame(cbind(y,x)) # create a data frame of all data (y and x)

# Create the model equations
eqSystem <- list() 
for(i in 1:N){
 eqSystem[[i]] <- as.formula(assign(paste0("eq",i), value=paste0("y",i," ~ x",i," + dummy"))) # define linear equations of SUR
}

# Estimate the model with `method="OLS"` and `method="WLS"`
m1a <- systemfit(formula=eqSystem, method="OLS", data=data1, control=systemfit.control(singleEqSigma=FALSE))
m1b <- systemfit(formula=eqSystem, method="OLS", data=data1, control=systemfit.control(singleEqSigma=TRUE ))
m2a <- systemfit(formula=eqSystem, method="WLS", data=data1, control=systemfit.control(residCovWeighted=TRUE ))
m2b <- systemfit(formula=eqSystem, method="WLS", data=data1, control=systemfit.control(residCovWeighted=FALSE))
#m3 <- systemfit(formula=eqSystem, method="SUR", data=data1)
summary(m1a, residCov=FALSE, equations=FALSE)
summary(m1b, residCov=FALSE, equations=FALSE)
summary(m2a, residCov=FALSE, equations=FALSE)
summary(m2b, residCov=FALSE, equations=FALSE)
#summary(m3, residCov=FALSE, equations=FALSE)

# Covariance matrix of the coefficients: OLS and WLS are equal
round(m1a$coefCov,4)
round(m1b$coefCov,4)
round(m2a$coefCov,4)
round(m2b$coefCov,4)
#round(m3$coefCov,4)
m1b$coefCov/m2a$coefCov # equal

# Covariance matrix of the residuals: all are equal
round(m1a$residCov,4)
round(m1b$residCov,4)
round(m2a$residCov,4)
round(m2b$residCov,4)
#round(m3$residCov,4)
m1b$residCov/m2a$residCov # equal

# Weighting matrix in the 2-step estimator: OLS is NULL, WLS and SUR differ
m1a$residCovEst
m1b$residCovEst
round(m2a$residCovEst,4)
round(m2b$residCovEst,4)
#round(m3$residCovEst,4)

#---------- Manual estimation under no restrictions

# Now generate x, dummy and y above, estimate the models above. Then:
# Manually construct the X matrix corresponding to equation (2) from the vignette
O_N3 <- matrix(0,nrow=n,ncol=N); 
X <- rbind( cbind( cbind(1,x[,c(1,N+1)]),O_N3                 ,O_N3                  ) , 
            cbind( O_N3                 ,cbind(1,x[,c(2,N+1)]),O_N3                  ) , 
            cbind( O_N3                 ,O_N3                 ,cbind(1,x[,c(3,N+1)]) ) )
#round(X,2)

# Obtain $\widehat{Cov}_{OLS}(\hat\beta)$ and compare to the output from `systemfit`
sigma2      <- var(unlist(resid(m1a)))
Omega_OLSa  <- sigma2 * diag(n*N)
inv_cov_OLSa<- t(X) %*% solve(Omega_OLSa) %*% X
cov_OLSa    <- solve(inv_cov_OLSa)
round(cov_OLSa,4)
round(m1a$coefCov,4)
#round(m1b$coefCov,4)
round(cov_OLSa/m1a$coefCov,4) 
#round(cov_OLSa/m1b$coefCov,4) 
# The ratio of estimated covariances between my solution and that of `systemfit` is constant. 
# The constant approaches 1 as the sample size increases (compare n=8 to n=100 to n=300). 
# I guess `systemfit` uses a small sample correction while I do not.

# Obtain $\widehat{Cov}_{WLS}(\hat\beta)$ and compare to the output from `systemfit`
Omega_WLSa  <- diag(diag(m2a$residCov)) %x% diag(n)
inv_cov_WLSa<- t(X) %*% solve(Omega_WLSa) %*% X
cov_WLSa    <- solve(inv_cov_WLSa)
round(cov_WLSa,4)
round(m2a$coefCov,4)
round(cov_WLSa/m2a$coefCov,4)
# The ratio of estimated covariances between my solution and that of `systemfit` is 1.

# Obtain $\widehat{Cov}_{WLS}(\hat\beta)$ and compare to the output from `systemfit`
Omega_WLSb  <- diag(diag(m2b$residCov)) %x% diag(n)
inv_cov_WLSb<- t(X) %*% solve(Omega_WLSb) %*% X
cov_WLSb    <- solve(inv_cov_WLSb)
round(cov_WLSb,4)
round(m2b$coefCov,4)
round(cov_WLSb/m2b$coefCov,4)
# The ratio of estimated covariances between my solution and that of `systemfit` is 1.

However, things change under restricted estimation. E.g. if we impose some of the coefficients to be equal across equations, systemfit reports different point estimates and estimated coefficient covariance matrices for OLS and WLS. This is puzzling, given the equations (26) and (27) from the package's vignette. There, the exact same equation is used for both OLS and WLS estimators. The only thing that could differ is the $\hat\Omega_{\text{OLS}}$ vs. $\hat\Omega_{\text{WLS}}$, but we have learned they are identical. So how come the actual numbers produced by systemfit are different for OLS and WLS?

# OLS and WLS turn out to be different when the model has linear restrictions on coefficients from different equations
Rmat0 <- matrix(0, nrow=N-1, ncol=N*3) # nrow ~ # of restrictions, ncol ~ # of unrestr. coeffs
for(i in 1:(N-1)) Rmat0[i, c(3*i,3*(i+1))] <- c(1,-1)
qvec0 <- rep(0,nrow(Rmat0))
m1r <- systemfit(formula=eqSystem, method="OLS", restrict.matrix=Rmat0, restrict.rhs=qvec0, data=data1) # restricted model
m2r <- systemfit(formula=eqSystem, method="WLS", restrict.matrix=Rmat0, restrict.rhs=qvec0, data=data1) # restricted model
summary(m1r, residCov=FALSE, equations=FALSE)
summary(m2r, residCov=FALSE, equations=FALSE)

Update 1: Perhaps the key to the answer is in footnote 3 of the vignette? The footnote says that $\hat\Omega$ may denote two different things: (1) the estimated error covariance matrix based on the first-stage OLS estimates and (2) the estimated error covariance matrix (?) based on the second-stage WLS estimates (in case of the WLS estimator only) - what is that?
We just have to figure out which is which. Now, to make sure follow, I am trying to replicate the OLS and WLS estimators manually (see the code below). I am able to replicate the point estimates of OLS. (I have not tried to replicating the estimated covariance matrix of coefficients yet.) I am also able to replicate the estimated covariance matrix of the coefficients up to a scalar multiple. The difference shrinks with sample size, so that suggests I am not using the small-sample correction that systemfit is using.
However, I cannot replicate WLS, as I do not understand what I should use for $\hat\Omega$. If I use the same as for OLS, then of course I am not getting the same results as by systemfit. How do I replicate restricted WLS? (See the answer right below.)

Update 2: I have found the weighting matrix $\hat\Omega$ that yields the restricted WLS result from systemfit. As an estimate of $\Sigma$, it takes the estimated error covariance matrix from the restricted OLS and sets the off-diagonal elements to zero. This is how to replicate the point estimate of restricted WLS. Similarly, I can replicate the estimated covariance matrix of the coefficients.

# Cross-equation restrictions followed by estimation by OLS and WLS
R <- matrix(0, nrow=N-1, ncol=N*3) # nrow ~ # of restrictions, ncol ~ # of unrestr. coeffs
for(i in 1:(N-1)) R[i, c(3*i,3*(i+1))] <- c(1,-1)
q <- rep(0,nrow(R))
m1ar <- systemfit(formula=eqSystem, method="OLS", restrict.matrix=R, restrict.rhs=q, data=data1, control=systemfit.control(singleEqSigma=FALSE))
m1br <- systemfit(formula=eqSystem, method="OLS", restrict.matrix=R, restrict.rhs=q, data=data1, control=systemfit.control(singleEqSigma=TRUE ))
m2ar <- systemfit(formula=eqSystem, method="WLS", restrict.matrix=R, restrict.rhs=q, data=data1, control=systemfit.control(residCovWeighted=TRUE ))
m2br <- systemfit(formula=eqSystem, method="WLS", restrict.matrix=R, restrict.rhs=q, data=data1, control=systemfit.control(residCovWeighted=FALSE))

# Covariance matrix of the coefficients: OLS (version b with different error variances across equations) and WLS are equal
round(m1ar$coefCov,4)
round(m1br$coefCov,4)
round(m2ar$coefCov,4)
round(m2br$coefCov,4)
#round(m3r$coefCov,4)
m1br$coefCov/m2ar$coefCov # equal

# Covariance matrix of the residuals: 1a=1b and 2a=2b but 1a!=2a
round(m1ar$residCov,4)
round(m1br$residCov,4)
round(m2ar$residCov,4)
round(m2br$residCov,4)
#round(m3r$residCov,4)
m1ar$residCov/m2ar$residCov # not equal

# Weighting matrix in the 2-step estimator: OLS is NULL, WLS and SUR differ
      m1ar$residCovEst
  m1br$residCovEst
round(m2ar$residCovEst,4)
round(m2br$residCovEst,4)
#round(m3r$residCovEst,4)

#---------- Manual estimation under cross-equation restrictions

# Choose only one of the following group of lines
Omega <- Omega_OLSa; inv_cov <- inv_cov_OLSa # OLS, replicates OLS from `systemfit`
#Omega <- Omega_WLSa; inv_cov <- inv_cov_WLSa # WLS
#Omega <- Omega_WLSb; #inv_cov <- inv_cov_WLSb # WLS, exact same results as the line above
Omega <- m2ar$residCovEst %x% diag(n); inv_cov <- t(X) %*% solve(Omega) %*% X # replicates WLS from `systemfit`

# Construct the restricted estimator: the first part
Mat0 <- matrix(0,nrow=nrow(R),ncol=nrow(R))
A1 <- cbind(inv_cov,t(R))
A2 <- cbind(R,Mat0)
A  <- rbind(A1,A2)
#round(A,3)

# Construct the restricted estimator: the second part
B1 <- t(X) %*% solve(Omega) %*% cbind(c(y))
B2 <- cbind(q)
B  <- rbind(B1,B2)
#round(B,3)

# Obtain the restricted estimator by combining the parts
beta_lambda <- solve(A) %*% B
#round(beta_lambda,3)
beta <- head(beta_lambda,-length(q))

# Compare the manually obtained restricted estimator with OLS and WLS from `systemfit`
table1 <- cbind(beta,coef(m1ar),coef(m1br),coef(m2ar),coef(m2br)); 
colnames(table1) <- c("Manual","OLSa","OLSb","WLSa","WLSb")
rownames(table1) <- c("eq1_(Intercept)","eq1_x1","eq1_dummy", "eq2_(Intercept)","eq2_x2","eq2_dummy", "eq3_(Intercept)","eq3_x3","eq3_dummy")
# OLS1 uses common error variance across equations; OLS2 uses different error variances in different equations (just like WLS)
round(table1,3)

# Covariance matrix of the coefficients
lenb=length(beta); cov_=solve(A)[1:lenb,1:lenb]; rownames(cov_)=colnames(cov_)=rownames(table1)
round(cov_,4)
round(m1ar$coefCov,4) # for cov_ obtained by manual OLS, this is a constant. It approaches 1 with sample size.
round(m1br$coefCov,4) # for cov_ obtained by manual WLS, this is 1, which is good.
round(m2ar$coefCov,4)
#round(m2br$coefCov,4)
round(cov_/m1ar$coefCov,4)
round(cov_/m1br$coefCov,4) 
round(cov_/m2ar$coefCov,4) 
#round(cov_/m2br$coefCov,4) 
# The ratio of 1 is good.
Improve this question
$\endgroup$
13
  • $\begingroup$ I see the problem, but it is a great mess and some work to disentangle that is not easily expressed in a simple post. As an answer I might create some overview. But, maybe you could explain first if there is something more of a specific background behind this problem? This is the third post on this systemfit package and your 'digging through it' seems to be hard (which makes sense that it creates multiple questions). I am wondering what the underlying final project/problem is that you are working on${\color{red}{^*}}$. $\endgroup$
    – Sextus Empiricus
    Commented Apr 18 at 13:41
  • $\begingroup$ @SextusEmpiricus, this is for a class I will be teaching in the future. I have a system of equations, and economic theory implies some restrictions (both zero restrictions and cross-equation equality restrictions). My examples posted here are the same as the real ones I wish to use (same variables), but I use simulated data here instead of real data. I want to demonstrate alternative estimation techniques, discussing assumptions behind them. It turns out I do not understand the techniques as well as I thought I did (when I studied the topic last time 11 years ago). $\endgroup$
    – Richard Hardy
    Commented Apr 18 at 13:47
  • $\begingroup$ ${\color{red}{^*}}$With that background in mind, an answer might be more directed towards your actual goals. To me these OLS and WLS methods in the systemfit package seem a bit superfluous and it the package appears to me to be more about the SUR situation. Although, while creating a solution to SUR the OLS and WLS designed on the way do have some usefull applications (by allowing fitting of restricted regression, but withoot restricted regression the package has not much use). $\endgroup$
    – Sextus Empiricus
    Commented Apr 18 at 13:48
  • $\begingroup$ My problem with SUR is that my sample is not very large but the number of equations is large. So estimation imprecision of the covariance matrix gets really wild; there are so many parameters and so few data points. Thus I think OLS and WLS may work better in practice than SUR would. Actually, I will go to system-GMM later, as I want HAC robust covariance matrix. That should amount to OLS with HAC robust covariance matrix. But for pedagogical purposes (and not least my own understanding) I think I should figure out OLS, WLS and SUR before system-GMM. $\endgroup$
    – Richard Hardy
    Commented Apr 18 at 13:51
  • $\begingroup$ So the goal seems to be to get from the basic understanding of the expression of a problem with multiple equations (which get connected either with constrains or with correlated errors), to the underlying algorithms that are used to solve such problem. $\endgroup$
    – Sextus Empiricus
    Commented Apr 18 at 13:55

1 Answer 1

Reset to default
2
$\begingroup$

About the different $\hat\Omega$

the definition of $\hat{\Omega}$ is driving me nuts. In my previous post, you helped me discover that it was identical for both OLS (allowing different variances in different equations) and WLS. The two were the exact same thing...

They are not exactly the exact same thing. For a single method, there are multiple $\hat\Omega$. The reason is because $\Omega$ is unknown in the beginning and the methods try to estimate this by repeatedly re-estimating and updating $\hat\Omega$ in several steps.

flow scheme

  • Method OLS

    This will not perform the loop and only does

    'GLS fitting'
   |
   V
'recomputing for estimating error'

  • Method WLS

    This will perform the extra loop once and does

    'GLS fitting'
   |
   V
'recomputing for new step'
   |   
   V
'GLS fitting (2nd time)'
   |
   V
'recomputing for estimating error'

The second 'GLS fitting' is using a different $\Omega$. OLS is like performing GLS with an identity matrix (which happens the first time) and WLS is like performing GLS with a diagonal matrix (which happens the second time).

The $\Omega$ are different for those 'GLS fitting' steps one and two. But... the $\Omega$'s after the step 'recomputing $\Omega$ for estimating error' will end up the same again if there are no restrictions (because then the estimates are the same).

Influence of restrictions

A simple model is when we have two equations with only an intercept to estimate. Below we see the fits for 2 x 2 variations:

  • On the left we see fitting without restrictions, which means that each equation is effectively fitted on it's own. (only when there are correlations between the error terms will there be an effect)

    The WLS makes no difference versus OLS. There are weights being applied, but within a single equation these weights are all the same. The data points that influence the intercept of equation 1 have all the same weight, and the data points that influence the intercept of equation 2 have all the same weight.

  • On the right we see fitting with restrictions. In this case we assume that the intercepts are the same (the restriction equation $a_1 + a_2 = 0$), which means that we effectively fit only one parameter/intercept.

    The WLS makes a difference versus OLS. When we apply weights, then the data points from equation 1 are made more important (since they have lower variance, which is estimated after the first fitting) and we see that the fitted line is slightly lower.

four different cases

So on the left, OLS and WLS result in the same fit

(But not the same $\Omega$ during the fitting, namely identity matrix versus diagonal matrix. Yet, the $\Omega$ recomputed for estimating the error will be the same again)

And on the right, OLS and WLS result in different fits.

Improve this answer
$\endgroup$
4
  • $\begingroup$ That looks promising, so +1! I will give it a go a bit later. (I am too exhausted to try it now, I would certainly make some mistakes.) But I can tell right away that I think I am already doing Step WLS a) as you described, and I fail to replicate the estimate from systemfit. (I have not attempted Step WLS b) yet, and my post says nothing about that part.) $\endgroup$
    – Richard Hardy
    Commented Apr 18 at 15:08
  • $\begingroup$ I am now able to replicate both the point estimate of WLS and the estimate of the covariance matrix of the coefficients for both OLS (almost) and WLS; see Update 2 to the OP. I do not think I follow Step WLS b), though. I am able to obtain the coefficient covariance matrix for WLS using only the $\hat\Omega'_\text{after OLS}$ from OLS. $\endgroup$
    – Richard Hardy
    Commented Apr 18 at 17:16
  • $\begingroup$ I have accepted your answer, as it was instrumental in solving the problem. I still think it is not entirely accurate, so you may still wish to updated it as indicated in my last comment. In any case, a big THANK YOU for all your help! You are one of the few users that I look up to, and you continue to impress me. $\endgroup$
    – Richard Hardy
    Commented Apr 23 at 6:51
  • $\begingroup$ With this update you knocked it out of the park! $\endgroup$
    – Richard Hardy
    Commented May 6 at 19:35

Your Answer

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

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