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`](http://cran.r-project.org/web/packages/spdep/index.html) 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?