Resolución Tarea 2

LAS RESPUESTAS SON REFERENCIALES. ES DECIR, SON LOS ELEMENTOS MÍNIMOS QUE DEBERÍAN CONTENER LAS RESPUESTAS

Incorporamos los paquetes necesarios para llevar a cabo el análisis. Asegúrese de que todos los paquetes estén instalados en R.

library(spatialreg)
library(spdep)
library(RColorBrewer)
library(leaflet)
library(dplyr)
library(ggplot2)
library(tmap)
library(tmaptools)
library(lmtest)
library(splm)
library(data.table)
library(ggfortify)
library(ade4)
library(adespatial)

Determinantes precios de vivienda

Exploración espacial precios de vivienda

Para realizar este ejercicio, cargamos el archivo shape que contiene los precios de vivienda. En este caso, las manzanas que contiene la base de datos no son contiguas (puede hacer zoom al mapa y ver que los polígonos están separados). Por tanto, para facilitar el análisis, usaré los centroides de las manzanas (el punto más central) para desarrollar el ejercicio. Esto también limita el uso de matrices de peso espacial, ya que ahora los datos son de punto (en realidad, dado que las manzanas están aisladas, deberíamos pensar en una base de datos de punto siempre) usando una matriz de distancia o una matriz de k-vecinos.

De esta forma, el archivo shape que contiene los precios sería el siguiente:

data.afta<-st_read("ANTOFAGASTA_CENTROIDES.shp")

## Reading layer `ANTOFAGASTA_CENTROIDES' from data source 
##   `/Users/yasnacortesgarriga/Library/CloudStorage/Dropbox/cursos/Cursos 2023/Econometría Espacial/Tareas Estudiantes/Tarea 2/ANTOFAGASTA_CENTROIDES.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3162 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -70.44201 ymin: -23.73881 xmax: -70.37203 ymax: -23.53374
## Geodetic CRS:  WGS 84

data<-data.frame(data.afta)
names(data)

##  [1] "COMUNA"     "MANZ_SII"   "NOM_COMUNA" "CMN_MZ"     "precio_viv"
##  [6] "unidad_viv" "ln_pr_viv"  "esc_priv"   "esc_pub"    "esc_parvu" 
## [11] "farmacias"  "buses"      "plazas"     "restau"     "salud"     
## [16] "super"      "x_coord"    "y_coord"    "geometry"

attach(data)

Construimos un mapa con los precios de vivienda (en logarítmo natural)

tm_shape(data.afta) +
  tm_dots(col="ln_pr_viv", style="quantile", palette="Reds", title = "Precios Vivienda (logaritmo)") + tm_layout(legend.outside.size = 0.2)

## Legend labels were too wide. The labels have been resized to 0.39, 0.39, 0.39, 0.39, 0.39. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Usamos el modo interactivo para corroborar la ubicación espacial de los datos:

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(data.afta) +
  tm_dots(col="ln_pr_viv", style="quantile", palette="Reds", title = "Precios Vivienda (logaritmo)") + tm_layout(legend.outside.size = 0.2)

También, usamos matrices de peso espacial para determinar si los precios de vivienda presentan autocorrelación espacial.

Matriz W basada en los 10 Vecinos más cercanos

coords<-cbind(data$x_coord, data$y_coord)
kn5<-nb2listw(knn2nb(knearneigh(coords,k=10)), style = "W")
kn5

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 3162 
## Number of nonzero links: 31620 
## Percentage nonzero weights: 0.3162555 
## Average number of links: 10 
## 3 disjoint connected subgraphs
## Non-symmetric neighbours list
## 
## Weights style: W 
## Weights constants summary:
##      n      nn   S0    S1       S2
## W 3162 9998244 3162 586.2 12841.94

Matriz W basada en Distancias (Distancia Euclideana)

dist.mat<-as.matrix(dist(coords, method="euclidean"))
dist.mat.inv<-1/dist.mat
diag(dist.mat.inv)<-0
dist.W<-mat2listw(dist.mat.inv, row.names = data.afta$MANZ_SII, style="W")

Moran Global

Estimamos un índice de moran para comprobar la existencia de dependencia espacial:

moran.viv.kn5<-moran.test(ln_pr_viv, kn5, alternative="two.sided"); moran.viv.kn5

