R/ModelB1.R

Defines functions ModelB1

Documented in ModelB1

#' Compute the maximum likelihood estimate of Model B1
#'
#'
#' @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 6-vector (r,h,a,b,c,d)
#' @export output a list with a llk variable storing the result of the log likelihood calculation
#' @examples
#' ModelB1(0.4,0.3,0.1,c(0.1,0.2,0.3,0.4,0.5,0.6))
#'
ModelB1 <- function(x,x0,del,param)
{

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

  r <- param[1]
  h <- param[2]
  a <- param[3]
  b <- param[4]
  c <- param[5]
  d <- param[6]

  Dv <- log(r) + log(x2) + (1/2)*log(1 - r^2)
  cm1 <- (3*h*(x1 - x10)*(x2 - x20)^2)/
  (4*(-r + r^3)*x20^(5/2)) - (x2 - x20)^3/(4*(-1 + r^2)*x20^2) -
  (h*(x1 - x10)^3*(x2 - x20)*(2*h^2 + x20 - r^2*x20)^2)/
  (24*r^3*(-1 + r^2)^2*x20^(11/2)) +
  ((x1 - x10)^4*(h^2 + x20 - r^2*x20)*(2*h^2 + x20 - r^2*x20)^2)/
  (96*r^4*(-1 + r^2)^2*x20^6) +
  ((x2 - x20)^4*(4*h^2 + 15*(-1 + r^2)*x20))/(96*(-1 + r^2)^2*x20^4) +
  ((x1 - x10)*(x2 - x20)^3*(-4*h^3 - 13*h*(-1 + r^2)*x20))/
  (24*r*(-1 + r^2)^2*x20^(9/2)) +
  ((x1 - x10)^2*(x2 - x20)^2*(12*h^4 + 16*h^2*(-1 + r^2)*x20 -
                              7*(-1 + r^2)^2*x20^2))/(48*r^2*(-1 + r^2)^2*x20^5) +
  ((x1 - x10)^2*(x2 - x20)*(2*h^2 + x20 - r^2*x20))/
  (4*r^2*x20^3 - 4*r^4*x20^3) +
  (h^2*(x1 - x10)^2 + (-((-1 + r^2)*x1^2) + 2*(-1 + r^2)*x1*x10 -
                       (-1 + r^2)*x10^2 + r^2*(x2 - x20)^2)*x20 +
    2*h*r*(x1 - x10)*sqrt(x20)*(-x2 + x20))/(2*r^2*(-1 + r^2)*x20^2)

  c0 <- ((x1 - x10)^2*
        (4*h^4 - 5*h^2*(-1 + r^2)*x20 + (-1 + r^2)^2*x20^2))/
  (24*r^2*(-1 + r^2)*x20^4) +
  ((x2 - x20)^2*(10*h^2*r + r*(1 + 12*c - r^2)*x20 -
                 6*h*sqrt(x20)*(3*a + b*x20)))/(24*r*(-1 + r^2)*x20^3) -
  ((x2 - x20)*(h^2*r - 4*h*sqrt(x20)*(a + b*x20) +
               4*r*x20*(c + d*x20)))/(4*r*(-1 + r^2)*x20^2) +
  ((x1 - x10)*(h^3*r - 4*h^2*sqrt(x20)*(a + b*x20) +
               4*(-1 + r^2)*x20^(3/2)*(a + b*x20) +
               h*r*x20*(-1 + 4*c + r^2 + 4*d*x20)))/(4*r^2*(-1 + r^2)*
                                                     x20^(5/2)) + ((x1 - x10)*(x2 - x20)*(-31*h^3*r -
                                                                                          24*a*(-1 + r^2)*x20^(3/2) + 24*h^2*sqrt(x20)*(2*a + b*x20) -
                                                                                          h*r*x20*(-5 + 36*c + 5*r^2 + 12*d*x20)))/
  (48*r^2*(-1 + r^2)*x20^(7/2))

  c1 <- (3*h^4*r^2 - 24*h^3*r*sqrt(x20)*
        (a + b*x20) - 24*h*r*x20^(3/2)*(a + b*x20)*
        (-1 + 4*c + r^2 + 4*d*x20) + h^2*x20*(48*a^2 + 13*r^4 +
                                              96*a*b*x20 + 48*b^2*x20^2 + r^2*(-13 + 24*c + 24*d*x20)) +
        8*x20^2*((-2 + 6*c)*r^4 + r^6 - 6*a^2*(-1 + r^2) -
                 12*a*b*(-1 + r^2)*x20 + 6*b^2*x20^2 +
                 r^2*(1 - 6*c + 6*c^2 + 12*c*d*x20 - 6*b^2*x20^2 + 6*d^2*x20^2)))/
  (96*r^2*(-1 + r^2)*x20^3)

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


  return(output)
  }
mfrdixon/MLEMVD documentation built on May 22, 2019, 8:52 p.m.