# Maximum likelihood estimation of dlmModReg

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. ## 3 Answers 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. • Thanks for your quick reply! Although it didn't contain the solution, it did provide the needed insight. May 12, 2011 at 19:24 Below is code which implements my solution and Paramonov's solution (a slight edit: I have changed dlmFilter(u,mod)$a in the orginally posted answer by 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]) }  When I run the above code, this is what I get: > 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  While the first set of estimates gives similar estimates for the second parameter, it occasionally gives values well off the mark for the first. I think the reason is that "tying" the state to its initial value with C0=matrix(0)  leads to numerical instability, but I am not sure. In any case, you may want to look at the issue. • I can reproduce the problem. The reason seems really to be numerical instability, as changing C0=matrix(0) to C0=matrix(1E-12) allowed me to get results very close to yours. May 15, 2011 at 14:14 After reading help for dlmFilter, I could come up with the following code: 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

• I see what the intent is and it seems corect to me; however, look in my next answer what happens (I enter a new answer as I cannot write code here). May 13, 2011 at 15:09