7 Programming IV and GMM Estimators in R

Chapter 6 developed the identification conditions, moment restrictions, consistency arguments, and limiting distributions of IV and GMM estimators for spatial models. This chapter turns those mathematical objects into transparent R programs.

The central principle is the same as in Chapter 5. Specialized software is used as a verification tool, not as a substitute for understanding. For each estimator, we will therefore proceed in the same order:

write the spatial model and identify the endogenous component;
construct the residual or innovation implied by a candidate parameter vector;
translate the population restrictions into sample moments;
assemble the estimator or numerical criterion;
construct the covariance matrix from the same objects;
inspect intermediate matrices and compare the result with established R software whenever a comparable implementation is available.

The chapter follows the order of Chapter 6:

spatial two-stage least squares for the SLM;
best spatial two-stage least squares;
GMM and optimally weighted GMM for the SLM;
feasible GLS for the SEM;
GMM and optimally weighted GMM for the SEM; and
feasible generalized spatial two-stage least squares for the SAC model.

Three objects that must not be confused

Throughout the chapter, keep the following distinction in mind:

the structural residual is obtained from the regression equation;
the innovation is obtained after applying the spatial-error transformation when an SEM or SAC model is present;
the moment vector combines the innovation with instruments or quadratic matrices.

Most programming errors in IV and GMM estimators arise from using the wrong one of these three objects at a particular step.

7.1 Packages Used in This Chapter

Four packages have distinct roles in this chapter.

spdep constructs neighbor lists and spatial-weight objects and converts them to the dense matrix representation used in our programs. spatialreg provides benchmark estimators such as stsls(), GMerrorsar(), and lagsarlm(). maxLik supplies the BFGS optimizer used for the SLM and SEM GMM criteria. Since maxLik maximizes rather than minimizes, our criterion functions return the negative of the GMM objective. Finally, numDeriv is used only to verify the analytical gradients numerically.

# Construction and manipulation of spatial-neighbor and weight objects.
library(spdep)

# Benchmark IV, GMM, and maximum-likelihood estimators.
library(spatialreg)

# Numerical maximization of the negative GMM criteria.
library(maxLik)

# Numerical verification of analytical gradients.
library(numDeriv)

The final comparison with sphet is conditional because that package is not needed by the estimators programmed here.

As in the maximum-likelihood programming chapter, dense matrices are used on purpose. They make projections, spatial transformations, quadratic forms, and matrix derivatives visible. For large samples, production implementations should exploit sparse matrices and avoid explicitly forming an \(n\times n\) projection matrix.

7.1.1 R operations used repeatedly

A few R expressions recur throughout the chapter:

A %*% B performs matrix multiplication;
crossprod(A, B) computes t(A) %*% B without explicitly constructing the transpose;
drop() converts a one-column matrix into a vector or a \(1\times1\) matrix into a scalar;
solve(A, b) solves \(\mathbf A\mathbf x=\mathbf b\) and is preferable to constructing solve(A) %*% b when only the solution is needed;
diag(n) creates \(\mathbf I_n\), whereas diag(v) places the elements of a vector on the diagonal of a matrix;
lm.fit() performs the numerical least-squares calculation without the additional formula-processing work of lm();
nlminb() minimizes a scalar criterion subject to box constraints; and
maxLik() maximizes a scalar objective and can use a gradient stored as the "gradient" attribute of that objective.

The user-facing estimators also share a small set of arguments:

Argument	Role
`formula`	specifies the dependent variable and exogenous regressors
`data`	supplies the variables appearing in `formula`
`W`	weights matrix for the spatial lag of the dependent variable
`M`	weights matrix for the spatial autoregressive error process
`instruments`	highest order used in \(\mathbf W\mathbf X,\ldots,\mathbf W^l\mathbf X\)
`estimator`	selects one-step GMM or feasible optimally weighted GMM
`gradient`	chooses analytical or numerical derivatives in `maxLik()`
`verbose`	prints values of a nonlinear criterion during optimization

7.2 Common Helper Functions

Several estimators reuse the same operations. We define them once so that the later functions can focus on the econometric step being implemented.

The most important helper is make.H(). Given an \(n\times n\) spatial weights matrix \(\mathbf W\), an \(n\times k\) regressor matrix \(\mathbf X\), and a positive integer \(l\), it constructs \[ \left[ \mathbf W\mathbf X_{-0}, \mathbf W^2\mathbf X_{-0}, \ldots, \mathbf W^l\mathbf X_{-0} \right], \] where \(\mathbf X_{-0}\) excludes the intercept. With row-standardized weights, \(\mathbf W\boldsymbol\iota_n=\boldsymbol\iota_n\), so spatially lagging the intercept would merely duplicate a column already present in \(\mathbf X\).

The helper section also performs four tasks:

it prevents missing observations from being silently removed without the corresponding rows and columns being removed from the weights matrix;
it checks matrix dimensions and the rank of the design matrix;
it computes an admissible interval for a spatial autoregressive parameter when a matrix inverse is required;
it protects starting values based on sample correlations from being missing or outside that interval.

# Validate a spatial matrix before it enters an estimator.
validate_spatial_matrix <- function(A, n, matrix_name = "W",
                                    diagonal_tolerance = 1e-12) {
  A <- as.matrix(A)
  if (!is.numeric(A)) stop(matrix_name, " must be numeric.", call. = FALSE)
  if (nrow(A) != n || ncol(A) != n) {
    stop(matrix_name, 
         " must be an n by n matrix conformable with the model data.",
         call. = FALSE
    )
  }
  if (anyNA(A) || any(!is.finite(A))) {
    stop(matrix_name, " must contain only finite, nonmissing values.", call. = FALSE)
  }
  if (max(abs(diag(A))) > diagonal_tolerance) {
    stop(matrix_name, " must have a zero diagonal.", call. = FALSE)
  }

  A
}

# Construct and validate the response, design matrix, and one spatial matrix.
spatial_model_components <- function(formula, data, A, matrix_name = "W") {
  model_frame <- model.frame(formula = formula, data = data, 
                             na.action = na.fail)
  y <- model.response(model_frame)
  if (is.matrix(y) && ncol(y) != 1L) {
    stop("The response must be univariate.", call. = FALSE)
  }
  if (!is.numeric(y)) {
    stop("The response must be numeric.", call. = FALSE)
  }
  y <- as.numeric(y)
  X <- model.matrix(formula, data = model_frame)
  n <- length(y)
  k <- ncol(X)
  A <- validate_spatial_matrix(A, n = n, matrix_name = matrix_name)

  if (anyNA(y) || anyNA(X) || any(!is.finite(y)) || any(!is.finite(X))) {
    stop("The response and regressors must be finite and nonmissing.", call. = FALSE)
  }
  if (qr(X)$rank < k) {
    stop("The model matrix X must have full column rank.", call. = FALSE)
  }

  list(
    model.frame = model_frame,
    y = y,
    X = X,
    A = A,
    n = n,
    k = k
  )
}

# Construct WX, W^2 X, ..., W^l X after omitting the intercept column.
make.H <- function(W, X, l = 3L) {
  W <- as.matrix(W)
  X <- as.matrix(X)

  if (!is.numeric(W) || !is.numeric(X)) {
    stop("W and X must be numeric.", call. = FALSE)
  }
  if (nrow(W) != ncol(W)) {
    stop("W must be a square matrix.", call. = FALSE)
  }
  if (nrow(W) != nrow(X)) {
    stop("W and X must contain the same number of spatial units.", call. = FALSE)
  }
  if (length(l) != 1L || !is.finite(l) || l < 1 || l != as.integer(l)) {
    stop("l must be a positive integer.", call. = FALSE)
  }
  
  l <- as.integer(l)
  x_names <- colnames(X)
  if (is.null(x_names)) {
    x_names <- paste0("x", seq_len(ncol(X)))
    colnames(X) <- x_names
  }

  intercept <- which(x_names == "(Intercept)")
  if (length(intercept) > 0L) {
    X_lag <- X[, -intercept, drop = FALSE]
    lag_names <- x_names[-intercept]
  } else {
    X_lag <- X
    lag_names <- x_names
  }

  if (ncol(X_lag) == 0L) {
    return(matrix(numeric(0), nrow = nrow(X), ncol = 0L))
  }

  H <- matrix(NA_real_, nrow = nrow(X), ncol = ncol(X_lag) * l)
  H_names <- character(ncol(H))
  current <- X_lag

  for (order in seq_len(l)) {
    columns <- ((order - 1L) * ncol(X_lag) + 1L):(order * ncol(X_lag))
    current <- W %*% current
    H[, columns] <- current
    prefix <- if (order == 1L) "W*" else paste0("W^", order, "*")
    H_names[columns] <- paste0(prefix, lag_names)
  }
  colnames(H) <- H_names
  H
}

# Retain a maximal linearly independent subset while preserving column order.
independent_columns <- function(H, tolerance = 1e-10) {
  H <- as.matrix(H)
  if (ncol(H) == 0L) return(H)
  qr_H <- qr(H, tol = tolerance)
  selected <- sort(qr_H$pivot[seq_len(qr_H$rank)])
  H[, selected, drop = FALSE]
}

# Matrix trace used by the quadratic moment functions.
tr <- function(A) sum(diag(A))

# Compute an interior interval between the nearest singular values of I - aA.
spatial_parameter_interval <- function(A, tolerance = 1e-10,
                                       margin = sqrt(.Machine$double.eps)) {
  eigenvalues <- eigen(A, only.values = TRUE)$values
  if (max(abs(Im(eigenvalues))) > tolerance) {
    stop(
      "This implementation requires a spatial matrix with a numerically real spectrum.",
      call. = FALSE
    )
  }

  eigenvalues <- Re(eigenvalues)
  negative <- eigenvalues[eigenvalues < -tolerance]
  positive <- eigenvalues[eigenvalues > tolerance]

  if (length(negative) == 0L || length(positive) == 0L) {
    stop("The spectrum must contain negative and positive eigenvalues.",
         call. = FALSE
    )
  }

  raw_interval <- c(lower = 1 / min(negative), upper = 1 / max(positive))
  buffer <- margin * max(1, max(abs(raw_interval)))
  interval <- c(
    lower = unname(raw_interval["lower"]) + buffer,
    upper = unname(raw_interval["upper"]) - buffer
  )

  if (interval["lower"] >= interval["upper"]) {
    stop("The admissible spatial-parameter interval is empty.", call. = FALSE)
  }

  interval
}

# Replace an undefined correlation by zero.
safe_correlation <- function(x, y, default = 0) {
  value <- suppressWarnings(cor(x, y))
  if (length(value) != 1L || !is.finite(value)) default else value
}

# Move a scalar starting value strictly inside a numerical interval.
clamp_to_interval <- function(value, interval, 
                              margin = sqrt(.Machine$double.eps)) {
  lower <- unname(interval["lower"])
  upper <- unname(interval["upper"])
  buffer <- margin * max(1, upper - lower)
  lower_interior <- lower + buffer
  upper_interior <- upper - buffer

  if (lower_interior >= upper_interior) {
    stop("The numerical interval has no usable interior.", call. = FALSE)
  }

  unname(min(max(value, lower_interior), upper_interior))
}

# Invert an estimated covariance matrix only when it is numerically positive definite.
invert_moment_covariance <- function(A, 
                                     matrix_name = "The moment covariance matrix",
                                     tolerance = 1e-10) {
  A <- (as.matrix(A) + t(as.matrix(A))) / 2
  eigenvalues <- eigen(A, symmetric = TRUE, only.values = TRUE)$values
  threshold  <- tolerance * max(1, max(abs(eigenvalues)))

  if (min(eigenvalues) <= threshold) {
    stop(
      matrix_name,
      " is not numerically positive definite and cannot be inverted.",
      call. = FALSE
    )
  }
  chol2inv(chol(A))
}

# Encode lower and upper bounds for one parameter in maxLik's C theta + d >= 0 form.
spatial_parameter_constraints <- function(number_parameters, parameter_position,
                                          interval) {
  C <- matrix(0, nrow = 2L, ncol = number_parameters)
  C[1L, parameter_position] <- 1
  C[2L, parameter_position] <- -1
  list(ineqA = C, ineqB = c(-interval["lower"], interval["upper"]))
}

The helper spatial_parameter_interval() is intentionally conservative. It is used only in examples whose row-standardized contiguity matrices have a numerically real spectrum. The function does not claim that every spatial weights matrix has real eigenvalues.

7.3 Spatial Two-Stage Least Squares for the SLM

The S2SLS estimator and its limiting variance are developed in Section 6.2. We first reproduce the artificial SLM used in the original notes.

The spatial lag model can be written as \[ \mathbf y = \rho\mathbf W\mathbf y + \mathbf X\boldsymbol\beta + \boldsymbol\varepsilon. \]

The regressor \(\mathbf W\mathbf y\) is endogenous because it contains spatially propagated innovations. Define \[ \mathbf Z = \left[ \mathbf X, \mathbf W\mathbf y \right]. \]

S2SLS replaces the endogenous column of \(\mathbf Z\) by its projection on instruments constructed from \(\mathbf X\), \(\mathbf W\mathbf X\), and higher-order spatial lags. The artificial example below keeps the data-generating process small enough to inspect every object.

7.3.1 Artificial SLM Data

# Use a 23 by 23 regular lattice, giving n = 529 spatial units.
set.seed(1986)
side_slm <- 23L
n_slm <- side_slm^2

# Rook contiguity and row-standardized spatial weights.
rho_slm <- 0.6
neighbors_slm <- cell2nb(side_slm, side_slm, type = "rook")
W_slm <- nb2mat(neighbors_slm, style = "W")

# Simulate three exogenous regressors.
x1_slm <- rnorm(n_slm)
x2_slm <- rnorm(n_slm)
x3_slm <- rnorm(n_slm)

