Type:	Package
Title:	Estimation of Spatial Weight Matrices
Version:	0.1.0
Author:	Tamas Krisztin [aut, cre], Philipp Piribauer [aut]
Maintainer:	Tamas Krisztin <krisztin@iiasa.ac.at>
Description:	Bayesian estimation of spatial weight matrices in spatial econometric panel models. Allows for estimation of spatial autoregressive (SAR), spatial error (SEM), spatial Durbin (SDM), spatial error Durbin (SDEM) and spatially lagged explanatory variable (SLX) type specifications featuring an unknown spatial weight matrix. Methodological details are given in Krisztin and Piribauer (2022) <doi:10.1080/17421772.2022.2095426>.
License:	GPL (≥ 3)
Depends:	R (≥ 3.5.0)
Imports:	graphics, Matrix, matrixcalc, plot.matrix, R6, stats, utils
Encoding:	UTF-8
Language:	EN-US
LazyData:	true
RoxygenNote:	7.3.2
NeedsCompilation:	no
Packaged:	2025-05-13 17:22:31 UTC; tamaskrisztin
Repository:	CRAN
Date/Publication:	2025-05-13 18:10:06 UTC

Set prior specifications for the spatial weight matrix

Description

Set prior specifications for the n by n spatial weight matrix W=f(\Omega), where \Omega is an n by n unknown binary adjacency matrix (with zeros on the main diagonal), and f() denotes the (optional) row-standardization function

Usage

W_priors(
  n,
  W_prior = matrix(0.5, n, n),
  symmetric_prior = FALSE,
  row_standardized_prior = TRUE,
  nr_neighbors_prior = bbinompdf(0:(n - 1), nsize = n - 1, a = 1, b = 1, min_k = 0, max_k
    = n - 1)
)

Arguments

An R6 class for sampling the elements of W

Description

An R6 class for sampling the elements of W

An R6 class for sampling the elements of W

Format

An R6Class generator object

Details

This class samples the spatial weight matrix. Use the function W_priors class for setup.

The sampling procedure relies on conditional Bernoulli posteriors outlined in Krisztin and Piribauer (2022).

Public fields

Methods

Public methods

• W_sampler$new()

• W_sampler$set_rho()

• W_sampler$sample()

Method new()

Usage

W_sampler$new(W_prior, curr_rho = NULL, spatial_error = FALSE)

Arguments

Method set_rho()

If the spatial autoregressive parameter \rho is updated during the sampling procedure the log determinant, the spatial projection matrix I - \rho W and it's inverse must be updated. This function should be used for a consistent update. At least the new scalar value for \rho must be supplied.

Usage

W_sampler$set_rho(new_rho, newLogdet = NULL, newA = NULL, newAI = NULL)

Arguments

Method sample()

Usage

W_sampler$sample(Y, curr_sigma, mu, lag_mu = matrix(0, nrow(tY), ncol(tY)))

Arguments

References

Krisztin, T., and Piribauer, P. (2022) A Bayesian approach for the estimation of weight matrices in spatial autoregressive models. Spatial Economic Analysis, 1-20.

Probability density for a hierarchical prior setup for the elements of the adjacency matrix based on the beta binomial distribution

Description

A hierarchical prior setup can be used in W_priors to anchor the prior number of expected neighbors. Assuming a fixed prior inclusion probability \underline{p}=1/2 for the off-diagonal entries in the binary n by n adjacency matrix \Omega implies that the number of neighbors (i.e. the row sums of \Omega) follows a Binomial distribution with a prior expected number of neighbors for the n spatial observations of (n-1)\underline{p}. However, such a prior structure has the potential undesirable effect of promoting a relatively large number of neighbors. To put more prior weight on parsimonious neighborhood structures and promote sparsity in \Omega, the beta binomial prior accounts for the number of neighbors in each row of \Omega.

Usage

bbinompdf(x, nsize, a, b, min_k = 0, max_k = nsize)

Arguments

Details

The beta-binomial distribution is the result of treating the prior inclusion probability \underline{p} as random (rather than being fixed) by placing a hierarchical beta prior on it. For the number of neighbors x, the resulting prior on the elements of \Omega, \omega_{ij}, can be written as:

p(\omega_{ij} = 1 | x)\propto \Gamma\left(a+ x \right)\Gamma\left(b+(n-1)-x\right),

where \Gamma(\cdot ) is the Gamma function, and a and b are hyperparameters from the beta prior. In the case of a = b = 1, the prior takes the form of a discrete uniform distribution over the number of neighbors. By fixing a = 1 the prior can be anchored around the expected number of neighbors m through b=[(n-1)-m]/m (see Ley and Steel, 2009).

