#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.