Reproducing the results from package dlm

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Reproducing the results from package dlm

Ashim Kapoor
Dear all,

I have a query with regards to the package dlm. My query is : will dlm
return the same results if I give it the same data set ?

Here is a MWE ( created from :
https://sites.ualberta.ca/~sfossati/e509/files/other/dlm_ex.R )

library(dlm)

# simulate AR(1) process
# phi = .8, sig2 = .25
nobs <- 250
yt <- arima.sim(n=nobs,list(ar=.8,ma=0),sd=.5)

# estimate AR(1) for comparison
model10 <- Arima(yt,order=c(1,0,0),method="ML",include.mean=FALSE)
model10

# set parameter restrictions
parm_rest <- function(parm){
     return( c(exp(parm[1])/(1+exp(parm[1])),exp(parm[2])) )
    }

# set up SS model
ssm1 <- function(parm){
    parm <- parm_rest(parm)
    return( dlm(FF=1,V=0,GG=parm[1],W=parm[2],
                m0=0,C0=solve(1-parm[1]^2)*parm[2]) )
    }
# estimate parameters
fit1 <- dlmMLE(y=yt,parm=c(0,1),build=ssm1,hessian=T)

# get estimates
coef <- parm_rest(fit1$par)
# get standard errors using delta method
dg1 <- exp(fit1$par[1])/(1+exp(fit1$par[1]))^2
dg2 <- exp(fit1$par[2])
dg <- diag(c(dg1,dg2))
var <- dg%*%solve(fit1$hessian)%*%dg
# print results
coef; sqrt(diag(var))

# Here are 2 runs of the latter part of the code: -

> fit1 <- dlmMLE(y=yt,parm=c(0,1),build=ssm1,hessian=T)
>
> # get estimates
> coef <- parm_rest(fit1$par)
> # get standard errors using delta method
> dg1 <- exp(fit1$par[1])/(1+exp(fit1$par[1]))^2
> dg2 <- exp(fit1$par[2])
> dg <- diag(c(dg1,dg2))
> var <- dg%*%solve(fit1$hessian)%*%dg
> # print results
> coef; sqrt(diag(var))
[1] 0.8442885 0.2513462
[1] 0.03361245 0.02248196

> fit1 <- dlmMLE(y=yt,parm=c(0,1),build=ssm1,hessian=T)
>
> # get estimates
> coef <- parm_rest(fit1$par)
> # get standard errors using delta method
> dg1 <- exp(fit1$par[1])/(1+exp(fit1$par[1]))^2
> dg2 <- exp(fit1$par[2])
> dg <- diag(c(dg1,dg2))
> var <- dg%*%solve(fit1$hessian)%*%dg
> # print results
> coef; sqrt(diag(var))
[1] 0.8442885 0.2513462
[1] 0.03361245 0.02248196
>

Seems to me that even without set.seed, dlm reports the SAME results. Is
that true? I have written a big program where I have slightly different
results from dlm and I wonder if I have made a mistake. I have created a
MWE above to verify my doubt.

Do we need a set.seed for reproducing results from dlm ? Please clarify.

Best Regards,
Ashim

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.