Si bien, los coeficientes de los modelos espaciales estimados podrían ser obtenidos a través de OLS/MCO, la presencia de endogeneidad entre los errores y las variables rezagadas espacialmente, hace que el estimador de OLS/MCO pierda sus propiedades óptimas (BLUE). Para subsanar este problema, máxima verosimilitud se transforma en una alternativa para lidiar con la endogeneidad existente aprovechando sus propiedades asintóticas que aseguran que los estimadores de MV sean consistentes y distribuidos normales asintóticamente.

Sin embargo, en términos prácticos, la función de log-verosimilitud proveniente del procedimiento estándar de Máxima Verosimilitud (MV) no puede ser maximizada analíticamente debido a su alto grado de no-linealidad en parámetros. Como consecuencia, no es posible obtener los parámetros del modelo espacial de manera simultánea.
Por tanto, a través de este documento revisaremos cómo numéricamente es posible obtener los parámetros de un modelo de regresión espacial a través de lo que se conoce como una función concentrada de log-verosimilitud. Esto se realiza reemplazando los estimaciones análiticas derivadas de Máxima Verosimilitud, tanto de \(\beta\) y \(\sigma^{2}\), para encontrar el valor de \(\rho\). Si bien, la función de log-verosimilitud difiere en propiedades en comparación a la función concentrada, es posible obtener los mismos valores que maximizan la función de log-verosimilitud.

Para lograr este objetivo, plantearemos el siguiente ejemplo hipotético para estimar los coeficientes de un modelo SAR. Un procedimiento similar puede ser utilizado para obtener los coeficientes del modelo SEM.

Insumos

Matriz de Regresores. Esta matriz de variables explicativas contiene el intercepto (vector de unos)

X<-cbind(c(1,1,1,1,1),
         c(0.6,1.0,1.6,2.6,2.2))
X
##      [,1] [,2]
## [1,]    1  0.6
## [2,]    1  1.0
## [3,]    1  1.6
## [4,]    1  2.6
## [5,]    1  2.2

Variable Dependiente

Y<-cbind(c(0.4,0.6,0.9,1.1,1.2))
Y
##      [,1]
## [1,]  0.4
## [2,]  0.6
## [3,]  0.9
## [4,]  1.1
## [5,]  1.2

Matriz de peso espacial (estandarizada por filas)

W<-rbind(c(0, 1/2, 1/2, 0, 0),
         c(1/3, 0, 1/3, 1/3, 0), 
         c(1/3, 1/3, 0, 1/3, 0),
         c(0, 1/3, 1/3, 0, 1/3), 
         c(0, 0, 0, 1, 0))
W
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.0000000 0.5000000 0.5000000 0.0000000 0.0000000
## [2,] 0.3333333 0.0000000 0.3333333 0.3333333 0.0000000
## [3,] 0.3333333 0.3333333 0.0000000 0.3333333 0.0000000
## [4,] 0.0000000 0.3333333 0.3333333 0.0000000 0.3333333
## [5,] 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000

Procedimiento

Paso 1: Regresión de Y en X. También, obtenemos los residuos de la regresión

BETA_YOLS<-solve(t(X)%*%X)%*%t(X)%*%Y
BETA_YOLS
##           [,1]
## [1,] 0.2164706
## [2,] 0.3897059
RES_YOLS<-Y-X%*%BETA_YOLS
RES_YOLS
##              [,1]
## [1,] -0.050294118
## [2,] -0.006176471
## [3,]  0.060000000
## [4,] -0.129705882
## [5,]  0.126176471

Paso 2: Regresión de WY en X. Obtenemos los residuos

WY<-W%*%Y #REZAGO ESPACIAL DE LA VARIABLE Y
BETA_WYOLS<-solve(t(X)%*%X)%*%t(X)%*%WY
BETA_WYOLS
##           [,1]
## [1,] 0.6558824
## [2,] 0.1213235
RES_WYOLS<-WY-X%*%BETA_WYOLS
RES_WYOLS
##             [,1]
## [1,]  0.02132353
## [2,]  0.02279412
## [3,] -0.15000000
## [4,] -0.07132353
## [5,]  0.17720588

Paso 3: Maximizamos la función concentrada de verosimilitud usando los residuos para obtener \(\rho\).

set.seed(1000)
n<-5
e_0<-RES_YOLS
e_1<-RES_WYOLS
e_0t<-t(e_0)
e_1t<-t(e_1)
I<-diag(5)
W<-rbind(c(0, 1/2, 1/2, 0, 0),
          c(1/3, 0, 1/3, 1/3, 0), 
          c(1/3, 1/3, 0, 1/3, 0),
          c(0, 1/3, 1/3, 0, 1/3), 
          c(0, 0, 0, 1, 0))