## 
##  Moran I test under randomisation
## 
## data:  ln_pr_viv  
## weights: kn5    
## 
## Moran I statistic standard deviate = 97.121, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      7.414226e-01     -3.163556e-04      5.832783e-05

moran.viv.dist<-moran.test(ln_pr_viv, dist.W, alternative="two.sided"); moran.viv.dist

## 
##  Moran I test under randomisation
## 
## data:  ln_pr_viv  
## weights: dist.W    
## 
## Moran I statistic standard deviate = 305.81, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.661228e-01     -3.163556e-04      7.590972e-07

Moran local

Significant Moran Local: Precios de Vivienda

En términos generales, podemos indicar lo siguiente: - A través de los mapas, podemos dar cuenta de la existencia de diferencias espaciales en los precios de vivienda, propias de las dinámicas que ocurren en dicho mercado. Podemos ver que el sector sur de la ciudad cuenta con los precios más altos de vivienda, mientras que el sector norte muestra una mayor mixtura. - En adición, podemos ver una clara distinción entre los precios de vivienda ubicados desde la línea del tren hacia el mar y los precios de vivienda ubicados desde la línea del tren hacia el cerro. Los precios de vivienda tienden a ser más bajos desde la línea del tren hacia el cerro, concentrándose mayoritariamente en el sector norte. Los precios de vivienda desde la línea de tren hasta la costa tienden a ser más altos. - Usando el índice de moran global, podemos dar cuenta de la existencia de autocorrelación global de la variable. Considerando distintas tipologías de matrices de peso espacial podemos dar cuenta que el índice es significativo, por tanto, podemos rechazar la hipótesis nula de aleatorización. Las magnitudes del índice de moran son diferentes entre ambas matrices (esperable, sobretodo con una matriz de distancia donde la magnitud es menor). - Usando el índice de moran local, podemos observar que en el sector sur de la ciudad y el sector norte/costa se concentran los precios de vivienda altos que se encuentran rodeados de precios de vivienda altos. En contraste, los precios de vivienda bajos que también están rodeados de precios de vivienda bajos se encuentran localizados en el sector centro-norte hacia el cerro.

Amenidades consideradas en el análisis:

Las amenidades consideradas en el análisis son las siguientes:

Tabla 1: Descripción de variables bajo estudio

Variable	Descripción
esc_priv	distancia inversa a la escuela privada más cercana.
esc_pub	distancia inversa a la escuela pública más cercana.
esc_parvu	distancia inversa a la escuela parvularia (nivel pre-básico) más cercana.
farmacias	distancia inversa a la farmacia más cercana.
buses	distancia inversa a la parada de autobus más cercana.
plazas	distancia inversa a la área verde más cercana.
restau	distancia inversa al restaurant más cercano.
salud	distancia inversa al centro de salud más cercano.
super	distancia inversa al supermercado más cercano.

Usamos mapas para reflejar el comportamiento espacial de las amenidades:

tmap_mode("view")

## tmap mode set to interactive viewing

am1<-tm_shape(data.afta) +
  tm_dots(col="esc_priv", style="quantile", palette="Reds", title = "Acceso a Escuelas Privadas") + tm_layout(legend.outside.size = 0.2)

am2<-tm_shape(data.afta) +
  tm_dots(col="esc_pub", style="quantile", palette="Reds", title = "Acceso a Escuelas Públicas") + tm_layout(legend.outside.size = 0.2)

am3<-tm_shape(data.afta) +
  tm_dots(col="esc_parvu", style="quantile", palette="Reds", title = "Acceso a Escuelas Parvularias") + tm_layout(legend.outside.size = 0.2)

am4<-tm_shape(data.afta) +
  tm_dots(col="farmacias", style="quantile", palette="Reds", title = "Acceso a Farmacias") + tm_layout(legend.outside.size = 0.2)

am5<-tm_shape(data.afta) +
  tm_dots(col="buses", style="quantile", palette="Reds", title = "Acceso a Paradas de Autobuses") + tm_layout(legend.outside.size = 0.2)

am6<-tm_shape(data.afta) +
  tm_dots(col="plazas", style="quantile", palette="Reds", title = "Acceso a Plazas") + tm_layout(legend.outside.size = 0.2)

am7<-tm_shape(data.afta) +
  tm_dots(col="restau", style="quantile", palette="Reds", title = "Acceso a Restaurantes") + tm_layout(legend.outside.size = 0.2)