The prior can be truncated by setting a minimum (min_k) and/or a maximum number of neighbors (max_k). Values outside this range have zero prior support.

Value

Prior density evaluated at x.

References

Ley, E., & Steel, M. F. (2009). On the effect of prior assumptions in Bayesian model averaging with applications to growth regression. Journal of Applied Econometrics, 24(4). doi:10.1002/jae.1057.

Set prior specifications for the slope parameters

Description

This function allows the user to specify custom values for Gaussian priors on the slope parameters.

Usage

beta_priors(
  k,
  beta_mean_prior = matrix(0, k, 1),
  beta_var_prior = diag(k) * 100
)

Arguments

Details

For the slope parameters \beta the package uses common Normal prior specifications. Specifically, p(\beta)\sim\mathcal{N}(\underline{\mu}_\beta,\underline{V}_\beta).

This function allows the user to specify custom values for the prior hyperparameters \underline{\mu}_\beta and \underline{V}_\beta. The default values correspond to weakly informative Gaussian priors with mean zero and a diagonal prior variance-covariance matrix with 100 on the main diagonal.

Value

A list with the prior mean vector (beta_mean_prior), the prior variance matrix (beta_var_prior) and the inverse of the prior variance matrix (beta_var_prior_inv).

An R6 class for sampling slope parameters

Description

This class samples slope parameters with a Gaussian prior from the conditional posterior. Use the beta_priors class for setup.

Format

An R6Class generator object

Public fields

Methods

Public methods

• beta_sampler$new()

• beta_sampler$sample()

Method new()

Usage

beta_sampler$new(beta_prior)

Arguments

Method sample()

Usage

beta_sampler$sample(Y, X, curr_sigma)

Arguments

The four-parameter Beta probability density function

Description

A four-parameter Beta specification as the prior for the spatial autoregressive parameter \rho, as proposed by LeSage and Parent (2007) .

Usage

betapdf(rho, a = 1, b = 1, rmin = 0, rmax = 1)

Arguments

Details

The prior density is given by:

p(\rho) \sim \frac{1}{Beta(a,b)} \frac{(\rho - \underline{\rho}_{min})^{(a-1)} (\underline{\rho}_{max} - \rho)^{(b-1)} }{2^{a + b - 1}}

where Beta(a, b) (a,b > 0) represents the Beta function, Beta(a, b)= \int_{0}^{1} t^{a-1} (1-t)^{b-1} dt.

Value

Density value evaluated at rho.

References

LeSage, J. P., and Parent, O. (2007) Bayesian model averaging for spatial econometric models. Geographical Analysis, 39(3), 241-267.

Covid incidences data

Description

COVID-19 data set provided by Johns Hopkins University (Dong et al., 2020). The database contains information on (official) daily infections for a large panel of countries around the globe in the very beginning of the outbreak from 17 February to 20 April 2020.

Usage

covid

Format

A data.frame object.

Details

Data is provided for countries: Australia (AUS), Bahrain (BHR), Belgium (BEL), Canada (CAN), China (CHN), Finland (FIN), France (FRA), Germany (DEU), Iran (IRN), Iraq (IRQ), Israel (ISR), Italy (ITA), Japan (JPN), Kuwait (KWT), Lebanon (LBN), Malaysia (MYS), Oman (OMN), Republic of Korea (KOR), Russian Federation (RUS), Singapore (SGP), Spain (ESP), Sweden (SWE), Thailand (THA), United Arab Emirates (ARE), United Kingdom (GBR), United States of America (USA), and Viet Nam (VNM).

The dataset includes daily data on the country specific maximum measured temperature (Temperature) and precipitation levels (Precipitation) as additional covariates (source: Dark Sky API). The stringency index (Stringency) put forward by Hale et al. (2020), which summarizes country-specific governmental policy measures to contain the spread of the virus. We use the biweekly average of the reported stringency index.

References

Dong, E., Du, H., and Gardner, L. (2020). An interactive web-based dashboard to track COVID-19 in real time. The Lancet Infectious Diseases, 20(5), 533–534. doi:10.1016/S1473-3099(20)30120-1.

Hale, T., Petherick, A., Phillips, T., and Webster, S. (2020). Variation in government responses to COVID-19. Blavatnik School of Government Working Paper, 31, 2020–2011. doi:10.1038/s41562-021-01079-8.

Krisztin, T., and Piribauer, P. (2022). A Bayesian approach for the estimation of weight matrices in spatial autoregressive models, Spatial Economic Analysis, 1-20. doi:10.1080/17421772.2022.2095426.

Krisztin, T., Piribauer, P., and Wögerer, M. (2020). The spatial econometrics of the coronavirus pandemic. Letters in Spatial and Resource Sciences, 13 (3), 209-218. doi:10.1007/s12076-020-00254-1.

