Estimación por máxima verosimilitud de dlmModReg

Question

Estoy estudiando el paquete dlm de R. Hasta ahora parece un paquete muy potente y flexible, con buenas interfaces de programación y buena documentación.

He podido utilizar con éxito dlmMLE y dlmModARMA para estimar los parámetros del proceso AR(1):

u <- arima.sim(list(ar = 0.3), 100)
fit <- dlmMLE(u, parm = c(0.5, sd(u)),
              build = function(x)
                dlmModARMA(ar = x[1], sigma2 = x[2]^2))
fit$par

Ahora estoy intentando utilizar un código similar para estimar los parámetros de un modelo de regresión lineal simple:

r <- rnorm(100)
u <- -1*r + 0.5*rnorm(100)
fit <- dlmMLE(u, parm = c(0, 1),
              build = function(x)
                dlmModReg(x[1]*r, FALSE, dV = x[2]^2))
fit$par

Espero que fit$par esté cerca de c(-1, 0.5), pero sigo obteniendo algo como

[1] -0.0002118851  0.4884367070

El coeficiente -1 no se estima correctamente. Sin embargo, lo extraño es que la varianza del ruido se devuelve correctamente.

Entiendo que la estimación de máxima verosimilitud puede fallar dados los malos valores iniciales, pero he observado que la función de verosimilitud devuelta por dlmLL es muy plana en la primera coordenada.

Así que me pregunto: ¿se puede estimar este modelo con dlm? Creo que el modelo es "no singular", sin embargo no estoy seguro de cómo se calcula la función de verosimilitud dentro del dlm.

Se agradece cualquier sugerencia.

Answer 1

Creo que tu configuración no es correcta. Intenta esto:

set.seed(1234)
r <- rnorm(100)
X <- r
u <- -1*X + 0.5*rnorm(100)
MyModel <- function(x)  dlmModReg(X, FALSE, dV = x[1]^2)
fit <- dlmMLE(u, parm = c(0.3), build = MyModel)
mod <- MyModel(fit$par)
dlmFilter(u,mod)$a

Se recupera la estimación de la varianza de la observación a partir del único elemento de fit$par:

> fit
$par
[1] 0.4431803

$value
[1] -20.69313

$counts
function gradient 
      17       17 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

mientras que su estimación del coeficiente (debería ser de alrededor de -1 en su caso) puede obtenerse como el último elemento de dlmFilter(u,mod)$a que da los valores del estado a medida que se procesan nuevas observaciones:

    > dlmFilter(u,mod)$m
  [1]  0.0000000 -1.1486921 -1.2123431 -1.1172783 -1.1231454 -1.1170222
  [7] -1.0974931 -1.1377114 -1.0378758 -1.0927136 -1.0955372 -1.0120210
 [13] -0.9874791 -1.0036429 -1.0765513 -1.0678725 -1.0795124 -1.1568597
 [19] -1.2044821 -1.2056687 -1.2102896 -1.2938958 -1.2922945 -1.2670604
 [25] -1.1789594 -1.1570172 -1.1601590 -1.1417200 -1.1585501 -1.1608675
 [31] -1.1616278 -1.1744861 -1.1717561 -1.1715025 -1.1568086 -1.1451311
 [37] -1.1520867 -1.1379211 -1.1270897 -1.1048035 -1.1015793 -1.1054597
 [43] -1.0621750 -1.0621218 -1.0696813 -1.0807651 -1.0816893 -1.0647963
 [49] -1.0643440 -1.0667282 -1.0626404 -1.0623697 -1.0586265 -1.0571205
 [55] -1.0569135 -1.0579224 -1.0607623 -1.0582257 -1.0495232 -1.0494288
 [61] -1.0539632 -1.0555427 -1.0553468 -1.0491239 -1.0488604 -1.0491036
 [67] -1.0510551 -1.0576294 -1.0611296 -1.0628612 -1.0626451 -1.0573650
 [73] -1.0629577 -1.0647724 -1.0658052 -1.0823839 -1.0753808 -1.0747229
 [79] -1.0747762 -1.0615243 -1.0630352 -1.0697431 -1.0666448 -1.0617227
 [85] -1.0585460 -1.0583981 -1.0563544 -1.0567715 -1.0544349 -1.0573228
 [91] -1.0588404 -1.0639155 -1.0625845 -1.0578004 -1.0571034 -1.0602645
 [97] -1.0604838 -1.0586019 -1.0580891 -1.0587096 -1.0577559  

Espero que esto ayude.

Answer 2

A continuación se muestra el código que implementa mi solución y la solución de Paramonov (una ligera edición: he cambiado dlmFilter(u,mod)$a en la respuesta publicada originalmente por dlmFilter(u,mod)$m ).

library(dlm)
set.seed(1234)
reps      <- 100
MyEstimates <- YourEstimates <- matrix(0,reps,2)
for (i in (1:reps) ) {
X <- r <- rnorm(100)
u <- -1*r + 0.5*rnorm(100)
#
fit <- dlmMLE(u, parm = c(1, sd(u)),
              build = function(x)
                dlmModReg(r, FALSE, dV = x[2]^2,
                          m0 = x[1], C0 = matrix(0)))
YourEstimates[i,] <- fit$par
#
MyModel <- function(x)  dlmModReg(X, FALSE, dV = x[1]^2)
fit <- dlmMLE(u, parm = c(0.3), build = MyModel)
mod <- MyModel(fit$par)
MyEstimates[i,] <- c(dlmFilter(u,mod)$m[101],fit$par[1])
}

Cuando ejecuto el código anterior, esto es lo que obtengo:

> summary(YourEstimates)
       V1                V2         
 Min.   :-9.5284   Min.   :-0.5747  
 1st Qu.:-1.4280   1st Qu.: 0.4710  
 Median :-0.9795   Median : 0.4937  
 Mean   :-0.9737   Mean   : 0.4369  
 3rd Qu.:-0.5636   3rd Qu.: 0.5215  
 Max.   : 4.5222   Max.   : 0.5980  
> summary(MyEstimates)
       V1                V2         
 Min.   :-1.1099   Min.   :-0.6010  
 1st Qu.:-1.0266   1st Qu.: 0.4736  
 Median :-0.9974   Median : 0.4961  
 Mean   :-0.9938   Mean   : 0.4469  
 3rd Qu.:-0.9635   3rd Qu.: 0.5158  
 Max.   :-0.8390   Max.   : 0.5776  

Mientras que el primer conjunto de estimaciones da estimaciones similares para el segundo parámetro, en ocasiones da valores muy alejados de la realidad para el primero. Creo que la razón es que "atar" el estado a su valor inicial con

C0=matrix(0)

conduce a la inestabilidad numérica, pero no estoy seguro. En cualquier caso, es posible que desee examinar la cuestión.

Answer 3

Después de leer la ayuda de dlmFilter, he podido llegar al siguiente código:

r <- rnorm(100)
u <- -1*r + 0.5*rnorm(100)
fit <- dlmMLE(u, parm = c(1, sd(u)),
              build = function(x)
                dlmModReg(r, FALSE, dV = x[2]^2,
                          m0 = x[1], C0 = matrix(0)))
fit$par

[1] -1.1330088  0.4788357