Dear all,
I need to simulate data which fit to a double poisson time series model
with certain parameters. Then, check whether the estimated parameter close
to the true parameter by using maximum likelihood estimation.
Simulation:
set.seed(10)
library("rmutil")
a0 = 1.5; a1 = 0.4; b1 = 0.3; g1= 0.7 ; n=100
#a0, a1 and b1 are parameter where n is size.
nu = h = rep(0, n)
h[1] = a0/(1-a1-b1)
nu[1] = rdoublepois(1,h[1],g1)
for (i in 2:n) {
h[i] = a0 + a1 * nu[i-1] + b1 * h[i-1]
nu[i] = rdoublepois(1,h[i],g1)
}
Maximum likelihood by double poisson density function from rmutil:
logl3 <- function(par) {
h.new <- vector()
a0 <- par[1]
a1 <- par[2]
b1 <- par[3]
g1 <- par[4]
for (i in 2:n) {
h.new[i] = a0 + a1 * nu[i-1] + b1 * h.new[i-1]
}
-sum(ddoublepois(nu, m=h.new[2:n], s=g1,log=TRUE))
}
nlminb(start = c(0.01,0.01,0.01,0.01), lower = 1e-3, upper = Inf,
logl3)$par
But there is error.
>Error in if (any(m <= 0)) stop("m must be positive") :
>missing value where TRUE/FALSE needed
So I further check the maximum likelihood stop at which iteration:
logl <- function(par,debug=FALSE) {
h.new <- vector()
a0 <- par[1]
a1 <- par[2]
b1 <- par[3]
g1 <- par[4]
h.new[1] <- a0/(1-a1-b1)
for (i in 2:n) {
h.new[i] <- a0 + a1 * nu[i-1] + b1 * h.new[i-1]
}
if(debug) cat(a0,a1,b1,g1,"")
r <- -sum(ddoublepois(nu, m=h.new[1:n], s=g1,log=TRUE))
if(debug) cat(r,"\n")
return(r)
}
nlminb(start = c(0.1,0.1,0.1,0.1), lower = -Inf,
upper = Inf, logl,debug=TRUE)
The results:
>0.1 0.1 0.1 0.1 1517.164
>0.1 0.1 0.1 0.1 1517.164
>0.1 0.1 0.1 0.1 1517.164
>0.1 0.1 0.1 0.1 1517.164
>0.1 0.1 0.1 0.1 1517.164
>0.218201 0.6245722 0.1792584 -0.739387
>Error in ddoublepois(nu, m = h.new[1:n], s = g1, log = TRUE) :
>s must be positive
I am guessing that there is problem at the part of initial value, lower
bound and upper bound of my optimization function.
Please give me some guides.
Thanks in advance.
[[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.