Dong, E., Du, H., and Gardner, L. (2020). An interactive web-based dashboard to track COVID-19 in real time. The Lancet Infectious Diseases, 20(5), 533–534. doi:10.1016/S1473-3099(20)30120-1.

Efficient update of the log-determinant and the matrix inverse

Description

While updating the elements of the spatial weight matrix in SAR and SDM type spatial models in a MCMC sampler, the log-determinant term has to be regularly updated, too. When the binary elements of the adjacency matrix are treated unknown, the Matrix Determinant Lemma and the Sherman-Morrison formula are used for computationally efficient updates.

Usage

logdetAinvUpdate(ch_ind, diff, AI, logdet)

Arguments

Details

Let A = (I_n - \rho W) be an invertible n by n matrix. v is an n by 1 column vector of real numbers and u is a binary vector containing a single one and zeros otherwise. Then the Matrix Determinant Lemma states that:

A + uv' = (1 + v'A^{-1}u)det(A)

.

This provides an update to the determinant, but the inverse of A has to be updated as well. The Sherman-Morrison formula proves useful:

(A + uv')^{-1} = A^{-1} \frac{A^{-1}uv'A^{-1}}{1 + v'A^{-1}u}

.

Using these two formulas, an efficient update of the spatial projection matrix determinant can be achieved.

Value

A list containing the updated n by n matrix A^{-1}, as well as the updated log determinant of A

References

Sherman, J., and Morrison, W. J. (1950) Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics, 21(1), 124-127.

Harville, D. A. (1998) Matrix algebra from a statistician's perspective. Taylor & Francis.

Pace and Barry's log determinant approximation

Description

Bayesian estimates of parameters of SAR and SDM type spatial models require the computation of the log-determinant of positive-definite spatial projection matrices of the form (I_n - \rho W), where W is a n by n spatial weight matrix. However, direct computation of the log-determinant is computationally expensive.

Usage

logdetPaceBarry(W, length.out = 200, rmin = -1, rmax = 1)

Arguments

Details

This function wraps the log-determinant approximation by Barry and Pace (1999), which can be used to pre-compute the log-determinants over a grid of \rho values.

Value

numeric length.out by 2 matrix; the first column contains the approximated log-determinants the second column the \rho values ranging between rmin and rmax.

References

Barry, R. P., and Pace, R. K. (1999) Monte Carlo estimates of the log determinant of large sparse matrices. Linear Algebra and its applications, 289(1-3), 41-54.

A Markov Chain Monte Carlo (MCMC) sampler for a linear panel model

Description

The sampler uses independent Normal-inverse-Gamma priors to estimate a linear panel data model. The function is used for an illustration on using the beta_sampler and sigma_sampler classes.

Usage

normalgamma(
  Y,
  tt,
  X = matrix(1, nrow(Y), 1),
  niter = 200,
  nretain = 100,
  beta_prior = beta_priors(k = ncol(X)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered model takes the form:

Y_t = X_t \beta + \varepsilon_t,

with \varepsilon_t \sim N(0,I_n \sigma^2).

Y_t (n \times 1) collects the n cross-sectional observations for time t=1,...,T. X_t (n \times k_1) is a matrix of explanatory variables. \beta (k_1 \times 1) is an unknown slope parameter matrix.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), X=[X_1',...,X_T']' (N \times k), with N=nT, the final model can be expressed as:

Y = X \beta + \varepsilon,

where \varepsilon \sim N(0,I_N \sigma^2). Note that the input data matrices have to be ordered first by the cross-sectional (spatial) units and then stacked by time.

Examples

n = 20; tt = 10; k = 3
X = matrix(stats::rnorm(n*tt*k),n*tt,k)
Y = X %*% c(1,0,-1) + stats::rnorm(n*tt,0,.5)
res = normalgamma(Y,tt,X)

Regional growth data

Description

Regional growth data set contains information on annual growth rates of GVA per worker (labor productivity) for 90 European NUTS-1 regions for the period 1999-2019.

Usage

nuts1growth

Format

A data.frame object.

Details

The dataset contains annual regional data for 26 European Union countries, disaggregated at the NUTS-1 level. The countries covered are: Austria, Belgium, Bulgaria, Cyprus, Czechia, Denmark, Estonia, Finland, France, Germany, Hungary, Ireland, Italy, Lithuania, Luxembourg, Latvia, Malta, Netherlands, Poland, Portugal, Romania, Sweden, Slovenia, and Slovakia. The NUTS-1 codes identify the first-level administrative regions within these countries.

The dataset includes the following variables: annual output per worker growth (in percent), the natural logarithm of initial gross value added (GVA) per worker, and the share of the population with low and high levels of education based on the International Standard Classification of Education (ISCED). This structure allows for analyzing regional economic performance in relation to educational attainment across the European Union.

References

Piribauer, P., Glocker, C., & Krisztin, T. (2023). Beyond distance: The spatial relationships of European regional economic growth. Journal of Economic Dynamics and Control, 155, 104735.

Graphical summary of the estimated adjacency matrix \Omega

Description

Graphical plot of the posterior probabilities of the estimated adjacency matrix \Omega.

Usage

## S3 method for class 'estimateW'
plot(
  x,
  cols = c("white", "lightgrey", "black"),
  breaks = c(0, 0.5, 0.75, 1),
  ...
)

Arguments

Graphical summary of a generated spatial weight matrix

Description

Graphical summary of a generated spatial weight matrix

Usage

## S3 method for class 'sim_dgp'
plot(x, ...)

Arguments

Specify prior for the spatial autoregressive parameter and sampling settings

Description

Specify prior for the spatial autoregressive parameter and sampling settings

Usage

rho_priors(
  rho_a_prior = 1,
  rho_b_prior = 1,
  rho_min = 0,
  rho_max = 1,
  init_rho_scale = 1,
  griddy_n = 60,
  use_griddy_gibbs = TRUE,
  mh_tune_low = 0.4,
  mh_tune_high = 0.6,
  mh_tune_scale = 0.1
)

Arguments

An R6 class for sampling the spatial autoregressive parameter \rho

Description

An R6 class for sampling the spatial autoregressive parameter \rho

An R6 class for sampling the spatial autoregressive parameter \rho

Format

An R6Class generator object

Details

This class samples the spatial autoregressive parameter using either a tuned random-walk Metropolis-Hastings or a griddy Gibbs step. Use the rho_priors class for setup.

For the griddy Gibbs algorithm see Ritter and Tanner (1992).

Public fields

Methods

Public methods

• rho_sampler$new()

• rho_sampler$stopMHtune()

• rho_sampler$setW()

• rho_sampler$sample()

• rho_sampler$sample_Griddy()

• rho_sampler$sample_MH()

Method new()

Usage

rho_sampler$new(rho_prior, W = NULL)

Arguments

Method stopMHtune()

Function to stop the tuning of the Metropolis-Hastings step. The tuning of the Metropolis-Hastings step is usually carried out until half of the burn-in phase. Call this function to turn it off.

Usage

rho_sampler$stopMHtune()

Method setW()

Usage

rho_sampler$setW(newW, newLogdet = NULL, newA = NULL, newAI = NULL)

Arguments

Method sample()

Usage

rho_sampler$sample(Y, mu, sigma)

Arguments

Method sample_Griddy()

Usage

rho_sampler$sample_Griddy(Y, mu, sigma)

Arguments

Method sample_MH()

Usage

rho_sampler$sample_MH(Y, mu, sigma)

Arguments

References

Ritter, C., and Tanner, M. A. (1992). Facilitating the Gibbs sampler: The Gibbs stopper and the griddy-Gibbs sampler. Journal of the American Statistical Association, 87(419), 861-868.

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial autoregressive model (SAR) with exogenous spatial weight matrix.

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter beta prior for the spatial autoregressive parameter \rho. The function is used as an illustration on using the beta_sampler, sigma_sampler, and rho_sampler classes.

Usage

sar(
  Y,
  tt,
  W,
  Z = matrix(1, nrow(Y), 1),
  niter = 200,
  nretain = 100,
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial autoregressive model (SAR) takes the form:

Y_t = \rho W Y_t + Z_t \beta + \varepsilon_t,

with \varepsilon_t \sim N(0,I_n \sigma^2). The row-stochastic n by n spatial weight matrix W is non-negative and has zeros on the main diagonal. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. Z_t (n \times k_3) is a matrix of explanatory variables. \beta (k_3 \times 1) is an unknown slope parameter matrix.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), Z=[Z_1',...,Z_T']' (N \times k_3), with N=nT, the final model can be expressed as:

Y = \rho \tilde{W}Y + Z \beta + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,I_N \sigma^2). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time. This is a wrapper function calling sdm with no spatially lagged dependent variables.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n =n, tt = tt, rho = .5, beta3 = c(1,.5), sigma2 = .5)
res = sar(Y = dgp_dat$Y,tt = tt, W = dgp_dat$W,
          Z = dgp_dat$Z,niter = 100,nretain = 50)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial autoregressive model (SAR) with unknown spatial weight matrix

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter beta prior for the spatial autoregressive parameter \rho. This is a wrapper function calling sdmw with no spatially lagged exogenous variables.

