# How to bootstrap a (spatial) lag operator?

In relation with this question of mine, which goes beyond what I am asking now.

Say one has the following data generating process (DGP) (known as the SAR ou spatial lag model in the literature):

$\boldsymbol{y} = \rho \boldsymbol{W}\boldsymbol{y} + \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{u}$

where $\boldsymbol{y}$ is a $n \times 1$ endogenous vector, $\boldsymbol{X}$ an exogenous $n \times k$ matrix, $\boldsymbol{\beta}$ a $k \times 1$ vector of coefficients, $\boldsymbol{W}$ a $n \times n$ spatial lag operator, $\rho$ a scalar coefficient of autoregression and $\boldsymbol{u}$ a vector of errors, wished iid normal.

How would you bootstrap $\boldsymbol{W}$ ? In the case of time series, assume that one has

$$\boldsymbol{y} = \begin{pmatrix} 2013 \\ 2014 \\ 2015 \\ 2016 \end{pmatrix}$$

and the following backshift (linear) operator

$$\boldsymbol{W} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{pmatrix}$$

Let's denote the (block) bootstrapping resampling process by an operator, say, $\mathcal{B}(.)$. What does one do in the times series case ?

$\mathcal{B}(\boldsymbol{W}\boldsymbol{y})$

OR

$\mathcal{B}(\boldsymbol{W})\boldsymbol{y}$ (equivalent to the previous one if one bootstraps only over the rows of $\boldsymbol{W}$)

OR

$\boldsymbol{W}\mathcal{B}(\boldsymbol{y})$

OR none of those

?

I would say $\mathcal{B}(\boldsymbol{W}\boldsymbol{y})$ but this has weird implications in space when computing the likelihood (cf my question). But let us forget these weird implications for a moment. What do we do in the time series case ?

...One of these weird implications is that $\boldsymbol{W}$ may not be hollow anymore if one only resamples its rows...

A python toy-example

import numpy as np

W = np.zeros((4,4))
W[1,0] = 1.
W[2,1] = 1.
W[3,2] = 1.

y = np.arange(2013,2017).reshape((4,1))
Wy= np.dot(W,y)

def B(obj):
"""
Operator which only resamples rows
In the case in which rows and columns are resampled
we would do:
return obj[idxs].T[idxs].T # (right ?)
"""
idxs = [2,3,0,1]
return obj[idxs]

print B(Wy)
print np.dot(B(W),y)
print np.dot(W,B(y))


python toy-example which shows that $\mathcal{B}(\boldsymbol{W})$ is not hollow

import numpy as np
from random import random as rd

W = np.array([[0.  ,0.,0.,0.],
[rd(),0.  ,0.,0.],
[rd(),rd(),0.  ,0.],
[rd(),rd(),rd(),0.]])
W += W.T                       #to make W be symmetric
print B(W)                     #not hollow !!!
print B(B(W).T).T              #while B(B(W).T).T is !


### Update

Optimal block size for variance estimation by a spatial block bootstrap is very likely to help me out.