am8<-tm_shape(data.afta) +
  tm_dots(col="salud", style="quantile", palette="Reds", title = "Acceso a Centros de Salud") + tm_layout(legend.outside.size = 0.2)

am9<-tm_shape(data.afta) +
  tm_dots(col="super", style="quantile", palette="Reds", title = "Acceso a Supermercados") + tm_layout(legend.outside.size = 0.2)

tmap_arrange(am1, am2, am3, am4, am5)

tmap_arrange(am6, am7, am8, am9)

Estimación de un modelo hedónico para los precios de vivienda

Planteamos el siguiente modelo hedónico para los precios de vivienda de Antofagasta:

\[\ln(precios vivienda) = X\beta + \epsilon\] Donde la variable dependiente es el logaritmo de los precios de vivienda. La matriz X contiene todas las amenidades consideradas en el estudio. Cabe resaltar que este modelo no considera características intrínsecas de la vivienda propiamente tal. Por tanto, estas variables no observadas podrían estar en el error.

Estimamos el modelo propuesto vía OLS.

model<-ln_pr_viv ~ esc_pub+esc_priv+esc_parvu+farmacias+buses+plazas+restau+salud+super
VI_OLS<-lm(model, data=data)
summary(VI_OLS)

## 
## Call:
## lm(formula = model, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3188 -0.4959 -0.1211  0.4092  6.2834 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  16.77484    0.03068 546.681  < 2e-16 ***
## esc_pub     -15.04706    4.18650  -3.594  0.00033 ***
## esc_priv      2.40796    4.05141   0.594  0.55232    
## esc_parvu   -36.69778    3.08187 -11.908  < 2e-16 ***
## farmacias    48.24788    5.26994   9.155  < 2e-16 ***
## buses         7.54040    1.81776   4.148 3.44e-05 ***
## plazas        3.83297    0.88319   4.340 1.47e-05 ***
## restau        0.77003    0.49182   1.566  0.11752    
## salud       -46.02260    4.99091  -9.221  < 2e-16 ***
## super        -4.78792    1.56102  -3.067  0.00218 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7298 on 3152 degrees of freedom
## Multiple R-squared:  0.1213, Adjusted R-squared:  0.1188 
## F-statistic: 48.36 on 9 and 3152 DF,  p-value: < 2.2e-16

autoplot(VI_OLS)

Con esta información, comprobamos la existencia de autocorrelación espacial en los residuos de la regresión.

lm.morantest(VI_OLS, listw=kn5, alternative = "two.sided")

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = model, data = data)
## weights: kn5
## 
## Moran I statistic standard deviate = 86.589, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.655987543     -0.001497303      0.000057656

lm.morantest(VI_OLS, listw=dist.W, alternative = "two.sided")

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = model, data = data)
## weights: dist.W
## 
## Moran I statistic standard deviate = 259.49, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     2.170902e-01    -5.152640e-04     7.032476e-07

Ya que existe autocorrelación espacial en los residuos, podríamos decir que en el error podrían existir variables no observadas que podrían tener un comportamiento espacial. Es decir, podríamos incorporar el rezago espacial de la variable dependiente como variable explicativa en el modelo y/o también, incorporar una estructura espacial en los residuos para modelarlos considerando la dependencia espacial. Sin embargo, antes de tomar esa decisión podríamos realizar un contraste de hipótesis para determinar y/o orientar nuestra búsqueda del modelo más adecuado para nuestra estructura de datos:

Para reducir el tiempo de cálculo usaremos la matriz basada en 10 vecinos.

lag.mul<-lm.LMtests(VI_OLS, listw = kn5, test = "all"); lag.mul

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = model, data = data)
## weights: kn5
## 
## LMerr = 7339.5, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = model, data = data)
## weights: kn5
## 
## LMlag = 7440.2, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = model, data = data)
## weights: kn5
## 
## RLMerr = 216.49, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = model, data = data)
## weights: kn5
## 
## RLMlag = 317.17, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = model, data = data)
## weights: kn5
## 
## SARMA = 7656.7, df = 2, p-value < 2.2e-16

De acuerdo a estos resultados, podemos señalar lo siguiente:

• Tanto el modelo SAR y SEM, serían modelos adecuados para explicar los precios de vivienda en la ciudad de Antofagasta. De acuerdo con los contrastes convencionales, ambos son significativos. Debemos revisar los robustos.
• Los contrastres robustos sugieren que tanto el modelo SAR y SEM serían adecuados para explicar los precios de vivienda. En este caso, deberíamos escoger el contraste que tiene el valor más alto (máxima confianza), el cual corresponde al modelo SAR. Esto también se condice con la teoría.
• Un modelo de tipo SARMA (considerando un coeficiente autorregresivo en la variable dependiente rezagada espacialmente y un coeficiente autorregresivo en el error) podría ser una versión alternativa a estimar.

Por tanto, estimamos los diferentes modelos solicitados. Por motivos de reducción del tiempo de cálculo, usaremos la matriz de 10 vecinos más cercanos:

Matriz basada en 10 vecinos más cercanos

AFTA.ERROR<-errorsarlm(model, listw=kn5, data=data)
summary(AFTA.ERROR)

## 
## Call:errorsarlm(formula = model, data = data, listw = kn5)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.977798 -0.164223 -0.024955  0.133071  5.342297 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value Pr(>|z|)
## (Intercept) 16.781927   0.073937 226.9755  < 2e-16
## esc_pub      7.700564   3.556117   2.1654  0.03035
## esc_priv    -1.806902   3.457213  -0.5226  0.60122
## esc_parvu   -0.558294   2.635002  -0.2119  0.83220
## farmacias    9.963236   4.552375   2.1886  0.02863
## buses        0.162093   1.293954   0.1253  0.90031
## plazas       0.200974   0.501946   0.4004  0.68887
## restau      -0.140445   0.265927  -0.5281  0.59741
## salud       -6.619183   4.248806  -1.5579  0.11926
## super       -1.395191   0.945173  -1.4761  0.13991
## 
## Lambda: 0.8986, LR test value: 3390, p-value: < 2.22e-16
## Asymptotic standard error: 0.0088205
##     z-value: 101.88, p-value: < 2.22e-16
## Wald statistic: 10379, p-value: < 2.22e-16
## 
## Log likelihood: -1790.77 for error model
## ML residual variance (sigma squared): 0.15838, (sigma: 0.39797)
## Number of observations: 3162 
## Number of parameters estimated: 12 
## AIC: 3605.5, (AIC for lm: 6993.5)

AFTA.SAR<-lagsarlm(model, listw=kn5, data=data)
summary(AFTA.SAR)

## 
## Call:lagsarlm(formula = model, data = data, listw = kn5)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.968823 -0.162010 -0.021154  0.128850  5.383248 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  1.9161241  0.1588838 12.0599 < 2.2e-16
## esc_pub      1.7673117  2.2859625  0.7731  0.439454
## esc_priv    -0.5243095  2.2119331 -0.2370  0.812628
## esc_parvu   -4.9113378  1.6915855 -2.9034  0.003691
## farmacias    8.3440289  2.8991293  2.8781  0.004001
## buses        0.7573518  0.9927308  0.7629  0.445525
## plazas       0.5843402  0.4822268  1.2118  0.225607
## restau       0.0095239  0.2684885  0.0355  0.971703
## salud       -7.5277941  2.7407481 -2.7466  0.006021
## super       -1.1420667  0.8525112 -1.3397  0.180359
## 
## Rho: 0.88629, LR test value: 3405.1, p-value: < 2.22e-16
## Asymptotic standard error: 0.009408
##     z-value: 94.206, p-value: < 2.22e-16
## Wald statistic: 8874.7, p-value: < 2.22e-16
## 
## Log likelihood: -1783.207 for lag model
## ML residual variance (sigma squared): 0.15873, (sigma: 0.39841)
## Number of observations: 3162 
## Number of parameters estimated: 12 
## AIC: NA (not available for weighted model), (AIC for lm: 6993.5)
## LM test for residual autocorrelation
## test value: 18.306, p-value: 1.8809e-05

AFTA.SAC<-sacsarlm(model, listw=kn5, data=data, type = "SAC")
summary(AFTA.SAC)