Usage

sarw(
  Y,
  tt,
  Z,
  niter = 100,
  nretain = 50,
  W_prior = W_priors(n = nrow(Y)/tt),
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial autoregressive model (SAR) with unknown (n by n) spatial weight matrix W takes the form:

Y_t = \rho W Y_t + Z \beta + \varepsilon_t,

with \varepsilon_t \sim N(0,I_n \sigma^2) and W = f(\Omega). The n by n matrix \Omega is an unknown binary adjacency matrix with zeros on the main diagonal and f(\cdot) is the (optional) row-standardization function. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. Z_t (n \times k_3) is a matrix of explanatory variables. \beta (k_3 \times 1) is an unknown slope parameter vector.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), and Z=[Z_1', ..., Z_T']' (N \times k_3), with N=nT. The final model can be expressed as:

Y = \rho \tilde{W}Y + Z \beta + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,I_N \sigma^2). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Estimation usually even works well in cases of n >> T. However, note that for applications with n > 200 the estimation process becomes computationally demanding and slow. Consider in this case reducing niter and nretain and carefully check whether the posterior chains have converged.

Value

List with posterior samples for the slope parameters, \rho, \sigma^2, W, and average direct, indirect, and total effects.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n =n, tt = tt, rho = .5, beta3 = c(.5,1),
            sigma2 = .05,n_neighbor = 3,intercept = TRUE)
