A Generalized Method of Moments Estimator for a Spatial Panel Model with an Endogenous Spatial Lag and Spatial Moving Average Errors

Abstract

This paper proposes a new generalized method of moments (GMM) estimator for spatial panel models with spatial moving average errors combined with a spatially autoregressive dependent variable. Monte Carlo results are given suggesting that the GMM estimator is consistent. The estimator is applied to English real estate price data.

Abstract

Une méthode généralisée d'estimateur de moments pour un modèle de panel spatial avec un décalage endogène spatial et des erreurs spatiales de type moyenne mobile

Cette étude propose un nouvel estimateur GMM pour des modèles de panel spatial avec des erreurs spatiales de type moyenne mobile combiné à une variable de dépendance spatiale autorégressive. Les résultats de Monte Carlo fournis suggèrent que l'estimateur GMM est cohérent. L'estimateur s'applique à des données sur des prix d'immobilier anglais.

Abstract

Un estimador de método generalizado de momentos para un modelo de panel espacial con retardo espacial endógeno y errores espaciales de media móvil

Este estudio propone un nuevo estimador GMM para modelos de panel espacial con errores espaciales de media móvil combinado con una variable dependiente autorregresiva. Se indican los resultados de Monte Carlo que revelan la coherencia del estimador GMM. El estimador se aplica a los datos de precios en inmobiliarias inglesas.

Keywords:

• Moving averages
• GMM
• real estate
• spatial econometrics
• panel data

Keywords:

• C21
• R12
• R31

Introduction

There is a growing literature dedicated to the analysis of panel data with spatial dependence, with various different approaches suggested. Probably the most useful starting point in the spatial econometrics literature is Anselin (1988), and among some of the more recent contributions, such as Conley (1999), Chen & Conley (2001), Baltagi et al. (2003), Druska & Horrace (2003), Elhorst (2003), Baltagi (2005), Baltagi & Li (2006) and Pinkse et al. (2006), we highlight the work of Kapoor et al. (2007), which generalizes the generalized moments estimators of Kelejian & Prucha (1999) to a panel data model with spatially and temporally correlated error components, and which provides a feasible generalized least squares procedure for the regression parameters, and formal large sample results for their estimators.

This paper draws on their contribution, which provides the necessary theoretical, computational, and mathematical background for the present paper. Given this context, the specific innovatory aspects of the current paper are:

1. the extension of the generalized moments estimators (GMM) estimation procedure to allow a spatial moving average (MA) error process rather than the spatial autoregressive process that has been the focus of attention thus far in the literature;

2. the extension of the methodology to incorporate an endogenous spatial lag, so that spatial dependence is not solely via the error process;

3. application of the method to real panel data involving real estate prices in England.

To summarize, the paper extends the scope of the approach suggested by Kapoor et al. (2007) by allowing different forms of spatial interaction for panel data.

2 The Model

Consider the N location cross-sectional time t regression specification

1

in which Y(t) is an N×1 vector of observations of the dependent variable, H(t) is the N×k matrix of regressors with full column rank, γ is a k×1 vector of parameters, and u(t) is an N×1 vector produced by a random error process. Also W ^E is an N×N matrix of non-stochastic time constant weights which defines the interdependence of Y(t) across areas, so that W ^E Y is an N×1 vector commonly referred to as an endogenous spatial lag and λ is a scalar parameter. Following Kapoor et al. (2007) and Kelejian & Prucha (1998), all of the diagonal elements of W ^E are zero, and is non-singular. Also W ^E is uniformly bounded in absolute value, meaning that a constant c exists such that and . Likewise, the elements of H(t) are uniformly bounded in absolute value.

The most widely used approach to modelling spatial error dependence involving N locations is to assume that in each period , in which u(t) is a vector of errors at time t, ρ is a parameter, W is also an N×N matrix of non-stochastic weights which defines the error interaction across areas and ξ(t) is an N×1 vector of time t innovations. All the diagonal elements of W are zero, (I–ρW) is non-singular and W is also uniformly bounded in absolute value. This is referred to as a spatial autoregressive (AR) process and implies complex interdependence between locations, so that a shock at location j is transmitted to all other locations, as indicated by the expansion of

2

which, assuming ∣ρ∣<1and a row-standardized W matrix with row sums equal to 1, is

3