# True structural parameters and homoskedastic innovations.
beta_slm <- c("(Intercept)" = 0, "x1" = -1, "x2" = 0, "x3" = 1)
sigma2_slm <- 2
epsilon_slm <- rnorm(n_slm, mean = 0, sd = sqrt(sigma2_slm))

# Construct the design matrix explicitly so its columns match beta_slm.
X_slm <- cbind(
  "(Intercept)" = 1,
  "x1" = x1_slm,
  "x2" = x2_slm,
  "x3" = x3_slm
)

# Solve (I - rho W)y = X beta + epsilon for y.
S_slm <- diag(n_slm) - rho_slm * W_slm
y_slm <- drop(solve(S_slm, X_slm %*% beta_slm + epsilon_slm))

data_slm <- data.frame(
  y = y_slm,
  x1 = x1_slm,
  x2 = x2_slm,
  x3 = x3_slm
)

The following calculation displays the instruments through order three.

H_example <- make.H(W = W_slm, X = X_slm, l = 3L)

# The first three columns are WX, the next three are W^2X,
# and the final three are W^3X.
dim(H_example)

[1] 529   9

head(H_example)

           W*x1       W*x2        W*x3      W^2*x1      W^2*x2      W^2*x3
[1,] -0.1838084  1.4871684  1.06237375 -0.19093708 -0.66721937  0.11419352
[2,] -0.1734224 -0.4801317  0.00889001 -0.18346010  1.16724914  0.84907535
[3,] -0.1376933  0.7858546  0.65742663 -0.07600407 -0.06559888 -0.33095264
[4,]  0.2516524  0.5410386 -0.66094372 -0.17572619  0.04624316  0.12822625
[5,] -0.2524966 -0.7357899 -0.27611007  0.53392070  0.46568956 -0.40883400
[6,]  0.7127103  0.4896574 -0.25597660 -0.08935028 -0.65274040 -0.05759884
         W^3*x1     W^3*x2      W^3*x3
[1,] -0.2364339  1.1506372  0.80406622
[2,] -0.1724420 -0.4360736 -0.07942842
[3,] -0.1916439  0.6032002  0.46356720
[4,]  0.2365006  0.1309099 -0.37543431
[5,] -0.1367966 -0.4243089 -0.02981086
[6,]  0.4918535  0.2858139 -0.35198454

The argument l in make.H() determines the highest spatial order. In this example, l = 3L produces nine excluded instrument columns because there are three nonconstant regressors. The full instrument matrix used by S2SLS also contains the original columns of \(\mathbf X\).

Before writing a wrapper, it is useful to map the estimator to the objects that will appear in the code:

Mathematical object	R object	Dimension
\(\mathbf y\)	`y`	\(n\times1\)
\(\mathbf X\)	`X`	\(n\times k\)
\(\mathbf Z=[\mathbf X,\mathbf W\mathbf y]\)	`Z`	\(n\times(k+1)\)
\(\mathbf H\)	`H`	\(n\times p\)
\(\mathbf P_H\)	`P_H`	\(n\times n\)
\(\widehat{\mathbf Z}=\mathbf P_H\mathbf Z\)	`Z_hat`	\(n\times(k+1)\)

7.3.2 Programming the S2SLS Estimator

The function below implements the estimator in Equation 6.14. It stores the matrices needed by the variance estimators and by subsequent procedures.

# S2SLS Estimator
slm.2sls <- function(formula, data, W, instruments = 2L) {
  call <- match.call(expand.dots = TRUE)
  components <- spatial_model_components(formula = formula, data = data, A = W,
                                         matrix_name = "W")

  y <- components$y
  X <- components$X
  W <- components$A
  n <- components$n

  # The only endogenous regressor is Wy.
  Wy <- drop(W %*% y)
  Z <- cbind(X, "Wy" = Wy)

  if (n <= ncol(Z)) {
    stop("The sample must contain more observations than structural parameters.")
  }

  # H contains X and the excluded spatial lags of the nonconstant regressors.
  H <- cbind(X, make.H(W = W, X = X, l = instruments))
  H <- independent_columns(H)

  if (ncol(H) < ncol(Z)) {
    stop("The model is underidentified: H has fewer columns than Z.", call. = FALSE)
  }
  if (qr(crossprod(H, Z))$rank < ncol(Z)) {
    stop("The sample first-stage rank condition fails.", call. = FALSE)
  }

  # Projection onto the column space of H.
  HH    <- crossprod(H)
  P_H   <- H %*% solve(HH, t(H))
  Z_hat <- P_H %*% Z

  # S2SLS = (Z_hat' Z_hat)^(-1) Z_hat' y.
  coefficients <- drop(solve(crossprod(Z_hat), crossprod(Z_hat, y)))
  names(coefficients) <- colnames(Z)
  fitted    <- drop(Z %*% coefficients)
  residuals <- drop(y - fitted)

  structure(
    list(
      coefficients = coefficients,
      call = call,
      X = X,
      H = H,
      Z = Z,
      y = y,
      W = W,
      P_H = P_H,
      fitted.values = fitted,
      residuals = residuals,
      first.stage.rank = qr(crossprod(H, Z))$rank
    ),
    class = "mys2sls"
  )
}

The function follows the two stages literally. First, P_H %*% Z constructs the fitted values from regressing every column of \(\mathbf Z\) on \(\mathbf H\). The exogenous columns of \(\mathbf X\) reproduce themselves because they are already included in \(\mathbf H\); the substantive first-stage operation is the projection of \(\mathbf W\mathbf y\).

The rank checks distinguish two different failures. Having fewer instruments than regressors is underidentification by construction. Having enough columns but a rank-deficient \(\mathbf H^\top\mathbf Z\) is a sample identification failure caused by redundant or uninformative instruments.

Now, we create the S3 method vcov() for our class mys2sls. The homoskedastic covariance estimator corresponds to Equation 6.25. The heteroskedasticity-robust version uses Equation 6.26 and Equation 6.24.

# vcov function 
vcov.mys2sls <- function(object, type = c("homoskedastic", "robust"), ...) {
  type <- match.arg(type)
  Z    <- object$Z
  H    <- object$H
  residuals <- object$residuals

  n <- nrow(Z)
  k <- ncol(Z)
  degrees_freedom <- n - k

  Q_HZ <- crossprod(H, Z) / n
  Q_HH <- crossprod(H) / n
  Q_HH_inverse <- solve(Q_HH)

  bread <- solve(t(Q_HZ) %*% Q_HH_inverse %*% Q_HZ)
  if (type == "homoskedastic") {
    sigma2_hat <- sum(residuals^2) / degrees_freedom
    return(sigma2_hat * bread / n)
  }

  H_residuals<- H * residuals
  Delta_hat  <- crossprod(H_residuals) / n
  middle     <- t(Q_HZ) %*% Q_HH_inverse %*% Delta_hat %*% Q_HH_inverse %*% Q_HZ

  bread %*% middle %*% bread / n
}

# S3 Methods for my2sls
summary.mys2sls <- function(object, type = c("homoskedastic", "robust"),
                            digits = max(3L, getOption("digits") - 3L), ...) {
  type            <- match.arg(type)
  estimates       <- drop(object$coefficients)
  standard_errors <- sqrt(diag(vcov(object, type = type)))
  statistics      <- estimates / standard_errors
  p_values        <- 2 * pnorm(-abs(statistics))

  coefficient_table <- cbind(
    "Estimate" = estimates,
    "Std. Error" = standard_errors,
    "z value" = statistics,
    "Pr(>|z|)" = p_values
  )

  structure(
    list(
      coefficients = coefficient_table,
      type = type,
      digits = digits,
      call = object$call
    ),
    class = "summary.mys2sls"
  )
}

print.summary.mys2sls <- function(x, digits = x$digits, ...) {
  cat("\nCall:\n")
  cat(paste(deparse(x$call), collapse = "\n"), "\n", sep = "")
  cat("\nVariance estimator:", x$type, "\n")
  cat("\nCoefficients:\n")
  printCoefmat(
    x$coefficients,
    digits = digits,
    P.values = TRUE,
    has.Pvalue = TRUE
  )
  invisible(x)
}

The covariance method reuses the same normalized cross-products that appear in Section 6.2. Under homoskedasticity, the only additional object is the residual variance. Under heteroskedasticity, the matrix \[ \widehat{\boldsymbol\Delta} = \frac{1}{n} \sum_{i=1}^n \widehat\varepsilon_i^2 \mathbf h_i\mathbf h_i^\top \] is computed by multiplying each row of H by the corresponding residual and then taking a cross-product.

The summary() method reports asymptotic normal tests. These are labeled as z value, rather than finite-sample \(t\) tests, because the spatial IV results in Chapter 6 are asymptotic.

7.3.3 Estimation and Package Verification

Now, we use our function fit_s2sls() and compare the results using spatialreg::stsls()

# Cache revision: recompute after the estimator and DGP audit.

# Estimate the model using our function
fit_s2sls <- slm.2sls(y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
                      instruments = 2L)
summary(fit_s2sls)


Call:
slm.2sls(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L)

Variance estimator: homoskedastic 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.06983    0.06656  -1.049    0.294    
x1          -1.00574    0.06347 -15.846   <2e-16 ***
x2          -0.04353    0.06276  -0.694    0.488    
x3           1.01485    0.06720  15.102   <2e-16 ***
Wy           0.64731    0.06222  10.404   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_s2sls, type = "robust")


Call:
slm.2sls(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L)

Variance estimator: robust 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.06983    0.06682  -1.045    0.296    
x1          -1.00574    0.06587 -15.269   <2e-16 ***
x2          -0.04353    0.06290  -0.692    0.489    
x3           1.01485    0.06669  15.217   <2e-16 ***
Wy           0.64731    0.05892  10.986   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We verify the implementation using spatialreg::stsls().

listw_slm <- mat2listw(W_slm, style = "W")

fit_stsls_homoskedastic <- stsls(y ~ x1 + x2 + x3, data = data_slm, 
                                 listw = listw_slm)
fit_stsls_robust <- stsls(y ~ x1 + x2 + x3, data = data_slm, listw = listw_slm,
                          robust = TRUE, HC = "HC0")

# Extract coefficients and standard errors by name rather than by position.
extract_stsls_result <- function(object, x_names) {
  estimates <- coef(object)
  rho_name <- intersect(c("Rho", "rho", "Rho_Wy"), names(estimates))

  if (length(rho_name) != 1L) {
    stop("Could not identify the spatial-lag coefficient in the stsls object.")
  }

  standard_errors <- sqrt(diag(object$var))
  list(
    coefficients = c(estimates[x_names], "rho" = estimates[rho_name]),
    standard.errors = c(standard_errors[x_names], "rho" = standard_errors[rho_name])
  )
}

benchmark_homoskedastic <- extract_stsls_result(
  fit_stsls_homoskedastic,
  x_names = colnames(X_slm)
)
benchmark_robust <- extract_stsls_result(
  fit_stsls_robust,
  x_names = colnames(X_slm)
)

cbind(
  "spatialreg" = benchmark_homoskedastic$coefficients,
  "custom" = fit_s2sls$coefficients
)

             spatialreg      custom
(Intercept) -0.06983166 -0.06983166
x1          -1.00573883 -1.00573883
x2          -0.04352912 -0.04352912
x3           1.01485440  1.01485440
rho.Rho      0.64730746  0.64730746

cbind(
  "spatialreg-homoskedastic" = benchmark_homoskedastic$standard.errors,
  "custom-homoskedastic" = sqrt(diag(vcov(fit_s2sls))),
  "spatialreg-robust" = benchmark_robust$standard.errors,
  "custom-robust" = sqrt(diag(vcov(fit_s2sls, type = "robust")))
)

            spatialreg-homoskedastic custom-homoskedastic spatialreg-robust
(Intercept)               0.06655890           0.06655890        0.06681918
x1                        0.06346892           0.06346892        0.06586903
x2                        0.06276194           0.06276194        0.06290024
x3                        0.06720175           0.06720175        0.06669244
rho.Rho                   0.06221908           0.06221908        0.05892332
            custom-robust
(Intercept)    0.06681918
x1             0.06586903
x2             0.06290024
x3             0.06669244
rho.Rho        0.05892332

7.4 Best Spatial Two-Stage Least Squares

The BS2SLS estimator uses the estimated optimal instrument described in Section 6.2.5 and Equation 6.27. The initial estimator is the S2SLS estimator programmed above.

The best instrument replaces the generic spatial lags by an estimate of \[ \mathbb E(\mathbf W\mathbf y\mid\mathbf X) = \mathbf W (\mathbf I-\rho\mathbf W)^{-1} \mathbf X\boldsymbol\beta. \] This produces a just-identified instrument matrix \[ \mathbf H^* = \left[ \mathbf X, \widehat{\mathbb E}(\mathbf W\mathbf y\mid\mathbf X) \right]. \]

The program therefore has two estimation rounds: ordinary S2SLS supplies consistent preliminary estimates, and those estimates construct the final instrument.

The inverse \((\mathbf I-\widetilde\rho\mathbf W)^{-1}\) requires an admissible preliminary spatial parameter. A finite-sample S2SLS estimate can fall outside the invertible interval even when it is consistent asymptotically. The function below therefore moves such an estimate to the nearest interior boundary and issues a warning. This safeguard affects only the construction of the feasible instrument.