res = sarw(Y = dgp_dat$Y,tt = tt,Z = dgp_dat$Z,niter = 20,nretain = 10)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial Durbin error model (SDEM) with exogenous spatial weight matrix.

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter prior for the spatial autoregressive parameter \rho. The function is used as an illustration on using the beta_sampler, sigma_sampler, and rho_sampler classes.

Usage

sdem(
  Y,
  tt,
  W,
  X = matrix(0, nrow(Y), 0),
  Z = matrix(1, nrow(Y), 1),
  niter = 200,
  nretain = 100,
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(X) * 2 + ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial Durbin error model (SDEM) takes the form:

Y_t = X_t \beta_1 + W X_t \beta_2 + Z \beta_3 + \varepsilon_t,

with \varepsilon_t \sim N(0,(I_n - \rho W) \sigma^2). The row-stochastic n by n spatial weight matrix W is non-negative and has zeros on the main diagonal. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. X_t (n \times k_1) and Z_t (n \times k_2) are matrices of explanatory variables, where the former will also be spatially lagged. \beta_1 (k_1 \times 1), \beta_2 (k_1 \times 1) and \beta_3 (k_2 \times 1) are unknown slope parameter vectors.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), X=[X_1',...,X_T']' (N \times k_1) and Z=[Z_1', ..., Z_T']' (N \times k_2), with N=nT, the final model can be expressed as:

Y = X \beta_1 + \tilde{W} X \beta_2 + Z \beta_3 + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,\sigma^2 (I_n \otimes (I_n - \rho W) ). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n = n, tt = tt, rho = .5, beta1 = c(.5,1), beta2 = c(-1,.5),
                   beta3 = c(1.5), sigma2 = .5, spatial_error = TRUE)
res = sdem(Y = dgp_dat$Y, tt = tt,  W = dgp_dat$W, X = dgp_dat$X,
           Z = dgp_dat$Z, niter = 100, nretain = 50)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial Durbin error model (SDEM) with unknown spatial weight matrix

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter beta prior for the spatial autoregressive parameter \rho. It is a wrapper around W_sampler.

Usage

sdemw(
  Y,
  tt,
  X = matrix(0, nrow(Y), 0),
  Z = matrix(1, nrow(Y), 1),
  niter = 100,
  nretain = 50,
  W_prior = W_priors(n = nrow(Y)/tt),
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(X) * 2 + ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial Durbin error model (SDEM) with unknown (n by n) spatial weight matrix W takes the form:

Y_t = X_t \beta_1 + W X_t \beta_2 + Z \beta_3 + \varepsilon_t,

with \varepsilon_t \sim N(0,(I_n - \rho W) \sigma^2) and W = f(\Omega). The n by n matrix \Omega is an unknown binary adjacency matrix with zeros on the main diagonal and f(\cdot) is the (optional) row-standardization function. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. X_t (n \times k_1) and Z_t (n \times k_2) are matrices of explanatory variables, where the former will also be spatially lagged. \beta_1 (k_1 \times 1), \beta_2 (k_1 \times 1) and \beta_3 (k_2 \times 1) are unknown slope parameter vectors.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), X=[X_1',...,X_T']' (N \times k_1) and Z=[Z_1', ..., Z_T']' (N \times k_2), with N=nT. The final model can be expressed as:

Y = X \beta_1 + \tilde{W} X \beta_2 + Z \beta_3 + \varepsilon,

