How to estimate spatial models with the use of the likelihood function
André Braz Golgher1
Abstract
This paper is part of a series that discusses introductory concepts of spatial econometrics 2. The
texts were written in Portuguese and intend to present this field of study to students at upper
undergraduate to graduate levels in Economy and in Regional Sciences. The main objective of
this third text is to discuss how to apply the likelihood function to estimate different spatial
models. First, I present this function in a general perspective and obtain the maximum likelihood
estimator for the normal distribution. Then, I show some of the proprieties of this estimator that
makes it very popular as an estimation method, in particular for spatial models. Afterwards, I
develop the mathematical expressions of the maximum likelihood estimators for the OLS, spatial
error and spatial lag models. Next, I present some illustrative simulations using Matlab that
address the concepts discussed in the text. Finally, I obtain the covariance matrix for the spatial
lag model as an example for the other spatial models.
Key words: spatial models; likelihood function; maximum likelihood estimators.
1
Associate Professor at the Economics Department at the Cedeplar/FACE/UFMG, visiting scholar at the Regional
Research Institute (RRI) at the West Virginia University (WVU) and visiting scholar at the Carolina Population
Center (CPC) at the University of North Carolina (UNC) in Chapel Hill.
2
A large part of the references selected for these texts comes from the readings assigned by professor Donald
Lacombe in his course ARE 693L Spatial Econometrics (spring 2012) at the WVU
(http://community.wvu.edu/~djl041/teaching.html).
1
Estimando modelos espaciais com a função de verossimilhança
André Braz Golgher3
Resumo
Esse texto faz parte de uma série que apresenta pontos introdutórios da econometria espacial. O
objetivo principal deste terceiro texto é discutir como utilizar a função de verossimilhança para
estimar diferentes modelos espaciais. Inicialmente, apresenta-se essa função e discute-se como
obter o estimador de máxima verossimilhança utilizando a distribuição normal como exemplo.
Em seguida, apresentam-se algumas das propriedades do estimador, que o tornam muito popular
em diversos tipos de estimação, em particular na estimação dos modelos espaciais. Depois,
desenvolvem-se as expressões matemáticas que são utilizadas na obtenção dos estimadores de
máxima verossimilhança para os modelos de mínimos quadrados ordinários, de lag espacial e de
erro espacial. Posteriormente, apresentam-se simulações feitas no Matlab que ilustram a
aplicação dos conceitos discutidos no texto. Por fim, apresenta-se a matriz de covariância para o
modelo de lag espacial como exemplo para os demais modelos espaciais.
Palavras chave: modelos espaciais; função de verossimilhança; estimador de máxima
verossimilhança.
3
Professor do Cedeplar/FACE/UFMG, pesquisador visitante do Regional Research Institute (RRI) da West Virginia
University (WVU) e pesquisador visitante do Carolina Population Center (CPC) na University of North Carolina
(UNC) em Chapel Hill.
2
1 - Introdução
Esse texto faz parte de uma série que apresenta a econometria espacial em pontos introdutórios4.
No primeiro texto, “Introdução à Econometria Espacial”, foram discutidos alguns conceitos
introdutórios sobre a econometria espacial, onde foram apresentados alguns dos modelos
espaciais, incluindo motivações teóricas para o uso desses modelos. No segundo texto da série,
“Interpretando os coeficientes dos modelos espaciais”, discutiu-se como interpretar os
coeficientes obtidos nos diversos modelos, e também como calcular os efeitos diretos, indiretos e
totais de diversas ordens.
Neste terceiro texto, discute-se como utilizar a função de verossimilhança para estimar diferentes
modelos espaciais. Inicialmente, apresenta-se a função de verossimilhança e discute-se o
estimador de máxima verossimilhança, utilizando a distribuição normal como exemplo. Em
seguida, apresentam-se algumas das propriedades do estimador, que o tornam muito popular em
diversos tipos de estimação. Depois, desenvolve-se passo-a-passo as expressões que são
utilizadas na obtenção dos estimadores de máxima verossimilhança de diferentes modelos
espaciais. Uma vez discutidos os conceitos teóricos envolvidos na estimação dos modelos
espaciais com estimadores de máxima verossimilhança, apresenta-se simulações feitas no Matlab
que exemplificam a aplicação desses conceitos. Por fim, apresenta-se a matriz de covariância
para o modelo de lag espacial como exemplo para os demais modelos espaciais.
Segundo Elhorst (2010), são três os principais métodos desenvolvidos para a estimação de
modelos espaciais: método da máxima verossimilhança (ML), variáveis instrumentais ou o
método de momentos generalizados (IV/GMM), e a abordagem Bayesiana com Monte Carlo via
cadeias de Markov (MCMC).
Segundo Anselin (1988), o uso da função de máxima verossimilhança era a abordagem mais
familiar para a estimação e teste de hipóteses de modelos espaciais na época de publicação de
seu influente livro. Greene (2003) afirmou que os modelos Bayesianos estavam se tornando
muito populares em aplicações da econometria, mas que o método da máxima verossimilhança
permanecia sendo o preferido na maioria dos campos de estudo. Neste texto discutimos o método
4
Grande parte do material citado aqui foi selecionado da ementa do curso ARE 693L Spatial Econometrics (spring
2012) (http://community.wvu.edu/~djl041/teaching.html) ministrado por Donald Lacombe do RRI da WVU.
3
mais usado que é o de máxima verossimilhança. Para uma discussão sobre o método Bayesiano
ver LeSage e Pace (2009).
Uma das razões principais do desenvolvimento do método com variáveis instrumentais ou o
método de momentos generalizados (IV/GMM) foi para facilitar os cálculos computacionais na
estimação dos modelos, que podem ser demandantes devido a matrizes de ordem n x n presentes
no processo de obtenção dos estimadores. Segundo Elhorst (2010), LeSage e Pace (2009)
mostram evidências conclusivas que essas limitações computacionais são coisas do passado. Para
uma discussão sobre esse método ver a extensiva série de publicações de Kelejian e Prucha
(Kelejian, 2011) .
Este texto foi dividido em dez seções, incluindo essa introdução. Na seção 2 é apresentado o
estimador de máxima verossimilhança, utilizando como exemplo a distribuição normal. Em
seguida, discutem-se algumas das propriedades desse estimador. A seção 4 apresenta os modelos
espaciais e discute algumas das similaridades entre eles que facilitam a estimação dos modelos.
Nas seções 5, 6 e 7 são desenvolvidos respectivamente os estimadores de máxima
verossimilhança para o MQO, para o modelo de lag espacial e para o modelo de erro espacial. A
seção 8 mostra simulações ilustrativas no Matlab sobre os conceitos discutidos nas seções
anteriores. Na seção seguinte, obtém-se a matriz de covariância do modelo de lag espacial, onde,
como veremos, os cálculos são bastante extensos. A última seção conclui o texto. Além dessas
seções são incluídos três apêndices econométricos no fim do texto, como forma de tornar a
discussão dos conceitos do corpo do texto mais completa.
2-O estimador de máxima verossimilhança
Nesta seção apresentamos alguns pontos relacionados com a função de verossimilhança e com a
utilização desta para a obtenção de estimadores. Inicialmente, apresenta-se a função de densidade
de probabilidade (pdf), utilizando a distribuição normal como exemplo. Em seguida, obtém-se a
função de verossimilhança e derivam-se os estimadores de máxima verossimilhança para essa
distribuição. O objetivo é dar ao leitor uma ideia inicial desse método de estimação, que será
posteriormente utilizada na obtenção das funções de verossimilhança dos modelos espaciais.
4
As funções de densidades de probabilidade, f (x) , de distribuições contínuas, como a normal,
b
têm as seguintes propriedades: f ( x)  0 , Pr ob(a  x  b)   f ( x)dx, e 0  Pr ob(a  x  b)  1
a
(Greene, 2003).
No caso particular da distribuição normal com média  e variância  2 , N ( ,  2 ) , essa função
tem o seguinte formato: f ( y;  , ) 

1
2
(2 2 )
1
e
1
2 2
 y   2
.
2
Se tivermos duas observações independentes, y1 e y 2 , obtemos a densidade conjunta destas
multiplicando as respectivas expressões:
1
1

 2  y1   2 
 2  y2   2 
1
1



.
2
2
F ( y1 , y2 ;  , ) 
e
e
 (2 2 ) 12
 (2 2 ) 12




2
Note que para distribuição normal temos dois parâmetros,  e  2 . Para simplificar e generalizar
a notação, definimos   ( , 2 ) , como o vetor composto por esses parâmetros e y  ( y1 , y2 )
como o vetor com as observações. Reescrevemos a expressão acima como:
1

 2  yi   2 
1

.
2
F ( y; )  
e
1


2
i 1 ( 2 ) 2


2
Seguindo esse raciocínio, a densidade conjunta de n observações independentes e identicamente
distribuídas, que no caso acima é a distribuição normal, é dada pelo produto das densidades
individuais, e é denominada função de verossimilhança:
n
L( y )   f ( yi  ).
i 1
Ao contrário da expressão da pdf, f ( yi  ), onde temos os dados condicionados aos parâmetros
( yi  ) , note que na função de verossimilhança, L( y), temos ( y ) , ou seja, os parâmetros são
condicionados aos dados. Assim, o objetivo do uso da função de verossimilhança (FV) eh, a
5
partir dos dados empíricos, estimar os parâmetros da função. Estima-se os parâmetros que
maximizam essa função, isto eh, maximizam a probabilidade de ocorrência dos dados empíricos.
Seguindo a discussão com a distribuição normal, a FV para n observações é escrita como:
n
L(  ,  y )  
i 1

1
2
(2 )
2
1
e
1
2
2
 yi   2

n
1
(2 )
2
2
n
e

1
2
2
 yi   2
2 i 1


1
(2 )
2
n
e
1
2 2
n
  yi  2
i 1
2
Para acharmos o máximo dessa função, podemos derivar parcialmente essa expressão com
relação aos dois parâmetros,  e  2 , igualando cada uma das derivadas a zero. Entretanto, em
geral, é mais simples trabalhar com o logaritmo dessa função. Note que o máximo da FV é o
mesmo ponto que o máximo do logaritmo dessa função, pois o logaritmo é uma função
monotonicamente crescente.
Vejamos.
1

 2   yi   2 
n
n
1
2 i 1
2
2
  ln  (2 2 )  2     1   yi   2
ln L(  ,  y )  l (  ,  y )  ln 
e
n
2

  2 i 1
 (2 2 ) 2



1 n
 n
 n

 yi   2 .
   ln 2      ln  2   
2 
 2
 2
 2 i 1

n



Uma vez obtido o logaritmo da FV, derivamos a expressão, inicialmente com relação a  ,
obtemos o estimador de máxima verossimilhança (MV) para esse parâmetro:
l (  ,  2 y )




n
 n 
 n
 1 n
2
 yi   2    12   yi     0
  2 ln 2     2  ln     2 2 i




 1

   i 1


  yi     0
n
i 1
ny  n  0
 MV  y
O estimador de MV para  MV é a média entre as observações.
Em seguida, derivamos com relação à  2 , e obtemos o estimador para a variância:
6
l (  , 2 y )
 2

  n 
 n
 1 n
2
2





ln
2



ln






  2   yi     
2 
  2 
 2
 2 i 1


 n   1
   2   
 2   2  2


n
n


  yi   2    1   n    1   yi  y 2   0
2 
2
2 

  2  
  i 1

i 1
 
 n   
2
 MV
 
1 n
2
  yi  y   0
2 
i

1
 
1 n
2
     yi  y 
 n i 1
 
Como discutido, a FV é muito utilizada na obtenção de estimadores em diversos campos de
aplicação e, em particular, também são aplicados aos modelos espaciais. Esse fato é devido as
interessantes propriedades desses estimadores (ver Greene 2003 para uma discussão mais
detalhada). A próxima seção discute algumas delas.
3 – Propriedades do estimador de MV
Os estimadores de MV apresentam algumas características que os tornam muito populares em
diversos tipos de aplicações, inclusive na estimação dos modelos espaciais. Os estimadores de
MV são consistentes e assintoticamente eficientes, com matriz de covariância assintótica
1
   2 ln L 
 . Além disso, os
definida no limite inferior de Cramér-Rao: Var (ˆ)  [ I ( )]1   E 




'



estimadores são assintoticamente normalmente distribuídos: ˆ ~ a ~ N[ , I ( ) 1 ] (ver apêndice
econométrico 1 para maiores detalhes sobre esses conceitos).
Para que esses estimadores tenham essas propriedades, existem certas condições formais que
devem ser satisfeitas. Para o caso particular dos modelos espaciais, essas foram descritas em
estudos realizados na década de 80 (Anselin, 1988). Dentre outras, essas condições formais são:
o logaritmo da FV deve existir para os valores dos parâmetros sendo considerados; essa função
deve ser continuamente derivável até a segunda ou a terceira ordens; e as derivadas parciais
devem assumir valores finitos.
7
Como vimos na seção anterior quando foram obtidos os estimadores de MV da distribuição
normal, todas essas condições são satisfeitas para essa distribuição. Uma vez apresentadas as
propriedades desses estimadores, continua-se esse exemplo, obtendo a matriz de covariância dos
estimadores, que deve ser positivamente definida.
Obtivemos os estimadores na seção anterior com o uso das derivadas de primeira ordem:
l  1  n

  yi   
   2 i 1
 1 n
l
 1 
2








n

y

y





i
  2 i 1
 2  2 2  



 
As covariâncias dos estimadores da distribuição normal são obtidas a partir do limite de CramérRao, com a seguinte matriz:
  2 ln L

2
  2 ln L 
   E  2 
I ( )   E 
  ln L
  ' 
 2
  
 2 ln L 

 2 
 2 ln L 

( 2 ) 2 
Para determinar matriz de covariância, calculamos as derivadas de segunda ordem:
 2l
  1 n
 n 

 2   yi      2 .
2

   i 1
 
 1 n

 2l
 2l
  l 
  1  n
 2 2   yi   .







y









i

 2  2   2     2    2 i 1

 ( ) i 1
 2l


2 2
 ( )
 2
 1  

 1 n
2 
 2   n    2   yi  y    
 2  
  i 1
 
 


 1 n
1 
 1    1  n
2
2 












n

y

y


y

y


  2    2 2  i
i
2 2 
  2 i
 

 1
 2( )  
  2    
i 1
 
 

 n   1  n
2



 

y

y

.

i
3
2 2 

 2( )    2 i 1
 
8
Em seguida, obtemos os valores esperados dessas expressões, lembrando que E yi    e que


E ( yi   ) 2  n 2 :
  2l 
  n 
 n 
E  2   E   2    2 
 
   
  
  1 n

  2l 
 1  n
E
 E   2 2   yi      4  E   yi     0
2

   i 1
  
  ( ) i 1

 n   1  n


 n
  2l 
  yi  y 2    n    1  E   yi  y 2 

E
 E 
3
3
2 2
2 2 
2
2



 2( )    2 i 1
  2( )    2  i 1
 ( ) 
 
 
 n   1 
 n 


 
(n 2 )  
3
2 2 
2 2 
 2( )    2 
 2( ) 
 
Obtemos assim a matriz de covariância com o limite de inferior de Cramér-Rao:
  2 ln L

 2
[ I ( )]   E  2
  ln L
 2
  
  2 


  n 
1
[ I ( )]  
 0


 2 ln L    n 

  2 
0 
2
     


 2 ln L  
 n 
 4 
  0
( 2 ) 2  
 2  



4 
 2 

 

 n 
0
Além dessas propriedades discutidas, outra propriedade do estimador de MV é que ele é
invariante a uma mudança por uma função injetiva e continuamente derivável. Ou seja, se temos
um estimador ˆ para um parâmetro  e temos uma função injetiva e continuamente derivável
f (x), como uma reta f ( x)  ax  b, então o estimador de f ( ) será f (ˆ) . Essa propriedade será
utilizada na obtenção dos estimadores para os modelos econométricos.
Uma vez obtidos os estimadores de MV para a distribuição normal, utilizaremos esses mesmos
conceitos para discutir os estimados de MV dos modelos espaciais. A próxima seção apresenta
9
esses modelos e discute algumas regularidades apresentadas por eles, que simplificarão a
discussão subsequente.
4-Modelos espaciais
Antes de discutirmos como estimar cada um dos modelos espaciais, eles são analisados em pares,
pois isso vai simplificar em muito a discussão futura. A figura 1, similar a figure 1 em Elhorst
(2010), mostra a relação entre diversos modelos espaciais.
Como mostra essa figura, o 7 - modelo de lag de X (SLX) é descrito pela seguinte equação:
Y  X  WX   ,  ~ N (0, 2 I n ).
Note que podemos reescrever esse modelo com Z  [ X ,WX ] e   [ , ] da seguinte forma:
y  Z   ,  ~ N (0, 2 I n ).
Essa expressão é exatamente igual ao 8 - MQO. Assim, os modelos são estimados de forma
similar depois de feitas essas mudanças no modelo SLX e, assim, não abordaremos este ultimo
em separado o modelo SLX.
De forma similar, o 4 - modelo de erro espacial de Durbin (SDEM) é descrito pela seguinte
equação:
Y  X  WX  u,
u  Wu   ,  ~ N (0,  2 I n )
Fazendo a mesma transformação acima, obtemos a equação do 6 - modelo de erro espacial.
Assim, também não abordaremos o modelo SDEM em separado.
Por sua vez, o 2 - modelo espacial de Durbin (SDM) pode ser escrito como o 5 - modelo de lag
espacial com essas mesmas transformações (LeSage e Pace, 2009) e esse primeiro também não é
discutido em separado. Fazendo uso dessa mesma transformação, 8 - modelo de Manski pode ser
escrito como o 3 - modelo de Kelejian-Prucha. Assim, esse primeiro também não é discutido.
10
Em todos os modelos mostrados na figura 1 assumimos que  ~ N (0, 2 I n ) . Assim, utilizamos a
distribuição normal como base para toda a discussão subsequente. Seguimos a discussão de
como estimar os modelos fazendo uso da FV inicialmente com o MQO, e em seguida com o
modelo de lag espacial e o modelo de erro espacial.
O modelo de Kelejian-Prucha engloba esses dois últimos modelos. Assim, pode-se estender a
discussão apresentada sobre a metodologia de estimação dos modelos de lag espacial e de erro
espacial na obtenção das expressões dos estimadores de MV desse primeiro. Assim, para não
estender em demasiado este texto, esse modelo não é discutido em separado.
11
Figura 1 – Modelos espaciais
1 - Modelo de Manski
Se
Se
Se
 0
 0
 0
4 - Modelo erro espacial de Durbin
(SDEM)
3 - Modelo de Kelejian-Prucha (SAC)
2 - Modelo espacial de Durbin (SDM)
 0
Se
 0
Se
Se
 0
Se 
Se
 0
 0
7 – Modelo de lag de X (SLX)
0
6 – Modelo de erro espacial (SEM)
5 – Modelo de lag espacial (SAR)
Se
Se 
Se
0
 0
 0
8 - MQO
12
5-Estimadores de MV para o MQO
Essa seção apresenta os estimadores de MV para o MQO, que, apesar de não ser um modelo
espacial, serve de base de comparação para os modelos espaciais, em particular para os modelos
de lag espacial e de erro espacial. Além disso, como vimos na seção anterior, os estimadores do
modelo de lag de X (SLX) são similares aos do modelo de MQO. O apêndice econométrico 3
apresenta algumas das hipóteses do MQO, além de mostrar o processo de obtenção do estimador
de mínimos quadrados ordinários. Essa apresentação é incluída neste texto para que seja feita
uma comparação entre este estimador e o estimador de MV.
O MQO tem como equação a seguinte expressão:
yi  X i    i ,  i ~ N (0, 2 )
Note que o erro tem distribuição normal, com média zero e variância  2 . Assim, a função de
densidade de probabilidade para o erro tem o seguinte formato:
f 

1
(2 2 )
1
e
1
 x   2
2 2

2

1
(2 2 )
1
e
1
 i 0 2
2 2
2


1
(2 2 )
1
e
1
 i 2
2 2
2
A função de verossimilhança é escrita como:
L

1
(2 2 )
n
e
2
n
1
2
2
  i 2
i 1


1
(2 2 )
n
e
1
 '
2 2
2
Entretanto, no modelo MQO desejamos obter o estimador para  e  2 a partir das variáveis
dependentes e independentes, obtidas empiricamente. Assim, reescrevemos a expressão acima
com:
 i  yi  X i 
O Jacobiano da transformação de n variáveis por m funções (ver apêndice econométrico 4) é
dado por:
13
x1
y1
J  abs ...
xn
y1
x1
y n
... ... .
xn
...
y n
...
Assim, para o caso particular do MQO, o Jacobiano da transformação de  i para yi tem o
seguinte valor: J 
 i
 1.
yi
Incorporando essas expressões na FV do erro, obtemos a FV do MQO:
L
1e
1
(2 2 )
n

1
2 2
( y  X )'( y  X )

2

1
(2 2 )
n
e
1
2 2
( y  X )'( y  X )
2
Tomando o logaritmo da função, utilizando a relação para matrizes ( A  B)'  A' B' , e a relação
para vetores X 'Y  Y ' X , temos:


 n
 n
 1 
l    ln 2      ln  2   2 ( y  X )' ( y  X )
 2
 2
 2 
 n
 n
 1 
   ln 2      ln  2   2  y' y  ( X )' y  y' X  ( X )' ( X )
 2
 2
 2 




 n
 n
 1 
   ln 2      ln  2   2  y' y  2( X )' y  ( X )' ( X ).
 2
 2
 2 
Derivando com relação ao vetor  , obtemos o estimador para esses parâmetros:
l
 1 
  2  2 X ' y  2 X ' X   0

 2 
X ' y  X ' X
ˆ
 ( X ' X ) 1 X ' y
MQO
Note que esse estimador é exatamente igual ao obtido para o estimador de mínimos quadrados
ordinários, discutido no apêndice econométrico 3.
14
Derivando com relação a  2 , obtemos o estimador para a variância:
l
 n   1 
( y  X )' ( y  X )  0
   2   
2
2 2 

 2   2( ) 
1
n
2
ˆ MQO
  ( y  X )' ( y  X ) 
e' e
.
n
Por fim, obtemos a matriz de covariância utilizando o limite de Cramér-Rao.

  l     1 
 1 
 X'X 

 
   2  2 X ' y  2 X ' X    2 2 X ' X    2 
   '     2 
 2 
  

  2l 
 X'X 
E
  2 

  
  ' 
  l 
   1 

 
 2 X ' y  2 X ' X    12 2  2 X ' y  2 X ' X 

2 
2  
2 
   '     2 
  2( ) 
 1 
  1 
  2l 
E
 E  2 2  X ' y  X ' X    2 2  E X ' y  X ' X   0
2
  
 ( ) 
  ( ) 

  l 
 
n   1 
 
( y  X )' ( y  X ) 

 
2 
2 
2 
2 
2 2 
       2   2( ) 

 n  
 n  
1 
1 
    2 3 ( y  X )' ( y  X )  
    2 3  ' 
 
2 2 
2 2 
 2( )   ( ) 
 2( )   ( ) 
 n  
  2l 
1    n  
1 
    2 3  '    
    2 3  E '  
E
 E 
2 2
2 2 
2 2 
 ( ) 
 2( )   ( )    2( )   ( ) 
 n  
 n 
1 
    2 3 (n 2 )  

 
2 2 
2 2 
 2( )   ( ) 
 2( ) 
15
   2 ln L
 
 2
1
[ I ( )]   E  2
   ln L
   2  '
 
 2 ln L 

 2 
 2 ln L 

 ( 2 ) 2 
1
  2

 X ' X
 
 0








4 
 2 

 

 n 
0
Como discutido na seção anterior, esse procedimento descrito aqui também é utilizado para o
modelo SLX.
6 – Estimadores de MV para o modelo de lag espacial
Segue a discussão de como obter os estimadores para o modelo de lag espacial e, como vimos,
também para o modelo espacial de Durbin. A matriz de covariância será obtida posteriormente
na seção nove do texto, devido a extensão dos cálculos. A discussão abaixo é baseada em
Doreian (1981).
Como vimos, o modelo de lag espacial é expresso pela seguinte equação:
Y  WY  X  
Seguindo os mesmos procedimentos da seção anterior, isolamos o termo do erro, obtemos o
Jacobiano da transformação de  i para yi , e substituindo os resultados na FV do erro:
  y  Wy  X  ( I  W ) y  X
J 
 i
 I  W .
yi
Assim, a FV do modelo de lag espacial é dada por:
 ( y  Wy  X )' ( y  Wy  X ) 
L(  ,  ,  ; y, X )  (2 2 ) n / 2 I  W exp  

2 2


Tomando o logaritmo da função e simplificando a expressão, temos:
16
 ( y  Wy  X )' ( y  Wy  X ) 
l (  ,  2 ,  ; y, X )  ln{(2 2 ) n / 2 I  W exp  
}
2 2


 ( y  Wy  X )' ( y  Wy  X ) 
 ln( 2 2 ) n / 2  ln I  W  ln[exp  
]
2 2





n
 1 
( y  Wy  X )' ( y  Wy  X )
ln( 2 )  ln( 2 )  ln I  W   
2 
2
 2 
Derivando com relação a  2 , obtemos o estimador da variância:
 1 
l
n
( y  Wy  X )' ( y  Wy  X )   0
  2  
2
2 2 

2
 2( ) 

 1 
 1 
 2   n    2 ( y  Wy  X )' ( y  Wy  X )   0
 2 
 

1
n
2
ˆ SAR
  ( y  Wy  X )' ( y  Wy  X )
Entretanto, não podemos estimar a variância, pois não sabemos os valores de  nem de  , que
estão presentes na expressão. Assim, devemos utilizar alguma expressão para estimarmos esses
parâmetros para depois estimarmos ˆ 2 .
Assim, derivando o log da verossimilhança com relação a  , temos:
l  1 
 1 
   2  X ' ( y  Wy  X )     2  X ' (( I  W ) y  X )   0
   
  
Para prosseguir o raciocínio, fazemos uma mudança de variável com z  ( I  W ) y. Daí
ficamos com:
 X ' ( z  X )  0
Manipulando essa expressão, obtemos o estimador do modelo:
ˆSAR  ( X ' X ) 1 X ' z
17
Note que para estimarmos ̂ SAR devemos estimar primeiro o parâmetro  , para depois
substituirmos em z.
Substituindo a mudança de variável z  ( I  W ) y e o resultado para ̂ SAR na expressão do
estimador da variância, lembrando que ( AB)'  B' A' , ficamos com:
1
n
2
ˆ SAR
  ( z  X )' ( z  X )
1
n
2
ˆ SAR
  ( z  X ( X ' X ) 1 X ' z )' ( z  X ( X ' X ) 1 X ' z )
1
  ( z ' z ' ( X ( X ' X ) 1 X ' )' )( z  X ( X ' X ) 1 X ' z )
n
1
  ( z ' z  2 z ' ( X ( X ' X ) 1 X ' ) z  z ' X ( X ' X ) 1 X ' X ( X ' X ) 1 X ' z )
n
1
  ( z ' z  2 z ' ( X ( X ' X ) 1 X ' ) z )  z ' X ( X ' X ) 1 IX ' z )
n
1
1
   z ' ( I  ( X ( X ' X ) 1 X ' ) z    z ' Mz,
n
n
onde M  I  ( X ( X ' X ) 1 X '.
Note que agora o estimador da variância depende somente de  . Ou seja, neste ponto da
2
(  ). Mas
discussão os dois estimadores já discutidos são função desse parâmetro: ˆSAR (  ) e ˆ SAR
ainda não sabemos o valor de  . Assim, devemos estimar o parâmetro  para depois
estimarmos os demais via essas expressões.
Retomamos a expressão do log da FV.



n
 1 
l ( ˆ (  ),ˆ 2 (  ),  ; y, X )   ln( 2 )  ln(ˆ 2 )  ln I  W   
( y  Wy  Xˆ )' ( y  Wy  Xˆ )
2 
ˆ
2
2



1
 
2
Sabendo que ˆ SAR
  ( y  Wy  Xˆ )' ( y  Wy  Xˆ ) , temos:
n
18


 
n
ln( 2 )  ln ˆ 2 )  ln I  W
2
n
 1
    ( y  Wy  Xˆ )' ( y  Wy  Xˆ )
ˆ
( y  Wy  X )' ( y  Wy  Xˆ )
 2
l



n
1  ln(2 )  n ln   1 ( y  Wy  Xˆ )' ( y  Wy  Xˆ )   ln I  W
2
2  n 

Um ponto chave na obtenção do estimador de MV é o logaritmo do Jacobiano. Para simplificar
essa relação devemos reescrevê-lo. Note que podemos calcular o determinante de uma matriz A
como o produto de seus autovalores: A   i .
i
Daí, temos:
W   i , onde i são os autovalores da matriz de peso.
i
I  W   (  i ) ,
i
I  W   (1  i ) .
i
Tomando o logaritmo, temos:
ln I  W  ln  (1  i )   ln(1  i ).
i
i
Note que o logaritmo só é definido para valores positivos. Assumindo que os autovalores são
reais, dessa relação, temos para todo i:
(1  i )  0
1  i
Se i é positivo:  
1
i
. Como essa relação é valida para todo i, daí temos que  
1
max
, onde
max é o maior autovalor real positivo.
19
Se i é negativo:  

1
min
1
i
. Como essa relação também é valida para todo i, dai temos que
, onde min é o menor autovalor real negativo.
 1
1
,
Assim obtemos o intervalo para os valores de  , onde a FV é definida:   
 min max

 .

Se a matriz de peso for normalizada temos os seguintes valores: min  1 e max  1 .
Esse mesmo raciocínio é valido para o parâmetro de correlação espacial do modelo de erro
espacial.
Substituindo a relação acima na função de log de FV e substituindo K  
n
1  ln(2 )
2
constante, temos:
l
ˆ
ˆ 

n
1  ln(2 )  n ln ( y  Wy  X )' ( y  Wy  X )   ln I  W
2
2 
n

 
n
l  K  ln ˆ 2   ln(1  i )
2
i
Sabemos que:
1
1
1
n
n
n
1
1
  ( y ' y 'W ' ) M ( y  Wy )   ( y ' My  2 y ' MWy   2 y 'W ' MWy)
n
n
ˆ 2    z ' Mz   (( I  W ) y )' M ( I  W ) y   ( y  Wy)' M ( y  Wy)
Substituindo essa expressão na função do log da FV, obtemos a seguinte equação:
 
n
l  K    ln ˆ 2   ln(1  i )
2
i

 n   1 
 K    ln   ( y ' My  2 y ' MWy   2 y 'W ' MWy)    ln(1  i )
 2   n 
 i
20
Simplificando essa expressão, suprimindo os valores constantes e reescrevendo, temos:


2
g (  ; y, X )  ln y' My  2 y' MWy   2 y'W ' MWy    ln(1  i )
n i
Devemos obter o mínimo dessa função para obtermos o valor do parâmetro  que maximiza a
FV. Em vez de resolvermos essa expressão analiticamente, ela é minimizada por procedimentos
numéricos com a utilização de métodos computacionais. Isto é, variam-se os valores de  passoa-passo em pequenos degraus em um intervalo pré-definido e obtém-se o mínimo da função.
Normalmente se utiliza o intervalo  [0,1). Os exemplos de correlação positiva são muito mais
numerosos que os de correlação negativa, o que justifica essa escolha. Valores negativos para os
parâmetros espaciais são observados em alguns trabalhos empíricos, como na competição entre
empresas em um mesmo espaço urbano ou na disputa por espaço geográfico entre espécies
vegetais ou animais (Griffth e Arbia, 2010) .
Ou seja, na prática os parâmetros  e  não são incluídos na otimização da função de
verossimilhança. Isso facilita em muito os procedimentos de maximização. Para uma discussão
mais aprofundada ver LeSage e Pace (2009).
Uma vez obtido o parâmetro  que minimiza a expressão acima, são estimados os demais
parâmetros, ˆ (  ) e ˆ 2 (  ), e o processo de obtenção dos parâmetros se encerra.
A seção seguinte discute o modelo de erro espacial de forma similar ao apresentado aqui.
7- Estimadores de MV do modelo de erro espacial
Seguindo um raciocínio similar ao descrito acima obtemos os estimadores do modelo de erro
espacial (SEM). Como vimos, o modelo tem a seguinte equação:
Y  X  u
u  Wu   ,  ~ N (0, 2 I n ).
21
Em seguida, isolamos  :
Y  X  ( I  W ) 1 
  ( I  W )Y  ( I  W ) X
  ( I  W )(Y  X ).
A partir dessa expressão, calculamos o Jacobiano da transformação de  para Y :
J 

 ( I  W ).
Y
De posse dessa informação, obtemos a FV do modelo:
 [( I  W )(Y  X )]'[( I  W )(Y  X )] 
L(  , ,  ; y, X )  (2 2 ) n / 2 I  W exp  

2 2


 [(Y  X )' ( I  W )' ][( I  W )(Y  X )] 
L(  , ,  ; y, X )  (2 2 ) n / 2 I  W exp  

2 2


Tomando o logaritmo dessa função, temos:
1 
n

l (  , ,  ; y, X )    ln( 2 2 )  ln I  W   
(Y  X )' ( I  W )' ( I  W )(Y  X )
2 
2
 2 
Derivando com relação a  , obtemos o estimador ̂ SEM :
l 
1 
 
(2)( X )' ( I  W )' ( I  W )(Y  X )  0
  2 2 
X ' ( I  W )' ( I  W )(Y  X )  0
[( I  W ) X ]' ( I  W )(Y  X )  0
[( I  W ) X ]' ( I  W )Y  [( I  W ) X ]' ( I  W ) X
[( I  W ) X ]'[( I  W )Y ]  [( I  W ) X ]'[( I  W ) X ]
Fazemos as seguintes mudanças de variáveis:
22
Y *  ( I  W )Y
X *  ( I  W )Y
Ficamos com:
( X *)' Y *  ( X *)' ( X *)
ˆ
 (( X *)' ( X *))1 ( X *)' Y *
SEM
Note que o estimador acima, assim como observado no modelo de lag espacial, depende do valor
da variável de correlação espacial, que no caso do modelo de erro espacial é o valor de  , que
ainda não conhecemos.
Prosseguindo, derivamos com relação à  2 e obtemos o estimador da variância:
l
 n   1 
(Y  X )' ( I  W )' ( I  W )(Y  X )  0
  2   
2
2 2 

 2   2( ) 

 1   1 
 2  n   2 (Y  X )' ( I  W )' ( I  W )(Y  X )   0
 2    

1
n
ˆ SEM   (( I  W )(Y  Xˆ ))' ( I  W )(Y  Xˆ )
2
Note que esse estimador também depende do valor de  . Ou seja, ambos estimadores, ̂ SEM e
ˆ SEM , dependem de  . Assim devemos estima-lo e em seguida obtemos os demais.
2
Retornamos a expressão do log da FV.
n
 1 
l (ˆ 2 ( ),  ; y, X )    ln( 2ˆ 2 )  ln I  W   
(Y  X )' ( I  W )' ( I  W )(Y  X )
2 
2
 2ˆ 
Para simplificar essa relação devemos reescrever o logaritmo do determinante como:
ln I  W   ln(1   i ), onde  i são os autovalores da matriz de peso (utiliza-se essa
i
nomenclatura para autovalores para não confusão com o parâmetro espacial  ).
23
Substituindo a expressão do estimador da variância, ficamos com:
n
 1 
l ( ˆ ,ˆ 2 , ; y, X )    ln( 2ˆ 2 )   ln(1   i )   
(Y  Xˆ )' ( I  W )' ( I  W )(Y  Xˆ )
2 
ˆ
2
2

 


i
n
n
   ln( 2 )  ln ˆ 2 )   ln(1   i )  

2
 2 
i


Simplificando, suprimindo as constantes e reescrevendo, obtemos a função que será minimizada
por procedimentos numéricos:


2
g ( )  ln ˆ 2    ln(1   i )
n i
Na próxima seção todo o procedimento discutido para os modelos MQO, de lag espacial e de
erro espacial é implementado em simulações ilustrativas no Matlab.
8 – Simulação ilustrativa
Este seção apresenta simulações ilustrativas feitas no Matlab que trata dos conceitos discutidos
nas seções anteriores. Como ponto de partida, temos o DGP do modelo de Kelejian-Prucha:
Y  WY  X  u
u  Wu   ,  ~ N (0, 2 I n )
Nas simulações estipulamos valores arbitrários para os termos desse DGP. Esses valores são
baseados em Florax et al (2003) e em Anselin et al (1996). Maiores detalhes das simulações
discutidas por esses autores são apresentados no próximo texto dessa serie.
As simulações foram feitas com 100 observações, como sugerido em Mur e Angulo (2009) como
número mínimo para uma estimativa confiável. As observações foram definidas como regiões,
que foram dispostas em uma estrutura quadrada regular, tipo tabuleiro de xadrez, com cada lado
de dimensão 10. Assim temos 10 x 10 = 100 observações.
Um ponto importante de toda essa análise é referente a relação entre o tamanho da amostra e
dificuldades computacionais de estimação. Note que para uma amostra com n observações, a
24
matriz de peso é n x n. Assim, uma amostra maior aumenta a acurácia da estimação, mas implica
em maior dificuldade computacional, pois as matrizes a serem manipuladas no processo de
estimação tem n2 elementos. Note, porém, que, de forma geral, as matrizes de peso têm muitos
elementos nulos, sendo que a proporção destes tende a aumentar quando o tamanho da amostra
aumenta. O uso de técnicas de matrizes esparsas, não utilizada aqui, pode ser usado para facilitar
a estimação dos parâmetros (Pace, 1997).
A matriz de peso foi definida como de contiguidade por torre, pelos lados, como uma torre se
movimento no xadrez, como mostra o diagrama abaixo. Note que a contiguidade em uma
estrutura regular pode ser definida de outras maneiras, como por exemplo, pelo movimento da
rainha no xadrez, que incluem também as regiões que se tocam apenas no vértice. A matriz foi
normalizada na linha.
Uma vez definida a matriz de peso, o próximo passo é obter a matriz X das variáveis
independentes. Assume-se que temos apenas uma variável exógena, sendo que o valor de cada
uma das 100 observações foi obtido a partir de uma distribuição uniforme com valores entre 0 e
10: X ~ U[0,10] . Os valores são mantidos os mesmos em todas as simulações. Para os
parâmetros temos  0  1 , referente ao intercepto, e 1  1 , referente a variável exógena X.
Esses valores foram incorporados ao DGP do modelo de Kelejian-Prucha. Foram analisados três
casos particulares desse modelo: o modelo MQO com   0 e   0 ; o modelo de lag espacial
com  
1
1
e   0 ; e o modelo de erro espacial com   0 e   .
2
2
25
Em seguida, determinam-se os erros estocásticos. Esses foram obtidos a partir de uma
distribuição normal com média zero, e assumindo  2  1 . Foram gerados 100 valores a partir
dessa distribuição para cada uma das simulações. No total são 100 simulações para o MQO, 20
para o modelo de lag espacial e 20 para o modelo de erro espacial. Foram assim gerados 1040
vetores com 100 valores aleatórios para o erro estocástico.
De posse de toda essa informação, obtêm-se com o DGP do modelo de Kelejian-Prucha, os
valores da variável dependente para esses modelos, também em 1040 vetores com 100 valores
cada um. Ou seja, sabemos quais são os parâmetros que geraram os valores da variável
dependente para cada um dos modelos. Os resultados são discutidos em separado para o MQO, o
modelo de lag espacial e o modelo de erro espacial.
8.1 – MQO
Com a informação descrita acima, geramos os 100 valores da variável dependente e os 100
valores da variável independente para cada uma das 1000 simulações. De posse dessa
informação, estima-se o modelo MQO 1000 vezes utilizando as expressões:
ˆMQO  ( X ' X ) 1 X 'Y
2
ˆ MQO

e' e
.
n
Compara-se, assim, o modelo real que gerou os dados com o modelo estimado.
Como vimos, o valor real para os parâmetros são  0  1 , 1  1 e  2  1 . O valor médio das
1000 estimações para esses parâmetros foi respectivamente de 1,0103, 0,9977 e 0,9749. Ou seja,
os valores estimados são muito próximos dos reais.
Os histogramas a seguir mostram a distribuição de cada um dos parâmetros, todos
aproximadamente normais, com a média descrita acima. Note que a dispersão dos valores de 1 ,
que é o coeficiente de maior interesse empírico, era menor que para os demais parâmetros. Na
26
grande maioria das simulações, os valores ficaram entre 0,95 e 1,05, variações relativamente
pequenas em torno do valor real. Porém, note que os valores absolutos dessas dispersões
dependem dos valores arbitrários incluídos nas estimações.
Histograma 1 – Dispersão dos valores de
0
250
200
150
100
50
0
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Histograma 2 – Dispersão dos valores de
2
1
300
250
200
150
100
50
0
0.85
0.9
0.95
1
1.05
1.1
1.15
1.2
27
Histograma 3 – Dispersão dos valores de
2
300
250
200
150
100
50
0
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
8.2 – Modelo de lag espacial
Também com a informação descrita acima, foram obtidos os 100 valores da variável dependente
e da variável independente para o modelo de lag espacial. Estimou-se o modelo empírico 20
vezes. Note que para cada uma dessas estimativas devemos fazer um procedimento que demanda
mais do computador do que o MQO. Essa é a razão desse número menor de simulações.
O primeiro passo é obter o valor de  que minimiza a expressão abaixo:
2
g (  ;Y , X )  ln Y ' MY  2 Y ' MWY   2Y 'W ' MWY     ln(1  i ),
n i
onde M  I  ( X ( X ' X ) 1 X '.
Assim, inicialmente, obtém-se a matriz M, e estimam-se os autovalores da matriz de peso.
Substituem-se esses resultados em conjunto com o vetor Y e a matriz W na expressão acima.
Esses valores serão constantes para todos os valores de  .
Qual valor de  minimiza essa expressão? Isso é obtido numericamente. No caso da simulação
aqui, variou-se o parâmetro de 0,05 até 0,95 em acréscimos de 0,05, em um total de 19 valores.
Note que poderíamos ter feito de 0 a 0,95 em acréscimos de 0,01. A única diferença seria o
tempo maior do computador. Com cada um dos valores de  , calculou-se o valor da função g.
28
Esse mesmo procedimento foi feito 20 vezes. Note que dependendo do erro gerado na simulação,
temos um valor de variável dependente distinto, o que leva a estimações de  distintas.
Lembrando que o valor real foi estipulado como  
1
, o diagrama abaixo mostra os resultados
2
obtidos nas 20 simulações. Na maioria das simulações, o valor de  que minimizava a expressão
acima ficou em torno desse valor real. O valor médio das 20 simulações foi 0,485.
Gráfico 1 – Valores de
g (  ; Y , X ) para diferentes valores de 
6.5
6
5.5
5
4.5
4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Uma vez obtido o valor de  que minimiza a expressão acima, estima-se os demais parâmetros
do modelo a partir das expressões abaixo:
ˆSAR  ( X ' X ) 1 X ' ( I  W )Y
1
n
2
ˆ SAR
  ( I  W )Y ' M ( I  W )Y
A média dos valores das 20 replicações para  0 , 1 e  2 foram respectivamente 1,0318, 0,9981
e 0,9729, também muito próximos dos valores reais.
29
8.3 – Modelo de erro espacial
De forma similar, de posse da informação descrita acima, foram obtidos os valores para a
variável dependente e independente do modelo de erro espacial e estimaram-se os modelos
empíricos. O primeiro passo foi estimar os 100 autovalores da matriz de peso, W, que são
expressos aqui como  i . Depois, deve-se obter o valor de  que minimiza a expressão abaixo,
sendo n o número de observações:


2
g ( )  ln ˆ 2    ln(1   i ),
n i
2
1
onde ˆ   (( I  W )(Y  Xˆ ))' ( I  W )(Y  Xˆ ),
n
'
ˆ  (( X *) ( X *))1 ( X *)' Y *,
X *  ( I  W ) X ,
e Y *  ( I  W )Y
Inicialmente, substituem-se os valores encontrados para os autovalores, e os vetores X, Y e a
matriz W nesta expressão. Assim, o único parâmetro a determinar é  . Varia-se o valor do
parâmetro numericamente e calcula-se o valor da função g para cada valor de  . No caso da
simulação aqui, variou-se o parâmetro de 0,05 até 0,95 em acréscimos de 0,05, em um total de 19
valores.
O diagrama abaixo mostra os resultados, onde na maioria das simulações, o valor ficou em torno
de  
1
, o valor real. O valor médio das 20 simulações foi de 0,482.
2
30
Gráfico 2 – Valores de
g (; Y , X ) para diferentes valores de 
6
5.8
5.6
5.4
5.2
5
4.8
4.6
4.4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
De posse do valor de  , utilizamos as expressões acima para obter os valores dos demais
parâmetros. As médias dos valores das 20 replicações para  0 , 1 e  2 foram respectivamente
0,9970, 0,9954 e 1,0009, também muito próximos dos valores reais.
9 - Matriz de covariância
Nas seções 3 e 5 foram obtidas as matrizes de covariância respectivamente para a distribuição
normal e para o MQO. Devido à extensão dos cálculos para a obtenção dessa matriz para os
modelos espaciais, a obtenção da matriz de covariância desses modelos será feita apenas para o
modelo de lag espacial. As passagens são feitas passo-a-passo e, assim, procedimentos análogos
podem ser adaptados para os demais modelos.
Como vimos, o estimador de máxima verossimilhança tem com propriedade ser assintoticamente
eficiente com covariância assintótica positivamente definida no limite inferior de Cramér-Rao:
1
   2 ln L 
 .
Var (ˆ)  [ I ( )]1   E 




'



31
Para o modelo de lag espacial, essa matriz de covariância tem o seguinte formato:
   2 ln L
  2 2
  2  '
   ln L
[ I ( )]1   E 
 2 '
  2
   ln L
   2 '
 2 ln L
 2  '
 2 ln L
 '
 2 ln L
 '
 2 ln L 

 2  ' 
 2 ln L 

 ' 

 2 ln L 
 ' 
1
Então partimos do log da FV do modelo:
l


n
 1 
( y  Wy  X )' ( y  Wy  X ),
ln( 2 )  ln( 2 )   ln(1  i )   
2 
2
 2 
i
e obtemos as derivadas de primeira ordem. As duas primeiras já foram feitas anteriormente:
l
 1 
  2  2 X ' y  2 X ' X 

 2 
 1 
l
n
( y  Wy  X )' ( y  Wy  X ) 
  2  
2
2 2 

2
 2( ) 
A terceira é realizada aqui pela primeira vez:

  n
1 

( y  Wy  X )' ( y  Wy  X ) 
l
 ln( 2 )  ln( 2 )   ln(1  i )   

2 

  2
 2 
i


 
i

i
 1 
  2 (Wy)' ( y  Wy  X ) 
(1  i )   
Em seguida, obtemos as derivadas de segunda ordem de todas as expressões acima, lembrando
que a matriz é simétrica. Além disso, estimamos os valores esperados de cada uma das
expressões.

 1 
 2 ln L
  l 
 
n

( y  Wy  X )' ( y  Wy  X ) 


 
2
2
2 
2 
2 
2
2 2 
  '       2
 2( ) 

 2 
n
( y  Wy  X )' ( y  Wy  X ) 

 
2 2
2( )  2( 2 ) 3 
32
Obtemos o valor esperado dessa expressão:
 n

 n

 2 
 2 








 E

(
y


Wy

X

)'
(
y


Wy

X

)


E


'




2 2
2 2
 2( 2 ) 3 
 2( 2 ) 3 




 2( )

 2( )

 n

 n

 2 
 2 
n
 E  '     
 n 2  
 
 
 
2 2
2 3 
4
2 3 
4
 2( )  2( ) 

 2( )  2( ) 
 2( )



 1 
 2 ln L
  l   
n









(
y


Wy

X

)'
(
y


Wy

X

)

2
2 
2
2
2



        2
 2( ) 

 1 
 1 
  2 2 (Wy )' ( y  Wy  X )    2 2  ' Wy 
 ( ) 
 ( ) 
Para simplificar a notação na obtenção do valor esperado definimos a matriz B  W I  W  .
1
Utilizamos as seguintes expressões: E[ X '  ]  E[ X ' ]  0; para uma matriz 1 x 1, tr (a11)  a11;
e temos a seguinte relação para o traço de uma matriz, tr ( AB)  tr ( BA). O valor esperado da
derivada acima tem o seguinte valor:

 1 

 1 
 1 
1
 E  2 2  ' Wy    2 2  E  ' Wy    2 2  E  ' W I  W  ( X   )
 ( ) 
 ( ) 
 ( ) 


 1 
 1 
 1 
  2 2  E ' B ( X   ))   2 2  E( ' B  )   2 2  E tr ( ' B )
 ( ) 
 ( ) 
 ( ) 
 1 
 1 
 1 
  2 2 tr BE  '     2 2 tr BI 2   2 tr B 
 
 ( ) 
 ( ) 



 1 
 2 ln L
  l   
n










(
y


Wy

X

)'
(
y


Wy

X

)



 2  '    2    2 2  2( 2 ) 2 

 1 
 1 
  2 2 ( X )' ( y  Wy  X )    2 2  X '  
 ( ) 
 ( ) 
33
Obtemos o valor esperado dessa expressão:
  1 

 E   2 2  X '    0
  ( ) 

Utilizando a notação   
i
(i ) 2
para simplificar, obtemos a próxima expressão:
(1  i ) 2

i
 2 ln L   l   
 1 
  

  2 (Wy )' ( y  Wy  X ) 
 
 '       i (1  i )   

(i ) 2
 1 
 1 
   2 (Wy )' (Wy )       2 ( y 'W 'Wy 
2
  
  
i (1  i )
 1 
     2  (( I  W ) 1 ( X   ))'W ' ( B)( X   ))
  
 1 
     2 (  ' X ' ' ) B' B( X   ) 
  
 1 
     2 [(  ' X ' B' BX  2 ' B' BX   ' B' B ]
  



 1 
 E{    2 [  ' X ' B ' BX  2 ' B' BX   ' B ' B ]}
  
 1 
    2  E[  ' X ' B' BX  2 ' B ' BX   ' B' B ]
 
 1 
    2 [  ' X ' B' BX  E[2 ' B' BX ]  E[ ' B' B ]]
 
 1 
    2 [  ' X ' B' BX  0  E[tr  ' B' B ]]
 
 1 
    2 [  ' X ' B' BX  trB ' BE[ '  ]]
 
 1 
    2 [  ' X ' B' BX  trB ' B ( I 2 )]
 
 1 
    2 [  ' X ' B' BX ]  trB ' B
 
34
  1 
 2 ln L   l     1 

 

   2 ( X )' ( y  Wy  X )    2  X ' (Wy) 
 '         
  
 1 

 1 
 1 
 E  2  X ' Wy    2  EX ' Wy    2  EX ' B( X   ) 
 
 
  

 1 
  2  X ' BX
 
Finalmente:
  1 
 2 ln L   l     1 

 

   2 ( X )' ( y  Wy  X )    2  X ' X 
 '         
  
 1 

 X'X 
 E  2  X ' X    2 
  
  

Substituindo todas essas expressões, obtemos a matriz de covariância:


n
 1 
0


 2 tr B 
4
 
 2( )



1
1
1






[ I ( )]1    2 tr B      2 [(  ' X ' B' BX ]  trB' B  2  X ' BX 
 
 
  

 1 
 X'X  

0
 2  X ' BX
 2  

 
   

1
10 – Conclusão
Esse texto faz parte de uma série que apresenta a econometria espacial em pontos introdutórios.
No primeiro texto da serie foram discutidos alguns conceitos introdutórios sobre a econometria
espacial. No segundo texto, discutiu-se como interpretar os coeficientes obtidos nos diversos
modelos. Este terceiro texto da serie, discute como utilizar a função de verossimilhança para
estimar alguns dos modelos espaciais comumente utilizados em análises empíricas. Conceitos
básicos da função de verossimilhança foram discutidos, e os estimadores de máxima
35
verossimilhança para os modelos de lag espacial e de erro espacial foram derivados. A
metodologia de estimação de vários dos demais modelos, a exceção do modelo de KelejianPrucha, é similar ao desses dois primeiros. Simulações no Matlab exemplificam a aplicação dos
conceitos discutidos no texto.
Como vimos, o modelo de Kelejian-Prucha engloba os modelos de lag espacial e de erro
espacial. Assim, pode-se obter a metodologia de estimação desse primeiro diretamente a partir
do apresentado para esses dois últimos. Assim, para que este texto não ficasse demasiadamente
extenso, esse modelo não foi discutido em separado.
Porém, existe uma característica na estimação desse modelo que deve ser mencionada, uma vez
que ela não aparece nos demais modelos espaciais citados. Como vimos na estimação do modelo
de lag espacial, o valor de  é obtido por procedimentos numéricos. Ou seja, varia-se o valor
desse parâmetro em pequenos acréscimos com o uso do computador, e verifica-se qual valor
minimiza uma função em particular. De forma similar, o valor de  também é obtido por
procedimentos numéricos na estimação dos modelos de erro espacial. Nos dois casos temos
análises unidimensionais: apenas um parâmetro,  ou  , varia de cada vez e temos um vetor
com os valores de uma função especifica que está sendo minimizada. No modelo de KelejianPrucha devemos variar esses dois parâmetros de forma simultânea. Ou seja, em vez de um vetor
com o valor da função sendo minimizada, obtemos uma matriz, e verificamos para qual par de
valores de  e  que a função é minimizada. Assim, conceitualmente as estimações dos três
modelos são similares, mas computacionalmente estimar o modelo de Kelejian-Prucha é mais
demandante.
O objetivo principal deste texto foi fornecer uma base para que o leitor possa entender de forma
mais sólida como aplicar o estimador de máxima verossimilhança em estudos com modelos
espaciais. Assim, caso seja necessário a utilização de modelos espaciais distintos dos
normalmente utilizados, o próprio leitor pode desenvolver a sintaxe de estimação.
36
Referências
Anselin, L. (1988) Spatial Econometrics: Methods and Models. Kluwer Academic, Dordrecht.
Anselin, L., Bera, A., Florax, R. e Yoon, M. (1996) Simple diagnostic tests for spatial
dependence. Regional Science and Urban Economics 26 (1), 77-104.
Doreian, P. (1981) Estimating Linear Models with Spatially Distributed Data. Sociological
Methodology 12, 359-388.
Elhorst, J. (2010) Applied spatial econometrics: raising the bar. Spatial Economic Analysis 5 (1).
Florax, R., Folmer, H. e Rey, S. (2003) Specification searches in spatial econometrics: the
relevance of Hendry's methodology. Regional Science and Urban Economics 33 (5), 557-579.
Greene, W. (2003) Econometric Analysis. Prentice Hall, 5º edition.
Griffith, D. e Arbia, G. (2010) Detecting negative spatial autocorrelation in georeferenced
random variables. International Journal of Geographical Information Science 24 (3), 417-437
Kelejian, H. (2011) Spatial models in Econometrics. Draft for the Spatial Econometric Advanced
Institute.
LeSage, J e Pace, R. (2009) Introduction to spatial econometrics. Taylor & Francis Group, Boca
Raton.
Mur, J. e Ângulo, A. (2009) Model selection strategies in a spatial setting: some additional
results. Regional Science and Urban Economics 39, 200-213
Pace, K. (1997) Performing Large Spatial Regressions and Autoregressions. Economic Letters 54
(3), 283–291.
Simon, C. e Blume, L.(2004) Matemática para Economistas. Bookman, 1ª edição.
37
Apêndice econométrico 1
Esse apêndice discute os seguintes pontos estimadores consistentes, matriz de covariância,
matriz positivamente definida, estimadores eficientes e limite inferior de Cramér-Rao. Para uma
discussão mais abrangente sobre esses temas ver Greene (2003) e Simon e Blume (2004).
Estimador consistente
Uma variável aleatória xn
converge em probabilidade para uma constante 
se,
lim Pr ob xn       0 para todo   0. Ou seja, a probabilidade de se observar uma
n
estimação muito diferente de  é muito pequena, tendendo para zero, quando o tamanho da
amostra tende a infinito.


No caso de um estimador não enviesado, com E[ˆn ]   , temos, lim Pr ob ˆn      0, ou
n
p lim ˆn   , onde n é o tamanho da amostra. Um estimador assim é consistente.
Matriz de covariância
Dada uma distribuição multivariada de n variáveis, xi , temos como valor esperado de cada uma
delas o seguinte vetor:
 1   E[ x1 ] 
  

  2   E[ x2 ] 
  
 E[ x].
...
... 
  

    E[ x ] 
n 
 n 
A matriz de covariância é dada por:
38
 ( x1  1 )( x1  1 ) ( x1  1 )( x2   2 )

...
 ( x   2 )( x1  1 )
E[( x   )( x   )' ]   2
...
...

 ( x   )( x   )
...
n
1
1
 n
... ( x1  1 )( xn   n ) 

...
...

.
...
...

... ( xn   n )( xn   n ) 
Ou de forma mais sintética:
  11  12

...

Var ( x)   21
... ...


 n1 ...
...  1n 

... ... 
 .
... ... 

...  nn 
Matriz positivamente definida
Uma matriz A é positivamente definida se z' Az  0 para todo z  0 .
Estimador eficiente
Tome dois estimadores não enviesados ˆ1 e ˆ2 de  : E[ˆ1 ]   e E[ˆ2 ]   .
O estimador ˆ1 é mais eficiente que ˆ2 se a variância do primeiro for menor que a do segundo,
Var (ˆ1 )  Var (ˆ2 ).
Quando temos um vetor de parâmetros, o estimador ˆ1 é mais eficiente que ˆ2 se a matriz obtida
pela subtração das respectivas matrizes de covariâncias, Var (ˆ2 )  Var (ˆ1 ) , for positivamente
definida.
39
Limite inferior de Cramér-Rao
Assumindo certas condições formais, a variância de um estimador não enviesado será pelo
1
1
   ln L( )  2 
   2 ln L( ) 
   E 
menos tão grande quanto: [ I ( )]   E 
  .
    
   ' 
1
Assim, um estimador não enviesado terá uma variância igual ou maior que esse limite. Um
estimador com variação definida por esse limite é, portanto, eficiente, como é o caso dos
estimadores de MV.
Apêndice econométrico 2 – Função injetiva
Uma função diz-se injetiva (ou injetora) se, para quaisquer valores de x1 e x2 , tais que x1  x2 ,
pertencentes ao domínio da função, tivermos f ( x1 )  f ( x2 ).
Apêndice econométrico 3 - Estimador mínimos quadrados ordinários
Aqui são incluídos alguns pontos introdutórios sobre a estimação de modelos MQO via mínimos
quadrados ordinários, que justamente dão nome ao modelo, que foram julgados úteis para o
entendimento das estimações utilizando a função de verossimilhança. A apresentação aqui se
baseia em Greene (2003).
O modelo clássico de regressão linear múltipla tem alguns pressupostos. Seguem alguns deles.
1) Existe uma linearidade na relação entre a variável dependente e as independentes, y  X   ;
2) A matriz das variáveis explicativas tem posto máximo, isto é, essas variáveis são
independentes uma das outras, e, assim, não podem ser escritas como uma combinação linear das
demais;
3) As variáveis explicativas são exógenas, ou seja, E[ X ]  0;
40
4) Os erros têm como características a homocedasticidade e a inexistência de autocorrelação,
E[ ' X ]   2 I ; e
5) Os erros são distribuídos de forma normal,  X ~ N (0,  2 I ) .
Nos modelos MQO estimasse os parâmetros minimizando a soma dos quadrados dos resíduos.
Para obtermos o estimador para os parâmetros  do modelo, partimos do pressuposto 1):
y  X   .
Isolamos o erro e manipulamos a expressão:
  y  X
 '   y  X '  y  X .
Lembrando que ( A  B)'  A' B' e que ( AB)'  B' A' , obtemos a seguinte expressão:
S ( )   '   y  X '  y  X    y'( X )' y  X   y' y  ( X )' y  y' ( X )  ( X )' ( X )
Lembrando-se da relação para vetores, X 'Y  Y ' X , temos:
S ( )  y' y  2( X )' y   ' X ' X  y' y  2 ' X ' y   ' X ' X
Derivando com relação a  , obtemos a condição de primeira ordem para minimizar a expressão:
S (  ) ( y' y  2 ' X ' y   ' X ' X )

 2 X ' y  2 X ' X  0


Daí, temos:
X ' X  X ' y
Como X tem posto máximo, como descrito no segundo pressuposto, a matriz X ' X tem inversa,
e obtemos o estimador de mínimos quadrados ordinários:
ˆMQO  ( X ' X ) 1 X ' y .
41
Apêndice econométrico 4 – O Jacobiano
Como vimos, o Jacobiano é incluído nas funções de verossimilhança. Esse apêndice
econométrico apresenta alguns conceitos referentes a este tema, inicialmente pela definição e, em
seguida, com a utilização na FV. A apresentação aqui se baseia em Greene (2003).
A4.1 - Definição de Jacobiano de uma transformação
A matriz Jacobiana é a matriz formada pelas derivadas parciais de primeira ordem de uma
função vetorial. O valor absoluto do determinante dessa matriz, como veremos, é utilizado
quando são feitas mudanças de função que define uma densidade de probabilidade conjunta,
como a função de verossimilhança.
Seguem exemplos. Partimos de uma função injetiva e de sua inversa:
f ( x)  y
x  f 1 ( y)
Derivando essa função inversa, temos:


dx d

f 1 ( y ) .
dy dy
Se f (x) é linear, por exemplo, y  ax  b , temos:
x
y b

a a
J
dx 1

dy a
A função acima foi definida como f :    . No caso de n funções e n variáveis, podemos
escrever o Jacobiano como:
42
 x1

y
x  1
J
  ...
y '  xn
 y
 1
x1 

yn 
... ... 
xn 
...
yn 

...
Segue um exemplo para um sistema de equações lineares n x n, onde podemos escrever o
sistema como:
Y  AX
Se a inversa da matriz existir, temos:
X  A1Y , onde J  A1
Exemplificando com um sistema 2 x 2, ficamos com:
 y1   a11 a12  x1 
   
 
 y 2   a21 a22  x2 
1
 x1   a11 a12   y1 
   
  
 x2   a21 a22   y2 
x  a11 a12 

J

y'  a21 a22 
1
O valor absoluto do determinante da matriz acima é o Jacobiano da transformação: abs J .
Uma vez definido o Jacobiano da transformação, segue uma explicação do por que do uso deste
na FV.
43
A4.2 – O uso do Jacobiano na função de verossimilhança
Tome uma densidade de probabilidade contínua, definida por f x (x) . Assim, temos:
b
Pr ob( x  b) 
f
( x)dx.
x

Seja y  g (x) uma transformação monotônica de x e, portanto, temos x  g 1 ( y) .
Qual é a densidade de probabilidade de y? Em outras palavras, qual é a função f y ( y ) , tal que:
a
Pr ob( y  a) 
f
y
( y)dy ?

Note que: f x ( x)  f x ( g 1 ( y)) e que
dx
 [ g 1 ( y )]'
dy
Substituindo esses termos na integral acima com g(a) = b, e, ainda, tomando apenas valores
positivos para não termos pdf com valores negativos, temos:
a
Pr ob( y  a) 
f
x
( g 1 ( y )) [ g 1 ( y )]' dy

Assim, f y ( y)  f x ( g 1 ( y)) [ g 1 ( y)]'  f x ( g 1 ( y)) J
Ou seja, no caso da densidade conjunta, temos: Ly ( ; y)  Lx ( g 1 ( ; y)) J
44
Download

How to estimate spatial models with the use of the likelihood