in which W ⁰=I, W ² is the matrix product of W and W, and W ⁱ is the matrix product of W ⁱ⁻¹ and W. The effect of shock at j is therefore felt directly at j, and there is an indirect effect due to ρWξ(t) which affects only those location pairs for which there is a non-zero element on the W matrix. If W were a contiguity matrix we might think of these as local effects. The global effect of a shock occurs because it is transmitted also to locations that are ‘neighbours of neighbours’ via the powers of W. Note that the effect rebounds. A shock to j affects the neighbours, and the neighbours of the neighbours, and eventually works its way back to j. In other words, the full effect of a shock to j is not simply the shock itself, but the initial shock plus the feedback from the other locations.

In contrast, the MA error process,1 which is the subject of this paper, is

4

so that a shock at location j will only affect the directly interacting locations as given by the non-zero elements in W. Hence shock-effects are local rather than global. Since we are considering a panel with T periods rather than purely cross-sectional data, we omit t to indicate that the observations are stacked. Hence

5

in which Y is a TN×1 vector of observations, X is a TN×(1+k) matrix of regressors, comprising the TN×1 vector (I _T ⊗W ^E )Y, which is the endogenous spatial lag, and H is a TN×k matrix of (exogenous) regressors. In addition, given ξ is an NT×1 vector of innovations, I _T is a T×T diagonal matrix with 1s on the main diagonal and 0s elsewhere, I _N is a similar N×N diagonal matrix, then I _TN =I _T ⊗I _N is a TN×TN diagonal matrix with 1s on the main diagonal and 0s elsewhere and the NT×1 vector u is given by the MA process

6

Regarding the error components in space-time, time dependency is introduced into the innovations ξ by specifying unobserved permanent unit-specific error components µ together with transient error components ν, where

7

8

and

9

Thus µ is the N×1 vector of errors specific to each area, ν is the NT×1 vector of errors specific to each area and time, ι _T is a T×1 matrix with 1s, and is a TN×N matrix equal to T stacked I _N matrices.

For ν, it is assumed that and we also make the standard assumption that the errors have finite fourth moments () to ensure a finite domain for estimation. Likewise for µ, and . Also, we assume that the error components are independent, hence , and each of the two error components µ and ν is subject to the same spatial moving average process, since

10

For areas i, j and times t, s:

11

For the purposes of estimation, it is useful to represent the TN×TN innovations variance–covariance matrix using the matrices Q ₀ and Q ₁ defined2 as follows

12

in which J _T is a T×T matrix of 1s, and

13

It follows that and

14

in which

15

3 The Moments Equations

Consider the TN×1 vector of residuals u

16

and

17

Pre-multiplying by Q ₀ gives

18

and

19

with

20

Likewise

21

and

22

Similarly

23

24

25

Also

26

27

28

To obtain the expectations of these variables, we know (see Kapoor et al., 2007) that

29

30

31

and using similar arguments it is also possible to show that

32

since

33

Likewise3

34

Also

35

By analogy, and following Kapoor et al. (2007),

36

37

38

39

40

41

Ignoring the expectations, we put these equations together using the 3×3 matrices Γ and , the 3×1 vectors φ and , and the 3×1 vectors γ and , using , , and so that

42

and

43

where

4 Estimation

The estimation procedure comprises three stages. At stage 1, because of the presence of the spatial lag, we obtain4 IV estimates of b and hence residuals . In stage 2 we use these IV residuals to obtain the estimates g and of γ and , and denoting Γ and by G and we have the sample counterpart of equation (42), which is

44

in which is a vector of residuals, and the non-linear least squares estimators are given by

From the estimated ρ and one can obtain the estimate of indirectly, using the fact that, for MA errors, , hence . Since , then

45

and therefore

46

In practice, to obtain a direct estimate of , we also use the sample counterpart of equation (43), which is

47

in which is a vector of residuals, and obtain as the minimum5 of F ₁+F ₂, where F ₁= and F ₂=.

In general, the variances associated with F ₁ and F ₂ differ, and Kapoor et al. (2007) suggest weighting to allow for this. However, in the Monte Carlo simulations that follow, for simplicity we have not introduced differential weighting. In the analogous situation examined by Kapoor et al. (2007) they note that giving equal weight to all six moments equations does give consistent estimates. While the small sample behaviour in the AR case is the worse of the alternative weighting schemes they examine, it seems appropriate commencing with MA errors to initially explore the behaviour of the simplest approach prior to more elaborate methods, which could be the subject of further research.

