R/ModelB2.R

In mfrdixon/MLEMVD: Maximum Likelihood Estimation for Diffusion Equations

#' Compute the maximum likelihood estimate of Model B2
#'
#'
#' @param x Observation of the state variable at time t
#' @param x0 Observation of the state variable at time t-1
#' @param del The time step between the current and previous observation
#' @param param The parameter 12-vector (a0,a1,a2,b0,b1,b2,c0,c1,c2,d0,d1,d2)
#' @export output a list with a llk variable storing the result of the log likelihood calculation
#' @examples
#' ModelB2(0.4,0.3,0.1,c(0.1,0.2,0.3,0.1,0.2,0.3,0.1,0.2,0.3,0.1,0.2,0.3))
#'
ModelB2 <- function (x,x0,del,param)
{

  x1 <- x[1]
  x2 <- x[2]
  x10 <- x0[1]
  x20 <- x0[2]

  a0 <- param[1]
  a1 <- param[2]
  a2 <- param[3]
  b0 <- param[4]
  b1 <- param[5]
  b2 <- param[4]
  c0 <- param[7]
  c1 <- param[8]
  c2 <- param[9]
  d0 <- param[10]
  d1 <- param[11]
  d2 <- param[12]

  Dv <- log(c0+c1 * x1+c2* x2) + log(d0+d1*x1+d2*x2)
  cmone <- (c1*(x1 - x10)^3)/(2*(c0 + c1*x10 + c2*x20)^3) + (c2*(x1 - x10)^2*(x2 - x20))/(2*(c0 + c1*x10 + c2*x20)^3) + (d1*(x1 - x10)*(x2 - x20)^2)/(2*(d0 + d1*x10 + d2*x20)^3) + (d2*(x2 - x20)^3)/(2*(d0 + d1*x10 + d2*x20)^3) + (1/144)*(x2 - x20)^4*((6*d1^2*(c0 + c1*x10 + c2*x20)^2)/(d0 + d1*x10 + d2*x20)^6 - (66*d2^2)/(d0 + d1*x10 + d2*x20)^4) + (1/2)*(-((x1 - x10)^2/(c0 + c1*x10 + c2*x20)^2) - (x2 - x20)^2/(d0 + d1*x10 + d2*x20)^2) + (1/144)*(x1 - x10)^4*(-((66*c1^2)/(c0 + c1*x10 + c2*x20)^4) + (6*c2^2*(d0 + d1*x10 + d2*x20)^2)/(c0 + c1*x10 + c2*x20)^6) - (c2*(x1 - x10)^3*(x2 - x20)*(d1*(c0 + c2*x20) + c1*(5*d0 + 6*d1*x10 + 5*d2*x20)))/(6*(c0 + c1*x10 + c2*x20)^4*(d0 + d1*x10 + d2*x20)) - (d1*(x1 - x10)*(x2 - x20)^3*(5*d2*(c0 + c1*x10) + c2*(d0 + d1*x10 + 6*d2*x20)))/(6*(c0 + c1*x10 + c2*x20)*(d0 + d1*x10 + d2*x20)^4) + (1/24)*(x1 - x10)^2*(x2 - x20)^2*(-((8*c2^2)/(c0 + c1*x10 + c2*x20)^4) - (2*c2*d2)/((c0 + c1*x10 + c2*x20)^3*(d0 + d1*x10 + d2*x20)) + (2*d1*(-4*d1 - (c1*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)))/(d0 + d1*x10 + d2*x20)^4)
  czero <- (x1 - x10)*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))) + (x2 - x20)*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))) + (1/2)*(x1 - x10)*(x2 - x20)*(-((c2*(a0 + a1*x10 + a2*x20))/(c0 + c1*x10 + c2*x20)^3) + a2/(c0 + c1*x10 + c2*x20)^2 + (c1*c2)/(c0 + c1*x10 + c2*x20)^2 - (d1*(b0 + b1*x10 + b2*x20))/(d0 + d1*x10 + d2*x20)^3 + b1/(d0 + d1*x10 + d2*x20)^2 + (d1*d2)/(d0 + d1*x10 + d2*x20)^2 - (c2*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))))/(c0 + c1*x10 + c2*x20) + (c2*(c1/(c0 + c1*x10 + c2*x20) + d1/(d0 + d1*x10 + d2*x20)))/(c0 + c1*x10 + c2*x20) - (d1*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))))/(d0 + d1*x10 + d2*x20) + (d1*(c2/(c0 + c1*x10 + c2*x20) + d2/(d0 + d1*x10 + d2*x20)))/(d0 + d1*x10 + d2*x20) - (c2*(d1*(c0 + c2*x20) + c1*(5*d0 + 6*d1*x10 + 5*d2*x20)))/ (2*(c0 + c1*x10 + c2*x20)^2*(d0 + d1*x10 + d2*x20)) - (d1*(5*d2*(c0 + c1*x10) + c2*(d0 + d1*x10 + 6*d2*x20)))/(2*(c0 + c1*x10 + c2*x20)*(d0 + d1*x10 + d2*x20)^2)) + (1/4)*(x2 - x20)^2*(c2^2/(c0 + c1*x10 + c2*x20)^2 - (d1*(a0 + a1*x10 + a2*x20))/(d0 + d1*x10 + d2*x20)^3 - (3*d2*(b0 + b1*x10 + b2*x20))/(d0 + d1*x10 + d2*x20)^3 + (2*c1*d1*(c0 + c1*x10 + c2*x20))/(d0 + d1*x10 + d2*x20)^3 + (2*b2)/(d0 + d1*x10 + d2*x20)^2 + (7*d2^2)/(d0 + d1*x10 + d2*x20)^2 + (1/12)*(d0 + d1*x10 + d2*x20)^2*((6*d1^2*(c0 + c1*x10 + c2*x20)^2)/(d0 + d1*x10 + d2*x20)^6 - (66*d2^2)/(d0 + d1*x10 + d2*x20)^4) + (d1*(c0 + c1*x10 + c2*x20)^2*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))))/(d0 + d1*x10 + d2*x20)^3 - (d2*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))))/(d0 + d1*x10 + d2*x20) - (3*d2*(c2/(c0 + c1*x10 + c2*x20) + d2/(d0 + d1*x10 + d2*x20)))/(d0 + d1*x10 + d2*x20) - (d1*(c0 + c1*x10 + c2*x20)*(d1*(c0 + c2*x20) + c1*(d0 + 2*d1*x10 + d2*x20)))/ (d0 + d1*x10 + d2*x20)^4 + (2*(d2^2*(c0 + c1*x10)^2 + 2*c2*d2*(c0 + c1*x10)*(d0 + d1*x10 + 2*d2*x20) - c2^2*(d0^2 + 2*d0*d1*x10 + d1^2*x10^2 - 2*d2^2*x20^2)))/ ((c0 + c1*x10 + c2*x20)^2*(d0 + d1*x10 + d2*x20)^2) + (1/12)*(c0 + c1*x10 + c2*x20)^2*(-((8*c2^2)/(c0 + c1*x10 + c2*x20)^4) - (2*c2*d2)/((c0 + c1*x10 + c2*x20)^3*(d0 + d1*x10 + d2*x20)) + (2*d1*(-4*d1 - (c1*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)))/(d0 + d1*x10 + d2*x20)^4)) + (1/4)*(x1 - x10)^2*(-((3*c1*(a0 + a1*x10 + a2*x20))/(c0 + c1*x10 + c2*x20)^3) - (c2*(b0 + b1*x10 + b2*x20))/(c0 + c1*x10 + c2*x20)^3 + (2*a1)/(c0 + c1*x10 + c2*x20)^2 + (7*c1^2)/(c0 + c1*x10 + c2*x20)^2 + d1^2/(d0 + d1*x10 + d2*x20)^2 + (2*c2*d2*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)^3 - (c1*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))))/(c0 + c1*x10 + c2*x20) - (3*c1*(c1/(c0 + c1*x10 + c2*x20) + d1/(d0 + d1*x10 + d2*x20)))/(c0 + c1*x10 + c2*x20) + (c2*(d0 + d1*x10 + d2*x20)^2*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))))/(c0 + c1*x10 + c2*x20)^3 + (1/12)*(c0 + c1*x10 + c2*x20)^2*(-((66*c1^2)/(c0 + c1*x10 + c2*x20)^4) + (6*c2^2*(d0 + d1*x10 + d2*x20)^2)/(c0 + c1*x10 + c2*x20)^6) - (c2*(d0 + d1*x10 + d2*x20)*(d2*(c0 + c1*x10) + c2*(d0 + d1*x10 + 2*d2*x20)))/(c0 + c1*x10 + c2*x20)^4 + (2*((-d1^2)*(c0 + c2*x20)^2 + 2*c1*d1*(c0 + c2*x20)*(d0 + d2*x20) + c1^2*(d0^2 + 4*d0*d1*x10 + 2*d1^2*x10^2 + 2*d0*d2*x20 + 4*d1*d2*x10*x20 + d2^2*x20^2)))/ ((c0 + c1*x10 + c2*x20)^2*(d0 + d1*x10 + d2*x20)^2) + (1/12)*(d0 + d1*x10 + d2*x20)^2*(-((8*c2^2)/(c0 + c1*x10 + c2*x20)^4) - (2*c2*d2)/((c0 + c1*x10 + c2*x20)^3*(d0 + d1*x10 + d2*x20)) + (2*d1*(-4*d1 - (c1*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)))/(d0 + d1*x10 + d2*x20)^4))
  cone <- (1/4)*(-4*a1 - 4*b2 + 4*c1^2 + 4*d2^2 - 2*(c0 + c1*x10 + c2*x20)^2*(-(c1^2/(c0 + c1*x10 + c2*x20)^2) - d1^2/(d0 + d1*x10 + d2*x20)^2) - 2*(d0 + d1*x10 + d2*x20)^2*(-(c2^2/(c0 + c1*x10 + c2*x20)^2) - d2^2/(d0 + d1*x10 + d2*x20)^2) - 4*(a0 + a1*x10 + a2*x20)*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))) + 8*c1*(c0 + c1*x10 + c2*x20)*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))) + 2*(c0 + c1*x10 + c2*x20)^2*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20)))^2 + 4*(a0 + a1*x10 + a2*x20)*(c1/(c0 + c1*x10 + c2*x20) + d1/(d0 + d1*x10 + d2*x20)) - 4*(c0 + c1*x10 + c2*x20)^2*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20)))*(c1/(c0 + c1*x10 + c2*x20) + d1/(d0 + d1*x10 + d2*x20)) - 4*(b0 + b1*x10 + b2*x20)*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))) + 8*d2*(d0 + d1*x10 + d2*x20)*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))) + 2*(d0 + d1*x10 + d2*x20)^2*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20)))^2 + 4*(b0 + b1*x10 + b2*x20)*(c2/(c0 + c1*x10 + c2*x20) + d2/(d0 + d1*x10 + d2*x20)) - 4*(d0 + d1*x10 + d2*x20)^2* (c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20)))*(c2/(c0 + c1*x10 + c2*x20) + d2/(d0 + d1*x10 + d2*x20)) - (8*c1*(d1*(c0 + c2*x20) + c1*(d0 + 2*d1*x10 + d2*x20)))/(d0 + d1*x10 + d2*x20) + (2*(d1*(c0 + c2*x20) + c1*(d0 + 2*d1*x10 + d2*x20))^2)/(d0 + d1*x10 + d2*x20)^2 - (8*d2*(d2*(c0 + c1*x10) + c2*(d0 + d1*x10 + 2*d2*x20)))/(c0 + c1*x10 + c2*x20) + (2*(d2*(c0 + c1*x10) + c2*(d0 + d1*x10 + 2*d2*x20))^2)/(c0 + c1*x10 + c2*x20)^2 + (d0 + d1*x10 + d2*x20)^2*(c2^2/(c0 + c1*x10 + c2*x20)^2 - (d1*(a0 + a1*x10 + a2*x20))/(d0 + d1*x10 + d2*x20)^3 - (3*d2*(b0 + b1*x10 + b2*x20))/(d0 + d1*x10 + d2*x20)^3 + (2*c1*d1*(c0 + c1*x10 + c2*x20))/(d0 + d1*x10 + d2*x20)^3 + (2*b2)/(d0 + d1*x10 + d2*x20)^2 + (7*d2^2)/(d0 + d1*x10 + d2*x20)^2 + (1/12)*(d0 + d1*x10 + d2*x20)^2*((6*d1^2*(c0 + c1*x10 + c2*x20)^2)/(d0 + d1*x10 + d2*x20)^6 - (66*d2^2)/(d0 + d1*x10 + d2*x20)^4) + (d1*(c0 + c1*x10 + c2*x20)^2*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))))/(d0 + d1*x10 + d2*x20)^3 - (d2*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))))/(d0 + d1*x10 + d2*x20) - (3*d2*(c2/(c0 + c1*x10 + c2*x20) + d2/(d0 + d1*x10 + d2*x20)))/(d0 + d1*x10 + d2*x20) - (d1*(c0 + c1*x10 + c2*x20)*(d1*(c0 + c2*x20) + c1*(d0 + 2*d1*x10 + d2*x20)))/ (d0 + d1*x10 + d2*x20)^4 + (2*(d2^2*(c0 + c1*x10)^2 + 2*c2*d2*(c0 + c1*x10)*(d0 + d1*x10 + 2*d2*x20) - c2^2*(d0^2 + 2*d0*d1*x10 + d1^2*x10^2 - 2*d2^2*x20^2)))/ ((c0 + c1*x10 + c2*x20)^2*(d0 + d1*x10 + d2*x20)^2) + (1/12)*(c0 + c1*x10 + c2*x20)^2*(-((8*c2^2)/(c0 + c1*x10 + c2*x20)^4) - (2*c2*d2)/((c0 + c1*x10 + c2*x20)^3*(d0 + d1*x10 + d2*x20)) + (2*d1*(-4*d1 - (c1*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)))/(d0 + d1*x10 + d2*x20)^4)) + (c0 + c1*x10 + c2*x20)^2*(-((3*c1*(a0 + a1*x10 + a2*x20))/(c0 + c1*x10 + c2*x20)^3) - (c2*(b0 + b1*x10 + b2*x20))/(c0 + c1*x10 + c2*x20)^3 + (2*a1)/(c0 + c1*x10 + c2*x20)^2 + (7*c1^2)/(c0 + c1*x10 + c2*x20)^2 + d1^2/(d0 + d1*x10 + d2*x20)^2 + (2*c2*d2*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)^3 - (c1*((a0 + a1*x10 + a2*x20)/(c0 + c1*x10 + c2*x20)^2 - c1/(2*(c0 + c1*x10 + c2*x20)) + d1/(2*(d0 + d1*x10 + d2*x20))))/(c0 + c1*x10 + c2*x20) - (3*c1*(c1/(c0 + c1*x10 + c2*x20) + d1/(d0 + d1*x10 + d2*x20)))/(c0 + c1*x10 + c2*x20) + (c2*(d0 + d1*x10 + d2*x20)^2*(c2/(2*(c0 + c1*x10 + c2*x20)) + (b0 + b1*x10 + b2*x20)/(d0 + d1*x10 + d2*x20)^2 - d2/(2*(d0 + d1*x10 + d2*x20))))/(c0 + c1*x10 + c2*x20)^3 + (1/12)*(c0 + c1*x10 + c2*x20)^2*(-((66*c1^2)/(c0 + c1*x10 + c2*x20)^4) + (6*c2^2*(d0 + d1*x10 + d2*x20)^2)/(c0 + c1*x10 + c2*x20)^6) - (c2*(d0 + d1*x10 + d2*x20)*(d2*(c0 + c1*x10) + c2*(d0 + d1*x10 + 2*d2*x20)))/(c0 + c1*x10 + c2*x20)^4 + (2*((-d1^2)*(c0 + c2*x20)^2 + 2*c1*d1*(c0 + c2*x20)*(d0 + d2*x20) + c1^2*(d0^2 + 4*d0*d1*x10 + 2*d1^2*x10^2 + 2*d0*d2*x20 + 4*d1*d2*x10*x20 + d2^2*x20^2)))/ ((c0 + c1*x10 + c2*x20)^2*(d0 + d1*x10 + d2*x20)^2) + (1/12)*(d0 + d1*x10 + d2*x20)^2*(-((8*c2^2)/(c0 + c1*x10 + c2*x20)^4) - (2*c2*d2)/((c0 + c1*x10 + c2*x20)^3*(d0 + d1*x10 + d2*x20)) + (2*d1*(-4*d1 - (c1*(d0 + d1*x10 + d2*x20))/(c0 + c1*x10 + c2*x20)))/(d0 + d1*x10 + d2*x20)^4)))

  output <- list()
  output$llk <- - log(2*pi*del) - Dv + cmone/del + czero + cone*del

  return(output)
}