# 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. - add comment ## 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. – Andrey Paramonov May 12 '11 at 19:24 add comment 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). –  F. Tusell May 13 '11 at 15:09

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. –  Andrey Paramonov May 15 '11 at 14:14