# BS2SLS estimator
slm.b2sls <- function(formula, data, W, instruments = 2L) {
  initial_fit <- slm.2sls(formula = formula, data = data,
                          W = W, instruments = instruments)
  X <- initial_fit$X
  Z <- initial_fit$Z
  y <- initial_fit$y
  W <- initial_fit$W
  n <- nrow(X)
  k <- ncol(X)

  initial_coefficients <- initial_fit$coefficients
  beta_initial <- initial_coefficients[seq_len(k)]
  rho_initial_raw <- unname(initial_coefficients[k + 1L])

  rho_interval <- spatial_parameter_interval(W)
  rho_initial  <- clamp_to_interval(rho_initial_raw, rho_interval)
  if (!isTRUE(all.equal(rho_initial, rho_initial_raw))) {
    warning("The preliminary S2SLS rho was moved inside the invertible interval.",
            call. = FALSE)
  }

  S_initial <- diag(n) - rho_initial * W
  optimal_spatial_instrument <- drop(W %*% solve(S_initial, X %*% beta_initial))
  H_star <- cbind(X, "E[Wy|X]" = optimal_spatial_instrument)

  HZ_star <- crossprod(H_star, Z)
  if (qr(HZ_star)$rank < ncol(Z)) {
    stop("The feasible optimal instrument matrix does not identify the model.")
  }

  coefficients        <- drop(solve(HZ_star, crossprod(H_star, y)))
  names(coefficients) <- colnames(Z)
  fitted              <- drop(Z %*% coefficients)
  residuals          <- drop(y - fitted)

  structure(
    list(
      coefficients = coefficients,
      coefficients.initial = initial_coefficients,
      call = match.call(),
      X = X,
      Z = Z,
      y = y,
      W = W,
      H.initial = initial_fit$H,
      H.star = H_star,
      rho.initial.raw = rho_initial_raw,
      rho.initial.instrument = rho_initial,
      rho.interval = rho_interval,
      residuals.initial = initial_fit$residuals,
      residuals = residuals,
      initial.fit = initial_fit
    ),
    class = "mybs2sls"
  )
}

The initial covariance matrix is the homoskedastic S2SLS covariance matrix. The final covariance matrix uses the estimated optimal instruments and the limiting variance in Equation 6.34.

# vcov function for mybs2sls
vcov.mybs2sls <- function(object, estimate = c("initial", "final"), ...) {
  estimate <- match.arg(estimate)
  if (estimate == "initial") return(vcov(object$initial.fit, type = "homoskedastic"))
  n <- nrow(object$X)
  number_parameters <- ncol(object$Z)
  degrees_freedom <- n - number_parameters
  sigma2_hat <- sum(object$residuals^2) / degrees_freedom
  drop(sigma2_hat) * solve(crossprod(object$H.star) / n) / n
}

# S3 methods for summary 
summary.mybs2sls <- function(object, estimate = c("initial", "final"),
                             digits = max(3L, getOption("digits") - 3L),
                             ...) {
  estimate <- match.arg(estimate)
  estimates <- if (estimate == "initial") {
    drop(object$coefficients.initial)
  } else {
    drop(object$coefficients)
  }

  standard_errors <- sqrt(diag(vcov(object, estimate = estimate)))
  statistics <- estimates / standard_errors
  p_values <- 2 * pnorm(-abs(statistics))

  coefficient_table <- cbind(
    "Estimate" = estimates,
    "Std. Error" = standard_errors,
    "z value" = statistics,
    "Pr(>|z|)" = p_values
  )

  structure(
    list(
      coefficients = coefficient_table,
      estimate = estimate,
      digits = digits,
      call = object$call
    ),
    class = "summary.mybs2sls"
  )
}

print.summary.mybs2sls <- function(x, digits = x$digits, ...) {
  cat("\nCall:\n")
  cat(paste(deparse(x$call), collapse = "\n"), "\n", sep = "")
  cat("\nEstimation round:", x$estimate, "\n")
  cat("\nCoefficients:\n")
  printCoefmat(
    x$coefficients,
    digits = digits,
    P.values = TRUE,
    has.Pvalue = TRUE
  )
  invisible(x)
}

Because \(\mathbf H^*\) has exactly as many columns as \(\mathbf Z\), the final estimator uses the just-identified IV formula rather than a projection matrix. The final covariance formula exploits the optimal-instrument result in Equation 6.34. It is not a generic heteroskedasticity-robust formula.

# Check our function
fit_bs2sls <- slm.b2sls(y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
                        instruments = 2L)
summary(fit_bs2sls, estimate = "initial")


Call:
slm.b2sls(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L)

Estimation round: initial 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.06983    0.06656  -1.049    0.294    
x1          -1.00574    0.06347 -15.846   <2e-16 ***
x2          -0.04353    0.06276  -0.694    0.488    
x3           1.01485    0.06720  15.102   <2e-16 ***
Wy           0.64731    0.06222  10.404   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_bs2sls, estimate = "final")


Call:
slm.b2sls(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L)

Estimation round: final 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.08156    0.06604  -1.235    0.217    
x1          -1.01266    0.06417 -15.781   <2e-16 ***
x2          -0.04195    0.06291  -0.667    0.505    
x3           1.02569    0.06655  15.413   <2e-16 ***
Wy           0.61460    0.05803  10.591   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The first summary reproduces the preliminary S2SLS round. The second summary uses the feasible optimal instrument. Reporting both rounds is useful because a surprising final estimate can then be traced either to the preliminary estimator or to the construction of \(\mathbf H^*\).

7.5 GMM and OGMM for the SLM

S2SLS uses only the linear restrictions \(\mathbb E(\mathbf H^\top\boldsymbol\varepsilon)=\mathbf0\). GMM adds quadratic restrictions that contain information about the spatial parameter. For a candidate vector \[ \boldsymbol\theta = \left( \rho, \boldsymbol\beta^\top \right)^\top, \] the structural innovation is \[ \boldsymbol\varepsilon(\boldsymbol\theta) = (\mathbf I-\rho\mathbf W)\mathbf y - \mathbf X\boldsymbol\beta. \]

We retain the two quadratic matrices, \[ \mathbf P_1=\mathbf W, \qquad \mathbf P_2=\mathbf W^2- \frac{\operatorname{tr}(\mathbf W^2)}{n}\mathbf I_n. \]

The second matrix is centered so that its trace is zero. Together with the linear instruments, the normalized sample moment vector is \[ \overline{\mathbf g}_n(\boldsymbol\theta) = \frac{1}{n} \begin{pmatrix} \boldsymbol\varepsilon^\top\mathbf P_1\boldsymbol\varepsilon\\ \boldsymbol\varepsilon^\top\mathbf P_2\boldsymbol\varepsilon\\ \mathbf H^\top\boldsymbol\varepsilon \end{pmatrix}. \]

The parameter order is important: the SLM GMM functions store rho first and the regression coefficients afterward.

7.5.1 Moment Function and Criterion

The function below evaluates the normalized sample moments and their analytical Jacobian. Recall from Chapter 6 that the of \(\partial\boldsymbol\varepsilon(\boldsymbol\theta)\) derivative with respect to the parameter vector is \[ \frac{\partial\boldsymbol\varepsilon(\boldsymbol\theta)} {\partial\boldsymbol\theta^\top} = \begin{pmatrix} -\mathbf W\mathbf y & -\mathbf X \end{pmatrix}. \]

For the quadratic moment \[ \boldsymbol\varepsilon(\boldsymbol\theta)^\top \mathbf P_j \boldsymbol\varepsilon(\boldsymbol\theta), \] the residual vector appears on both sides of \(\mathbf P_j\). Its derivative therefore depends on \[ \mathbf P_j^s = \mathbf P_j+\mathbf P_j^\top. \]

Using the same notation as in Chapter 6, the Jacobian of the complete moment vector is \[ \frac{\partial\mathbf g_n(\boldsymbol\theta)} {\partial\boldsymbol\theta^\top} = \begin{pmatrix} \boldsymbol\varepsilon(\boldsymbol\theta)^\top\mathbf P_1^s\\ \boldsymbol\varepsilon(\boldsymbol\theta)^\top\mathbf P_2^s\\ \mathbf H^\top \end{pmatrix} \begin{pmatrix} -\mathbf W\mathbf y & -\mathbf X \end{pmatrix}. \]

The objects P1_symmetric and P2_symmetric in the code correspond to \(\mathbf P_1^s\) and \(\mathbf P_2^s\). The object derivative_epsilon corresponds to \[ \frac{\partial\boldsymbol\varepsilon(\boldsymbol\theta)} {\partial\boldsymbol\theta^\top}. \]

Writing the code in this form makes the correspondence between the theoretical derivation and its implementation explicit.

# Function to generate moments and their derivatives
moments.lee2007 <- function(theta, y, X, H, W) {
  X <- as.matrix(X)
  H <- as.matrix(H)
  W <- as.matrix(W)
  y <- as.numeric(y)

  k <- ncol(X)
  n <- nrow(X)
  if (length(theta) != k + 1L) {
    stop("theta must contain rho followed by the k regression coefficients.")
  }

  rho <- theta[1L]
  beta <- theta[2L:(k + 1L)]

  S       <- diag(n) - rho * W
  epsilon <- drop(S %*% y - X %*% beta)

  P1 <- W
  W2 <- W %*% W
  P2 <- W2 - (tr(W2) / n) * diag(n)

  moments_linear     <- crossprod(H, epsilon)
  moment_quadratic_1 <- drop(crossprod(epsilon, P1 %*% epsilon))
  moment_quadratic_2 <- drop(crossprod(epsilon, P2 %*% epsilon))
  moments <- rbind(moment_quadratic_1, moment_quadratic_2, moments_linear)
  rownames(moments) <- c("q1", "q2", colnames(H))

  P1_symmetric <- P1 + t(P1)
  P2_symmetric <- P2 + t(P2)
  derivative_epsilon <- cbind("rho" = W %*% y, X)

  derivative <- -rbind(
    crossprod(epsilon, P1_symmetric),
    crossprod(epsilon, P2_symmetric),
    t(H)
  ) %*% derivative_epsilon

  list(
    g = moments / n,
    D = derivative / n,
    epsilon = epsilon,
    P1 = P1,
    P2 = P2
  )
}

# Function to be minimized
Qmin.slm <- function(theta, y, X, H, W, Psi, gradient = TRUE) {
  moment_fit <- moments.lee2007(theta, y = y, X = X, H = H, W = W)

  # maxLik maximizes, so return the negative GMM criterion.
  criterion <- -drop(crossprod(moment_fit$g, Psi %*% moment_fit$g))

  if (gradient) {
    criterion_gradient <- -2 * crossprod(moment_fit$D, Psi %*% moment_fit$g)
    attr(criterion, "gradient") <- drop(criterion_gradient)
  }

  criterion
}

The GMM criterion is \[ Q_n(\boldsymbol\theta) = \overline{\mathbf g}_n(\boldsymbol\theta)^\top \boldsymbol\Psi_n \overline{\mathbf g}_n(\boldsymbol\theta). \] maxLik maximizes its objective, so Qmin.slm() returns \(-Q_n\). The analytical gradient attached to that scalar is \[ -2\mathbf D_n(\boldsymbol\theta)^\top \boldsymbol\Psi_n \overline{\mathbf g}_n(\boldsymbol\theta). \] Keeping the sign convention explicit prevents a common bug: supplying the gradient of \(Q_n\) while asking the optimizer to maximize \(-Q_n\).

The next function estimates the covariance matrix of the normalized moments in Equation 6.46.

# vcov of moments
make.vmom <- function(coefficients, y, X, H, W) {
  X <- as.matrix(X)
  H <- as.matrix(H)
  W <- as.matrix(W)
  y <- as.numeric(y)

  k <- ncol(X)
  n <- nrow(X)
  number_instruments <- ncol(H)
  rho <- coefficients[1L]
  beta <- coefficients[2L:(k + 1L)]

  S <- diag(n) - rho * W
  epsilon <- drop(S %*% y - X %*% beta)
  sigma2 <- mean(epsilon^2)

  P1 <- W
  W2 <- W %*% W
  P2 <- W2 - (tr(W2) / n) * diag(n)
  P1_symmetric <- P1 + t(P1)
  P2_symmetric <- P2 + t(P2)

  Delta <- matrix(0, nrow = 2L, ncol = 2L)
  Delta[1L, 1L] <- tr(P1 %*% P1_symmetric)
  Delta[1L, 2L] <- tr(P1 %*% P2_symmetric)
  Delta[2L, 1L] <- tr(P2 %*% P1_symmetric)
  Delta[2L, 2L] <- tr(P2 %*% P2_symmetric)

  V <- matrix(0, nrow = number_instruments + 2L, ncol = number_instruments + 2L)
  V[1L:2L, 1L:2L] <- Delta
  V[3L:(number_instruments + 2L), 3L:(number_instruments + 2L)] <-
    crossprod(H) / sigma2
  V <- sigma2^2 * V

  omega   <- cbind(diag(P1), diag(P2))
  mu3_hat <- mean(epsilon^3)
  mu4_hat <- mean(epsilon^4)

  nonnormal <- matrix(0, nrow = number_instruments + 2L, 
                         ncol = number_instruments + 2L)
  nonnormal[1L:2L, 1L:2L] <-
    (mu4_hat - 3 * sigma2^2) * crossprod(omega)
  nonnormal[1L:2L, 3L:(number_instruments + 2L)] <-
    mu3_hat * crossprod(omega, H)
  nonnormal[3L:(number_instruments + 2L), 1L:2L] <-
    mu3_hat * crossprod(H, omega)

  Omega_hat <- (nonnormal + V) / n

  # Remove negligible numerical asymmetry before inversion.
  (Omega_hat + t(Omega_hat)) / 2
}

make.vmom() estimates the asymptotic covariance of \(\sqrt n\,\overline{\mathbf g}_n\). Its blocks distinguish three sources of variation: quadratic moments, linear moments, and their covariance when the innovation distribution is asymmetric. This is not an arbitrary heteroskedasticity-robust covariance matrix; it corresponds to the maintained independent-innovation assumptions in the theoretical chapter.

The final averaging with the transpose removes only floating-point asymmetry. It does not change the analytical formula.

