Package 'estimateW'

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.

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).

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

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.

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.

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.

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 precompute 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.

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)

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

Description

Graphical summary of a generated spatial weight matrix

Usage

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

Arguments

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

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.

Description

This function is intended to be called from the R6 class W_sampler.

Usage

sampleW_fast(
  Y,
  curr_sigma,
  mu,
  lag_mu,
  W_prior,
  curr_W,
  curr_w,
  curr_A,
  curr_AI,
  curr_logdet,
  curr_rho,
  nr_neighbors_prior,
  symmetric,
  spatial_error,
  row_standardized
)

Arguments

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)

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)

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)

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)

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)

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)