Maximum likelihood estimation of dlmModReg

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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