7.5.2 Estimation Function

# Main function
slm.gmm <- function(formula, data, W, instruments = 2L,
                    estimator = c("gmm", "ogmm"), gradient = TRUE) {
  call <- match.call(expand.dots = TRUE)
  estimator <- match.arg(estimator)
  components <- spatial_model_components(formula = formula, data = data, A = W,
                                         matrix_name = "W")
  y <- components$y
  X <- components$X
  W <- components$A
  n <- components$n
  k <- components$k
  
  H <- cbind(X, make.H(W = W, X = X, l = instruments))
  H <- independent_columns(H)

  # The compact parameter space keeps I - rho W invertible for inference.
  rho_interval <- spatial_parameter_interval(W)
  Wy <- drop(W %*% y)
  rho_start <- safe_correlation(Wy, y)
  rho_start <- clamp_to_interval(rho_start, rho_interval)

  beta_start <- lm.fit(x = X, y = y)$coefficients
  starting_values <- c("rho" = rho_start, beta_start)
  names(starting_values)[-1L] <- colnames(X)

  constraints <- spatial_parameter_constraints(
    number_parameters = k + 1L,
    parameter_position = 1L,
    interval = rho_interval
  )

  # First-step GMM uses the identity weighting matrix.
  Psi <- diag(ncol(H) + 2L)
  fit_first <- maxLik::maxLik(logLik = Qmin.slm, 
                              start = starting_values,
                              method = "BFGS",
                              constraints = constraints,
                              y = y,
                              X = X,
                              H = H,
                              W = W,
                              Psi = Psi,
                              gradient = gradient,
                              print.level = 0,
                              finalHessian = FALSE)
  fit <- fit_first
  optimizer_codes <- c(first = maxLik::returnCode(fit_first))

  # OGMM freezes the covariance estimate from the first step and reoptimizes.
  if (estimator == "ogmm") {
    Omega_hat <- make.vmom(coef(fit_first), y = y, X = X, H = H, W = W)
    Psi <- invert_moment_covariance(
      Omega_hat,
      matrix_name = "The estimated SLM moment covariance matrix"
    )

    fit_second <- maxLik::maxLik(
      logLik = Qmin.slm,
      start = coef(fit_first),
      method = "BFGS",
      constraints = constraints,
      y = y,
      X = X,
      H = H,
      W = W,
      Psi = Psi,
      gradient = gradient,
      print.level = 0,
      finalHessian = FALSE
    )
    fit <- fit_second
    optimizer_codes <- c(
      optimizer_codes,
      second = maxLik::returnCode(fit_second)
    )
  }

  if (any(optimizer_codes != 0L)) {
    warning("At least one SLM GMM optimization stage did not report convergence.")
  }

  final_derivative <- moments.lee2007(
    coef(fit),
    y = y,
    X = X,
    H = H,
    W = W
  )$D
  if (qr(final_derivative)$rank < k + 1L) {
    warning("The sample GMM derivative matrix is rank deficient at the estimate.")
  }

  structure(
    list(
      coefficients = coef(fit),
      call = call,
      X = X,
      H = H,
      y = y,
      W = W,
      Psi = Psi,
      rho.interval = rho_interval,
      estimator = estimator,
      optimizer.codes = optimizer_codes,
      optimizer = fit
    ),
    class = "gmm.slm"
  )
}

The wrapper makes the two-step logic explicit. estimator = "gmm" performs only the identity-weighted step. estimator = "ogmm" estimates the moment covariance at the first-step estimate, inverts it, and holds that matrix fixed during the second optimization.

The interval constraint is not a likelihood restriction. It implements the compact parameter space used by the GMM theory and guarantees that the inverse needed by the population derivative matrix exists. The regression coefficients remain unrestricted.

The population counterpart of the derivative matrix is used as an alternative to the sample gradient. It corresponds to Equation 6.45.

make.D <- function(rho, beta, y, X, H, W) {
  n <- nrow(X)
  k <- ncol(X)
  number_instruments <- ncol(H)

  S <- diag(n) - rho * W
  epsilon <- drop(S %*% y - X %*% beta)
  sigma2 <- mean(epsilon^2)

  P1 <- W
  W2 <- W %*% W
  P2 <- W2 - (tr(W2) / n) * diag(n)
  P1_symmetric <- P1 + t(P1)
  P2_symmetric <- P2 + t(P2)
  G <- W %*% solve(S)

  D <- matrix(
    0,
    nrow = number_instruments + 2L,
    ncol = k + 1L,
    dimnames = list(
      c("q1", "q2", colnames(H)),
      c("rho", colnames(X))
    )
  )

  D[1L, 1L] <- sigma2 * tr(P1_symmetric %*% G)
  D[2L, 1L] <- sigma2 * tr(P2_symmetric %*% G)
  D[3L:(number_instruments + 2L), 1L] <-
    drop(crossprod(H, G %*% X %*% beta))
  D[3L:(number_instruments + 2L), 2L:(k + 1L)] <- crossprod(H, X)

  D
}

vcov.gmm.slm <- function(object, D = c("population", "gradient"), ...) {
  D <- match.arg(D)

  X <- object$X
  H <- object$H
  y <- object$y
  W <- object$W
  n <- nrow(X)
  k <- ncol(X)
  coefficients <- object$coefficients
  rho <- coefficients[1L]
  beta <- coefficients[2L:(k + 1L)]

  derivative <- if (D == "population") {
    make.D(rho, beta, y = y, X = X, H = H, W = W) / n
  } else {
    moments.lee2007(coefficients, y = y, X = X, H = H, W = W)$D
  }

  if (object$estimator == "gmm") {
    Omega <- make.vmom(coefficients, y = y, X = X, H = H, W = W)
    inverse_DD <- solve(crossprod(derivative))
    return(
      inverse_DD %*% t(derivative) %*% Omega %*%
        derivative %*% inverse_DD / n
    )
  }

  solve(t(derivative) %*% object$Psi %*% derivative) / n
}

summary.gmm.slm <- function(
  object,
  D = c("population", "gradient"),
  digits = max(3L, getOption("digits") - 3L),
  ...
) {
  D <- match.arg(D)
  estimates <- drop(object$coefficients)
  standard_errors <- sqrt(diag(vcov(object, D = D)))
  statistics <- estimates / standard_errors
  p_values <- 2 * pnorm(-abs(statistics))

  coefficient_table <- cbind(
    "Estimate" = estimates,
    "Std. Error" = standard_errors,
    "z value" = statistics,
    "Pr(>|z|)" = p_values
  )

  structure(
    list(
      coefficients = coefficient_table,
      estimator = object$estimator,
      derivative = D,
      optimizer.codes = object$optimizer.codes,
      digits = digits,
      call = object$call
    ),
    class = "summary.gmm.slm"
  )
}

print.summary.gmm.slm <- function(x, digits = x$digits, ...) {
  cat("\nCall:\n")
  cat(paste(deparse(x$call), collapse = "\n"), "\n", sep = "")
  cat("\nEstimator:", toupper(x$estimator), "\n")
  cat("Optimizer code(s):", paste(x$optimizer.codes, collapse = ", "), "\n")
  cat("Derivative matrix:", x$derivative, "\n")
  cat("\nCoefficients:\n")
  printCoefmat(
    x$coefficients,
    digits = digits,
    P.values = TRUE,
    has.Pvalue = TRUE
  )
  invisible(x)
}

7.5.3 Verifying the SLM Criterion Gradient

Before estimating the model, we compare the analytical gradient attached to Qmin.slm() with a numerical derivative. The check is performed at an interior parameter vector rather than at the optimum, where every gradient should be close to zero and sign errors are harder to detect.

H_slm_check <- independent_columns(
  cbind(X_slm, make.H(W = W_slm, X = X_slm, l = 2L))
)
theta_slm_check <- c("rho" = 0.35, beta_slm + c(0.10, -0.05, 0.08, -0.04))
Psi_slm_check <- diag(ncol(H_slm_check) + 2L)

criterion_slm_check <- Qmin.slm(
  theta_slm_check,
  y = y_slm,
  X = X_slm,
  H = H_slm_check,
  W = W_slm,
  Psi = Psi_slm_check,
  gradient = TRUE
)

analytical_gradient_slm <- attr(criterion_slm_check, "gradient")
numerical_gradient_slm <- numDeriv::grad(
  func = function(theta) {
    Qmin.slm(
      theta,
      y = y_slm,
      X = X_slm,
      H = H_slm_check,
      W = W_slm,
      Psi = Psi_slm_check,
      gradient = FALSE
    )
  },
  x = theta_slm_check
)

cbind(
  analytical = analytical_gradient_slm,
  numerical = numerical_gradient_slm,
  difference = analytical_gradient_slm - numerical_gradient_slm
)

            analytical  numerical    difference
             1.7421351  1.7421351 -3.354783e-11
(Intercept) -0.9618739 -0.9618739  2.166676e-10
x1          -0.2084953 -0.2084953 -2.239070e-12
x2          -0.1679047 -0.1679047 -5.997092e-12
x3           0.4429805  0.4429805 -7.280843e-12

max(abs(analytical_gradient_slm - numerical_gradient_slm))

[1] 2.166676e-10

A small maximum discrepancy provides a direct diagnostic of the programmed Jacobian, the sign of the criterion, and the ordering of the parameters.

7.5.4 GMM and OGMM Estimates

# Cache revision: recompute after the estimator and DGP audit.

fit_gmm_slm <- slm.gmm(
  y ~ x1 + x2 + x3,
  data = data_slm,
  W = W_slm,
  instruments = 2L,
  estimator = "gmm"
)

summary(fit_gmm_slm, D = "population")


Call:
slm.gmm(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L, estimator = "gmm")

Estimator: GMM 
Optimizer code(s): 0 
Derivative matrix: population 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
rho          0.53223    0.04352  12.230   <2e-16 ***
(Intercept) -0.10692    0.06530  -1.637    0.102    
x1          -1.02829    0.06365 -16.155   <2e-16 ***
x2          -0.04131    0.06355  -0.650    0.516    
x3           1.05348    0.06586  15.997   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_gmm_slm, D = "gradient")


Call:
slm.gmm(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L, estimator = "gmm")

Estimator: GMM 
Optimizer code(s): 0 
Derivative matrix: gradient 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
rho          0.53223    0.04983  10.681   <2e-16 ***
(Intercept) -0.10692    0.06623  -1.614    0.106    
x1          -1.02829    0.06400 -16.067   <2e-16 ***
x2          -0.04131    0.06359  -0.650    0.516    
x3           1.05348    0.06914  15.237   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit_ogmm_slm <- slm.gmm(
  y ~ x1 + x2 + x3,
  data = data_slm,
  W = W_slm,
  instruments = 2L,
  estimator = "ogmm"
)

summary(fit_ogmm_slm, D = "population")


Call:
slm.gmm(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L, estimator = "ogmm")

Estimator: OGMM 
Optimizer code(s): 0, 0 
Derivative matrix: population 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
rho          0.56341    0.03739  15.068   <2e-16 ***
(Intercept) -0.09800    0.06469  -1.515    0.130    
x1          -1.03166    0.06328 -16.302   <2e-16 ***
x2          -0.03459    0.06335  -0.546    0.585    
x3           1.03197    0.06530  15.803   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_ogmm_slm, D = "gradient")


Call:
slm.gmm(formula = y ~ x1 + x2 + x3, data = data_slm, W = W_slm, 
    instruments = 2L, estimator = "ogmm")

Estimator: OGMM 
Optimizer code(s): 0, 0 
Derivative matrix: gradient 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
rho          0.56341    0.04321  13.040   <2e-16 ***
(Intercept) -0.09800    0.06533  -1.500    0.134    
x1          -1.03166    0.06324 -16.313   <2e-16 ***
x2          -0.03459    0.06324  -0.547    0.584    
x3           1.03197    0.06711  15.377   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The following table compares the estimators programmed in this chapter. The maximum-likelihood benchmark is obtained with spatialreg::lagsarlm() rather than relying on an external function defined outside this chapter.

fit_mle_slm <- lagsarlm(
  y ~ x1 + x2 + x3,
  data = data_slm,
  listw = listw_slm,
  method = "eigen"
)

parameter_names_slm <- c(colnames(X_slm), "rho")

# Extract each estimator in the common order
# (Intercept), x1, x2, x3, rho. Using the stored regression
# coefficients avoids relying on the package-specific ordering of coef().
mle_parameters_slm <- c(
  fit_mle_slm$coefficients[colnames(X_slm)],
  fit_mle_slm$rho
)

s2sls_parameters_slm <- setNames(
  as.numeric(fit_s2sls$coefficients),
  parameter_names_slm
)

bs2sls_parameters_slm <- setNames(
  as.numeric(fit_bs2sls$coefficients),
  parameter_names_slm
)

gmm_parameters_slm <- c(
  setNames(
    as.numeric(fit_gmm_slm$coefficients[-1L]),
    colnames(X_slm)
  ),
  "rho" = as.numeric(fit_gmm_slm$coefficients[1L])
)

ogmm_parameters_slm <- c(
  setNames(
    as.numeric(fit_ogmm_slm$coefficients[-1L]),
    colnames(X_slm)
  ),
  "rho" = as.numeric(fit_ogmm_slm$coefficients[1L])
)

true_parameters_slm <- c(beta_slm, "rho" = rho_slm)

coefficient_comparison_slm <- cbind(
  "True" = true_parameters_slm[parameter_names_slm],
  "MLE" = mle_parameters_slm[parameter_names_slm],
  "S2SLS" = s2sls_parameters_slm[parameter_names_slm],
  "BS2SLS" = bs2sls_parameters_slm[parameter_names_slm],
  "GMM" = gmm_parameters_slm[parameter_names_slm],
  "OGMM" = ogmm_parameters_slm[parameter_names_slm]
)