where \varepsilon \sim N(0,\sigma^2 (I_n \otimes (I_n - \rho W) ). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Estimation usually even works well in cases of n >> T. However, note that for applications with n > 200 the estimation process becomes computationally demanding and slow. Consider in this case reducing niter and nretain and carefully check whether the posterior chains have converged.

Value

List with posterior samples for the slope parameters, \rho, \sigma^2, W, and average direct, indirect, and total effects.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n =n, tt = tt, rho = .75, beta1 = c(.5,1),beta2 = c(-1,.5),
            beta3 = c(1.5), sigma2 = .05,n_neighbor = 3,intercept = TRUE, spatial_error = TRUE)
res = sdemw(Y = dgp_dat$Y,tt = tt,X = dgp_dat$X,Z = dgp_dat$Z,niter = 20,nretain = 10)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial Durbin model (SDM) with exogenous spatial weight matrix.

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter prior for the spatial autoregressive parameter \rho. The function is used as an illustration on using the beta_sampler, sigma_sampler, and rho_sampler classes.

Usage

sdm(
  Y,
  tt,
  W,
  X = matrix(0, nrow(Y), 0),
  Z = matrix(1, nrow(Y), 1),
  niter = 200,
  nretain = 100,
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(X) * 2 + ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial Durbin model (SDM) takes the form:

Y_t = \rho W Y_t + X_t \beta_1 + W X_t \beta_2 + Z \beta_3 + \varepsilon_t,

with \varepsilon_t \sim N(0,I_n \sigma^2). The row-stochastic n by n spatial weight matrix W is non-negative and has zeros on the main diagonal. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. X_t (n \times k_1) and Z_t (n \times k_2) are matrices of explanatory variables, where the former will also be spatially lagged. \beta_1 (k_1 \times 1), \beta_2 (k_1 \times 1) and \beta_3 (k_2 \times 1) are unknown slope parameter vectors.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), X=[X_1',...,X_T']' (N \times k_1) and Z=[Z_1', ..., Z_T']' (N \times k_2), with N=nT, the final model can be expressed as:

Y = \rho \tilde{W}Y + X \beta_1 + \tilde{W} X \beta_2 + Z \beta_3 + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,I_N \sigma^2). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n = n, tt = tt, rho = .5, beta1 = c(.5,1), beta2 = c(-1,.5),
                  beta3 = c(1.5), sigma2 = .5)
res = sdm(Y = dgp_dat$Y, tt = tt,  W = dgp_dat$W, X = dgp_dat$X,
          Z = dgp_dat$Z, niter = 100, nretain = 50)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial Durbin model (SDM) with unknown spatial weight matrix

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter beta prior for the spatial autoregressive parameter \rho. It is a wrapper around W_sampler.

Usage

sdmw(
  Y,
  tt,
  X = matrix(0, nrow(Y), 0),
  Z = matrix(1, nrow(Y), 1),
  niter = 100,
  nretain = 50,
  W_prior = W_priors(n = nrow(Y)/tt),
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(X) * 2 + ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial Durbin model (SDM) with unknown (n by n) spatial weight matrix W takes the form:

Y_t = \rho W Y_t + X_t \beta_1 + W X_t \beta_2 + Z \beta_3 + \varepsilon_t,

with \varepsilon_t \sim N(0,I_n \sigma^2) and W = f(\Omega). The n by n matrix \Omega is an unknown binary adjacency matrix with zeros on the main diagonal and f(\cdot) is the (optional) row-standardization function. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. X_t (n \times k_1) and Z_t (n \times k_2) are matrices of explanatory variables, where the former will also be spatially lagged. \beta_1 (k_1 \times 1), \beta_2 (k_1 \times 1) and \beta_3 (k_2 \times 1) are unknown slope parameter vectors.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), X=[X_1',...,X_T']' (N \times k_1) and Z=[Z_1', ..., Z_T']' (N \times k_2), with N=nT. The final model can be expressed as:

Y = \rho \tilde{W}Y + X \beta_1 + \tilde{W} X \beta_2 + Z \beta_3 + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,I_N \sigma^2). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Estimation usually even works well in cases of n >> T. However, note that for applications with n > 200 the estimation process becomes computationally demanding and slow. Consider in this case reducing niter and nretain and carefully check whether the posterior chains have converged.

Value

List with posterior samples for the slope parameters, \rho, \sigma^2, W, and average direct, indirect, and total effects.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n =n, tt = tt, rho = .75, beta1 = c(.5,1),beta2 = c(-1,.5),
            beta3 = c(1.5), sigma2 = .05,n_neighbor = 3,intercept = TRUE)