## 
## Call:sacsarlm(formula = model, data = data, listw = kn5, type = "SAC")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.859223 -0.161533 -0.017439  0.129174  5.233426 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     1.40758    0.17769  7.9217 2.442e-15
## esc_pub        10.03004    3.60939  2.7789  0.005455
## esc_priv       -2.32246    3.50241 -0.6631  0.507265
## esc_parvu       2.24363    2.67727  0.8380  0.402013
## farmacias       7.68834    4.60543  1.6694  0.095036
## buses          -0.58017    1.29733 -0.4472  0.654727
## plazas          0.17376    0.50646  0.3431  0.731532
## restau         -0.10941    0.27021 -0.4049  0.685554
## salud          -4.14934    4.29866 -0.9653  0.334413
## super          -1.49994    0.95119 -1.5769  0.114819
## lag.esc_pub   -13.11071    4.81677 -2.7219  0.006491
## lag.esc_priv    1.77525    4.25490  0.4172  0.676514
## lag.esc_parvu -10.04603    3.55995 -2.8220  0.004773
## lag.farmacias  -3.86094    5.74501 -0.6721  0.501551
## lag.buses       2.73158    1.88697  1.4476  0.147728
## lag.plazas      1.06436    0.98706  1.0783  0.280897
## lag.restau      1.31554    0.69245  1.8998  0.057455
## lag.salud      -3.08466    5.38133 -0.5732  0.566498
## lag.super       1.62707    1.56109  1.0423  0.297289
## 
## Rho: 0.91672
## Asymptotic standard error: 0.010483
##     z-value: 87.448, p-value: < 2.22e-16
## Lambda: -0.31491
## Asymptotic standard error: 0.066098
##     z-value: -4.7643, p-value: 1.8949e-06
## 
## LR test value: 3450.7, p-value: < 2.22e-16
## 
## Log likelihood: -1760.375 for sacmixed model
## ML residual variance (sigma squared): 0.15244, (sigma: 0.39043)
## Number of observations: 3162 
## Number of parameters estimated: 22 
## AIC: NA (not available for weighted model), (AIC for lm: 6993.5)

AFTA.SDM<-lagsarlm(model, listw=kn5, data=data, type = "mixed")
summary(AFTA.SDM)

## 
## Call:lagsarlm(formula = model, data = data, listw = kn5, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.884179 -0.163774 -0.019904  0.129459  5.285630 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     2.125858   0.172598 12.3168 < 2.2e-16
## esc_pub         9.479239   3.657638  2.5916  0.009552
## esc_priv       -2.218250   3.486852 -0.6362  0.524662
## esc_parvu       2.023521   2.702859  0.7487  0.454062
## farmacias       8.472337   4.616521  1.8352  0.066473
## buses          -0.142611   1.303310 -0.1094  0.912868
## plazas          0.240925   0.505032  0.4770  0.633327
## restau         -0.047457   0.270602 -0.1754  0.860785
## salud          -5.119583   4.303051 -1.1898  0.234142
## super          -1.478883   0.947207 -1.5613  0.118451
## lag.esc_pub   -13.358506   5.257356 -2.5409  0.011056
## lag.esc_priv    1.363123   4.566273  0.2985  0.765306
## lag.esc_parvu -12.530939   3.840943 -3.2625  0.001104
## lag.farmacias  -2.549145   6.201940 -0.4110  0.681055
## lag.buses       2.943420   2.112470  1.3934  0.163513
## lag.plazas      1.652262   1.166460  1.4165  0.156637
## lag.restau      1.440600   0.809284  1.7801  0.075061
## lag.salud      -4.880944   5.795484 -0.8422  0.399677
## lag.super       1.216484   1.798979  0.6762  0.498909
## 
## Rho: 0.87403, LR test value: 3035.4, p-value: < 2.22e-16
## Asymptotic standard error: 0.010155
##     z-value: 86.072, p-value: < 2.22e-16
## Wald statistic: 7408.4, p-value: < 2.22e-16
## 
## Log likelihood: -1769.906 for mixed model
## ML residual variance (sigma squared): 0.15842, (sigma: 0.39802)
## Number of observations: 3162 
## Number of parameters estimated: 21 
## AIC: NA (not available for weighted model), (AIC for lm: 6615.3)
## LM test for residual autocorrelation
## test value: 16.99, p-value: 3.7576e-05

También, podemos indicar que debemos seguir buscando el mejor modelo porque aun persiste la autocorrelación espacial en los residuos de los modelos. Por ejemplo, los residuos del modelo SAR todavía conserva autocorrelación espacial, por tanto, aun existen variables omitidas del modelo y necesita un mayor tratamiento agregando variables que corresponden a las características propias de la vivienda y/u otras variables que pueden ser relevantes a la hora de estudiar los precios de vivienda en la ciudad de Antofasgasta.