round(coefficient_comparison_slm, 4L)

            True     MLE   S2SLS  BS2SLS     GMM    OGMM
(Intercept)  0.0 -0.1086 -0.0698 -0.0816 -0.1069 -0.0980
x1          -1.0 -1.0286 -1.0057 -1.0127 -1.0283 -1.0317
x2           0.0 -0.0383 -0.0435 -0.0420 -0.0413 -0.0346
x3           1.0  1.0507  1.0149  1.0257  1.0535  1.0320
rho          0.6  0.5391  0.6473  0.6146  0.5322  0.5634

The comparison should not be read as a contest decided by one simulated sample. Its purpose is diagnostic. The true column reveals sampling error, while the agreement among independently programmed estimators helps detect parameter-order or extraction mistakes.

7.6 Feasible GLS for the SEM

The spatial error model is \[ \mathbf y = \mathbf X\boldsymbol\beta + \mathbf u, \qquad \mathbf u = \lambda\mathbf M\mathbf u + \boldsymbol\varepsilon. \] Here \(\mathbf u\) is the regression disturbance, whereas \(\boldsymbol\varepsilon=(\mathbf I-\lambda\mathbf M)\mathbf u\) is the innovation. This distinction is essential: OLS initially estimates \(\mathbf u\), the moment equations estimate \(\lambda\), and GLS is applied only after the innovation transformation has been estimated.

The feasible procedure has three steps:

estimate the nonspatial regression and retain the OLS residuals;
estimate \(\lambda\) and \(\sigma^2\) from the nonlinear moment equations in Section 6.4;
transform both \(\mathbf y\) and \(\mathbf X\) by \(\widehat{\mathbf R}=\mathbf I-\widehat\lambda\mathbf M\) and run OLS on the transformed model.

7.6.1 Moment Equations for the Spatial Error Parameter

Let \(\widehat{\mathbf u}\) denote the OLS residual vector, \(\widehat{\mathbf u}_L=\mathbf M\widehat{\mathbf u}\), and \(\widehat{\mathbf u}_{LL}=\mathbf M^2\widehat{\mathbf u}\). The sample system can be written as \[ \mathbf g \approx \mathbf G \begin{pmatrix} \lambda\\ \lambda^2\\ \sigma^2 \end{pmatrix}. \] mom.sem() constructs the three-vector \(\mathbf g\) and the \(3\times3\) matrix \(\mathbf G\). Qn.sem() imposes the nonlinear restriction that the second parameter is the square of the first instead of estimating three unrelated coefficients.

mom.sem <- function(u, M) {
  u <- as.numeric(u)
  M <- validate_spatial_matrix(M, n = length(u), matrix_name = "M")
  n <- length(u)

  u_lag <- drop(M %*% u)
  u_lag2 <- drop(M %*% u_lag)
  trace_MM <- tr(crossprod(M))

  G <- matrix(0, nrow = 3L, ncol = 3L)
  G[1L, 1L] <- 2 * crossprod(u, u_lag)
  G[2L, 1L] <- 2 * crossprod(u_lag2, u_lag)
  G[3L, 1L] <- crossprod(u, u_lag2) + crossprod(u_lag)
  G[1L, 2L] <- -crossprod(u_lag)
  G[2L, 2L] <- -crossprod(u_lag2)
  G[3L, 2L] <- -crossprod(u_lag, u_lag2)
  G[1L, 3L] <- n
  G[2L, 3L] <- trace_MM
  G <- G / n
  colnames(G) <- c("lambda", "lambda^2", "sigma2")

  g <- c(
    "u'u" = crossprod(u),
    "uL'uL" = crossprod(u_lag),
    "u'uL" = crossprod(u, u_lag)
  ) / n

  list(G = G, g = g)
}

Qn.sem <- function(parameters, moments, verbose = FALSE) {
  lambda <- parameters[1L]
  sigma2 <- parameters[2L]
  discrepancy <- moments$g - moments$G %*% c(lambda, lambda^2, sigma2)
  criterion <- drop(crossprod(discrepancy))

  if (verbose) {
    cat(
      "criterion:", criterion,
      "lambda:", lambda,
      "sigma2:", sigma2,
      "\n"
    )
  }

  criterion
}

7.6.2 FGLS Estimation Function

sem.sfgls <- function(formula, data, M, verbose = FALSE) {
  call <- match.call(expand.dots = TRUE)
  components <- spatial_model_components(
    formula = formula,
    data = data,
    A = M,
    matrix_name = "M"
  )

  y <- components$y
  X <- components$X
  M <- components$A
  n <- components$n

  # Step 1: OLS residuals estimate the spatially correlated disturbance u.
  initial_ols <- lm.fit(x = X, y = y)
  u_initial <- initial_ols$residuals
  moments <- mom.sem(u = u_initial, M = M)

  # Step 2: nonlinear minimum-distance estimates of lambda and sigma^2.
  lambda_interval <- spatial_parameter_interval(M)
  lambda_start <- safe_correlation(drop(M %*% u_initial), u_initial)
  lambda_start <- clamp_to_interval(lambda_start, lambda_interval)
  sigma2_start <- max(mean(u_initial^2), 10 * .Machine$double.eps)

  optimization <- nlminb(
    start = c("lambda" = lambda_start, "sigma2" = sigma2_start),
    objective = Qn.sem,
    moments = moments,
    verbose = verbose,
    lower = c(lambda_interval["lower"], .Machine$double.eps),
    upper = c(lambda_interval["upper"], Inf)
  )

  if (optimization$convergence != 0L) {
    warning("The SEM moment optimization did not report convergence.", call. = FALSE)
  }

  # Step 3: feasible GLS on the spatially transformed model.
  lambda_hat <- unname(optimization$par["lambda"])
  R_hat <- diag(n) - lambda_hat * M
  y_transformed <- drop(R_hat %*% y)
  X_transformed <- R_hat %*% X
  beta_hat <- drop(solve(
    crossprod(X_transformed),
    crossprod(X_transformed, y_transformed)
  ))
  names(beta_hat) <- colnames(X)
  epsilon_hat <- drop(y_transformed - X_transformed %*% beta_hat)

  structure(
    list(
      coefficients = c(beta_hat, "lambda" = lambda_hat),
      call = call,
      X = X,
      y = y,
      M = M,
      X.transformed = X_transformed,
      y.transformed = y_transformed,
      residuals = epsilon_hat,
      moments = moments,
      lambda.interval = lambda_interval,
      optimizer = optimization
    ),
    class = "myfgls.sem"
  )
}

vcov.myfgls.sem <- function(object, ...) {
  n <- nrow(object$X)
  sigma2_hat <- sum(object$residuals^2) / n
  sigma2_hat * solve(crossprod(object$X.transformed))
}

summary.myfgls.sem <- function(
  object,
  digits = max(3L, getOption("digits") - 3L),
  ...
) {
  k <- ncol(object$X)
  estimates <- object$coefficients[seq_len(k)]
  standard_errors <- sqrt(diag(vcov(object)))
  statistics <- estimates / standard_errors
  p_values <- 2 * pnorm(-abs(statistics))

  coefficient_table <- cbind(
    "Estimate" = estimates,
    "Std. Error" = standard_errors,
    "z value" = statistics,
    "Pr(>|z|)" = p_values
  )

  structure(
    list(
      coefficients = coefficient_table,
      lambda = unname(object$coefficients["lambda"]),
      convergence = object$optimizer$convergence,
      digits = digits,
      call = object$call
    ),
    class = "summary.myfgls.sem"
  )
}

print.summary.myfgls.sem <- function(x, digits = x$digits, ...) {
  cat("\nCall:\n")
  cat(paste(deparse(x$call), collapse = "\n"), "\n", sep = "")
  cat("\nOptimizer convergence code:", x$convergence, "\n")
  cat("Estimated spatial error parameter:", format(x$lambda, digits = digits), "\n")
  cat("\nRegression coefficients:\n")
  printCoefmat(
    x$coefficients,
    digits = digits,
    P.values = TRUE,
    has.Pvalue = TRUE
  )
  invisible(x)
}

The transformed regression has innovation \(\widehat{\boldsymbol\varepsilon}=\widehat{\mathbf R}(\mathbf y-\mathbf X\widehat{\boldsymbol\beta})\). The covariance method therefore uses the cross-product of \(\widehat{\mathbf R}\mathbf X\). The code now assigns the column names of \(\mathbf X\) explicitly to beta_hat; without that conversion, solve() returns a one-column matrix and c() may discard the regression-coefficient names.

The simple covariance method above reports inference for the regression coefficients. It does not construct a standard error for the first-stage moment estimator of the spatial error parameter.

7.6.3 Artificial SEM Data and Estimation

set.seed(1)
side_sem <- 23L
n_sem <- side_sem^2
lambda_sem <- 0.6

neighbors_sem <- cell2nb(side_sem, side_sem, type = "rook")
M_sem <- nb2mat(neighbors_sem, style = "W")

x1_sem <- rnorm(n_sem)
x2_sem <- rnorm(n_sem)
x3_sem <- rnorm(n_sem)
X_sem <- cbind(
  "(Intercept)" = 1,
  "x1" = x1_sem,
  "x2" = x2_sem,
  "x3" = x3_sem
)
beta_sem <- c("(Intercept)" = 0, "x1" = -1, "x2" = 0, "x3" = 1)
sigma2_sem <- 2
epsilon_sem <- rnorm(n_sem, mean = 0, sd = sqrt(sigma2_sem))

R_sem <- diag(n_sem) - lambda_sem * M_sem
u_sem <- drop(solve(R_sem, epsilon_sem))
y_sem <- drop(X_sem %*% beta_sem + u_sem)

data_sem <- data.frame(
  y = y_sem,
  x1 = x1_sem,
  x2 = x2_sem,
  x3 = x3_sem
)

# Verify both equations of the SEM data-generating process.
c(
  regression_equation = max(abs(y_sem - X_sem %*% beta_sem - u_sem)),
  error_equation = max(abs(R_sem %*% u_sem - epsilon_sem))
)

regression_equation      error_equation 
       8.881784e-16        3.552714e-15

# Cache revision: recompute after the estimator and DGP audit.

fit_fgls_sem <- sem.sfgls(
  y ~ x1 + x2 + x3,
  data = data_sem,
  M = M_sem,
  verbose = FALSE
)
summary(fit_fgls_sem)


Call:
sem.sfgls(formula = y ~ x1 + x2 + x3, data = data_sem, M = M_sem, 
    verbose = FALSE)

Optimizer convergence code: 0 
Estimated spatial error parameter: 0.556 

Regression coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.192429   0.149854  -1.284    0.199    
x1          -1.079338   0.063535 -16.988   <2e-16 ***
x2           0.006391   0.060054   0.106    0.915    
x3           0.975246   0.064009  15.236   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cbind(
  "True" = c(beta_sem, "lambda" = lambda_sem),
  "FGLS" = fit_fgls_sem$coefficients
)

            True         FGLS
(Intercept)  0.0 -0.192428998
x1          -1.0 -1.079338334
x2           0.0  0.006390999
x3           1.0  0.975245554
lambda       0.6  0.556015812

The comparison with the true parameters separates implementation errors from ordinary sampling variation before we turn to the package benchmark.

We compare the estimates with spatialreg::GMerrorsar().

# Cache revision: recompute after the estimator and DGP audit.

fit_gmerrorsar <- GMerrorsar(
  y ~ x1 + x2 + x3,
  data = data_sem,
  listw = mat2listw(M_sem, style = "W"),
  verbose = FALSE,
  legacy = TRUE
)
summary(fit_gmerrorsar)


Call:
GMerrorsar(formula = y ~ x1 + x2 + x3, data = data_sem, listw = mat2listw(M_sem, 
    style = "W"), verbose = FALSE, legacy = TRUE)

Residuals:
      Min        1Q    Median        3Q       Max 
-3.748647 -0.946972 -0.050107  1.076211  4.413975 

Type: GM SAR estimator
Coefficients: (GM standard errors) 
             Estimate Std. Error  z value Pr(>|z|)
(Intercept) -0.192429   0.149854  -1.2841   0.1991
x1          -1.079338   0.063535 -16.9880   <2e-16
x2           0.006391   0.060054   0.1064   0.9152
x3           0.975246   0.064009  15.2362   <2e-16

Lambda: 0.55602 (standard error): 0.10098 (z-value): 5.5061
Residual variance (sigma squared): 2.3414, (sigma: 1.5302)
GM argmin sigma squared: 2.3535
Number of observations: 529 
Number of parameters estimated: 6

7.7 GMM and OGMM for the SEM

The preceding feasible-GLS estimator estimates \(\lambda\) from three nonlinear moments and then estimates \(\boldsymbol\beta\). The SEM GMM estimator instead estimates \[ \boldsymbol\theta = \left( \boldsymbol\beta^\top, \lambda \right)^\top \] jointly from linear and quadratic moments.

For a candidate parameter vector, define \[ \mathbf u(\boldsymbol\beta) = \mathbf y-\mathbf X\boldsymbol\beta, \] and \[ \boldsymbol\varepsilon(\boldsymbol\theta) = (\mathbf I-\lambda\mathbf M)\mathbf u(\boldsymbol\beta). \] The parameter ordering now differs from the SLM GMM section: the regression coefficients come first and lambda comes last. Making that order explicit is important for both the gradient and the covariance matrix.

7.7.1 Moment Function and Moment Covariance

