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]]