Why does `systemfit` yield identical results for OLS and WLS?

Question

I am estimating a system of seemingly unrelated regressions (SUR) using the systemfit package in R. Each of the equations has one unique regressor and one common regressor. I use two alternative estimation methods: ordinary least squares, OLS (method="OLS") and weighted least squares, WLS (method="WLS"), as discussed in section 2.1 (p. 3-4) of the systemfit vignette. I get the same point estimates and standard errors from both. The covariance matrices of coefficients are also identical. That puzzles me. Shouldn't the standard errors and more generally, the covariance matrices of coefficients differ?

Question: Why does systemfit yield identical results for method="OLS" and method="WLS"?

library(systemfit)

# Generate and prepare the data
n <- 1000 # sample size
m <- 100  # length of the "second part" of the sample
N <- 3    # number of eqations
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"`
m1 <- systemfit(formula=eqSystem, method="OLS", data=data1)
m2 <- systemfit(formula=eqSystem, method="WLS", data=data1)
summary(m1, residCov=FALSE, equations=FALSE)
summary(m2, residCov=FALSE, equations=FALSE)
m1$coefCov # covariance matrix of coefficients
m2$coefCov # covariance matrix of coefficients

A follow-up question: "Why does systemfit yield different results for OLS and WLS under cross-equation restrictions?"

Answer 1

The weights only play a role when there are restrictions on the regression. Below is an example for

$$\begin{array}{} y_1 &=& a_1 + b_1 x_1 + \epsilon_1 \\ y_2& = &a_2 + b_2 x_2 + \epsilon_2 \\ \rlap{\text{with the condition $b_1 + b_2 = 0$}} \end{array}$$

the restriction makes that the estimation of $b_1$ and $b_2$ are coupled. effectively we are guessing a single parameter:

$$\begin{array}{} y_1 &=& a_1 + b x_1 + \epsilon_1 \\ y_2& = &a_2 - b x_2 + \epsilon_2 \\ \end{array}$$

and this makes that the potential differences in variance of $\epsilon_1$ and $\epsilon_2$ play a role.

Code example:

library(systemfit, quietly = TRUE)

# Generate and prepare the data
set.seed(1)
n <- 1000 # sample size
x1 = rnorm(n)
x2 = rnorm(n)
epsilon = rnorm(n)
y1 = 1 + x1 + epsilon
y2 = 1 - x2 + epsilon*2+rnorm(n)

### regressions
eqSystem <- list(eq1 = as.formula("y1~x1"),
                 eq2 = as.formula("y2~x2")) 


### perform fits 
restrict = matrix(c(0,1,0,1), ncol = 4)
rhs= c(0)


m1 <- systemfit(formula=eqSystem, 
                method="OLS", 
                control = systemfit.control(residCovWeighted = FALSE),
                restrict.matrix = restrict,
                restrict.rhs = rhs
)
m2 <- systemfit(formula=eqSystem, 
                method="WLS", 
                control = systemfit.control(residCovWeighted = TRUE),
                restrict.matrix = restrict,
                restrict.rhs = rhs
)

summary(m1, equations=FALSE)
summary(m2, equations=FALSE)

The covariance matrices of coefficients are also identical.

The difference between WLS and OLS methods is just that, after the OLS step, the WLS adds an extra weighted least squares step (equivalent to a GLS step with a diagonal covariance matrix). But those weights have no influence if the equations are without restrictions.

In the output you can see that the WLS has an extra covariance matrix

in my example this is:

The covariance matrix of the residuals used for estimation
        eq1     eq2
eq1 1.06489 0.00000
eq2 0.00000 5.38194

The covariance matrix of the residuals
        eq1     eq2
eq1 1.06201 2.13663
eq2 2.13663 5.38773 

This first diagonal matrix does not occur with OLS.

• The OLS method just computes a single OLS step and makes an estimate of the covariance matrix.
• Before computing the covariance matrix, the WLS method computes a weighted regression with a diagonal covariance matrix.

Without restrictions, the weights have only influence on the observations of a single equation, and they will be all the same weight for each row of observations within that single equation. So the extra step with the WLS method has effectively no influence on the fitted coefficients.

The regression of the multiple equations can be seen as a single equation that is stacked.

If the regressor matrices are like

$$\begin{bmatrix} 1 & x_{1,1} &\text{dummy_{1,1}} \\ 1 & x_{1,2} &\text{dummy_{1,2}}\\ 1 & x_{1,3} &\text{dummy_{1,3}} \\ 1 & x_{1,4} &\text{dummy_{1,4}}\\ 1 & x_{1,5} &\text{dummy_{1,5}} \end{bmatrix} \quad \begin{bmatrix} 1 & x_{2,1} &\text{dummy_{2,1}} \\ 1 & x_{2,2} &\text{dummy_{2,2}}\\ 1 & x_{2,3} &\text{dummy_{2,3}} \\ 1 & x_{2,4} &\text{dummy_{2,4}}\\ 1 & x_{2,5} &\text{dummy_{2,5}} \end{bmatrix}\quad \begin{bmatrix} 1 & x_{3,1} &\text{dummy_{3,1}} \\ 1 & x_{3,2} &\text{dummy_{3,2}}\\ 1 & x_{3,3} &\text{dummy_{3,3}} \\ 1 & x_{3,4} &\text{dummy_{3,4}}\\ 1 & x_{3,5} &\text{dummy_{3,5}} \end{bmatrix} $$