In the third stage, because the errors are not constant, the appropriate method is generalized least squares (GLS), estimated by IV to also allow for the presence of the endogenous spatial lag. The estimated error covariance matrix is obtained using from stage 2, but first is used to perform a Cochrane–Orcutt (C-O)-type transformation to account for the spatial dependence in the residuals.

Normally with C-O the assumption is an autoregressive error process, hence , in which case one pre-multiplies through by to obtain the innovations ξ. However, the MA error process requires pre-multiplication6 by the inverse7 to obtain ξ, thus

48

In both the first and third stages, to carry out the IV estimation, as instruments we employ a linearly independent subset of the exogenous variables, so that Z is a TN× matrix of instruments. Assume matrices X and Z are full column rank with , and following what is evidently a comparatively robust approach for IV estimation with non-spherical disturbances (Bowden & Turkington, 1984), calculate , which is a symmetric matrix ( is idempotent) and hence

49

The estimated variance–covariance matrix of the parameters is given by

50

Greene (2003) also gives the equivalent of (49) and (50) as generalized methods of moments (instrumental variables) estimators with non-spherical disturbances. The standard errors of the are given by the squares roots of the values on the main diagonal of , which allows ‘t-ratios’ to be calculated for purposes of inference.

5 Example 1: the Data-generating Process

In this first example the data are purely artificial, and correspond to model (6), which is repeated here for convenience:

51

52

In the MA error process, W is a contiguity matrix8 on a × square. Matrix W is standardized by dividing each row cell by its row total, so that the maximum and minimum eigenvalues are 1 and –1. Also in practice, for simplicity we assume that W ^E =W.

Matrix H had columns equal to the TN×1 vectors ι _TN , H ₁, H ₂, and H ₃, in which ι _TN is a TN×1 vector with 1s. To obtain each H, first generate time t=0, N×1 vectors H ₁(0), H ₂(0), H ₃(0) by sampling at random from a uniform (rectangular) distribution with minimum equal to 0 and maximum equal to 1. Then for t equal to 1 … T, , in which π ₁∼N(0,1), and likewise for H ₂(t), H ₃(t) using π ₂∼N(0,1) and π ₃∼N(0,1). Then stacking these N×1 vectors we obtain H ₁, H ₂ and H ₃. In his way the exogenous variables in H have some time dependency, as seems reasonable for panel data. Once generated, the variables H ₁, H ₂ and H ₃ remain fixed. Also, in practice T=2 and T=4 below.

Given the exogenous variables, we next obtain the innovations vector NT×1 vector ξ. The innovations vector depends on the N×1 vector µ obtained by sampling from an distribution and on the NT×1 vector ν obtained by sampling from an distribution, so that . This is repeated for each iteration k=1 … K, to obtain

53

Given Y, W and H, K estimates are obtained of the known parameters and using the three-stage method outlined above. This is achieved by using instruments Z for the endogenous spatial lag comprising the exogenous variables, H, together with the TN×1 vector comprising T stacked identical time ‘zero’ N×1 spatial lag vectors WY(0), which is assumed to be exogenous with respect to the endogenous lag .

5.1 Monte Carlo Results

Monte Carlo results are given both here and in more detail in Appendix B. Those given here are illustrative, while those in Appendix B provide more substantive empirical evidence of the consistency of the estimator. In this first example, the values ρ=−0.25, λ=0.75, γ ₀=1, γ₁=10, γ ₂=10, γ ₃=10, and are used to generate Y, and the three-stage estimation method employed K=100 times9 gives a set of K estimates and (k=1 … K) of these parameters. Table 1 summarizes the parameter estimate distributions. It is evident that the parameter estimate means are close to the true values, although it is shown later (see Appendix B) that there is evidence of small sample bias in the estimator of ρ, although it is apparently consistent. The distributions are relatively symmetrical, and on the whole have a degree of kurtosis consistent with the normal distribution. To formally test the null of normality, the K estimates are divided into groups with upper and lower bounds defined so that each group has approximately the same observed frequency (O _i ). These observed frequencies are then compared to expected frequencies (E _i ) calculated using the data to obtain maximum-likelihood estimates of the normal mean and variance.10 The test statistic is then referred to the distribution. It is apparent that none of the distributions differs significantly from normal, using the upper 5% point (14.07) of the distribution.

