R/ModelU13.R

In mfrdixon/MLEMVD: Maximum Likelihood Estimation for Diffusion Equations

#' Compute the maximum likelihood estimate of Model U13
#'
#' @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 7-vector (am1,a0,a1,a2,a3,sigma,rho)
#' @export output a list with a llk variable storing the result of the log likelihood calculation
#' @examples
#' ModelU13(0.4,0.3,0.1,c(0.1,0.3,0.2,0.1,0.2,1.3,0.5))
#'

ModelU13 <- function(x,x0,del,param)
{

  a0  <- param[1]
  am1  <- param[2]
  a1  <- param[3]
  a2  <- param[4]
  a3  <- param[5]
  sigma  <- param[6]
  rho  <- param[7]

  s_x  <- sigma*x^rho

  y  <- x^(1 - rho)/((-1 + rho)*sigma)
  y0  <- x0^(1 - rho)/((-1 + rho)*sigma)

  cm1  <- (-(1/2))*(y - y0)^2


  c0  <- (1/2)*(a1*(1 - rho)*y^2 - a1*(1 - rho)*y0^2 - (2*a2*(-1 + rho)^(1 + (-2 + rho)/(-1 + rho))*sigma^(1/(1 - rho))*       (y^(2 + 1/(1 - rho)) - y0^(2 + 1/(1 - rho))))/(-3 + 2*rho) - (a3*(-1 + rho)^(1 + (-3 + rho)/(-1 + rho))*(y^(2 - 2/(-1 + rho)) - y0^(2 - 2/(-1 + rho))))/      (sigma^(2/(-1 + rho))*(-2 + rho)) - (2*a0*(-1 + rho)^(1 + rho/(-1 + rho))*sigma^(1/(-1 + rho))*(y^(2 + 1/(-1 + rho)) - y0^(2 + 1/(-1 + rho))))/      (-1 + 2*rho) - (am1*(-1 + rho)^(1 + (1 + rho)/(-1 + rho))*sigma^(2/(-1 + rho))*(y^((2*rho)/(-1 + rho)) - y0^((2*rho)/(-1 + rho))))/rho +      (2*rho*log(y/y0))/(-2 + 2*rho))


  c1  <- (-(1/(y - y0)))*(a1*(1 - rho)*y + (a1*(1 - rho)*rho*y)/(-2 + 2*rho) -      (1/3)*a3*am1*(-1 + rho)^((-3 + rho)/(-1 + rho) + (1 + rho)/(-1 + rho))*y^3 + (rho*(1/y - 1/y0))/(-2 + 2*rho) + (rho^2*(-(1/y) + 1/y0))/(2*(-2 + 2*rho)^2) +      (rho^2*(-(1/y) + 1/y0))/(2*(-1 + rho)*(-2 + 2*rho)) + (a1*(1 - rho)*(y - y0))/(2*(-1 + rho)) - (a1*(1 - rho)*rho*(y - y0))/(2*(-1 + rho)) -      a1*(1 - rho)*y0 - (a1*(1 - rho)*rho*y0)/(-2 + 2*rho) + (1/3)*a3*am1*(-1 + rho)^((-3 + rho)/(-1 + rho) + (1 + rho)/(-1 + rho))*y0^3 +      ((rho*y)/(-2 + 2*rho) - (rho*y0)/(-2 + 2*rho))/(2*y*y0 - 2*rho*y*y0) + (1/6)*a1^2*(1 - rho)^2*(y^3 - y0^3) +      (2/3)*a3*am1*(-1 + rho)^((-3 + rho)/(-1 + rho) + (1 + rho)/(-1 + rho))*(y^3 - y0^3) + (1/3)*a0*a2*(-1 + rho)^((-2 + rho)/(-1 + rho) + rho/(-1 + rho))*      sigma^(1/(1 - rho) + 1/(-1 + rho))*(y^3 - y0^3) + (a2*(-1 + rho)^((-2 + rho)/(-1 + rho))*sigma^(1/(1 - rho))*(y^(1 + 1/(1 - rho)) - y0^(1 + 1/(1 - rho))))/      (-2 + rho) - (a2*(-1 + rho)^((-2 + rho)/(-1 + rho))*rho*sigma^(1/(1 - rho))*(y^(1 + 1/(1 - rho)) - y0^(1 + 1/(1 - rho))))/(2*(-2 + rho)) -      (a2*(-1 + rho)^(1 + (-2 + rho)/(-1 + rho))*rho*sigma^(1/(1 - rho))*(y^(1 + 1/(1 - rho)) - y0^(1 + 1/(1 - rho))))/((-2 + rho)*(-2 + 2*rho)) -
      (a1*a2*(1 - rho)*(-1 + rho)^(1 + (-2 + rho)/(-1 + rho))*sigma^(1/(1 - rho))*(y^(3 + 1/(1 - rho)) -
  y0^(3 + 1/(1 - rho))))/(-4 + 3*rho) +      (a0*a3*(-1 + rho)^(1 + (-3 + rho)/(-1 + rho) + rho/(-1 + rho))*(y^(3 + 1/(1 - rho)) - y0^(3 + 1/(1 - rho))))/(sigma^(1/(-1 + rho))*(-4 + 3*rho)) +      (a3^2*(-1 + rho)^(1 + (2*(-3 + rho))/(-1 + rho))*(y^(3 - 4/(-1 + rho)) - y0^(3 - 4/(-1 + rho))))/(sigma^(4/(-1 + rho))*(2*(-7 + 3*rho))) +      (a2*a3*(-1 + rho)^(1 + (-3 + rho)/(-1 + rho) + (-2 + rho)/(-1 + rho))*sigma^(1/(1 - rho) - 2/(-1 + rho))*(y^(3 - 3/(-1 + rho)) - y0^(3 - 3/(-1 + rho))))/      (3*(-2 + rho)) + (a2^2*(-1 + rho)^(1 + (2*(-2 + rho))/(-1 + rho))*sigma^(2/(1 - rho))*(y^(3 - 2/(-1 + rho)) - y0^(3 - 2/(-1 + rho))))/(2*(-5 + 3*rho)) -      (a1*a3*(1 - rho)*(-1 + rho)^(1 + (-3 + rho)/(-1 + rho))*(y^(3 - 2/(-1 + rho)) - y0^(3 - 2/(-1 + rho))))/(sigma^(2/(-1 + rho))*(-5 + 3*rho)) +      (a2*am1*(-1 + rho)^(1 + (-2 + rho)/(-1 + rho) + (1 + rho)/(-1 + rho))*sigma^(1/(1 - rho) + 2/(-1 + rho))*(y^(3 + 1/(-1 + rho)) - y0^(3 + 1/(-1 + rho))))/      (-2 + 3*rho) - (a0*a1*(1 - rho)*(-1 + rho)^(1 + rho/(-1 + rho))*sigma^(1/(-1 + rho))*(y^(3 + 1/(-1 + rho)) - y0^(3 + 1/(-1 + rho))))/(-2 + 3*rho) +
    (a0^2*(-1 + rho)^(1 + (2*rho)/(-1 + rho))*sigma^(2/(-1 + rho))*(y^(3 + 2/(-1 + rho)) - y0^(3 + 2/(-1 + rho))))/(-1 + 3*rho) -      (a1*am1*(1 - rho)*(-1 + rho)^(1 + (1 + rho)/(-1 + rho))*sigma^(2/(-1 + rho))*(y^(3 + 2/(-1 + rho)) - y0^(3 + 2/(-1 + rho))))/(-1 + 3*rho) -      (a0^2*(-1 + rho)^(1 + (2*rho)/(-1 + rho))*sigma^(2/(-1 + rho))*(y^(3 + 2/(-1 + rho)) - y0^(3 + 2/(-1 + rho))))/(-2 + 6*rho) +      (am1^2*(-1 + rho)^(1 + (2*(1 + rho))/(-1 + rho))*sigma^(4/(-1 + rho))*(y^(3 + 4/(-1 + rho)) - y0^(3 + 4/(-1 + rho))))/(1 + 3*rho) -      (am1^2*(-1 + rho)^(1 + (2*(1 + rho))/(-1 + rho))*sigma^(4/(-1 + rho))*(y^(3 + 4/(-1 + rho)) - y0^(3 + 4/(-1 + rho))))/(2 + 6*rho) -      (3*a3*(-1 + rho)^((-3 + rho)/(-1 + rho))*(y^((-3 + rho)/(-1 + rho)) - y0^((-3 + rho)/(-1 + rho))))/(sigma^(2/(-1 + rho))*(3 - rho)) -      (a3*(-1 + rho)^((-3 + rho)/(-1 + rho))*rho*(y^((-3 + rho)/(-1 + rho)) - y0^((-3 + rho)/(-1 + rho))))/(sigma^(2/(-1 + rho))*(6 - 2*rho)) -      (1/2)*a0*(-1 + rho)^(rho/(-1 + rho))*sigma^(1/(-1 + rho))*(y^(rho/(-1 + rho)) - y0^(rho/(-1 + rho))) -      (a0*(-1 + rho)^(1 + rho/(-1 + rho))*sigma^(1/(-1 + rho))*(y^(rho/(-1 + rho)) - y0^(rho/(-1 + rho))))/(-2 + 2*rho) -
    (3*a3*(-1 + rho)^((-3 + rho)/(-1 + rho))*(y^(rho/(-1 + rho))*y0^(3/(-1 + rho)) - y^(3/(-1 + rho))*y0^(rho/(-1 + rho))))/(sigma^(2/(-1 + rho))*(y*y0)^(3/(-1 + rho))*(2*(-3 + rho))) - (a3*(-1 + rho)^((-3 + rho)/(-1 + rho))*rho*       (y^(rho/(-1 + rho))*y0^(3/(-1 + rho)) - y^(3/(-1 + rho))*y0^(rho/(-1 + rho))))/(sigma^(2/(-1 + rho))*(y*y0)^(3/(-1 + rho))*(-3 + rho)) -      (a3*(-1 + rho)^(1 + (-3 + rho)/(-1 + rho))*rho*(y^(rho/(-1 + rho))*y0^(3/(-1 + rho)) - y^(3/(-1 + rho))*y0^(rho/(-1 + rho))))/      (sigma^(2/(-1 + rho))*(y*y0)^(3/(-1 + rho))*((-3 + rho)*(-2 + 2*rho))) + (a0*am1*(-1 + rho)^(1 + rho/(-1 + rho) + (1 + rho)/(-1 + rho))*       sigma^(3/(-1 + rho))*(y^((3*rho)/(-1 + rho)) - y0^((3*rho)/(-1 + rho))))/(3*rho) -      (am1*(-1 + rho)^((1 + rho)/(-1 + rho))*sigma^(2/(-1 + rho))*(y^((1 + rho)/(-1 + rho)) - y0^((1 + rho)/(-1 + rho))))/(2*(1 + rho)) -      (am1*(-1 + rho)^((1 + rho)/(-1 + rho))*rho*sigma^(2/(-1 + rho))*(y^((1 + rho)/(-1 + rho)) - y0^((1 + rho)/(-1 + rho))))/(2*(1 + rho)) -      (am1*(-1 + rho)^(1 + (1 + rho)/(-1 + rho))*rho*sigma^(2/(-1 + rho))*(y^((1 + rho)/(-1 + rho)) - y0^((1 + rho)/(-1 + rho))))/((1 + rho)*(-2 + 2*rho)))


  output <- list()
  output$llk  <- -log(2*pi*del)/2 - log(s_x) + cm1/del + c0 + c1*del

  return(output)
}