then the three regressions can be solved as a single regression

$$\begin{bmatrix} 1 & x_{1,1} &\text{dummy_{1,1}}&0&0&0&0&0&0 \\ 1 & x_{1,2} &\text{dummy_{1,2}}&0&0&0&0&0&0\\ 1 & x_{1,3} &\text{dummy_{1,3}}&0&0&0&0&0&0 \\ 1 & x_{1,4} &\text{dummy_{1,4}}&0&0&0&0&0&0\\ 1 & x_{1,5} &\text{dummy_{1,5}} &0&0&0&0&0&0\\ 0&0&0&1 & x_{2,1} &\text{dummy_{2,1}}&0&0&0 \\ 0&0&0&1 & x_{2,2} &\text{dummy_{2,2}}&0&0&0\\ 0&0&0&1 & x_{2,3} &\text{dummy_{2,3}}&0&0&0 \\ 0&0&0&1 & x_{2,4} &\text{dummy_{2,4}}&0&0&0\\ 0&0&0&1 & x_{2,5} &\text{dummy_{2,5}}&0&0&0\\ 0&0&0&0&0&0&1 & x_{3,1} &\text{dummy_{3,1}} \\ 0&0&0&0&0&0&1 & x_{3,2} &\text{dummy_{3,2}}\\ 0&0&0&0&0&0&1 & x_{3,3} &\text{dummy_{3,3}} \\ 0&0&0&0&0&0&1 & x_{3,4} &\text{dummy_{3,4}}\\ 0&0&0&0&0&0&1 & x_{3,5} &\text{dummy_{3,5}} \end{bmatrix}$$

Multiplying this with a diagonal weights matrix will only multiply the columns of the same regression, and it doesn't change the fit.

Answer 2

With a ton of patient guidance by Sextus Empiricus, I got convinced that the point estimates from OLS and WLS are not supposed to differ. (Note: I did not claim the opposite in the OP, but I was unsure and lost in detail.) However, I was still puzzled by the fact that the estimated covariance matrices of the coefficients do not differ in systemfit. Using the notation from the systemfit vignette, the estimated covariance matrix is $\widehat{\text{Cov}}(\hat\beta)=(X'\hat\Omega^{-1}X)^{-1}$ (equation 8 on p. 3).

• In the case of OLS, $$\hat\Omega_{\text{OLS}}=\hat\sigma^2 I_{TG}\color{red}{^*}$$ where $T$ is the number of observations (n in my code) and $G$ the number of equations (N in my code).
• In the case of WLS, $$\hat\Omega_{\text{OLS}}=\hat\Sigma \bigotimes I_{T}$$ where the elements $\hat\sigma_{ij}$ of $\hat\Sigma$ are such: $\hat\sigma_{ij}=0$ for $i\neq j$ and $\hat\sigma_{ii}=\sigma_i^2$, so the latter is the estimated variance of the disturbance terms in the $i$th equation.

I do not see how $$\widehat{\text{Cov}}_{\text{OLS}}(\hat\beta)=(X'\hat\Omega_{\text{OLS}}^{-1}X)^{-1}$$ is the same as $$\widehat{\text{Cov}}_{\text{WLS}}(\hat\beta)=(X'\hat\Omega_{\text{WLS}}^{-1}X)^{-1}.$$

The following R code illustrates manual calculations of $\widehat{\text{Cov}}_{\text{OLS}}(\hat\beta)$ and $\widehat{\text{Cov}}_{\text{WLS}}(\hat\beta)$.$\color{red}{^{**}}$ It requires running the code of the OP first.

# 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(m1)))
cov_inv_OLS= 1/sigma2 * t(X) %*% X
cov_OLS=solve(cov_inv_OLS)
round(cov_OLS,4)
round(m1$coefCov,4)
round(cov_OLS/m1$coefCov,4) 
# The ratio of estimated covariances between my solution and that of `systemfit` 
# is constant per equation. However, it differs across equations. 

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

I am able to replicate the systemfit result for WLS but not OLS. It seems to me systemfit is doing WLS instead of OLS - contrary to what the vignette claims at the bottom of p. 3 and the first sentence on p. 4. However, right after that it presents the alternative of defining the OLS estimator in exactly the same way as they define the WLS estimator a few lines below, and that seems to be the definition of OLS that they actually implement. (Compare lines 3-5 of p. 4 with lines 3-5 of the subsection Weighted Least Squares (WLS) on the same page.) Not sure what to conclude from all that, but at least we see some details...

$\color{red}{^*}$Actually, the vignette says $\hat\Omega_{\text{OLS}}^{-1}=I_{TG}$, so it misses the factor $\hat\sigma^2$. Correct me if I am wrong.

$\color{red}{^{**}}$Sorry about less than elegant code. Print some intermediate results to convince yourselves that it is doing what it is doing. Set n=8, m=2 and N=3 in the code of the OP for more transparent results.