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

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

• @SextusEmpiricus, sorry about the code. You can print out intermediary results such as eqSystem to check what they are. There may be a typo (please let me know where), but I got the same paradoxical results with real data. This is where the question originated. I only use simulated data for reproducibility. Commented Apr 16 at 20:29
• I thought there was a typo when you added y[,i] to y[,i], but it was the concise coding. You used y[,i] as the error to be added to y[,i] which was confusing. (also, I am reading this from a phone so scrolling the lines is nog so easy). Commented Apr 16 at 20:31
• Yes, I saved some memory by not creating a separate matrix of epsilons. Sorry about the confusion! I am hopeful to have you on this case! It is bed time for me now, but I will come back in the morning. Commented Apr 16 at 20:32
• Somewhat related: stats.stackexchange.com/questions/645175 Commented Apr 16 at 20:34

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.

• You found a case where OLS and WLS are different, and that does not surprise me, as OLS and WLS are not the same estimator. What surprises me is the case where they are the same. I take your answer to suggest that without restrictions, point estimates and standard errors are algebraically the same in OLS and WLS. Is that true? And then what is different betwee them? Which elements of the covariance matrix $(X'\hat\Omega^{-1}X)^{-1}$ are different? Here $\Omega:=\mathbb{E}(u u')=\Sigma\bigotimes I_T$ and $\hat\Omega$ must be different between OLS and WLS. (Notation from systemfit vignette.) Commented Apr 17 at 7:34
• @RichardHardy the weights are irrelevant when there are no restrictions. Each coefficient is estimated with the same weights. Without restrictions, the weights of equation x have only influence on estimation of the parameters in equation x. Commented Apr 17 at 7:36
• I believe they must be used for obtaining the system-wide covriance matrix $(X'\hat\Omega^{-1}X)^{-1}$. However, it could happen that the part of that matrix that we see in the summary output (namely, the standard errors of individual coefficients) is not affected by the use of weights. Commented Apr 17 at 7:37
• @RichardHardy the covariance matrix is estimated regardless of the weights (and later on used in a generalized least squares fit if you are using SUR method). Commented Apr 17 at 7:38
• The weights are built in, they determine $\hat\Omega$. Am I mistaken? Commented Apr 17 at 7:39

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.

• The OLS and WLS method parameter are related to the fitting, but not to the estimation of the error where the heteroscedasticity or homoscedasticity is regulated with the control parameter. With the setting control = systemfit.control(singleEqSigma = FALSE) the error estimation will not use a seperate sigma for each seperate equation. Commented Apr 17 at 17:15
• @SextusEmpiricus, thanks! Commented Apr 17 at 17:45