O artigo de Harrinson e Rubinfield (1878) é um dos primeiros estudos aplicados em Economia Urbana a avaliar como a qualidade do ar pode influenciar a precificação no mercado imobiliário. Usando dados do Censo norte-americano e especificando um modelo hedônico para precificação de imóveis residenciais, os autores estimaram a disposição a pagar (WTP - willingness to pay) pelo ar puro na cidade de Boston, a partir de informações sobre o valor de mercado da moradia em 506 áreas censitárias em 1970. A qualidade do ar mensurada pela emissão de óxidos de nitrogênio (NOX - nitrogen oxides) foi aferida por 122 estações TASSIM (Transportation and Air Shed Simulation Model), das quais 96 estavam localizadas dentro da amostra. Esses valores de NOX foram então atribuídos proporcionalmente a todas áreas censitários pertencentes a amostra utilizada no modelo hedônico estimado por OLS (Ordinary Least Square) por Harrinson e Rubinfield (1878). Os autores verificaram que os danos marginais da poluição do ar aumentaram com o nível de poluição do ar e com a renda das famílias. Os resultados mostram ainda que a precificação no mercado habitacional é relativamente mais sensível a mudança, do que a demanda pela qualidade do ar.
Apesar do esforço, estudos recentes mostram que o impacto da poluição do ar sobre a disposição a pargar pode ser bastante variável quando o modelo hedônico incorpora o efeito da autocorrelação espacial (Gilley e Pace, 1996; Pace e Gilley, 1997). Além da forte influência desse efeito, o problema de especificação do modelo WTP também deve considerar a forma funcional e a natureza desconhecida da heterocedásticidade espacial. Motivado por essas questões, Bivand (2017) revisita o modelo hedônico proposto por Harrinson e Rubinfield (1878), com o propósito de identificar o modelo espacial que melhor represente a disposição a pagar pelo ar puro. Outra questão interessante apontada por Bivand (2017), diz respeito a configuração amostral. Segundo o autor, a forte autocorrelação espacial pode ser causada pelas características específicas das áreas censitárias ou foi introduzida pela mudança nas unidades observacionais utilizadas para as diferentes variáveis.
Usando os dados de Harrinson e Rubinfield (1878), neste exercício iremos explorar apenas a sensibilidade do modelo WTP em diferentes especificações espaciais. Com base no método de Bivand (2017), testamos o poder de ajuste dos modelos SLX, SEM e SDEM; para tanto, testes de especificação da família LM são computados.
Dado que o objetivo é explorar diferentes formas funcionais para a hipótese de WTP pelo ar puro, três modelos espaciais são estimados por Máxima Verossimilhança.
Estimação do Spatial Lag in X (SLX), especificado por: \[\begin{eqnarray} \tag{1} y & = & X \beta + W X \delta + \epsilon \\ \tag{2} \epsilon &\sim& N(0, \sigma^2 I). \end{eqnarray}\]
Em que, é um vetor \((N \times 1)\) da variável dependente, \(X\) é uma matriz \((N \times K)\) de variáveis explicativas exógenas, \(W\) é uma matriz \((N \times N)\). Nesse modelo, \(\beta\) é um vetor de coeficientes angulares responsável pelo efeito interno das variáveis exógenas, enquanto \(\beta\) mede o impacto dos spillovers espaciais. O termo de erro \(\epsilon\) é assumido homocedástico e não autocorrelacionado.
Estimação do Spatial Error Model (SEM), especificado por: \[\begin{eqnarray} \tag{3} y & = & X \beta + u \\ \tag{4} u &=& \lambda W u + \epsilon. \end{eqnarray}\]
No SEM, \(u\) é um vetor \((N \times 1)\) de distúrbios espacialmente autocorrelacionados e \(\lambda\) é o coeficiente autorregressivo erro espacial. Os termos \(y\), \(X\), \(W\) e \(\beta\) seguem as mesmas descrições anteriores.
O modelo Spatial Durbin Error Model (SDEM), é especificado por: \[\begin{eqnarray} \tag{5} y & = & X \beta + W X \delta + u \\ \tag{6} u &=& \lambda W u + \epsilon. \end{eqnarray}\]
No SLX, os spillovers espaciais são operados pelas variáveis explicativas exógenas do modelo WTP. No modelo SEM o processo autorregressivo espacial é operado unicamente pelos erros. O modelo SDEM, por fim, é uma combinação dos dois primeiros, ou seja, há spillovers operados pelas variáveis exógenas e um processo autorregressivo espacial nos erros.
A base de dados é a boston.c (Corrected Boston
Housing Data), acessível pelo pacote spData. A base
original Harrinson e Rubinfield (1878) além de sofrer pequenos ajustes,
foi expandida com a inserção da longitude e latitude das observações. A
variável dependente do modelo WTP é a mediana do valor do imóvel (USD
1000) por área censitária, variando de USD 5,00 a USD 50,00. A poluição
do ar é medida pela concentração de óxidos de nitrogênio (NOX -
parts per 10 million), também por área censitária. As variáveis
exógenas contidas em X são especificadas abaixo.
CRIM é o número de crime per capita;
ZN é a proporção de terrenos residenciais zoneados para lotes com mais de 25.000 pés quadrados por cidade (constante para todas as áreas de Boston);
INDUS é a proporção de hectares comerciais não varejistas por cidade (constante para todas as áreas de Boston);
CHAS é uma dummy com valor 1 se a área faz fronteira com o Rio Charles; 0 caso contrário;
RM é o número médio de quartos por residência;
AGE é um vetor numérico com a proporção de residência ocupada e construída pelo proprietário antes de 1940;
DIS são as distâncias ponderadas para cinco centros de emprego em Boston;
RAD é um índice de acessibilidade a rodovias radiais;
TAX é um imposto sobre propriedade;
PTRATIO é a proporção professor-aluno por área censitária;
B é uma taxa proporcional da concentração de negros por área censitária;
LSTAT é uma proporção de pobres.
Carregando os pacotes e a base de dados
library(spatialreg)
library(spData)
library(spdep)
data(boston)Os testes de especificação LM comparam a estrutura dos modelos espaciais, com o modelo clássico não espacial. Por isso, computamos primeiramente o estimador OLS.
fm <- log(MEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT)
res0 <- lm(formula = fm, data = boston.c)
summary(res0)##
## Call:
## lm(formula = fm, data = boston.c)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71176 -0.09169 -0.00566 0.09895 0.79780
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.558e+00 1.544e-01 29.512 < 2e-16 ***
## CRIM -1.186e-02 1.245e-03 -9.532 < 2e-16 ***
## ZN 8.016e-05 5.056e-04 0.159 0.874105
## INDUS 2.395e-04 2.364e-03 0.101 0.919318
## CHAS1 9.140e-02 3.320e-02 2.753 0.006129 **
## I(NOX^2) -6.380e-01 1.131e-01 -5.639 2.88e-08 ***
## I(RM^2) 6.328e-03 1.312e-03 4.823 1.89e-06 ***
## AGE 9.074e-05 5.263e-04 0.172 0.863179
## log(DIS) -1.913e-01 3.339e-02 -5.727 1.78e-08 ***
## log(RAD) 9.571e-02 1.913e-02 5.002 7.91e-07 ***
## TAX -4.203e-04 1.227e-04 -3.426 0.000664 ***
## PTRATIO -3.112e-02 5.013e-03 -6.208 1.14e-09 ***
## B 3.637e-04 1.031e-04 3.527 0.000460 ***
## log(LSTAT) -3.712e-01 2.501e-02 -14.841 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1825 on 492 degrees of freedom
## Multiple R-squared: 0.8059, Adjusted R-squared: 0.8008
## F-statistic: 157.1 on 13 and 492 DF, p-value: < 2.2e-16Com base nos resíduos OLS, computamos as estatísticas LM para análise da autocorrelação espacial. As estatísticas significantes a menos de 5% de probabilidade de erro mostram que tanto o modelo SEM (Spatial Error Model) quanto o SAR (Spatial Autorregressive) oferecem melhor ajuste que o modelo clássico. Os testes não esclarecem, entretanto, a natureza do processo autorregressivo.
resLMtest <- lm.LMtests(model = res0, listw = nb2listw(boston.soi), test = c("LMerr","LMlag","RLMerr","RLMlag"))
summary(resLMtest)## Lagrange multiplier diagnostics for spatial dependence
## data:
## model: lm(formula = fm, data = boston.c)
## weights: nb2listw(boston.soi)
##
## statistic parameter p.value
## LMerr 186.570 1 < 2.2e-16 ***
## LMlag 190.713 1 < 2.2e-16 ***
## RLMerr 37.606 1 8.658e-10 ***
## RLMlag 41.749 1 1.038e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Estimamos então os modelos SLX, SEM e SDEM por Máxima Verossimilhança.
## Estimating the SLX (Spatial Lag in X)
res1 <- lmSLX(formula = fm, data = boston.c, listw = nb2listw(boston.soi), weights = NULL, Durbin = TRUE)
## Estimating the SEM (Spatial Error Model)
res2 <- errorsarlm(formula = fm, data = boston.c, listw = nb2listw(boston.soi), weights = NULL, Durbin = FALSE, method = 'Matrix')
## Estimating the SDM (Spatial Durbin Model)
res3 <- errorsarlm(formula = fm, data = boston.c, listw = nb2listw(boston.soi), Durbin = TRUE, method = 'Matrix')O objetivo é verificar se diferença da taxa log-verossimilhança extraída de modelo de mesma família (nested, nesting) é significante a 5% de probabilidade de erro. Os resultados são demonstrados a seguir.
lr1 <- LR.Sarlm(res0,res1)
lr2 <- LR.Sarlm(res0,res2)
lr3 <- LR.Sarlm(res0,res3)
lr4 <- LR.Sarlm(res1,res3)
lr5 <- LR.Sarlm(res2,res3)
mres <- matrix(data = NA, nrow = 5, ncol = 5)
mres[,1] <- c(lr1$estimate[1],lr2$estimate[1],lr3$estimate[1],lr4$estimate[1],lr5$estimate[1])
mres[,2] <- c(lr1$estimate[2],lr2$estimate[2],lr3$estimate[2],lr4$estimate[2],lr5$estimate[2])
mres[,3] <- c(lr1$statistic,lr2$statistic,lr3$statistic,lr4$statistic,lr5$statistic)
mres[,4] <- c(lr1$p.value,lr2$p.value,lr3$p.value,lr4$p.value,lr5$p.value)
mres[,5] <- c(lr1$parameter,lr2$parameter,lr3$parameter,lr4$parameter,lr5$parameter)
colnames(mres) <- c("Log-like nested","Log-like nesting","LR","p-value","df")
rownames(mres) <- c("OLS:SLX","OLS:SEM","OLS:SDEM","SLX:SDEM","SEM:SDEM")
mres## Log-like nested Log-like nesting LR p-value df
## OLS:SLX 149.9548 199.3880 -98.86649 2.775558e-15 13
## OLS:SEM 149.9548 255.8946 -211.87968 0.000000e+00 1
## OLS:SDEM 149.9548 282.4843 -265.05904 0.000000e+00 14
## SLX:SDEM 199.3880 282.4843 -166.19256 0.000000e+00 1
## SEM:SDEM 255.8946 282.4843 -53.17937 8.409851e-07 13Na tabela acima, todas as taxas log-verossimilhança (terceira coluna) são negativas e estatisticamente diferentes de zero ao nível de significância de 1%. Portanto, os testes indicam haver diferença sistemática entre os quatro modelos estimados. Nesse sentido, o SEM oferece melhor ajuste que o OLS ou SLX. O SDEM, por sua vez, oferece melhor ajuste que o SLX ou SEM.
Apresentamos agora os coeficientes estimados para os modelos SEM e SDEM.
summary(res2)##
## Call:errorsarlm(formula = fm, data = boston.c, listw = nb2listw(boston.soi),
## weights = NULL, Durbin = FALSE, method = "Matrix")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6476342 -0.0676007 0.0011091 0.0776939 0.6491629
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.85706017 0.16083868 23.9809 < 2.2e-16
## CRIM -0.00545832 0.00097262 -5.6120 2.000e-08
## ZN 0.00049195 0.00051835 0.9491 0.3425906
## INDUS 0.00019244 0.00282240 0.0682 0.9456389
## CHAS1 -0.03303429 0.02836929 -1.1644 0.2442464
## I(NOX^2) -0.23369331 0.16219195 -1.4408 0.1496287
## I(RM^2) 0.00800078 0.00106472 7.5145 5.707e-14
## AGE -0.00090974 0.00050116 -1.8153 0.0694827
## log(DIS) -0.10889419 0.04783715 -2.2764 0.0228249
## log(RAD) 0.07025730 0.02108181 3.3326 0.0008604
## TAX -0.00049870 0.00012072 -4.1311 3.611e-05
## PTRATIO -0.01907770 0.00564160 -3.3816 0.0007206
## B 0.00057442 0.00011101 5.1744 2.286e-07
## log(LSTAT) -0.27212780 0.02323159 -11.7137 < 2.2e-16
##
## Lambda: 0.70175, LR test value: 211.88, p-value: < 2.22e-16
## Asymptotic standard error: 0.032698
## z-value: 21.461, p-value: < 2.22e-16
## Wald statistic: 460.59, p-value: < 2.22e-16
##
## Log likelihood: 255.8946 for error model
## ML residual variance (sigma squared): 0.018098, (sigma: 0.13453)
## Number of observations: 506
## Number of parameters estimated: 16
## AIC: NA (not available for weighted model), (AIC for lm: -269.91)summary(res3)##
## Call:errorsarlm(formula = fm, data = boston.c, listw = nb2listw(boston.soi),
## Durbin = TRUE, method = "Matrix")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5096008 -0.0676523 -0.0054992 0.0655289 0.7493297
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.1222e+00 3.0704e-01 13.4256 < 2.2e-16
## CRIM -7.6205e-03 1.0270e-03 -7.4199 1.172e-13
## ZN 6.2747e-04 5.0475e-04 1.2431 0.2138229
## INDUS -2.5789e-04 2.8905e-03 -0.0892 0.9289079
## CHAS1 2.9755e-03 2.9198e-02 0.1019 0.9188291
## I(NOX^2) -1.0080e-01 1.8046e-01 -0.5586 0.5764683
## I(RM^2) 9.1072e-03 1.1277e-03 8.0763 6.661e-16
## AGE -1.1267e-03 4.8906e-04 -2.3038 0.0212366
## log(DIS) -1.5780e-01 8.3553e-02 -1.8886 0.0589484
## log(RAD) 7.0701e-02 2.1313e-02 3.3172 0.0009091
## TAX -4.7725e-04 1.1803e-04 -4.0434 5.269e-05
## PTRATIO -1.6751e-02 5.5828e-03 -3.0004 0.0026962
## B 4.9917e-04 1.0717e-04 4.6579 3.194e-06
## log(LSTAT) -2.5705e-01 2.2908e-02 -11.2211 < 2.2e-16
## lag.CRIM -1.1448e-02 2.2982e-03 -4.9813 6.315e-07
## lag.ZN -6.4948e-04 9.3189e-04 -0.6970 0.4858316
## lag.INDUS 6.3762e-04 4.6117e-03 0.1383 0.8900344
## lag.CHAS1 1.2942e-01 5.1856e-02 2.4959 0.0125656
## lag.I(NOX^2) -6.3595e-01 2.5636e-01 -2.4807 0.0131132
## lag.I(RM^2) 4.8033e-03 2.0105e-03 2.3891 0.0168884
## lag.AGE -4.0924e-04 9.3495e-04 -0.4377 0.6615957
## lag.log(DIS) -1.0050e-01 9.6207e-02 -1.0446 0.2961921
## lag.log(RAD) 3.7316e-02 3.7376e-02 0.9984 0.3180892
## lag.TAX 1.3743e-04 2.1371e-04 0.6431 0.5201737
## lag.PTRATIO -1.4181e-02 9.3931e-03 -1.5097 0.1311204
## lag.B -8.0535e-05 1.8920e-04 -0.4256 0.6703633
## lag.log(LSTAT) 1.0969e-02 4.0707e-02 0.2695 0.7875715
##
## Lambda: 0.62859, LR test value: 166.19, p-value: < 2.22e-16
## Asymptotic standard error: 0.037651
## z-value: 16.695, p-value: < 2.22e-16
## Wald statistic: 278.73, p-value: < 2.22e-16
##
## Log likelihood: 282.4843 for error model
## ML residual variance (sigma squared): 0.016945, (sigma: 0.13017)
## Number of observations: 506
## Number of parameters estimated: 29
## AIC: NA (not available for weighted model), (AIC for lm: -342.78)O coeficiente autorregreesivo espacial do SEM, \(\lambda=0.70175\), se aproxima do modelo
não poderado estimado por Bivand (2017), \(\lambda=0.73273\), disposto na Tabela A.10.
Analogamente, o coeficiente do modelo SDEM \(\lambda=0.62859\), também se aproxima da
estimativa produzida pelo autor, \(\lambda=0.65797\). Ambos são significantes
a menos de 5% de probabilidade de erro. É importante reforçar que a base
boston.c não é exatamente a mesma usada por Bivand
(2017).
No que diz respeito ao efeito direto da poluição do ar, os resultados
revelam um coeficiente angular não significante associado à variável
I(NOX^2), o que corrobora com as estimativas de Bivand
(2017), organizadas na Tabela A.11. Entretanto, o efeito indireto da
poluição do ar é estatisticamente significante
lag.I(NOX^2). Há, portanto, evidências de que a
precificação dos imóveis pode ser bastante variável em relação à
poluição gerada dentro das áreas censitárias. Esse resultado não é
verdadeiro, entretando, em termos da análise de vizinhança. O
coeficiente \(\delta_{NOX}=-0.63595\)
estimado significante no SDEM mostra que a poluição do ar afeta
sobretudo o bem-estar dos vizinhos mais próximos.
Bivand, Roger. Revisiting the Boston data set - changing the units of observation affects estimated willingness to play for clean air. REGION, 4(1):109, 05 2017.
Harrison, David and Rubinfield, Daniel. Hedonic housing prices and the demand for clean air. Journal of Environmental Economics and Management, 5(1):81-102, 03 1978.
Gilley, Otis. and Pace, Kelley. On the Harrison and Rubinfeld data. Journal of Environmental Economics and Management 31: 403-405, 1996.
Pace, Kelley and Gilley, Otis. Using the spatial configuration of the data to improve estimation. Journal of the Real Estate Finance and Economics 14: 333-340, 1997.