[R] Lower bound and upper bound in maximum likelihood

[Siow Yun Tan]

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