The normalized moment vector places the \(k\) linear moments first and the two quadratic moments last: \[ \overline{\mathbf g}_n(\boldsymbol\theta) = \frac{1}{n} \begin{pmatrix} \mathbf X^\top\boldsymbol\varepsilon\\ \boldsymbol\varepsilon^\top\mathbf P_1\boldsymbol\varepsilon\\ \boldsymbol\varepsilon^\top\mathbf P_2\boldsymbol\varepsilon \end{pmatrix}. \] Because \[ \frac{\partial\boldsymbol\varepsilon}{\partial\boldsymbol\theta^\top} = - \left[ (\mathbf I-\lambda\mathbf M)\mathbf X, \mathbf M\mathbf u \right], \] the same code can assemble every row of the Jacobian by matrix multiplication.

moments.lee2010 <- function(theta, y, X, M) {
  X <- as.matrix(X)
  M <- as.matrix(M)
  y <- as.numeric(y)

  k <- ncol(X)
  n <- nrow(X)
  if (length(theta) != k + 1L) {
    stop("theta must contain the k regression coefficients followed by lambda.")
  }

  beta <- theta[seq_len(k)]
  lambda <- theta[k + 1L]

  R <- diag(n) - lambda * M
  u <- drop(y - X %*% beta)
  epsilon <- drop(R %*% u)

  P1 <- M
  M2 <- M %*% M
  P2 <- M2 - (tr(M2) / n) * diag(n)

  moments_linear <- crossprod(X, epsilon)
  moment_quadratic_1 <- drop(crossprod(epsilon, P1 %*% epsilon))
  moment_quadratic_2 <- drop(crossprod(epsilon, P2 %*% epsilon))
  moments <- rbind(moments_linear, moment_quadratic_1, moment_quadratic_2)
  rownames(moments) <- c(colnames(X), "q1", "q2")

  P1_symmetric <- P1 + t(P1)
  P2_symmetric <- P2 + t(P2)
  derivative_epsilon <- cbind(R %*% X, "lambda" = M %*% u)

  derivative <- -rbind(
    t(X),
    crossprod(epsilon, P1_symmetric),
    crossprod(epsilon, P2_symmetric)
  ) %*% derivative_epsilon

  list(
    g = moments / n,
    D = derivative / n,
    epsilon = epsilon,
    u = u,
    P1 = P1,
    P2 = P2
  )
}

Q.sem <- function(theta, y, X, M, Psi, gradient = TRUE) {
  moment_fit <- moments.lee2010(theta, y = y, X = X, M = M)
  criterion <- -drop(crossprod(moment_fit$g, Psi %*% moment_fit$g))

  if (gradient) {
    criterion_gradient <- -2 * crossprod(moment_fit$D, Psi %*% moment_fit$g)
    attr(criterion, "gradient") <- drop(criterion_gradient)
  }

  criterion
}

make.vmom.sem <- function(coefficients, y, X, M) {
  X <- as.matrix(X)
  M <- as.matrix(M)
  y <- as.numeric(y)

  k <- ncol(X)
  n <- nrow(X)
  beta <- coefficients[seq_len(k)]
  lambda <- coefficients[k + 1L]

  R <- diag(n) - lambda * M
  u <- drop(y - X %*% beta)
  epsilon <- drop(R %*% u)
  sigma2 <- mean(epsilon^2)

  P1 <- M
  M2 <- M %*% M
  P2 <- M2 - (tr(M2) / n) * diag(n)
  P1_symmetric <- P1 + t(P1)
  P2_symmetric <- P2 + t(P2)

  Delta <- matrix(0, nrow = 2L, ncol = 2L)
  Delta[1L, 1L] <- tr(P1_symmetric %*% P1)
  Delta[1L, 2L] <- tr(P1_symmetric %*% P2)
  Delta[2L, 1L] <- tr(P2_symmetric %*% P1)
  Delta[2L, 2L] <- tr(P2_symmetric %*% P2)

  V <- matrix(0, nrow = k + 2L, ncol = k + 2L)
  V[seq_len(k), seq_len(k)] <- crossprod(X) / sigma2
  V[(k + 1L):(k + 2L), (k + 1L):(k + 2L)] <- Delta
  V <- sigma2^2 * V

  omega <- cbind(diag(P1), diag(P2))
  mu3_hat <- mean(epsilon^3)
  mu4_hat <- mean(epsilon^4)

  nonnormal <- matrix(0, nrow = k + 2L, ncol = k + 2L)
  nonnormal[(k + 1L):(k + 2L), (k + 1L):(k + 2L)] <-
    (mu4_hat - 3 * sigma2^2) * crossprod(omega)
  nonnormal[(k + 1L):(k + 2L), seq_len(k)] <-
    mu3_hat * crossprod(omega, X)
  nonnormal[seq_len(k), (k + 1L):(k + 2L)] <-
    mu3_hat * crossprod(X, omega)

  Omega_hat <- (nonnormal + V) / n
  (Omega_hat + t(Omega_hat)) / 2
}

The block ordering in make.vmom.sem() matches the moment ordering exactly: linear moments first, quadratic moments last. As in the SLM, the averaging with the transpose removes only numerical asymmetry.

7.7.2 SEM GMM Estimation Function

sem.gmm <- function(
  formula,
  data,
  M,
  estimator = c("gmm", "ogmm"),
  gradient = TRUE
) {
  call <- match.call(expand.dots = TRUE)
  estimator <- match.arg(estimator)

  components <- spatial_model_components(
    formula = formula,
    data = data,
    A = M,
    matrix_name = "M"
  )
  y <- components$y
  X <- components$X
  M <- components$A
  k <- components$k

  ols <- lm.fit(x = X, y = y)
  ols_residuals <- ols$residuals

  lambda_interval <- spatial_parameter_interval(M)
  lambda_start <- safe_correlation(drop(M %*% ols_residuals), ols_residuals)
  lambda_start <- clamp_to_interval(lambda_start, lambda_interval)

  starting_values <- c(ols$coefficients, "lambda" = lambda_start)
  names(starting_values)[seq_len(k)] <- colnames(X)

  constraints <- spatial_parameter_constraints(
    number_parameters = k + 1L,
    parameter_position = k + 1L,
    interval = lambda_interval
  )

  Psi <- diag(k + 2L)
  fit_first <- maxLik::maxLik(
    logLik = Q.sem,
    start = starting_values,
    method = "BFGS",
    constraints = constraints,
    y = y,
    X = X,
    M = M,
    Psi = Psi,
    gradient = gradient,
    print.level = 0,
    finalHessian = FALSE
  )
  fit <- fit_first
  optimizer_codes <- c(first = maxLik::returnCode(fit_first))

  if (estimator == "ogmm") {
    Omega_hat <- make.vmom.sem(coef(fit_first), y = y, X = X, M = M)
    Psi <- invert_moment_covariance(
      Omega_hat,
      matrix_name = "The estimated SEM moment covariance matrix"
    )

    fit_second <- maxLik::maxLik(
      logLik = Q.sem,
      start = coef(fit_first),
      method = "BFGS",
      constraints = constraints,
      y = y,
      X = X,
      M = M,
      Psi = Psi,
      gradient = gradient,
      print.level = 0,
      finalHessian = FALSE
    )
    fit <- fit_second
    optimizer_codes <- c(
      optimizer_codes,
      second = maxLik::returnCode(fit_second)
    )
  }

  if (any(optimizer_codes != 0L)) {
    warning("At least one SEM GMM optimization stage did not report convergence.")
  }

  final_derivative <- moments.lee2010(
    coef(fit),
    y = y,
    X = X,
    M = M
  )$D
  if (qr(final_derivative)$rank < k + 1L) {
    warning("The sample SEM GMM derivative matrix is rank deficient at the estimate.")
  }

  structure(
    list(
      coefficients = coef(fit),
      call = call,
      X = X,
      y = y,
      M = M,
      Psi = Psi,
      lambda.interval = lambda_interval,
      estimator = estimator,
      optimizer.codes = optimizer_codes,
      optimizer = fit
    ),
    class = "gmm.sem"
  )
}

The estimation wrapper repeats the same two-step logic used for the SLM. The only differences are the parameter ordering and the absence of an external instrument matrix: the linear moments use \(\mathbf X\) directly.

The population derivative below is the negative expectation of the sample Jacobian. Its zero blocks encode orthogonality results, not omitted calculations. For example, the expected derivative of the linear moments with respect to \(\lambda\) is zero at the true parameter under exogeneity.

make.D.sem <- function(coefficients, y, X, M) {
  k <- ncol(X)
  n <- nrow(X)
  beta <- coefficients[seq_len(k)]
  lambda <- coefficients[k + 1L]

  R <- diag(n) - lambda * M
  u <- drop(y - X %*% beta)
  epsilon <- drop(R %*% u)
  sigma2 <- mean(epsilon^2)

  P1 <- M
  M2 <- M %*% M
  P2 <- M2 - (tr(M2) / n) * diag(n)
  P1_symmetric <- P1 + t(P1)
  P2_symmetric <- P2 + t(P2)
  Q <- M %*% solve(R)

  D <- matrix(
    0,
    nrow = k + 2L,
    ncol = k + 1L,
    dimnames = list(
      c(colnames(X), "q1", "q2"),
      c(colnames(X), "lambda")
    )
  )

  D[seq_len(k), seq_len(k)] <- t(X) %*% R %*% X
  D[k + 1L, k + 1L] <- sigma2 * tr(P1_symmetric %*% Q)
  D[k + 2L, k + 1L] <- sigma2 * tr(P2_symmetric %*% Q)
  D
}

vcov.gmm.sem <- function(object, D = c("population", "gradient"), ...) {
  D <- match.arg(D)

  X <- object$X
  y <- object$y
  M <- object$M
  n <- nrow(X)
  coefficients <- object$coefficients

  derivative <- if (D == "population") {
    make.D.sem(coefficients, y = y, X = X, M = M) / n
  } else {
    moments.lee2010(coefficients, y = y, X = X, M = M)$D
  }

  if (object$estimator == "gmm") {
    Omega <- make.vmom.sem(coefficients, y = y, X = X, M = M)
    inverse_DD <- solve(crossprod(derivative))
    return(
      inverse_DD %*% t(derivative) %*% Omega %*%
        derivative %*% inverse_DD / n
    )
  }

  solve(t(derivative) %*% object$Psi %*% derivative) / n
}

summary.gmm.sem <- function(
  object,
  D = c("population", "gradient"),
  digits = max(3L, getOption("digits") - 3L),
  ...
) {
  D <- match.arg(D)
  estimates <- drop(object$coefficients)
  standard_errors <- sqrt(diag(vcov(object, D = D)))
  statistics <- estimates / standard_errors
  p_values <- 2 * pnorm(-abs(statistics))

  coefficient_table <- cbind(
    "Estimate" = estimates,
    "Std. Error" = standard_errors,
    "z value" = statistics,
    "Pr(>|z|)" = p_values
  )

  structure(
    list(
      coefficients = coefficient_table,
      estimator = object$estimator,
      derivative = D,
      optimizer.codes = object$optimizer.codes,
      digits = digits,
      call = object$call
    ),
    class = "summary.gmm.sem"
  )
}

print.summary.gmm.sem <- function(x, digits = x$digits, ...) {
  cat("\nCall:\n")
  cat(paste(deparse(x$call), collapse = "\n"), "\n", sep = "")
  cat("\nEstimator:", toupper(x$estimator), "\n")
  cat("Optimizer code(s):", paste(x$optimizer.codes, collapse = ", "), "\n")
  cat("Derivative matrix:", x$derivative, "\n")
  cat("\nCoefficients:\n")
  printCoefmat(
    x$coefficients,
    digits = digits,
    P.values = TRUE,
    has.Pvalue = TRUE
  )
  invisible(x)
}

7.7.3 Verifying the SEM Criterion Gradient

theta_sem_check <- c(beta_sem + c(0.10, -0.05, 0.08, -0.04), "lambda" = 0.35)
Psi_sem_check <- diag(ncol(X_sem) + 2L)

criterion_sem_check <- Q.sem(
  theta_sem_check,
  y = y_sem,
  X = X_sem,
  M = M_sem,
  Psi = Psi_sem_check,
  gradient = TRUE
)

analytical_gradient_sem <- attr(criterion_sem_check, "gradient")
numerical_gradient_sem <- numDeriv::grad(
  func = function(theta) {
    Q.sem(
      theta,
      y = y_sem,
      X = X_sem,
      M = M_sem,
      Psi = Psi_sem_check,
      gradient = FALSE
    )
  },
  x = theta_sem_check
)

cbind(
  analytical = analytical_gradient_sem,
  numerical = numerical_gradient_sem,
  difference = analytical_gradient_sem - numerical_gradient_sem
)

             analytical   numerical    difference
(Intercept) -0.46129555 -0.46129555 -1.711681e-11
x1           0.03073563  0.03073563  9.319250e-12
x2          -0.07150477 -0.07150477 -8.978487e-11
x3           0.08815837  0.08815837 -1.576013e-11
             1.53330627  1.53330627 -6.689227e-11

max(abs(analytical_gradient_sem - numerical_gradient_sem))

[1] 8.978487e-11

This check is especially useful here because the position of lambda differs from the SLM parameter vector.

7.7.4 SEM GMM and OGMM Estimates

# Cache revision: recompute after the estimator and DGP audit.

fit_gmm_sem <- sem.gmm(
  y ~ x1 + x2 + x3,
  data = data_sem,
  M = M_sem,
  estimator = "gmm"
)
summary(fit_gmm_sem, D = "population")


Call:
sem.gmm(formula = y ~ x1 + x2 + x3, data = data_sem, M = M_sem, 
    estimator = "gmm")

Estimator: GMM 
Optimizer code(s): 0 
Derivative matrix: population 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.19250    0.15049  -1.279    0.201    
x1          -1.06289    0.06607 -16.087   <2e-16 ***
x2           0.01368    0.06197   0.221    0.825    
x3           0.98552    0.06651  14.818   <2e-16 ***
lambda       0.55800    0.04690  11.898   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_gmm_sem, D = "gradient")


