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]] ______________________________________________ R-help@r-project.org 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.
