How does SYSTEMFIT from R work (from statistical perspective)?

In R, there is a package for estimating systems of equations with OLS, called systemfit. Initially, this package does not seem to add that much of a value, as I can in practice estimate each equation separately by the OLS (use the function lm() multiple times)... Except that one case, when I need to make restrictions across equations, which is something I could not do by estimating equations separately... Which brings me to the point of realization that I cannot even imagine how does this package work and how to estimate multiple equations as a system.

So, to start: We have two equations A and B.

\begin{align} \text{A)} \qquad y_1 & = \alpha_0 + \alpha_1 x_1 + \alpha_2 x_2 + \alpha_3 x_3 + u \\ \text{B)} \qquad y_1 & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + v \\ \end{align}

Right now, I can estimate both equations separately via $$\boldsymbol{\hat\beta} = \left(\mathbf{X}' \mathbf{X} \right)^{-1} \mathbf{X}' \mathbf{Y}$$ or as a system by the systemfit and get the same results.

Now, if I had a single-equation assumption that $$\alpha_2 = 1$$, I could simply solve this in two ways, using the separate estimation of both equations:

Case 1: I just substitute the right value in the equation, shift it accordingly, and then estimate via OLS the new equation:

$$y_1 - 1 \cdot x_2= \alpha_0 + \alpha_1 x_1 + \alpha_3 x_3 + u$$

Case 2: Alternatively, I can come up with a matrix $$\boldsymbol{R} = [0,0,1,0]$$ of restrictions, the right hand side vector $$r = [1]$$, and then estimate the new model parameters as $$\boldsymbol{\hat\beta}_p = \boldsymbol{\hat\beta} - \left(\mathbf{X}' \mathbf{X} \right)^{-1} \boldsymbol{R}'\boldsymbol{R} \left(\mathbf{X}' \mathbf{X} \right)^{-1} \boldsymbol{R}' ( \boldsymbol{R} \boldsymbol{\hat\beta} - r )$$

Nevertheless, if I wanted to apply this to an across-equations assumption that $$\alpha_2 = \beta_2$$, I do not know how I could do that. Both case 1 and case 2 are inapplicable in that scenario. I could code this only via systemfit, which has relatively simple syntax. But I do not see inside of the code.

So, the question is: How does systemfit estimate both equations as a system and how does it adjust it for restrictions; alternatively, how can we do it in general (which equation we should use)?

• This post mischaracterizes the problem--I suspect there's a typo. Don't you mean for the left hand side of the second equation to be $y_2$ rather than $y_1$? It also implicitly is too narrow: when such equation are written, readers usually understand all error terms ("$u$" and "$v$") to be independent, but the whole point of this package is to handle models where those terms are correlated between variables. Read Section 2, "Statistical background," of the vignette.
– whuber
Commented Jul 30 at 15:35

Actually you can view the code, that's the great thing about open source software like R!

Below I will explain how. The key take-aways are:

• All source code is available on the CRAN GitHub page.
• Mathematical details can often be found in an accompanying vignette.
• Most code is input checks, error handling, etc. Try to focus on the important parts.

The package can do quite a few things, including constrained least squares, two-stage least squares, weighted two-stage least squares and three-stage least squares. I will focus on what it appears you are trying to do: Add a (set of) restriction(s) on certain parameter estimates.

I'm not too familiar with constrained least squares, but I think this is the final equation for that:

$$\begin{pmatrix} \mathbf{X'X} & \mathbf{R'} \\ \mathbf{R} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\lambda} \end{pmatrix} = \begin{pmatrix} \mathbf{X'y} \\ \mathbf{q} \end{pmatrix}$$

Finding the relevant code to understand what is going on

Just Google "r github <package name>", search the R folder and there you go.

On l.450 of the systemfit.R function, the relevant part begins:

  if( method %in% c( "OLS", "WLS", "SUR" ) ) {
if(is.null(R.restr)) {
coef <- solve( crossprod( xMatAll ), crossprod( xMatAll, yVecAll ), tol=control$solvetol ) # estimated coefficients } else { W <- .prepareWmatrix( crossprod( xMatAll ), R.restr, useMatrix = control$useMatrix )
V <- c( as.numeric( crossprod( xMatAll, yVecAll ) ), q.restr )
if( method == "OLS" || control$residCovRestricted ){ coef <- solve( W, V, tol=control$solvetol )[ 1:ncol(xMatAll) ]
} else {
coef <- solve( crossprod( xMatAll ), crossprod( xMatAll, yVecAll ),
tol = control\$solvetol )
}
}
}


Apparently, without any restrictions, it computes (as expected): $$\boldsymbol{\hat\beta} = (\mathbf{X'X})^{-1}\mathbf{X'y}$$

If using a restriction matrix $$\mathbf{R}$$, two new matrices $$\mathbf{W}$$ and $$\mathbf{V}$$ are constructed and it solves:

$$\mathbf{W} \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\lambda} \end{pmatrix} = \mathbf{V}$$

The $$\mathbf{V}$$ matrix is apparently equal to

$$\mathbf{V} = \begin{pmatrix} \mathbf{X'y} \\ \mathbf{q} \end{pmatrix},$$

where $$\mathbf{q}$$ is given by q.restr.

To delve further, we need to know what .prepareWmatrix does and what q.restr is.

Creating a block matrix $$\mathbf{W}$$ with .prepareMatrix

This function can be found on l.25 of the stackMatList.R script:

.prepareWmatrix <- function( upperleft, R.restr, useMatrix = FALSE ){
if( nrow( R.restr ) == 1 ){
lowerRows <- c( R.restr, 0 )
} else {
lowerRows <- cbind2( R.restr,
matrix( 0, nrow( R.restr ), nrow( R.restr ) ) )
}
result <- rbind2( cbind2( as.matrix( upperleft ), t(R.restr) ), lowerRows )

if( useMatrix ){
result <- as( result, "denseMatrix")
}

return( result )
}


This function takes the design matrix $$\mathbf{X}$$ and the restriction matrix $$\mathbf{R}$$ and then combines them as follows:

$$\mathbf{W} = \begin{pmatrix} \mathbf{X'X} & \mathbf{R'} \\ \mathbf{R} & \mathbf{0} \end{pmatrix}$$

Creating $$\mathbf{V}$$ using the restraints vector $$\mathbf{q}$$ (q.restr)

The q.restr part can be found on l.375 of the systemfit.R script:

q.restr <- R.restr[ , ncol( R.restr ), drop = FALSE ]
R.restr <- R.restr[ , -ncol( R.restr ), drop = FALSE ]


All it does is take the last column of R.restr. This column is then removed from R.restr itself in the next line.

(If you didn't supply a restriction matrix $$\mathbf{R}$$, it is generated in l.371 using car::makeHypothesis.)

To answer the question then, this is the system of linear equations being solved:

$$\mathbf{W} \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\lambda} \end{pmatrix} = \mathbf{V},$$

where $$\boldsymbol{\lambda}$$ is a vector of Lagrange multipliers associated with the constraints. Plugging in $$\mathbf{W}$$ and $$\mathbf{V}$$:

$$\begin{pmatrix} \mathbf{X'X} & \mathbf{R'} \\ \mathbf{R} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\lambda} \end{pmatrix} = \begin{pmatrix} \mathbf{X'y} \\ \mathbf{q} \end{pmatrix}$$

So there we have the final equation. You could inspect some example restriction matrices to see what they look like, and maybe gain some insight into why it works. Or perhaps someone more knowledgeable on constrained least squares can chime in. The vignette also contains references to literature.

• The entire point of this package is to implement a specific weighted least squares version of this model.
– whuber
Commented Jul 30 at 15:37
• @whuber Thanks for your comment. I read the question as being primarily interested in how the constraints are applied. I added a note at the top to clarify that this is not the only thing the package does. Commented Jul 30 at 15:48