Spatial Durbin Model

mwe.R

library(rstan)
## options(mc.cores = parallel::detectCores())

## Data should be here:  <https://drive.google.com/open?id=1BUV8mChrZ0UkyW2G4EHB5wc1taUpMZUM>
load('mwe.Rdata')
## dataf: The data; response is `w_use`
## W_3nn:  row-normalized spatial weights matrix, based on 3 nearest neighbors, in triplet sparse form
## W_dnn:  row-normalized spatial weights matrix, distance-based, in triplet sparse form

parsed = stanc(file = 'sdm.stan', verbose = TRUE)
compiled = stan_model(stanc_ret = parsed, verbose = FALSE)

tictoc::tic()
## ~300 sec on my laptop
samples = sampling(compiled, data = list(N = nrow(dataf), 
                                         p = ncol(dataf) - 1,
                                         y = dataf$w_use, 
                                         X = dataf[-1],
                                         # W = as.matrix(W_3nn)
                                         Wn = length(W_3nn@i), 
                                         Wi = W_3nn@i+1, 
                                         Wj = W_3nn@j+1,
                                         Ww = W_3nn@x
                                         ),
                   chains = 2, iter = 1000, 
                   control = list(adapt_delta = 0.8))
tictoc::toc()

library(spdep)
freq = lagsarlm(w_use ~ hispanicP + densityE, data = dataf, 
                type = 'Durbin',
                listw = mat2listw(W_3nn))

## Compare estimates
summary(samples, pars = c('rho', 'alpha', 'beta', 'gamma'))$summary
freq

sdm.stan

// This gist contains Stan code for a Spatial Durbin Model
// The code is a bit slow, but appears to converge well enough and yields the same point estimates as spdep::lagsarlm(type = 'Durbin')
// As written, the spatial weights matrix W is passed as a triple sparse matrix; 
// in R, use the `dgTMatrix` class from the `Matrix` package for this. 
data {
    int<lower=0> N; // number of observations / locations
    int<lower=0> p; // number of covariates
    matrix[N, p] X; // covariate matrix, w/o intercept
    vector<lower=0>[N] y; // response
    // matrix[N, N] W; // spatial weights matrix (non-sparse version)
    int Wn; // num. non-zero entries in spatial weights
    int Wi[Wn]; // i index of spatial weights
    int Wj[Wn]; // j index of spatial weights
    vector[Wn] Ww; // spatial weight
    
}
parameters {
    real<lower = -1, upper = 1> rho; // spatial autoregression coef.
    real alpha; // intercept
    vector[p] beta; // covariate effects
    vector[p] gamma; // lagged covariate effects
    real<lower = 0> sigma; // sd of noise
}
transformed parameters {
    vector[N] y_hat;
    
    // Non-sparse version
    // y_hat = rho*W*y + alpha*rep_vector(1, N) + X*beta + W*X*gamma;
    // Sparse version
    {
        vector[N] y_hat_local;
        vector[N] y_hat_lags;
        
        y_hat_local = alpha*rep_vector(1, N) + X*beta;
        y_hat_lags = rep_vector(0, N);
        for (w_idx in 1:Wn) {
            y_hat_lags[Wi[w_idx]] += rho*Ww[w_idx]*y[Wj[w_idx]];
            for (col_idx in 1:p) {
            y_hat_lags[Wi[w_idx]] += Ww[w_idx]*X[Wj[w_idx], col_idx]*gamma[col_idx];
            }
        }
        
        y_hat = y_hat_local + y_hat_lags;
    }
}
model {
    y ~ normal(y_hat, sigma);
    
    // priors
    alpha ~ normal(0, 3);
    beta ~ normal(0, 3);
    gamma ~ normal(0, 3);
    
    sigma ~ cauchy(0, 3);
}