x<-seq(-1,1,0.0001)
Lc<-x
rho<- -0.9999
for (j in 1:20000) {
  if(j<10000) {
  d <- -(10000-j)/10000 
  }
  else {
  d <- (j-1-10000)/10000  
  }
  Lc[j] <- - n/2*log((e_0t-d*e_1t)%*%(e_0-d*e_1)/n) + log(det(I-d*W))
  if (j==1){
  Lcmax <- Lc[1] 
  }
  if (Lc[j]>Lcmax) {
  rho <- d 
  Lcmax <- Lc[j]
  }
}
Graficaremos la función log concentrada de Máxima Verosimilitud
plot(x,Lc, main="The Log Concentrated Function for SLM",
     ylab="ln L(rho)", xlab="rho",
     type="l", 
     col="blue")
grid(4, 5, lty = "solid")

Nos aseguramos que el valor de \(\rho\) sea aquel que maximiza la función de maxima verosimilitud. Podemos obtener este valor reemplazando \(\rho\) en la función log concentrada de máxima verosimilitud.

lnC<-- n/2*log((e_0t-rho*e_1t)%*%(e_0-rho*e_1)/n) + log(det(I-rho*W))
lnC
##          [,1]
## [1,] 12.57298

Agregamos el valor de \(\rho\) en el gráfico. Con ello, podemos que el valor obtenido máximiza la función log concentrada de máxima verosimilitud.

plot(x,Lc, main="The Log Concentrated Function for SLM",
     ylab="ln L(rho)", xlab="rho",
     type="l", 
     col="blue")
grid(4, 5, lty = "solid")
points(rho, lnC, col="red", pch=19)

Paso 4: Obtenemos a través del valor de \(\rho\) los valores de los coeficientes \(\beta\) y \(\sigma^2\).

#Para Beta
b_mle<-BETA_YOLS-rho*BETA_WYOLS
b_mle
##            [,1]
## [1,] 0.02810118
## [2,] 0.35486176
#ML estimator for the error variance
e_mle<-e_0-rho*e_1

#Para la Varianza 
sigma2_mle<-(1/n)*(e_0t-rho*e_1t)%*%(e_0-rho*e_1)
sigma2_mle
##             [,1]
## [1,] 0.006313458
Opcional: Usando el paquete “spatialreg” podemos obtener los mismos coeficientes.

Este comando, como todos existentes se basan en la lógica descrita anteriormente para obtener los coeficientes de regresión.

library(spdep) #cargamos paquete para convertir matrix de peso espacial en el lenguaje que necesita el comando.
## Loading required package: spData
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(spatialreg) #Este paquete contiene el modelo SAR.
## Loading required package: Matrix
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
#Matriz de peso espacial
SW<-mat2listw(W)
## Warning in mat2listw(W): style is M (missing); style should be set to a valid
## value
SW
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 5 
## Number of nonzero links: 12 
## Percentage nonzero weights: 48 
## Average number of links: 2.4 
## 
## Weights style: M 
## Weights constants summary:
##   n nn S0  S1       S2
## M 5 25  5 4.5 21.05556
#Variables 
Y<-cbind(c(0.4,0.6,0.9,1.1,1.2))
X2<-cbind(c(0.6,1.0,1.6,2.6,2.2))

#Estimamos el modelo SAR. 
MOD.SAR<-lagsarlm(Y~X2, listw = SW)
summary(MOD.SAR)
## 
## Call:lagsarlm(formula = Y ~ X2, listw = SW)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.109224 -0.056418 -0.012722  0.075287  0.103076 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 0.028118   0.164543  0.1709    0.8643
## X2          0.354865   0.056629  6.2664 3.694e-10
## 
## Rho: 0.28717, LR test value: 0.86648, p-value: 0.35193
## Asymptotic standard error: 0.21686
##     z-value: 1.3242, p-value: 0.18543
## Wald statistic: 1.7536, p-value: 0.18543
## 
## Log likelihood: 5.478291 for lag model
## ML residual variance (sigma squared): 0.0063135, (sigma: 0.079458)
## Number of observations: 5 
## Number of parameters estimated: 4 
## AIC: -2.9566, (AIC for lm: -4.0901)
## LM test for residual autocorrelation
## test value: 3.9562, p-value: 0.046698

Podemos ver que los valores obtenidos son los mismos en comparación al procedimiento realizado.