res = sdmw(Y = dgp_dat$Y,tt = tt,X = dgp_dat$X,Z = dgp_dat$Z,niter = 20,nretain = 10)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial error model (SEM) with exogenous spatial weight matrix.

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter prior for the spatial autoregressive parameter \rho. The function is used as an illustration on using the beta_sampler, sigma_sampler, and rho_sampler classes.

Usage

sem(
  Y,
  tt,
  W,
  Z = matrix(1, nrow(Y), 1),
  niter = 200,
  nretain = 100,
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial error model (SDEM) takes the form:

Y_t = Z \beta + \varepsilon_t,

with \varepsilon_t \sim N(0,(I_n - \rho W) \sigma^2). The row-stochastic n by n spatial weight matrix W is non-negative and has zeros on the main diagonal. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. Z_t (n \times k_2) is a matrix of explanatory variables, where the former will also be spatially lagged. \beta (k_3 \times 1) is an unknown slope parameter vector.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1) and Z=[Z_1', ..., Z_T']' (N \times k_2), with N=nT, the final model can be expressed as:

Y = Z \beta_3 + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,\sigma^2 (I_n \otimes (I_n - \rho W) ). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n =n, tt = tt, rho = .5, beta3 = c(.5,1),
            sigma2 = .05,n_neighbor = 3,intercept = TRUE,spatial_error = TRUE)
res = sem(Y = dgp_dat$Y, tt = tt,  W = dgp_dat$W,
           Z = dgp_dat$Z, niter = 100, nretain = 50)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial error model (SEM) with unknown spatial weight matrix

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters, as well as a four-parameter beta prior for the spatial autoregressive parameter \rho. This is a wrapper function calling sdemw with no spatially lagged exogenous variables.

Usage

