I would like to estimate a model of the following form:

$$ y = \sigma G y + \beta X + \delta G^* X + \epsilon $$

where $G$ and $G^*$ are quadratic adjacency matrices, $y$ is a vector of a dependent variable subject to peer effects and $X$ is a matrix of controls/exogenous characteristics. $G$ is a binary matrix with entries equal to 1 if two individuals are peers of each other; as such, $G$ is also symmetric with a 0-diagonal. $G^*$ is the row-normalized version of $G$ (in fact, it's a stacked matrix with the row-normalized $G$ along the diagonal an 0 in the off-diagonals because $X$ is a matrix).

Econometricians will immediately recognize the similarity to spatial simultaneous autoregressive lag models, whose general form is

$$ y = \rho W y + \beta X + \epsilon $$

(with $W = G$). Because of the high similarity, I tried Roger Bivand's spdep package for R. The regression command is lagsarlm(), but when I am not mistaken, it only estimates

$$ y = \sigma G y + \beta X + \epsilon $$

using the following command:

lagsarlm(y ~ x1 + x2, data=reg_data, listw=G, method="eigen",
    quiet=FALSE, zero.policy = FALSE, tol.solve=1e-14)

which works and also returns results, but obviously an entire term is missing. How can I incorporate the $ \delta G^* X$ term? Are there packages the serve my needs better?