Call:
sem.gmm(formula = y ~ x1 + x2 + x3, data = data_sem, M = M_sem, 
    estimator = "gmm")

Estimator: GMM 
Optimizer code(s): 0 
Derivative matrix: gradient 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.19250    0.15049  -1.279    0.201    
x1          -1.06289    0.06622 -16.051   <2e-16 ***
x2           0.01368    0.06186   0.221    0.825    
x3           0.98552    0.06657  14.804   <2e-16 ***
lambda       0.55800    0.04923  11.334   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit_ogmm_sem <- sem.gmm(
  y ~ x1 + x2 + x3,
  data = data_sem,
  M = M_sem,
  estimator = "ogmm"
)
summary(fit_ogmm_sem, D = "population")


Call:
sem.gmm(formula = y ~ x1 + x2 + x3, data = data_sem, M = M_sem, 
    estimator = "ogmm")

Estimator: OGMM 
Optimizer code(s): 0, 0 
Derivative matrix: population 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.19255    0.15183  -1.268    0.205    
x1          -1.06348    0.06607 -16.096   <2e-16 ***
x2           0.01109    0.06196   0.179    0.858    
x3           0.98598    0.06652  14.823   <2e-16 ***
lambda       0.56190    0.04561  12.319   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_ogmm_sem, D = "gradient")


Call:
sem.gmm(formula = y ~ x1 + x2 + x3, data = data_sem, M = M_sem, 
    estimator = "ogmm")

Estimator: OGMM 
Optimizer code(s): 0, 0 
Derivative matrix: gradient 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.19255    0.15183  -1.268    0.205    
x1          -1.06348    0.06621 -16.061   <2e-16 ***
x2           0.01109    0.06175   0.180    0.857    
x3           0.98598    0.06658  14.810   <2e-16 ***
lambda       0.56190    0.04909  11.445   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

parameter_names_sem <- c(colnames(X_sem), "lambda")
cbind(
  "True" = c(beta_sem, "lambda" = lambda_sem)[parameter_names_sem],
  "FGLS" = fit_fgls_sem$coefficients[parameter_names_sem],
  "GMM" = fit_gmm_sem$coefficients[parameter_names_sem],
  "OGMM" = fit_ogmm_sem$coefficients[parameter_names_sem]
)

            True         FGLS         GMM        OGMM
(Intercept)  0.0 -0.192428998 -0.19250015 -0.19255343
x1          -1.0 -1.079338334 -1.06289181 -1.06348122
x2           0.0  0.006390999  0.01368492  0.01109481
x3           1.0  0.975245554  0.98551854  0.98597574
lambda       0.6  0.556015812  0.55800375  0.56189516

The table places the sequential FGLS estimator beside the joint GMM estimators. Differences are expected in one sample because the procedures use different moment information and weighting matrices.

7.8 GS2SLS for the SAC Model

The SAC model combines an endogenous spatial lag with spatially autoregressive errors: \[ \begin{aligned} \mathbf y &= \rho\mathbf W\mathbf y + \mathbf X\boldsymbol\beta + \mathbf u,\\ \mathbf u &= \lambda\mathbf M\mathbf u + \boldsymbol\varepsilon. \end{aligned} \] The matrices \(\mathbf W\) and \(\mathbf M\) may be different. The estimator must therefore address two distinct problems: instrument \(\mathbf W\mathbf y\) and estimate the spatial-error transformation.

The programmed estimator follows four visible stages:

S2SLS: estimate \((\boldsymbol\beta,\rho)\) and obtain structural residuals;
initial GMM: estimate \(\lambda\) with an identity weighting matrix;
first efficient GMM and GS2SLS: update the weighting matrix, estimate \(\lambda\), spatially transform the equation, and re-estimate \((\boldsymbol\beta,\rho)\);
second efficient GMM: update \(\lambda\) using the GS2SLS residuals.

This section is longer than the preceding estimators because no single optimization contains all four stages.

7.8.1 Artificial SAC Data

set.seed(666)
side_sac <- 23L
n_sac <- side_sac^2
rho_sac <- 0.6
lambda_sac <- 0.6

# W uses rook contiguity; M uses queen contiguity.
neighbors_W_sac <- cell2nb(side_sac, side_sac, type = "rook")
W_sac <- nb2mat(neighbors_W_sac, style = "W")

neighbors_M_sac <- cell2nb(side_sac, side_sac, type = "queen")
M_sac <- nb2mat(neighbors_M_sac, style = "W")

x1_sac <- rnorm(n_sac)
x2_sac <- rnorm(n_sac)
x3_sac <- rnorm(n_sac)
X_sac <- cbind(
  "(Intercept)" = 1,
  "x1" = x1_sac,
  "x2" = x2_sac,
  "x3" = x3_sac
)
beta_sac <- c("(Intercept)" = 0, "x1" = -1, "x2" = 0, "x3" = 1)
sigma2_sac <- 1
epsilon_sac <- rnorm(n_sac, mean = 0, sd = sqrt(sigma2_sac))

R_sac <- diag(n_sac) - lambda_sac * M_sac
S_sac <- diag(n_sac) - rho_sac * W_sac
u_sac <- drop(solve(R_sac, epsilon_sac))
y_sac <- drop(solve(S_sac, X_sac %*% beta_sac + u_sac))

data_sac <- data.frame(
  y = y_sac,
  x1 = x1_sac,
  x2 = x2_sac,
  x3 = x3_sac
)

c(
  structural_equation = max(abs(S_sac %*% y_sac - X_sac %*% beta_sac - u_sac)),
  error_equation = max(abs(R_sac %*% u_sac - epsilon_sac))
)

structural_equation      error_equation 
       3.552714e-15        1.776357e-15

7.8.2 Quadratic Moments and Weighting Matrix

For a candidate \(\lambda\), the innovation implied by a structural residual vector \(\mathbf u\) is \[ \boldsymbol\varepsilon(\lambda) = (\mathbf I-\lambda\mathbf M)\mathbf u. \] The two normalized quadratic moments are \[ \overline{\mathbf m}_n(\lambda) = \frac{1}{n} \begin{pmatrix} \boldsymbol\varepsilon(\lambda)^\top\mathbf A_1 \boldsymbol\varepsilon(\lambda)\\ \boldsymbol\varepsilon(\lambda)^\top\mathbf A_2 \boldsymbol\varepsilon(\lambda) \end{pmatrix}. \] The original implementation divided the final quadratic criterion by \(n\) while leaving the moments unnormalized. That scaling does not change the minimizer, but normalizing the moments directly makes the code agree with the notation and keeps the criterion on a more stable scale as \(n\) changes.

moments.kp <- function(lambda, u, M, A1, A2) {
  n <- length(u)
  R <- diag(n) - lambda * M
  epsilon <- drop(R %*% u)

  rbind(
    "q1" = drop(crossprod(epsilon, A1 %*% epsilon)),
    "q2" = drop(crossprod(epsilon, A2 %*% epsilon))
  ) / n
}

Q.min.sac <- function(lambda, u, M, A1, A2, Upsilon, verbose = FALSE) {
  moments <- moments.kp(lambda, u = u, M = M, A1 = A1, A2 = A2)
  criterion <- drop(crossprod(moments, Upsilon %*% moments))

  if (verbose) {
    cat("criterion:", criterion, "lambda:", lambda, "\n")
  }

  criterion
}

make.Upsilon <- function(
  lambda,
  HH.inverse,
  H,
  Z,
  A1,
  A2,
  u,
  M,
  step = c("first", "second")
) {
  step <- match.arg(step)
  n <- nrow(H)

  A1_symmetric <- A1 + t(A1)
  A2_symmetric <- A2 + t(A2)
  R <- diag(n) - lambda * M
  Z_transformed <- R %*% Z
  u_transformed <- drop(R %*% u)

  alpha1 <- -drop(crossprod(Z_transformed, A1_symmetric %*% u_transformed)) / n
  alpha2 <- -drop(crossprod(Z_transformed, A2_symmetric %*% u_transformed)) / n

  if (step == "first") {
    HZ <- crossprod(H, Z) / n
    P <- HH.inverse %*% HZ %*% solve(t(HZ) %*% HH.inverse %*% HZ)
    R_inverse_transpose <- solve(t(R))
    a1 <- drop(R_inverse_transpose %*% H %*% P %*% alpha1)
    a2 <- drop(R_inverse_transpose %*% H %*% P %*% alpha2)
  } else {
    HZ <- crossprod(H, Z_transformed) / n
    P <- HH.inverse %*% HZ %*% solve(t(HZ) %*% HH.inverse %*% HZ)
    a1 <- drop(H %*% P %*% alpha1)
    a2 <- drop(H %*% P %*% alpha2)
  }

  epsilon <- drop(R %*% u)
  Sigma <- diag(epsilon^2)

  Psi11 <- tr(A1_symmetric %*% Sigma %*% A1_symmetric %*% Sigma) / (2 * n) +
    drop(crossprod(a1, Sigma %*% a1)) / n
  Psi12 <- tr(A1_symmetric %*% Sigma %*% A2_symmetric %*% Sigma) / (2 * n) +
    drop(crossprod(a1, Sigma %*% a2)) / n
  Psi22 <- tr(A2_symmetric %*% Sigma %*% A2_symmetric %*% Sigma) / (2 * n) +
    drop(crossprod(a2, Sigma %*% a2)) / n

  Upsilon <- matrix(c(Psi11, Psi12, Psi12, Psi22), nrow = 2L)

  list(
    Upsilon = Upsilon,
    a1 = a1,
    a2 = a2,
    epsilon = epsilon,
    P = P,
    Sigma = Sigma
  )
}

make.G <- function(u, M, A1, A2) {
  n <- length(u)
  u_lag <- drop(M %*% u)

  G <- matrix(0, nrow = 2L, ncol = 2L)
  # A1 is symmetric, so u_lag'(A1 + A1')u = 2 u_lag'A1u.
  # Do not multiply that symmetric expression by an additional factor of two.
  G[1L, 1L] <- crossprod(u_lag, (A1 + t(A1)) %*% u)
  G[1L, 2L] <- -crossprod(u_lag, A1 %*% u_lag)
  G[2L, 1L] <- crossprod(u_lag, (A2 + t(A2)) %*% u)
  G[2L, 2L] <- -crossprod(u_lag, A2 %*% u_lag)

  G / n
}

make.Upsilon() estimates the covariance matrix of the two quadratic moments while accounting for the fact that the structural residuals themselves contain estimated regression parameters. The vectors a1 and a2 carry that first-stage estimation effect into the moment covariance. In the first update, the formula requires \[ (\mathbf I-\lambda\mathbf M^\top)^{-1}\mathbf H\mathbf P\boldsymbol\alpha_r = \mathbf R^{-\top}\mathbf H\mathbf P\boldsymbol\alpha_r. \] The original code used crossprod(R_inverse_transpose, H), which equals \(\mathbf R^{-1}\mathbf H\), not \(\mathbf R^{-\top}\mathbf H\). This transpose error can affect the first efficient estimate of \(\lambda\), the subsequent GS2SLS estimates, and the final covariance matrix. The corrected code multiplies R_inverse_transpose %*% H directly.

make.G() constructs the derivative terms used for inference on \(\lambda\). For the first moment, \(\mathbf A_1\) is symmetric by construction. Therefore, \[ \mathbf u_L^\top(\mathbf A_1+\mathbf A_1^\top)\mathbf u = 2\mathbf u_L^\top\mathbf A_1\mathbf u. \] The original code multiplied the left-hand side by another factor of two. That did not change the point estimates, because make.G() is used only by the covariance method, but it distorted the standard error of \(\lambda\). The code above removes the duplicated factor.

7.8.3 GS2SLS Estimation Function