semw(
  Y,
  tt,
  Z,
  niter = 100,
  nretain = 50,
  W_prior = W_priors(n = nrow(Y)/tt),
  rho_prior = rho_priors(),
  beta_prior = beta_priors(k = ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered panel spatial error model (SEM) with unknown (n by n) spatial weight matrix W takes the form:

Y_t = Z \beta + \varepsilon_t,

with \varepsilon_t \sim N(0,(I_n - \rho W) \sigma^2) and W = f(\Omega). The n by n matrix \Omega is an unknown binary adjacency matrix with zeros on the main diagonal and f(\cdot) is the (optional) row-standardization function. \rho is a scalar spatial autoregressive parameter.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. Z_t (n \times k_3) is a matrix of explanatory variables. \beta (k_3 \times 1) is an unknown slope parameter vector.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), and Z=[Z_1', ..., Z_T']' (N \times k_3), with N=nT. The final model can be expressed as:

Y = Z \beta + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,\sigma^2 (I_n \otimes (I_n - \rho W) ). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Estimation usually even works well in cases of n >> T. However, note that for applications with n > 200 the estimation process becomes computationally demanding and slow. Consider in this case reducing niter and nretain and carefully check whether the posterior chains have converged.

Value

List with posterior samples for the slope parameters, \rho, \sigma^2, W, and average direct, indirect, and total effects.

Examples

n = 20; tt = 10
dgp_dat = sim_dgp(n =n, tt = tt, rho = .5, beta3 = c(.5,1),
            sigma2 = .05,n_neighbor = 3,intercept = TRUE,spatial_error = TRUE)
res = semw(Y = dgp_dat$Y,tt = tt,Z = dgp_dat$Z,niter = 20,nretain = 10)

Set prior specification for the error variance using an inverse Gamma distribution

Description

Set prior specification for the error variance using an inverse Gamma distribution

Usage

sigma_priors(sigma_rate_prior = 0.001, sigma_shape_prior = 0.001)

Arguments

An R6 class for sampling for sampling \sigma^2

Description

This class samples nuisance parameter which an inverse Gamma prior from the conditional posterior. Use the sigma_priors class for setup.

Format

An R6Class generator object

Public fields

Methods

Public methods

• sigma_sampler$new()

• sigma_sampler$sample()

Method new()

Usage

sigma_sampler$new(sigma_prior)

Arguments

Method sample()

Usage

sigma_sampler$sample(Y, mu)

Arguments

Simulating from a data generating process

Description

This function can be used to generate data from a data generating process for SDM, SAR, SEM, and SLX type models.

Usage

sim_dgp(
  n,
  tt,
  rho,
  beta1 = c(),
  beta2 = c(),
  beta3 = c(),
  sigma2,
  n_neighbor = 4,
  W = NULL,
  do_symmetric = FALSE,
  intercept = FALSE,
  spatial_error = FALSE
)

Arguments

Details

For the SDM, SAR, and SLX models the generated spatial panel model takes the form

Y = \rho W Y + X \beta_1 + W X \beta_2 + Z \beta_3 + \epsilon,

with \epsilon \sim N(0,I_n\sigma^2).

For the SEM model the generated spatial panel model takes the form

Y = X \beta_1 + W X \beta_2 + Z \beta_3 + \epsilon,

with \epsilon \sim N(0,(I_n - \rho W)\sigma^2).

The function generates the N \times 1 vector Y. The elements of the explanatory variable matrices X (N \times k_1) and Z (N \times k_2) are randomly generated from a Gaussian distribution with zero mean and unity variance (N(0,1)).

The non-negative, row-stochastic n by n matrix W is constructed using a k-nearest neighbor specification based on a randomly generated spatial location pattern, with coordinates sampled from a standard normal distribution.

Values for the parameters \beta_1, \beta_2, and \beta_3, as well as \rho and \sigma^2 have to be provided by the user. The length of \beta_1 and \beta_2 have to be equal.

• A spatial Durbin model (SDM) is constructed if \rho is not equal to zero, spatial_error is FALSE, and \beta_1, \beta_2, and \beta_3 are all supplied by the user.

• A spatial autoregressive model is constructed if \rho is not equal to zero, spatial_error is FALSE, and only \beta_3 is supplied by the user.

• An SLX type model is constructed if \rho is equal to zero, spatial_error is FALSE, and \beta_1, \beta_2 are supplied by the user.

• An SEM type model is constructed if spatial_error is TRUE and either only \beta_3 or \beta_1, \beta_2, and \beta_3 are supplied by the user.

Value

A list with the generated X, Y and W and a list of parameters.

Examples

# SDM data generating process
dgp_dat = sim_dgp(n =20, tt = 10, rho = .5, beta1 = c(1,-1),
                  beta2 = c(0,.5),beta3 = c(.2),sigma2 = .5)

A Markov Chain Monte Carlo (MCMC) sampler for the panel spatial SLX model with unknown spatial weight matrix

Description

The sampler uses independent Normal-inverse-Gamma priors for the slope and variance parameters. It is a wrapper around W_sampler.

Usage

slxw(
  Y,
  tt,
  X = matrix(0, nrow(Y), 0),
  Z = matrix(1, nrow(Y), 1),
  niter = 100,
  nretain = 50,
  W_prior = W_priors(n = nrow(Y)/tt),
  beta_prior = beta_priors(k = ncol(X) * 2 + ncol(Z)),
  sigma_prior = sigma_priors()
)

Arguments

Details

The considered spatial panel SLX model with unknown (n by n) spatial weight matrix W takes the form:

Y_t = X_t \beta_1 + W X_t \beta_2 + Z \beta_3 + \varepsilon_t,

with \varepsilon_t \sim N(0,I_n \sigma^2) and W = f(\Omega). The n by n matrix \Omega is an unknown binary adjacency matrix with zeros on the main diagonal and f(\cdot) is the (optional) row-standardization function.

Y_t (n \times 1) collects the n cross-sectional (spatial) observations for time t=1,...,T. X_t (n \times k_1) and Z_t (n \times k_2) are matrices of explanatory variables, where the former will also be spatially lagged. \beta_1 (k_1 \times 1), \beta_2 (k_1 \times 1) and \beta_3 (k_2 \times 1) are unknown slope parameter vectors.

After vertically staking the T cross-sections Y=[Y_1',...,Y_T']' (N \times 1), X=[X_1',...,X_T']' (N \times k_1) and Z=[Z_1', ..., Z_T']' (N \times k_2), with N=nT. The final model can be expressed as:

Y = X \beta_1 + \tilde{W} X \beta_2 + Z \beta_3 + \varepsilon,

where \tilde{W}=I_T \otimes W and \varepsilon \sim N(0,I_N \sigma^2). Note that the input data matrices have to be ordered first by the cross-sectional spatial units and then stacked by time.

Estimation usually even works well in cases of n >> T. However, note that for applications with n > 200 the estimation process becomes computationally demanding and slow. Consider in this case reducing niter and nretain and carefully check whether the posterior chains have converged.

Value

List with posterior samples for the slope parameters, \sigma^2, W, and average direct, indirect, and total effects.

Examples

set.seed(123)
n = 20; tt = 10
dgp_dat = sim_dgp(n = 20, tt = 10, rho = 0, beta1 = c(1,-1),
                  beta2 = c(3,-2.5), beta3 = c(.2), sigma2 = .05,
                  n_neighbor = 3,intercept = TRUE)
res = slxw(Y = dgp_dat$Y, tt = tt, X = dgp_dat$X, Z = dgp_dat$Z,
                  niter = 20, nretain = 10)