Solved – Maximum likelihood estimation of dlmModReg

dlmmaximum likelihoodrregression

I'm studying R package dlm. So far it seems very powerful and flexible package, with nice programming interfaces and good documentation.

I've been able to successfully use dlmMLE and dlmModARMA to estimate the parameters of AR(1) process:

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

Now I'm trying to use similar code to estimate the parameters of simple linear regression model:

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

I expect fit$par to be close to c(-1, 0.5), but I keep getting something like

[1] -0.0002118851  0.4884367070

The coefficient -1 is not estimated correctly. However, the strange thing is that the variance of the noise is returned correctly.

I understand that max-likelihood estimation might fail given bad initial values, but I observed that the likelihood function returned by dlmLL is very flat in the first coordinate.

So I wonder: can such model be estimated at all using dlm? I believe the model is "non-singular", however I'm not sure how the likelihood function is calculated inside the dlm.

Any hint greatly appreciated.

Best Answer

I think your setup is not correct. Try this:

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

You recover the estimate of the observation variance from the only element of 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"

while your estimate of the coefficient (should be around -1 in your case) can be obtained as the last element of dlmFilter(u,mod)$a, which gives the values of the state as new observations are processed:

    > 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

Hope this helps.

Related Question