gstsls.sac <- function(
  formula,
  data,
  W,
  M = NULL,
  instruments = 2L,
  verbose = FALSE
) {
  call <- match.call(expand.dots = TRUE)
  components <- spatial_model_components(
    formula = formula,
    data = data,
    A = W,
    matrix_name = "W"
  )

  y <- components$y
  X <- components$X
  W <- components$A
  n <- components$n

  if (is.null(M)) {
    M <- W
  } else {
    M <- validate_spatial_matrix(M, n = n, matrix_name = "M")
  }

  different_error_weights <- !isTRUE(all.equal(M, W))
  lambda_interval <- spatial_parameter_interval(M)

  # Step 1a: S2SLS for delta = (beta', rho)'.
  WX <- make.H(W = W, X = X, l = instruments)
  H <- cbind(X, WX)

  if (different_error_weights) {
    intercept <- which(colnames(X) == "(Intercept)")
    X_without_intercept <- if (length(intercept) > 0L) {
      X[, -intercept, drop = FALSE]
    } else {
      X
    }
    H <- cbind(H, M %*% X_without_intercept, M %*% WX)
  }

  H <- independent_columns(H)
  Wy <- drop(W %*% y)
  Z <- cbind(X, "Wy" = Wy)

  if (n <= ncol(Z)) {
    stop("The sample must contain more observations than structural parameters.")
  }
  if (ncol(H) < ncol(Z) || qr(crossprod(H, Z))$rank < ncol(Z)) {
    stop("The SAC first-stage rank condition fails.", call. = FALSE)
  }

  HH <- crossprod(H)
  HH_inverse <- solve(HH)
  P_H <- H %*% solve(HH, t(H))
  Z_hat <- P_H %*% Z
  delta_s2sls <- drop(solve(crossprod(Z_hat), crossprod(Z_hat, y)))
  names(delta_s2sls) <- colnames(Z)
  u_s2sls <- drop(y - Z %*% delta_s2sls)

  # Step 1b: initial identity-weighted GMM estimate of lambda.
  lambda_start <- drop(
    crossprod(M %*% u_s2sls, u_s2sls) / crossprod(u_s2sls)
  )
  lambda_start <- clamp_to_interval(lambda_start, lambda_interval)

  MM <- crossprod(M)
  A1 <- MM - diag(diag(MM))
  A2 <- M

  fit_initial_lambda <- nlminb(
    start = lambda_start,
    objective = Q.min.sac,
    u = u_s2sls,
    M = M,
    A1 = A1,
    A2 = A2,
    Upsilon = diag(2L),
    verbose = verbose,
    lower = lambda_interval["lower"],
    upper = lambda_interval["upper"]
  )

  # Step 1c: efficient GMM based on the S2SLS residuals.
  lambda_step1_initial <- unname(fit_initial_lambda$par)
  Upsilon_step1 <- make.Upsilon(
    lambda = lambda_step1_initial,
    HH.inverse = HH_inverse,
    H = H,
    Z = Z,
    A1 = A1,
    A2 = A2,
    u = u_s2sls,
    M = M,
    step = "first"
  )$Upsilon

  fit_efficient_step1 <- nlminb(
    start = lambda_step1_initial,
    objective = Q.min.sac,
    u = u_s2sls,
    M = M,
    A1 = A1,
    A2 = A2,
    Upsilon = invert_moment_covariance(
      Upsilon_step1,
      matrix_name = "The first SAC moment covariance matrix"
    ),
    verbose = verbose,
    lower = lambda_interval["lower"],
    upper = lambda_interval["upper"]
  )

  # Step 2a: GS2SLS after spatially filtering y and Z.
  lambda_hat_step1 <- unname(fit_efficient_step1$par)
  R_hat_step1 <- diag(n) - lambda_hat_step1 * M
  y_transformed <- drop(R_hat_step1 %*% y)
  Z_transformed <- R_hat_step1 %*% Z
  Z_transformed_hat <- P_H %*% Z_transformed
  delta_gs2sls <- drop(solve(
    crossprod(Z_transformed_hat),
    crossprod(Z_transformed_hat, y_transformed)
  ))
  names(delta_gs2sls) <- colnames(Z)
  u_gs2sls <- drop(y - Z %*% delta_gs2sls)

  # Step 2b: efficient GMM based on the GS2SLS residuals.
  Upsilon_step2 <- make.Upsilon(
    lambda = lambda_hat_step1,
    HH.inverse = HH_inverse,
    H = H,
    Z = Z,
    A1 = A1,
    A2 = A2,
    u = u_gs2sls,
    M = M,
    step = "second"
  )$Upsilon

  fit_efficient_step2 <- nlminb(
    start = lambda_hat_step1,
    objective = Q.min.sac,
    u = u_gs2sls,
    M = M,
    A1 = A1,
    A2 = A2,
    Upsilon = invert_moment_covariance(
      Upsilon_step2,
      matrix_name = "The second SAC moment covariance matrix"
    ),
    verbose = verbose,
    lower = lambda_interval["lower"],
    upper = lambda_interval["upper"]
  )

  convergence_codes <- c(
    initial = fit_initial_lambda$convergence,
    efficient.first = fit_efficient_step1$convergence,
    efficient.second = fit_efficient_step2$convergence
  )
  if (any(convergence_codes != 0L)) {
    warning("At least one SAC optimization stage did not report convergence.")
  }

  lambda_final <- unname(fit_efficient_step2$par)
  coefficients <- c(delta_gs2sls, "lambda" = lambda_final)

  structure(
    list(
      coefficients = coefficients,
      call = call,
      X = X,
      H = H,
      HH.inverse = HH_inverse,
      Z = Z,
      y = y,
      W = W,
      M = M,
      A1 = A1,
      A2 = A2,
      G = make.G(u = u_gs2sls, M = M, A1 = A1, A2 = A2),
      residuals = u_gs2sls,
      lambda.interval = lambda_interval,
      steps = list(
        delta.s2sls = delta_s2sls,
        lambda.initial = unname(fit_initial_lambda$par),
        lambda.efficient.first = lambda_hat_step1,
        delta.gs2sls = delta_gs2sls,
        lambda.efficient.second = lambda_final
      ),
      optimizers = list(
        initial = fit_initial_lambda,
        efficient.first = fit_efficient_step1,
        efficient.second = fit_efficient_step2
      )
    ),
    class = "mygs2sls"
  )
}

The returned object stores every intermediate estimate in steps. This is not needed to calculate the final coefficient vector, but it makes the multistep algorithm auditable. In particular, a problematic final estimate of \(\lambda\) can be traced to the identity-weighted step, the first efficient update, or the residuals produced by GS2SLS.

The bounds now come from the spectrum of \(\mathbf M\) rather than the hard-coded interval \((-0.9,0.9)\). The old interval happened to contain the true parameter in this simulation, but it was not tied to the matrix supplied by the user.

The covariance estimator implements the block expression in Equation 6.82.

vcov.mygs2sls <- function(object, ...) {
  coefficients <- object$coefficients
  Z <- object$Z
  H <- object$H
  n <- nrow(Z)
  k <- ncol(Z)
  p <- ncol(H)

  delta_hat <- coefficients[seq_len(k)]
  lambda_hat <- coefficients[k + 1L]
  u_gs2sls <- drop(object$y - Z %*% delta_hat)

  fit_lambda <- make.Upsilon(
    lambda = lambda_hat,
    HH.inverse = object$HH.inverse,
    H = H,
    Z = Z,
    A1 = object$A1,
    A2 = object$A2,
    u = u_gs2sls,
    M = object$M,
    step = "second"
  )

  Psi_delta_delta <- crossprod(H, fit_lambda$Sigma %*% H) / n
  Psi_delta_lambda <- crossprod(
    H,
    fit_lambda$Sigma %*% cbind(fit_lambda$a1, fit_lambda$a2)
  ) / n

  Psi_joint <- rbind(
    cbind(Psi_delta_delta, Psi_delta_lambda),
    cbind(t(Psi_delta_lambda), fit_lambda$Upsilon)
  )
  Psi_joint <- (Psi_joint + t(Psi_joint)) / 2

  Psi_lambda_inverse <- invert_moment_covariance(
    fit_lambda$Upsilon,
    matrix_name = "The final SAC moment covariance matrix"
  )
  J <- object$G %*% c(1, 2 * lambda_hat)
  lower_bread <- solve(
    crossprod(J, Psi_lambda_inverse %*% J)
  ) %*% t(J) %*% Psi_lambda_inverse

  bread <- matrix(0, nrow = k + 1L, ncol = p + 2L)
  bread[seq_len(k), seq_len(p)] <- t(fit_lambda$P)
  bread[k + 1L, (p + 1L):(p + 2L)] <- lower_bread

  bread %*% Psi_joint %*% t(bread) / n
}

summary.mygs2sls <- function(
  object,
  digits = max(3L, getOption("digits") - 3L),
  ...
) {
  estimates <- drop(object$coefficients)
  standard_errors <- sqrt(diag(vcov(object)))
  statistics <- estimates / standard_errors
  p_values <- 2 * pnorm(-abs(statistics))

  coefficient_table <- cbind(
    "Estimate" = estimates,
    "Std. Error" = standard_errors,
    "z value" = statistics,
    "Pr(>|z|)" = p_values
  )

  structure(
    list(
      coefficients = coefficient_table,
      digits = digits,
      call = object$call
    ),
    class = "summary.mygs2sls"
  )
}

print.summary.mygs2sls <- function(x, digits = x$digits, ...) {
  cat("\nCall:\n")
  cat(paste(deparse(x$call), collapse = "\n"), "\n", sep = "")
  cat("\nCoefficients:\n")
  printCoefmat(
    x$coefficients,
    digits = digits,
    P.values = TRUE,
    has.Pvalue = TRUE
  )
  invisible(x)
}

7.8.4 Estimation and Package Verification

# Cache revision: recompute after the estimator and DGP audit.

fit_gs2sls_sac <- gstsls.sac(
  y ~ x1 + x2 + x3,
  data = data_sac,
  W = W_sac,
  M = M_sac
)
summary(fit_gs2sls_sac)


Call:
gstsls.sac(formula = y ~ x1 + x2 + x3, data = data_sac, W = W_sac, 
    M = M_sac)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.05118    0.06762   0.757 0.449139    
x1          -1.00304    0.04398 -22.805  < 2e-16 ***
x2           0.07097    0.04275   1.660 0.096940 .  
x3           0.98948    0.04216  23.468  < 2e-16 ***
Wy           0.71231    0.04321  16.485  < 2e-16 ***
lambda       0.35028    0.10244   3.419 0.000628 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sac_step_comparison <- rbind(
  "True" = c(beta_sac, "rho" = rho_sac, "lambda" = lambda_sac),
  "Initial S2SLS" = c(fit_gs2sls_sac$steps$delta.s2sls, "lambda" = NA_real_),
  "GS2SLS + final GMM" = fit_gs2sls_sac$coefficients
)
round(sac_step_comparison, 4L)

                   (Intercept)     x1     x2     x3    rho lambda
True                    0.0000 -1.000 0.0000 1.0000 0.6000 0.6000
Initial S2SLS           0.0534 -0.988 0.0809 1.0040 0.7162     NA
GS2SLS + final GMM      0.0512 -1.003 0.0710 0.9895 0.7123 0.3503

c(
  lambda.initial = fit_gs2sls_sac$steps$lambda.initial,
  lambda.efficient.first = fit_gs2sls_sac$steps$lambda.efficient.first,
  lambda.efficient.second = fit_gs2sls_sac$steps$lambda.efficient.second
)

         lambda.initial  lambda.efficient.first lambda.efficient.second 
              0.3385977               0.3400930               0.3502774

The first table shows how spatially filtering the equation changes the estimate of \((\boldsymbol\beta,\rho)\). The second isolates the three estimates of \(\lambda\) produced by the successive GMM updates.

The original example compares the estimates with sphet::spreg() using step1.c = TRUE. Some versions of sphet fail internally when assigning names to the coefficient vector for this specification. The code therefore preserves the original call and, only if that call fails, retries the same model with step1.c = FALSE. If both calls fail, the benchmark is skipped without stopping the render of the chapter.

# Cache revision: recompute after the estimator and DGP audit.

if (requireNamespace("sphet", quietly = TRUE)) {
  listw_W_sac <- spdep::mat2listw(W_sac, style = "W")
  listw_M_sac <- spdep::mat2listw(M_sac, style = "W")

  fit_sphet_sac <- tryCatch(
    sphet::spreg(
      y ~ x1 + x2 + x3,
      data = data_sac,
      listw = listw_W_sac,
      listw2 = listw_M_sac,
      model = "sarar",
      het = TRUE,
      step1.c = TRUE
    ),
    error = function(error_step1c) {
      message(
        "sphet::spreg() failed with step1.c = TRUE: ",
        conditionMessage(error_step1c),
        " Retrying with step1.c = FALSE."
      )

      tryCatch(
        sphet::spreg(
          y ~ x1 + x2 + x3,
          data = data_sac,
          listw = listw_W_sac,
          listw2 = listw_M_sac,
          model = "sarar",
          het = TRUE,
          step1.c = FALSE
        ),
        error = function(error_default) {
          message(
            "The sphet benchmark was skipped because the fallback call also failed: ",
            conditionMessage(error_default)
          )
          NULL
        }
      )
    }
  )

  if (!is.null(fit_sphet_sac)) {
    summary(fit_sphet_sac)
  }
} else {
  message("Install the sphet package to run the SAC benchmark.")
}

sphet::spreg() failed with step1.c = TRUE: length of 'dimnames' [1] not equal to array extent Retrying with step1.c = FALSE.

====================================================
====================================================
               GS2SLS SARAR Model
====================================================
====================================================

Call:
sphet::spreg(formula = y ~ x1 + x2 + x3, data = data_sac, listw = listw_W_sac, 
    listw2 = listw_M_sac, model = "sarar", het = TRUE, step1.c = FALSE)

Residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-2.92546 -0.74662  0.04502  0.00348  0.67160  3.13620 

Coefficients:
             Estimate Std. Error  t-value  Pr(>|t|)    
(Intercept)  0.051194   0.067618   0.7571 0.4489913    
x1          -1.002991   0.043984 -22.8035 < 2.2e-16 ***
x2           0.071008   0.042755   1.6608 0.0967509 .  
x3           0.989534   0.042163  23.4693 < 2.2e-16 ***
lambda       0.712337   0.043208  16.4863 < 2.2e-16 ***
rho          0.350223   0.102443   3.4187 0.0006292 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald test that rho and lambda are both zero:
 Statistics: 190.26 p-val: 2.7909e-43

7.9 What We Have Learned

The estimators in this chapter differ in their details, but their programs share a common architecture.

Construct the correct residual object. For the SLM this is the structural innovation \((\mathbf I-\rho\mathbf W)\mathbf y-\mathbf X\boldsymbol\beta\). For the SEM and SAC models, the structural residual must also be transformed by \(\mathbf I-\lambda\mathbf M\).
Keep moment ordering and parameter ordering synchronized. The SLM GMM code stores rho first; the SEM GMM code stores lambda last. The rows and columns of every derivative and covariance matrix must follow those orders.
Separate estimation from inference. The SAC factor-of-two error in make.G() illustrates why correct point estimates do not guarantee correct standard errors.
Verify analytical derivatives numerically. A gradient check away from the optimum tests the residual formula, the Jacobian, the parameter order, and the sign of the optimized criterion simultaneously.
Treat package estimates as benchmarks. Agreement is valuable evidence, but it does not replace checking the instruments, moment assumptions, parameter space, and covariance formula.

Reading this chapter together with the theory

The functions expose the matrix calculations behind the estimators. Their validity depends on the moment, identification, and regularity conditions stated in Chapter 6. Numerical optimization cannot repair invalid instruments, failed identification, or moment conditions that do not hold for the maintained disturbance process.