Table 1 . Estimated parameter distributions

T=2		Mean	SD	W ^E =W (Rook's case) Skewness	H 0∼U[0,1] Kurtosis	K=100 Normality
λ=0.75		0.7486	0.0071	−0.04	−0.21	4.44
γ 0 =1		1.0960	0.4594	0.19	−0.11	8.58
γ 1 =10		10.0013	0.0667	−0.16	−0.39	13.57
γ 2 =10		10.0089	0.0680	0.08	−0.47	3.28
γ 3 =10		9.9992	0.0651	0.13	0.09	8.05
ρ=−0.25		−0.2373	0.0902	−0.04	−0.75	11.22
=1		1.0167	0.0953	0.73	1.76	11.15
=3		2.9282	0.2878	0.53	−0.18	4.08
Notes: T=time periods, N=sample size, K=number of replications, W is a standardized N×N contiguity matrix.
H 0 denotes time 0 distribution of exogenous variables H.
U[0,1] denotes a uniform (rectangular) distribution with minimum equal to 0 and maximum equal to 1.
Skewness is calculated as Σ(x i – m)^3/(n – 1)s ^3.
Kurtosis is calculated as Σ{(x i – m)^4/(n – 1)s ^4} – 3 in which m=Σx i /n.
The goodness of fit to the normal distribution is indicated by the residual deviance which has an asymptotic chi-squared distribution with the specified degrees of freedom. The table is formed by dividing the data into groups of approximately equal observed frequency. The degrees of freedom are c – p – 1, where c is the number of cells in the table of fitted values and p is the number of parameters (2) estimated in the model. Here there are 7 df.

Table 2 summarizes the estimate distributions obtained with T=4 time periods and assuming a different set of parameter values. In this case the true values are ρ=−0.5, λ=0.25, γ ₀=1, γ ₁=2, γ ₂=4, γ ₃=6, and . Once again the distribution means are all quite close to expectation, in most cases with low levels of skewness and kurtosis and acceptable approximations to normality, although the estimates are clearly furthest from normality.

Table 2 . Estimated parameter distributions

T=4		∼N(0,) Mean	∼N(0,) SD	W ^E =W (Rook's case) Skewness	H 0∼U[0,1] Kurtosis	K=100 Normality
λ=0.25		0.2481	0.0291	−0.37	−0.15	9.03
γ 0=1		1.0254	0.2868	0.35	−0.37	10.68
γ 1=2		1.9996	0.0250	−0.36	0.11	10.05
γ 2=4		3.9950	0.0229	0.02	0.35	16.17
γ 3=6		6.0001	0.0212	0.08	0.61	8.69
ρ=−0.5		−0.4879	0.0590	0.06	−0.46	5.28
=1		1.0227	0.0678	0.16	0.64	11.54
=1.4		1.5026	0.2832	2.16	7.42	20.59

The more detailed Monte Carlo results given in Appendix B (Tables A2 to A8) use both Rook's (edge touching) and Queen's (edge and corner touching) definitions of contiguity on a lattice, and also a torus, with opposite sides of the lattice treated as being contiguous so as to eliminate edges. Measures of bias and of a single indicator combining both precision (variance) and accuracy (bias), as given by a variant of the RMSE statistic (see Appendix B), are calculated from 1,000 Monte Carlo replications of equation (53). Summarizing the outcomes obtained under the various alternative assumptions detailed in Tables A2 to A8, it is evident that there is a small sample bias in . Attention is focused on positive dependence (negative ρ), which gives positive bias, and these outcomes are mirrored in the case of negative dependence,11 which gives negative bias, so that in both cases the estimated parameter is closer to zero than the true value. The bias is increasing in , but the most significant result is the clear evidence that as the sample size (N) increases, the bias in diminishes and the RMSE falls, suggesting consistency.

These indications of consistency are precisely what one might anticipate on the basis of the theoretical results given by Kapoor et al. (2007) (see also Kelejian & Prucha, 1998, 1999). An essential difference between their analysis and what is done here is, of course, that here we are assuming a spatial moving average process rather than spatially autoregressive errors. In addition, in this paper we also introduce an endogenous spatial lag in the panel context, a feature absent from their analysis. Consistency of the generalized moments estimators of ρ and σ ² is maintained by utilizing IV estimates of b leading to consistent disturbances. Thus, although the formal proofs given by Kapoor et al. (2007) are in the context of exogenous regressors (no spatial lag) and autoregressive rather than moving average errors, it is clear that their results carry through to the present set-up. Finally, it seems that although there is a small sample positive bias in the estimator , in many applied situations will be effectively unbiased. One advantage of GMM estimation is its comparative simplicity and computational efficaciousness12 in applications in which the number of locations is far in excess of those subject to Monte Carlo exploration in this paper. It is clear from the results presented here that as the number of locations rises into the thousands, as for example with the 3,000 plus counties in the USA, small sample bias in the estimator should be minimal.

6 Example 2: Real Estate Prices

In this example the GMM estimator is applied to a panel of average house prices in N=353 small areas13 of England in the T=2 years 2000, 2001, denoted by the NT×1 vector p. If the price at j is comparatively high, then demand may be displaced to nearby location k. On the other hand, supply may be displaced from k to j as investors in property seek higher returns. We therefore assume that price in area j interacts contemporaneously with price in area k, and model this interaction by the presence of an endogenous spatial lag Wp. In this case we again use the row normalized contiguity matrix W for both the endogenous spatial lag and the MA error process. The other explanatory variables14 are income from local jobs (wE), equal to the local wage rate (w) times the local employment level (E), and income from wages and employments within commuting distance (w ^c E ^c). In order to be able to treat these as exogenous variables, 1 year lags are introduced, so that year 2001 prices are a function of income in 2000, and year 2000 prices are a function of income in 1999.

There are many other variables that one might wish to introduce were panel data available, such as air quality, the quality of local schooling, the size of the existing stock of properties, demand coming from non-wage earners such as the retired and students, and the effects of criminality, social quality of the neighbourhood, amenity, local taxes, the nature of the housing stock, planning and building regulations, vulnerability to flooding and therefore the additional insurance premiums for areas on flood plains, and various other social, demographic, labour market, environmental and cultural differences. These omitted variables are likely to be spatially autocorrelated, the net effect of which is to induce an organized residual pattern (Dubin, 1988). We model these omitted variables by the spatial MA error process. While displaced demand or supply may causes price interactions that cascade outwards in an autoregressive process, I assume no such chain reaction for these variables, so that a shock, on its own, has a limited spatial extent, which is a property of the spatial MA process.

In order to obtain the estimates given in Table 3, the exogenous variables wE and w ^c E ^c and their first spatial lags, obtained by pre-multiplying these vectors by W, were used as instruments for the endogenous spatial lag in the first stage of the three-stage estimation process. This then provided estimated IV residuals which facilitated the second stage, enabling ρ, and estimates to be obtained. The third-stage estimates are given in Table 3, showing that there is a significant endogenous spatial lag effect, so that prices are directly positively related to contemporaneous prices in contiguous areas, and there are significant effects due to income from local jobs and jobs within commuting distance.

Table 3 . GMM estimates for the real estate price panel data with spatial moving average errors

	Parameter	Par. est.	SE	t-ratio
Constant	γ 0	44.4949 (43.6422)	12.3936 (12.5396)	3.59 (3.48)
Wp	λ	0.3941 (0.4011)	0.1174 (0.1196)	3.36 (3.35)
wE	γ 1	0.2221 (0.2494)	0.0976 (0.0978)	2.27 (2.55)
w ^c E ^c	γ 2	0.0507 (0.0493)	0.0066 (0.0067)	7.68 (7.34)
MA	ρ	−0.49279 (0.34998)
		42.9822 (37.4959)
		2722.83133 (2566.9143)
Note: Estimates using AR errors are given in parentheses.

7 Conclusion

This paper considers panel data in which spatial interaction comes from the effects of an endogenous lag and also from the MA error process. Monte Carlo results are given suggesting that the GMM estimator is consistent. It appears that this is the first paper to consider panel analysis with spatial MA errors, and also to jointly consider an endogenous lag together with spatial and temporal correlation in the error components, although much of this has been presaged in the earlier spatial econometrics literature (Anselin, 1988), and also in the time series context (Harvey, 1990). Indeed, in the conclusion to their paper, Kapoor et al. (2007) suggest that they would like to extend their results to models containing spatially lagged dependent variables. The present paper raises many issues which should be the subject of further study, such as the choice of appropriate instruments, the most efficient optimization method, and the small sample properties of the estimator, but the evidence presented here does suggest that there is scope for the practical implementation via GMM of panel data models with an endogenous spatial lag and spatial error processes.

Acknowledgements

The author is grateful to Michael Pfaffermayr and the other participants at the International Workshop on Spatial Econometrics and Statistics, Rome, 25–27 May 2006, and participants at the 13th International Conference on Panel Data, University of Cambridge, 7–9 July 2006, for their contributions to the discussion of this paper.

Notes

1. An early detailed account of the MA spatial process is given by Haining (1978).

2. Pre-multiplication of a TN×1 vector θ by Q ₀ creates a TN×1 vector of deviations from the mean, where the mean is obtained by averaging θ over time. Pre-multiplication of a TN×1 vector θ by Q ₁ creates a TN×1 vector, comprising N across time area-specific means stacked for each T.

3. Note that Tr(W′WW)=0 for the Rook's case contiguity matrix.

4. So that we can use equation (49) in both stage 1 and stage 3, it is assumed that and (so that Ω _ξ is a diagonal matrix of 1s) and that ρ=0. The result is that at stage 1 we simply obtain IV estimates.

5. Using unconstrained non-linear least squares estimation. The method is a modified Newton–Raphson method which is suitable for minimizing any non-linear function, and which depends on numerical differences rather than derivatives.

6. In contrast, at stage 1, ρ is assumed to equal 0, so that in that case Y ^* =Y, X ^* =X.

7. Moore–Penrose generalized inverses are used to avoid singularities.

8. In the main body of the text W is a Rook's case contiguity matrix based on a 15×15 lattice. In the Monte Carlo simulations described in Appendix B, the lattice size is varied and an irregular spatial partitioning is also considered. Also, alternative contiguity definitions, namely the Queen's case and torus, are also implemented.

9. Appendix B gives the results obtained using K=1,000 replications.

10. These procedures were carried out using the DISTRIBUTION directive of the programming language GENSTAT. The DISTRIBUTION directive is used to fit an observed sample of data to a theoretical distribution function, in order to obtain maximum-likelihood estimates of the parameters of the distribution and test the goodness of fit.

11. To save space these results have not been reported here.

12. With the MA error process the C-O transform involves the inverse. See Smirnov & Anselin (2001) for a discussion of the use of the power expansion to approximate the matrix inverse with large matrices.

13. Unitary authority and local authority districts, or UALADs.

14. Appendix A gives details of the sources and construction of these variables.

15. Small administrative areas, with median area equal to 250.77 km².

16. Available on the NOMIS website (the ONS online labour market statistics database).

17. 1991 Census of Population—Special Workplace Statistics, available from NOMIS.

18. Total employees and self-employed with a workplace coded, tabulated by residents in each zone (10% sample).

19. Minimum of the sum of the squared deviations of the observed proportions in each distance band up to 40 km and the proportions of the sum of the function exp(−δ _i d _ij ) calculated using the upper limit of each distance band.

Appendix A

The dependent variable p is the mean transaction price (all types of residential property) by area for the period July–September 2000, and July–September 2001 for the N=353 English Unitary Authority and Local Authority Districts15 (UALADs). The data were provided by the Land Registry. The wage rate (w) is the gross weekly pay for all occupations and both males and females taken from the Office for National Statistics’16 (ONS) New Earnings Survey. The employment level for the years 1999 and 2000 is based on the annual business enquiry employee analysis, also carried out by the ONS and available on the NOMIS database.

Total earnings in an area is the product of the average wage rate (w) in 1999 and 2000 and the total level of employment in 1999 and 2000, denoted by wE. These are assumed to be predetermined with respect to year 2000 and 2001 price levels.The vector w ^c E ^c denotes total earnings within commuting distance of a UALAD. This is equal to the matrix product of the n×n matrix C and the n×1 vector wE. Matrix C is defined as follows:

Cell (i, j) of the C matrix is a function of the (straight line) distance (d _ij ) between areas i and j and an area-specific coefficient δ _i . This allows for the different levels of transport infrastructure and commuting in different areas, with the choice of exponent δ _i based on empirical comparisons with observed census data17 on travel-to-work patterns. Table A1 shows the overall proportion of workers18 living in England and Wales travelling various distances from home to work. Given observed travel percentages comparable to Table A1 for each area, the exponent δ _i for each area was chosen by iterating the function exp(–δ _i d _ij ) through a range of values to obtain the value giving the closest fit19 to each area's commuting data.

Table A1 . Commuting distances to work in England and Wales

	Distance (km)
	<2	2–4	5–9	10–19	20–29	30–39	>40
%	26.63	25.28	20.93	15.90	4.96	2.05	4.25	100

Appendix B: Appendix B: Monte Carlo Investigation

Bias=median – true parameter value

where IQ is the interquantile range, equal to the difference between the 0.75 and 0.25 quantile. While this approximation is based on IQ rather than the variance, under normality the median is equal to the mean and, apart from slight rounding, IQ/1.35 is the standard deviation, so this measure reduces to the standard RMSE statistic (see Kapoor et al., 2007).

In all cases W is normalized to row totals equal to 1 and the bias is based on 1,000 Monte Carlo replications.

Table A2 . Bias and RMSE with increasing lattice size, positive dependence, Rook's case. ρ=−0.25, λ=0.75, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10, , . Matrix W is a Rook's case contiguity matrix. There are two layers (T=2)

Bias
T=2				W ^E=W Rook's case	H∼U[0,10]	K=1,000
Lattice size	5×5	7×7	9×9	11×11	13×13	15×15
γ 0	0.31131	0.25106	0.10097	0.007300	–0.09312	0.010786
λ	–0.00031	–0.00048	–0.00006	0.000041	0.00012	–0.000092
γ 1	0.00009	–0.00078	–0.00121	0.000398	0.00051	0.000472
γ 2	–0.00385	0.00209	–0.00120	–0.000090	0.00043	0.001027
γ 3	0.00647	–0.00226	–0.00037	0.002743	0.00251	–0.000860
ρ	0.09523	0.07720	0.05598	0.039236	0.02294	0.012190
RMSE
T=2				H∼U[0,10]	K=1,000	W ^E=W Rook's case
Lattice size	5×5	7×7	9×9	11×11	13×13	15×15
γ 0	5.489	3.035	2.519	1.265	1.665	1.368
λ	0.008289	0.004939	0.004075	0.002079	0.002743	0.002337
γ 1	0.07910	0.06366	0.04657	0.03653	0.02868	0.02700
γ 2	0.06890	0.06400	0.05137	0.03187	0.02646	0.02800
γ 3	0.09746	0.05135	0.04204	0.03859	0.03179	0.02603
ρ	0.2180	0.1744	0.1391	0.1192	0.1029	0.08351

Table A3 . Bias and RMSE with increasing lattice size, with negative dependence, Rook's case. ρ=0.5, λ=0.75, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10 , . Matrix W is a Rook's case contiguity matrix. There are two layers (T=2)

Bias
T=2			W ^E=W Rook's case	H∼U[0,10]	K=1,000
Lattice size	5×5	7×7	9×9	11×11	13×13
γ 0	−0.05185	0.09947	−0.004240	−0.007804	0.04937
λ	0.00007	−0.00014	0.000012	−0.000016	−0.00009
γ 1	0.00082	0.00102	0.000088	−0.000737	0.00235
γ 2	−0.00586	−0.00128	−0.000204	−0.001198	0.00009
γ 3	−0.01132	−0.00224	−0.001816	−0.001393	0.00060
ρ	−0.04482	−0.05278	−0.019962	−0.022022	−0.01821
RMSE
			W ^E=W Rook's case	H∼U[0,10]	K=1,000
Lattice size	5×5	7×7	9×9	11×11	13×13
γ 0	2.644	2.202	0.8189	0.4571	0.8922
λ	0.005466	0.004334	0.001814	0.000999	0.001852
γ 1	0.09958	0.05722	0.05082	0.03614	0.03045
γ 2	0.08222	0.06473	0.04433	0.03452	0.02953
γ 3	0.1147	0.05843	0.04015	0.03460	0.03073
ρ	0.2690	0.1859	0.1483	0.1337	0.1168

Table A4 . Bias with increases in the time dimension, positive dependence, Rook's case. ρ=−0.25, λ=0.75, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10, , . Matrix W is a Rook's case contiguity matrix for a 7×7 square, hence the number of locations is N=49, and W is of dimension 49×49. The number of layers is from 3 up to 10

			K=1,000	W ^E=W Rook's case
T	3	4	5	10
γ 0	−0.11550	−0.03176	−0.15351	−0.04148
λ	0.00017	0.00006	0.00029	0.00013
γ 1	0.00238	0.00174	0.00289	−0.00148
γ 2	−0.00289	0.00141	0.00103	0.00096
γ 3	0.00143	−0.00157	0.00088	0.00352
ρ	0.05291	0.05349	0.05032	0.07608

Table A5 . Bias with increasing positive dependence, Rook's case. λ=0.75, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10, , , T=2. Matrix W is a Rook's case contiguity matrix for a 7×7 square, hence the number of locations is N=49, and W is of dimension 49×49. There are two layers (T=2)

	ρ=0	ρ=−0.25	ρ=−0.5	ρ=−0.75	ρ=−0.95
γ 0	−0.005857	−0.08888	0.02602	0.06492	0.01550
λ	0.000164	0.00010	0.00010	0.00003	0.00000
γ 1	0.000023	0.00168	0.00042	−0.00056	0.00125
γ 2	0.000085	−0.00157	−0.00111	−0.00005	−0.00230
γ 3	−0.005534	−0.00041	−0.00445	−0.00084	−0.00028
ρ	0.044455	0.06901	0.12872	0.19019	0.25428

Table A6 . Bias with increasing positive dependence, Queen's case. λ=0.75, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10, , , T=2. Matrix W is a Queen's case contiguity matrix for a 9×9 square, hence the number of locations is N=81, and W is of dimension 81×81. There are two layers (T=2)

	ρ=0	ρ=−0.25	ρ=−0.5	ρ=−0.75	ρ=−0.95
γ 0	−0.13457	−0.10588	−0.17422	−0.17873	−0.26565
λ	0.00019	0.00025	0.00026	0.00025	0.00029
γ 1	−0.00075	−0.00077	−0.00063	−0.00051	−0.00073
γ 2	−0.00024	0.00001	0.00029	0.00056	0.00021
γ 3	−0.00053	0.00020	−0.00024	−0.00065	−0.00036
ρ	0.07173	0.10511	0.13368	0.17192	0.19795

Table A7 . Bias and RMSE with increasing lattice size, with positive dependence, Queen's case. ρ=−0.25, λ=0.75, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10, , . Matrix W is a Queen's case contiguity matrix. There are two layers (T=2)

Bias
T=2				W ^E=W Queen's case	H∼U[0,10]	K=1,000
Lattice size	5×5	7×7	9×9	11×11	13×13	15×15
γ 0	0.15583	0.40579	0.05725	0.00685	−0.07037	0.07598
λ	−0.00002	−0.00077	−0.00001	−0.00001	0.00012	−0.00018
γ 1	0.00077	−0.00088	−0.00048	−0.00058	0.00105	0.00057
γ 2	−0.00188	0.00266	−0.00157	0.00240	0.00067	0.00131
γ 3	0.00056	0.00049	−0.00163	0.00119	0.00364	−0.00119
ρ	0.27036	0.14180	0.06174	0.06889	0.04701	0.03344
RMSE
				H∼U[0,10]	K=1,000	W ^E=W Queen's case
Lattice size	5×5	7×7	9×9	11×11	13×13	15×15
γ 0	5.117	4.296	4.441	1.678	2.345	1.794
λ	0.008873	0.007006	0.007649	0.002720	0.003938	0.003026
γ 1	0.06928	0.06651	0.04646	0.03442	0.02895	0.02723
γ 2	0.1033	0.06679	0.04325	0.03177	0.02564	0.02739
γ 3	0.1116	0.05538	0.04431	0.03356	0.03217	0.02701
ρ	0.3914	0.2835	0.2160	0.1789	0.1529	0.1297

Table A8 . Bias with increasing lattice size, with positive dependence, torus. ρ=−0.75, λ=0.25, γ ₀=1, γ ₁=10, γ ₂=10, γ ₃=10, , . Matrix W is a torus, in other words a Rook's case contiguity matrix with no edges. There are three layers (T=3)

T=3		H∼U[0,10]	K=1,000	W ^E=W Rook's case
Lattice size	5×5	7×7	9×9	11×11
γ 0	−0.01686	−0.04101	−0.02184	−0.01271
λ	−0.00008	0.00018	0.00022	−0.00020
γ 1	−0.00450	0.00127	−0.00189	0.00103
γ 2	−0.00061	0.00089	−0.00009	−0.00149
γ 3	−0.00013	0.00046	0.00175	−0.00044
ρ	0.20835	0.14658	0.09758	0.08673