R/ModelB6.R
In mfrdixon/MLEMVD: Maximum Likelihood Estimation for Diffusion Equations
Defines functions
ModelB6
Documented in
ModelB6
#' Compute the maximum likelihood estimate of Model B6
#'
#' x1 is the log of the stock price, x2 the instantaneous variance.
#'
#' dx1 <- (m - x2/2) dt + sqrt(x2) dW1
#' dx2 <- (a - b x2) dt + s sqrt(1-r^2) Sqrt(x2) dW1 + s r sqrt(x2) dW2
#'
#' @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 5-vector (m, a, b, s, r):
#' @export output a list with a llk variable storing the result of the log likelihood calculation
#' @examples
#' ModelB6(0.4,0.3,0.1,c(0.1,0.2,0.3,0.3,-0.2))
#'
ModelB6<-function(x,x0,del,param){
#
#
# lnpX(del,x1,x2,x10,x20) = - log(2*del*pi) - Dv(x1,x2) + c(x1,x2,x10,x20,-1)/del + c(x1,x2,x10,x20,0) + del*c(x1,x2,x10,x20,1) + (1/2)*del^2*c(x1,x2,x10,x20,2);
# Dv(x1, x2) = (1/2)*log(1 - r^2) + log(s) + log(x2);
# c(x1, x2, x10, x20, -1) = (7*s^4*(x1 - x10)^6)/(11520*(-1 + r^2)^3*x20^5) - (7*r*s^3*(x1 - x10)^5*(x2 - x20))/(1920*(-1 + r^2)^3*x20^5) + ((-193 + 298*r^2)*s^2*(x1 - x10)^4*(x2 - x20)^2)/(11520*(-1 + r^2)^3*x20^5) + (r*(193 - 228*r^2)*s*(x1 - x10)^3*(x2 - x20)^3)/(2880*(-1 + r^2)^3*x20^5) + ((745 - 2648*r^2 + 2008*r^4)*(x1 - x10)^2*(x2 - x20)^4)/(11520*(-1 + r^2)^3*x20^5) + ((-745*r + 1876*r^3 - 1152*r^5)*(x1 - x10)*(x2 - x20)^5)/(5760*(-1 + r^2)^3*s*x20^5) + ((945 - 2090*r^2 + 1152*r^4)*(x2 - x20)^6)/(11520*(-1 + r^2)^3*s^2*x20^5) - (s^2*(x1 - x10)^4*(x2 - x20))/(64*(-1 + r^2)^2*x20^4) + (r*s*(x1 - x10)^3*(x2 - x20)^2)/(16*(-1 + r^2)^2*x20^4) + ((3 - 6*r^2)*(x1 - x10)^2*(x2 - x20)^3)/(32*(-1 + r^2)^2*x20^4) + (r*(-3 + 4*r^2)*(x1 - x10)*(x2 - x20)^4)/(16*(-1 + r^2)^2*s*x20^4) + ((7 - 8*r^2)*(x2 - x20)^5)/(64*(-1 + r^2)^2*s^2*x20^4) + (s^2*(x1 - x10)^4)/(96*(-1 + r^2)^2*x20^3) - (r*s*(x1 - x10)^3*(x2 - x20))/ (24*(-1 + r^2)^2*x20^3) + ((-7 + 10*r^2)*(x1 - x10)^2*(x2 - x20)^2)/ (48*(-1 + r^2)^2*x20^3) + ((7*r - 8*r^3)*(x1 - x10)*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s*x20^3) + ((-15 + 16*r^2)*(x2 - x20)^4)/ (96*(-1 + r^2)^2*s^2*x20^3) - ((x1 - x10)^2*(x2 - x20))/(4*(-1 + r^2)*x20^2) + (r*(x1 - x10)*(x2 - x20)^2)/(2*(-1 + r^2)*s*x20^2) - (x2 - x20)^3/(4*(-1 + r^2)*s^2*x20^2) + (s^2*(x1 - x10)^2 - 2*r*s*(x1 - x10)*(x2 - x20) + (x2 - x20)^2)/ (2*(-1 + r^2)*s^2*x20);
# c(x1, x2, x10, x20, 0) = (7*s^4*(x1 - x10)^4)/(1920*(-1 + r^2)^2*x20^4) - (s*(30*a*r + s*(-30*m + 7*r*s))*(x1 - x10)^3*(x2 - x20))/ (480*(-1 + r^2)^2*x20^4) + ((540*a*r^2 + s*(-540*m*r + (-97 + 160*r^2)*s))* (x1 - x10)^2*(x2 - x20)^2)/(2880*(-1 + r^2)^2*x20^4) + ((-270*a*r*(-1 + 2*r^2) + s*(270*m*(-1 + 2*r^2) + r*(97 - 118*r^2)*s))* (x1 - x10)*(x2 - x20)^3)/(1440*(-1 + r^2)^2*s*x20^4) + ((360*a*(-4 + 5*r^2) + s*(-360*m*r*(-3 + 4*r^2) + (-215 + 236*r^2)*s))* (x2 - x20)^4)/(5760*(-1 + r^2)^2*s^2*x20^4) + (s*(a*r - m*s)*(x1 - x10)^3)/ (24*(-1 + r^2)^2*x20^3) + ((-3*a*r^2 + s*(3*m*r + s - r^2*s))* (x1 - x10)^2*(x2 - x20))/(24*(-1 + r^2)^2*x20^3) + ((a*r*(-7 + 10*r^2) + s*(m*(7 - 10*r^2) + 2*r*(-1 + r^2)*s))*(x1 - x10)* (x2 - x20)^2)/(24*(-1 + r^2)^2*s*x20^3) + ((a*(8 - 9*r^2) + s*(m*r*(-7 + 8*r^2) + s - r^2*s))*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s^2*x20^3) + (s^2*(x1 - x10)^2)/(24*(-1 + r^2)*x20^2) + ((12*a + s*(-12*m*r + s))*(x2 - x20)^2)/(24*(-1 + r^2)*s^2*x20^2) - ((x1 - x10)*(2*a*r - 2*m*s - 2*b*r*x20 + s*x20))/ (2*s*x20 - 2*r^2*s*x20) + ((x2 - x20)*(2*a - 2*m*r*s - 2*b*x20 + r*s*x20))/(2*s^2*x20 - 2*r^2*s^2*x20) + ((6*a*r - 6*m*s + r*s^2)*(x1 - x10)*(x2 - x20))/ (12*s*x20^2 - 12*r^2*s*x20^2);
# c(x1, x2, x10, x20, 1) = (s*(a*r - m*s)*(x1 - x10))/(12*(-1 + r^2)*x20^2) + ((x1 - x10)^2*(60*a^2*(1 + 2*r^2) + 180*m^2*s^2 + 2*s^4 - 2*r^2*s^4 - 60*a*s*(6*m*r + s - r^2*s) - 60*b^2*x20^2 + 60*b*r*s*x20^2 - 15*s^2*x20^2))/(2880*(-1 + r^2)^2*x20^3) + ((x2 - x20)*(-12*a^2 - 12*m^2*s^2 + 4*m*r*s^3 - 2*s^4 + 2*r^2*s^4 + 4*a*s*(6*m*r + (3 - 4*r^2)*s) + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2))/(48*(-1 + r^2)*s^2*x20^2) + (1/(2880*(-1 + r^2)^2*s^2*x20^3))* ((x2 - x20)^2*(180*a^2*(-3 + 4*r^2) + 60*m^2*(-7 + 10*r^2)*s^2 - 240*m*r*(-1 + r^2)*s^3 - 94*s^4 + 190*r^2*s^4 - 96*r^4*s^4 + 60*a*s*(m*(14*r - 20*r^3) + (9 - 23*r^2 + 14*r^4)*s) + 60*b^2*x20^2 - 120*b^2*r^2*x20^2 - 60*b*r*s*x20^2 + 120*b*r^3*s*x20^2 + 15*s^2*x20^2 - 30*r^2*s^2*x20^2)) + (1/(24*(-1 + r^2)*s^2*x20))* (12*a^2 + 12*m^2*s^2 + 2*s^4 - 2*r^2*s^4 + 12*m*(2*b*r - s)*s* x20 + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2 - 12*a*(2*m*r*s + s^2 - r^2*s^2 + 2*b*x20 - r*s*x20)) + (1/(1440*(-1 + r^2)^2*s*x20^3))*((x1 - x10)*(x2 - x20)*(-60*a^2*(r + 2*r^3) - 180*m^2*r*s^2 + 120*m*(-1 + r^2)*s^3 + 180*a*r*s* (2*m*r + s - r^2*s) + r*(2*(-1 + r^2)*s^4 + 60*b^2*x20^2 - 60*b*r*s*x20^2 + 15*s^2*x20^2)));
# c(x1, x2, x10, x20, 2) = -((60*a^2*(-2 + r^2) - 60*m^2*s^2 + 23*(-1 + r^2)*s^4 + 120*a*s*(m*r + s - r^2*s))/(720*(-1 + r^2)*x20^2));
x1 <- x[1]
x2 <- x[2]
x10 <- x0[1]
x20 <- x0[2]
m <- param[1]
a <- param[2]
b <- param[3]
s <- param[4]
r <- param[5]
# old
# cm1 <- -((7*s^4*(x1 - x10)^6)/(11520*r^6*x20^5)) + (7*sqrt(1 - r^2)*s^3*(x1 - x10)^5*(x2 - x20))/(1920*r^6*x20^5) + ((-105 + 298*r^2)*s^2*(x1 - x10)^4*(x2 - x20)^2)/(11520*r^6*x20^5) + ((35 - 228*r^2)*sqrt(1 - r^2)*s*(x1 - x10)^3*(x2 - x20)^3)/(2880*r^6*x20^5) - ((105 - 1368*r^2 + 2008*r^4)*(x1 - x10)^2*(x2 - x20)^4)/(11520*r^6*x20^5) + (sqrt(1 - r^2)*(21 - 428*r^2 + 1152*r^4)*(x1 - x10)*(x2 - x20)^5)/(5760*r^6*s*x20^5) - ((7 - 214*r^2 + 1152*r^4)*(x2 - x20)^6)/(11520*r^6*s^2*x20^5) - (s^2*(x1 - x10)^4*(x2 - x20))/(64*r^4*x20^4) + (sqrt(1 - r^2)*s*(x1 - x10)^3*(x2 - x20)^2)/(16*r^4*x20^4) + (3*(-1 + 2*r^2)*(x1 - x10)^2*(x2 - x20)^3)/(32*r^4*x20^4) + ((1 - 4*r^2)*sqrt(1 - r^2)*(x1 - x10)*(x2 - x20)^4)/(16*r^4*s*x20^4) + ((-1 + 8*r^2)*(x2 - x20)^5)/(64*r^4*s^2*x20^4) + (s^2*(x1 - x10)^4)/(96*r^4*x20^3) - (sqrt(1 - r^2)*s*(x1 - x10)^3*(x2 - x20))/(24*r^4*x20^3) + ((3 - 10*r^2)*(x1 - x10)^2*(x2 - x20)^2)/(48*r^4*x20^3) + (sqrt(1 - r^2)*(-1 + 8*r^2)*(x1 - x10)*(x2 - x20)^3)/(24*r^4*s*x20^3) + ((1 - 16*r^2)*(x2 - x20)^4)/(96*r^4*s^2*x20^3) + ((x1 - x10)^2*(x2 - x20))/(4*r^2*x20^2) - (sqrt(1 - r^2)*(x1 - x10)*(x2 - x20)^2)/(2*r^2*s*x20^2) + (x2 - x20)^3/(4*r^2*s^2*x20^2) - (s^2*(x1 - x10)^2 - 2*sqrt(1 - r^2)*s*(x1 - x10)*(x2 - x20) + (x2 - x20)^2)/(2*r^2*s^2*x20)
# c0 <- (7*s^4*(x1-x10)^4)/(1920*r^4*x20^4)-(1/(480*r^4*x20^4))*(s*(30*a*sqrt(1-r^2)+s*(-30*m+7*sqrt(1-r^2)*s))*(x1-x10)^3*(x2-x20))+(1/(2880*r^4*x20^4))*((-540*a*(-1+r^2)+s*(-540*m*sqrt(1-r^2)+(63-160*r^2)*s))*(x1-x10)^2*(x2-x20)^2)+(1/(1440*r^4*s*x20^4))*((270*a*sqrt(1-r^2)*(-1+2*r^2)+s*(-270*m*(-1+2*r^2)+sqrt(1-r^2)*(-21+118*r^2)*s))*(x1-x10)*(x2-x20)^3)+(1/(5760*r^4*s^2*x20^4))*((-360*a*(-1+5*r^2)+s*(360*m*sqrt(1-r^2)*(-1+4*r^2)+(21-236*r^2)*s))*(x2-x20)^4)+(s*(a*sqrt(1-r^2)-m*s)*(x1-x10)^3)/(24*r^4*x20^3)+((3*a*(-1+r^2)+s*(3*m*sqrt(1-r^2)+r^2*s))*(x1-x10)^2*(x2-x20))/(24*r^4*x20^3)+(1/(24*r^4*s*x20^3))*((a*(3-10*r^2)*sqrt(1-r^2)+s*(m*(-3+10*r^2)-2*r^2*sqrt(1-r^2)*s))*(x1-x10)*(x2-x20)^2)+((a*(-1+9*r^2)+s*(m*(1-8*r^2)*sqrt(1-r^2)+r^2*s))*(x2-x20)^3)/(24*r^4*s^2*x20^3)-(s^2*(x1-x10)^2)/(24*r^2*x20^2)+((6*a*sqrt(1-r^2)-6*m*s+sqrt(1-r^2)*s^2)*(x1-x10)*(x2-x20))/(12*r^2*s*x20^2)-((12*a-12*m*sqrt(1-r^2)*s+s^2)*(x2-x20)^2)/(24*r^2*s^2*x20^2)-((x1-x10)*(2*a*sqrt(1-r^2)-2*m*s-2*b*sqrt(1-r^2)*x20+s*x20))/(2*r^2*s*x20)+((x2-x20)*(2*a-2*m*sqrt(1-r^2)*s-2*b*x20+sqrt(1-r^2)*s*x20))/(2*r^2*s^2*x20)
# c1 <- (s*((-a)*sqrt(1-r^2)+m*s)*(x1-x10))/(12*r^2*x20^2)+(1/(2880*r^4*x20^3))*((x1-x10)^2*(-60*a^2*(-3+2*r^2)+180*m^2*s^2+2*r^2*s^4-60*a*s*(6*m*sqrt(1-r^2)+r^2*s)-60*b^2*x20^2+60*b*sqrt(1-r^2)*s*x20^2-15*s^2*x20^2))+(1/(48*r^2*s^2*x20^2))*((x2-x20)*(12*a^2+12*m^2*s^2-4*m*sqrt(1-r^2)*s^3+2*r^2*s^4-4*a*s*(6*m*sqrt(1-r^2)+(-1+4*r^2)*s)-12*b^2*x20^2+12*b*sqrt(1-r^2)*s*x20^2-3*s^2*x20^2))+(1/(2880*r^4*s^2*x20^3))*((x2-x20)^2*(-180*a^2*(-1+4*r^2)-60*m^2*(-3+10*r^2)*s^2+240*m*r^2*sqrt(1-r^2)*s^3+2*r^2*s^4-96*r^4*s^4+60*a*s*(2*m*sqrt(1-r^2)*(-3+10*r^2)+r^2*(-5+14*r^2)*s)-60*b^2*x20^2+120*b^2*r^2*x20^2+60*b*sqrt(1-r^2)*s*x20^2-120*b*r^2*sqrt(1-r^2)*s*x20^2-15*s^2*x20^2+30*r^2*s^2*x20^2))+(1/(1440*r^4*s*x20^3))*((x1-x10)*(x2-x20)*(60*a^2*sqrt(1-r^2)*(-3+2*r^2)-180*m^2*sqrt(1-r^2)*s^2-120*m*r^2*s^3-2*r^2*sqrt(1-r^2)*s^4+180*a*s*(-2*m*(-1+r^2)+r^2*sqrt(1-r^2)*s)+60*b^2*sqrt(1-r^2)*x20^2-60*b*s*x20^2+60*b*r^2*s*x20^2+15*sqrt(1-r^2)*s^2*x20^2))-(1/(24*r^2*s^2*x20))*(12*a^2+12*m^2*s^2+2*r^2*s^4-12*m*s*(-2*b*sqrt(1-r^2)+s)*x20+12*b^2*x20^2-12*b*sqrt(1-r^2)*s*x20^2+3*s^2*x20^2-12*a*(2*m*sqrt(1-r^2)*s+r^2*s^2+2*b*x20-sqrt(1-r^2)*s*x20))
# c2 <- -((60*a^2*(1+r^2)+60*m^2*s^2+23*r^2*s^4-120*a*s*(m*sqrt(1-r^2)+r^2*s))/(720*r^2*x20^2))
Dv <- (1/2)*log(1 - r^2) + log(s) + log(x2)
cm1 <- (7*s^4*(x1 - x10)^6)/(11520*(-1 + r^2)^3*x20^5) - (7*r*s^3*(x1 - x10)^5*(x2 - x20))/(1920*(-1 + r^2)^3*x20^5) + ((-193 + 298*r^2)*s^2*(x1 - x10)^4*(x2 - x20)^2)/(11520*(-1 + r^2)^3*x20^5) + (r*(193 - 228*r^2)*s*(x1 - x10)^3*(x2 - x20)^3)/(2880*(-1 + r^2)^3*x20^5) + ((745 - 2648*r^2 + 2008*r^4)*(x1 - x10)^2*(x2 - x20)^4)/(11520*(-1 + r^2)^3*x20^5) + ((-745*r + 1876*r^3 - 1152*r^5)*(x1 - x10)*(x2 - x20)^5)/(5760*(-1 + r^2)^3*s*x20^5) + ((945 - 2090*r^2 + 1152*r^4)*(x2 - x20)^6)/(11520*(-1 + r^2)^3*s^2*x20^5) - (s^2*(x1 - x10)^4*(x2 - x20))/(64*(-1 + r^2)^2*x20^4) + (r*s*(x1 - x10)^3*(x2 - x20)^2)/(16*(-1 + r^2)^2*x20^4) + ((3 - 6*r^2)*(x1 - x10)^2*(x2 - x20)^3)/(32*(-1 + r^2)^2*x20^4) + (r*(-3 + 4*r^2)*(x1 - x10)*(x2 - x20)^4)/(16*(-1 + r^2)^2*s*x20^4) + ((7 - 8*r^2)*(x2 - x20)^5)/(64*(-1 + r^2)^2*s^2*x20^4) + (s^2*(x1 - x10)^4)/(96*(-1 + r^2)^2*x20^3) - (r*s*(x1 - x10)^3*(x2 - x20))/ (24*(-1 + r^2)^2*x20^3) + ((-7 + 10*r^2)*(x1 - x10)^2*(x2 - x20)^2)/ (48*(-1 + r^2)^2*x20^3) + ((7*r - 8*r^3)*(x1 - x10)*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s*x20^3) + ((-15 + 16*r^2)*(x2 - x20)^4)/ (96*(-1 + r^2)^2*s^2*x20^3) - ((x1 - x10)^2*(x2 - x20))/(4*(-1 + r^2)*x20^2) + (r*(x1 - x10)*(x2 - x20)^2)/(2*(-1 + r^2)*s*x20^2) - (x2 - x20)^3/(4*(-1 + r^2)*s^2*x20^2) + (s^2*(x1 - x10)^2 - 2*r*s*(x1 - x10)*(x2 - x20) + (x2 - x20)^2)/ (2*(-1 + r^2)*s^2*x20)
c0 <- (7*s^4*(x1 - x10)^4)/(1920*(-1 + r^2)^2*x20^4) - (s*(30*a*r + s*(-30*m + 7*r*s))*(x1 - x10)^3*(x2 - x20))/ (480*(-1 + r^2)^2*x20^4) + ((540*a*r^2 + s*(-540*m*r + (-97 + 160*r^2)*s))* (x1 - x10)^2*(x2 - x20)^2)/(2880*(-1 + r^2)^2*x20^4) + ((-270*a*r*(-1 + 2*r^2) + s*(270*m*(-1 + 2*r^2) + r*(97 - 118*r^2)*s))* (x1 - x10)*(x2 - x20)^3)/(1440*(-1 + r^2)^2*s*x20^4) + ((360*a*(-4 + 5*r^2) + s*(-360*m*r*(-3 + 4*r^2) + (-215 + 236*r^2)*s))* (x2 - x20)^4)/(5760*(-1 + r^2)^2*s^2*x20^4) + (s*(a*r - m*s)*(x1 - x10)^3)/ (24*(-1 + r^2)^2*x20^3) + ((-3*a*r^2 + s*(3*m*r + s - r^2*s))* (x1 - x10)^2*(x2 - x20))/(24*(-1 + r^2)^2*x20^3) + ((a*r*(-7 + 10*r^2) + s*(m*(7 - 10*r^2) + 2*r*(-1 + r^2)*s))*(x1 - x10)* (x2 - x20)^2)/(24*(-1 + r^2)^2*s*x20^3) + ((a*(8 - 9*r^2) + s*(m*r*(-7 + 8*r^2) + s - r^2*s))*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s^2*x20^3) + (s^2*(x1 - x10)^2)/(24*(-1 + r^2)*x20^2) + ((12*a + s*(-12*m*r + s))*(x2 - x20)^2)/(24*(-1 + r^2)*s^2*x20^2) - ((x1 - x10)*(2*a*r - 2*m*s - 2*b*r*x20 + s*x20))/ (2*s*x20 - 2*r^2*s*x20) + ((x2 - x20)*(2*a - 2*m*r*s - 2*b*x20 + r*s*x20))/(2*s^2*x20 - 2*r^2*s^2*x20) + ((6*a*r - 6*m*s + r*s^2)*(x1 - x10)*(x2 - x20))/ (12*s*x20^2 - 12*r^2*s*x20^2)
c1 <- (s*(a*r - m*s)*(x1 - x10))/(12*(-1 + r^2)*x20^2) + ((x1 - x10)^2*(60*a^2*(1 + 2*r^2) + 180*m^2*s^2 + 2*s^4 - 2*r^2*s^4 - 60*a*s*(6*m*r + s - r^2*s) - 60*b^2*x20^2 + 60*b*r*s*x20^2 - 15*s^2*x20^2))/(2880*(-1 + r^2)^2*x20^3) + ((x2 - x20)*(-12*a^2 - 12*m^2*s^2 + 4*m*r*s^3 - 2*s^4 + 2*r^2*s^4 + 4*a*s*(6*m*r + (3 - 4*r^2)*s) + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2))/(48*(-1 + r^2)*s^2*x20^2) + (1/(2880*(-1 + r^2)^2*s^2*x20^3))* ((x2 - x20)^2*(180*a^2*(-3 + 4*r^2) + 60*m^2*(-7 + 10*r^2)*s^2 - 240*m*r*(-1 + r^2)*s^3 - 94*s^4 + 190*r^2*s^4 - 96*r^4*s^4 + 60*a*s*(m*(14*r - 20*r^3) + (9 - 23*r^2 + 14*r^4)*s) + 60*b^2*x20^2 - 120*b^2*r^2*x20^2 - 60*b*r*s*x20^2 + 120*b*r^3*s*x20^2 + 15*s^2*x20^2 - 30*r^2*s^2*x20^2)) + (1/(24*(-1 + r^2)*s^2*x20))* (12*a^2 + 12*m^2*s^2 + 2*s^4 - 2*r^2*s^4 + 12*m*(2*b*r - s)*s* x20 + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2 - 12*a*(2*m*r*s + s^2 - r^2*s^2 + 2*b*x20 - r*s*x20)) + (1/(1440*(-1 + r^2)^2*s*x20^3))*((x1 - x10)*(x2 - x20)*(-60*a^2*(r + 2*r^3) - 180*m^2*r*s^2 + 120*m*(-1 + r^2)*s^3 + 180*a*r*s* (2*m*r + s - r^2*s) + r*(2*(-1 + r^2)*s^4 + 60*b^2*x20^2 - 60*b*r*s*x20^2 + 15*s^2*x20^2)))
c2 <- -((60*a^2*(-2 + r^2) - 60*m^2*s^2 + 23*(-1 + r^2)*s^4 + 120*a*s*(m*r + s - r^2*s))/(720*(-1 + r^2)*x20^2))
output <- list()
output$llk <- -log(2*del*pi) - Dv + cm1/del + c0 + del*c1 + (1/2)*del^2*c2
# print(paste('r',r,sep='='))
# print(paste('s',s,sep='='))
# print(paste('x2',x2,sep='='))
# a1 <- param[1]
# b1 <- -0.5
# a2 <- param[2]
# b2 <- param[3]
# s <- param[4]
# r <- param[5]
#
# Dv <- 0.5*log((1-r^2)*s*x2)
# c1 <- (r*s*a2*b1+r*s*a1*b2-s^2*a1*b1-a2*b2)/((1-r^2)*s^2) + (2*r*s*b1*b2-s^2*b1^2-b2^2)*x20/(2*(1-r^2)*s^2) - (s^4-r^2*s^4+6*s^2*a1^2-6*s^2*a2+6*r^2*s^2*a2 - 12*r*s*a1*a2 +6*a2^2)/(12*(1-r^2)*s^2*x20)
# c0 <- (x2-x20)*(x20*b2-r*s*x20*b1 -r*s*a1+a2)/((1-r^2)*s^2*x20) -(x1-x10)*(r*x20*b2-s*x20*b1-s*a1+r*a2)/((1-r^2)*s*x20) -s^2*(x1-x10)^2/(24*(1-r^2)*x20^2) - (x2-x20)^2*(s*(s-12*r*a1) + 12*a2)/(24*(1-r^2)*s^2*x20^2) + (x2-x20)*(x1-x10)*(r*s^2-6*s*a1 + 6*r*a2)/(12*s*x20^2*(1-r^2))
# cm1 <- -((x2-x20)^2 - 2*r*s*(x2-x20)*(x1-x10) + s^2*(x1-x10)^2)/(2*(1-r^2)*s^2*x20) + (x2-x20)^3/(4*(1-r^2)*s^2*x20^2) - r*(x2-x20)^2*(x1-x10)/(2*(1-r^2)*s*x20^2) + (x2-x20)*(x1-x10)^2/(4*(1-r^2)*x20^2) + (7*r-8*r^3)*(x2-x20)^3*(x1-x10)/(24*(1-r^2)^2*s*x20^3) - (7-10*r^2)*(x2-x20)^2*(x1-x10)^2/(48*(1-r^2)^2*x20^3) - r*s*(x2-x20)*(x1-x10)^3/(24*(1-r^2)^2*x20^3) + (s^2*(x1-x10)^4)/(96*(1-r^2)^2*x20^3) - (15-16*r^2)*(x2-x20)^4/(96*(1-r^2)^2*s^2*x20^3)
# output <- list()
# output$llk <- -log(2*del*pi) - Dv + cm1/del + c0 + del*c1
# print(paste('Dv',Dv,sep='='))
# print(paste('c1',c1,sep='='))
# print(paste('c0',c0,sep='='))
# print(paste('cm1',cm1,sep='='))
return(output)
}
mfrdixon/MLEMVD documentation built on May 22, 2019, 8:52 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ModelB6
Documented in ModelB6
#' Compute the maximum likelihood estimate of Model B6
#'
#' x1 is the log of the stock price, x2 the instantaneous variance.
#'
#' dx1 <- (m - x2/2) dt + sqrt(x2) dW1
#' dx2 <- (a - b x2) dt + s sqrt(1-r^2) Sqrt(x2) dW1 + s r sqrt(x2) dW2
#'
#' @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 5-vector (m, a, b, s, r):
#' @export output a list with a llk variable storing the result of the log likelihood calculation
#' @examples
#' ModelB6(0.4,0.3,0.1,c(0.1,0.2,0.3,0.3,-0.2))
#'
ModelB6<-function(x,x0,del,param){
#
#
# lnpX(del,x1,x2,x10,x20) = - log(2*del*pi) - Dv(x1,x2) + c(x1,x2,x10,x20,-1)/del + c(x1,x2,x10,x20,0) + del*c(x1,x2,x10,x20,1) + (1/2)*del^2*c(x1,x2,x10,x20,2);
# Dv(x1, x2) = (1/2)*log(1 - r^2) + log(s) + log(x2);
# c(x1, x2, x10, x20, -1) = (7*s^4*(x1 - x10)^6)/(11520*(-1 + r^2)^3*x20^5) - (7*r*s^3*(x1 - x10)^5*(x2 - x20))/(1920*(-1 + r^2)^3*x20^5) + ((-193 + 298*r^2)*s^2*(x1 - x10)^4*(x2 - x20)^2)/(11520*(-1 + r^2)^3*x20^5) + (r*(193 - 228*r^2)*s*(x1 - x10)^3*(x2 - x20)^3)/(2880*(-1 + r^2)^3*x20^5) + ((745 - 2648*r^2 + 2008*r^4)*(x1 - x10)^2*(x2 - x20)^4)/(11520*(-1 + r^2)^3*x20^5) + ((-745*r + 1876*r^3 - 1152*r^5)*(x1 - x10)*(x2 - x20)^5)/(5760*(-1 + r^2)^3*s*x20^5) + ((945 - 2090*r^2 + 1152*r^4)*(x2 - x20)^6)/(11520*(-1 + r^2)^3*s^2*x20^5) - (s^2*(x1 - x10)^4*(x2 - x20))/(64*(-1 + r^2)^2*x20^4) + (r*s*(x1 - x10)^3*(x2 - x20)^2)/(16*(-1 + r^2)^2*x20^4) + ((3 - 6*r^2)*(x1 - x10)^2*(x2 - x20)^3)/(32*(-1 + r^2)^2*x20^4) + (r*(-3 + 4*r^2)*(x1 - x10)*(x2 - x20)^4)/(16*(-1 + r^2)^2*s*x20^4) + ((7 - 8*r^2)*(x2 - x20)^5)/(64*(-1 + r^2)^2*s^2*x20^4) + (s^2*(x1 - x10)^4)/(96*(-1 + r^2)^2*x20^3) - (r*s*(x1 - x10)^3*(x2 - x20))/ (24*(-1 + r^2)^2*x20^3) + ((-7 + 10*r^2)*(x1 - x10)^2*(x2 - x20)^2)/ (48*(-1 + r^2)^2*x20^3) + ((7*r - 8*r^3)*(x1 - x10)*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s*x20^3) + ((-15 + 16*r^2)*(x2 - x20)^4)/ (96*(-1 + r^2)^2*s^2*x20^3) - ((x1 - x10)^2*(x2 - x20))/(4*(-1 + r^2)*x20^2) + (r*(x1 - x10)*(x2 - x20)^2)/(2*(-1 + r^2)*s*x20^2) - (x2 - x20)^3/(4*(-1 + r^2)*s^2*x20^2) + (s^2*(x1 - x10)^2 - 2*r*s*(x1 - x10)*(x2 - x20) + (x2 - x20)^2)/ (2*(-1 + r^2)*s^2*x20);
# c(x1, x2, x10, x20, 0) = (7*s^4*(x1 - x10)^4)/(1920*(-1 + r^2)^2*x20^4) - (s*(30*a*r + s*(-30*m + 7*r*s))*(x1 - x10)^3*(x2 - x20))/ (480*(-1 + r^2)^2*x20^4) + ((540*a*r^2 + s*(-540*m*r + (-97 + 160*r^2)*s))* (x1 - x10)^2*(x2 - x20)^2)/(2880*(-1 + r^2)^2*x20^4) + ((-270*a*r*(-1 + 2*r^2) + s*(270*m*(-1 + 2*r^2) + r*(97 - 118*r^2)*s))* (x1 - x10)*(x2 - x20)^3)/(1440*(-1 + r^2)^2*s*x20^4) + ((360*a*(-4 + 5*r^2) + s*(-360*m*r*(-3 + 4*r^2) + (-215 + 236*r^2)*s))* (x2 - x20)^4)/(5760*(-1 + r^2)^2*s^2*x20^4) + (s*(a*r - m*s)*(x1 - x10)^3)/ (24*(-1 + r^2)^2*x20^3) + ((-3*a*r^2 + s*(3*m*r + s - r^2*s))* (x1 - x10)^2*(x2 - x20))/(24*(-1 + r^2)^2*x20^3) + ((a*r*(-7 + 10*r^2) + s*(m*(7 - 10*r^2) + 2*r*(-1 + r^2)*s))*(x1 - x10)* (x2 - x20)^2)/(24*(-1 + r^2)^2*s*x20^3) + ((a*(8 - 9*r^2) + s*(m*r*(-7 + 8*r^2) + s - r^2*s))*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s^2*x20^3) + (s^2*(x1 - x10)^2)/(24*(-1 + r^2)*x20^2) + ((12*a + s*(-12*m*r + s))*(x2 - x20)^2)/(24*(-1 + r^2)*s^2*x20^2) - ((x1 - x10)*(2*a*r - 2*m*s - 2*b*r*x20 + s*x20))/ (2*s*x20 - 2*r^2*s*x20) + ((x2 - x20)*(2*a - 2*m*r*s - 2*b*x20 + r*s*x20))/(2*s^2*x20 - 2*r^2*s^2*x20) + ((6*a*r - 6*m*s + r*s^2)*(x1 - x10)*(x2 - x20))/ (12*s*x20^2 - 12*r^2*s*x20^2);
# c(x1, x2, x10, x20, 1) = (s*(a*r - m*s)*(x1 - x10))/(12*(-1 + r^2)*x20^2) + ((x1 - x10)^2*(60*a^2*(1 + 2*r^2) + 180*m^2*s^2 + 2*s^4 - 2*r^2*s^4 - 60*a*s*(6*m*r + s - r^2*s) - 60*b^2*x20^2 + 60*b*r*s*x20^2 - 15*s^2*x20^2))/(2880*(-1 + r^2)^2*x20^3) + ((x2 - x20)*(-12*a^2 - 12*m^2*s^2 + 4*m*r*s^3 - 2*s^4 + 2*r^2*s^4 + 4*a*s*(6*m*r + (3 - 4*r^2)*s) + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2))/(48*(-1 + r^2)*s^2*x20^2) + (1/(2880*(-1 + r^2)^2*s^2*x20^3))* ((x2 - x20)^2*(180*a^2*(-3 + 4*r^2) + 60*m^2*(-7 + 10*r^2)*s^2 - 240*m*r*(-1 + r^2)*s^3 - 94*s^4 + 190*r^2*s^4 - 96*r^4*s^4 + 60*a*s*(m*(14*r - 20*r^3) + (9 - 23*r^2 + 14*r^4)*s) + 60*b^2*x20^2 - 120*b^2*r^2*x20^2 - 60*b*r*s*x20^2 + 120*b*r^3*s*x20^2 + 15*s^2*x20^2 - 30*r^2*s^2*x20^2)) + (1/(24*(-1 + r^2)*s^2*x20))* (12*a^2 + 12*m^2*s^2 + 2*s^4 - 2*r^2*s^4 + 12*m*(2*b*r - s)*s* x20 + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2 - 12*a*(2*m*r*s + s^2 - r^2*s^2 + 2*b*x20 - r*s*x20)) + (1/(1440*(-1 + r^2)^2*s*x20^3))*((x1 - x10)*(x2 - x20)*(-60*a^2*(r + 2*r^3) - 180*m^2*r*s^2 + 120*m*(-1 + r^2)*s^3 + 180*a*r*s* (2*m*r + s - r^2*s) + r*(2*(-1 + r^2)*s^4 + 60*b^2*x20^2 - 60*b*r*s*x20^2 + 15*s^2*x20^2)));
# c(x1, x2, x10, x20, 2) = -((60*a^2*(-2 + r^2) - 60*m^2*s^2 + 23*(-1 + r^2)*s^4 + 120*a*s*(m*r + s - r^2*s))/(720*(-1 + r^2)*x20^2));
x1 <- x[1]
x2 <- x[2]
x10 <- x0[1]
x20 <- x0[2]
m <- param[1]
a <- param[2]
b <- param[3]
s <- param[4]
r <- param[5]
# old
# cm1 <- -((7*s^4*(x1 - x10)^6)/(11520*r^6*x20^5)) + (7*sqrt(1 - r^2)*s^3*(x1 - x10)^5*(x2 - x20))/(1920*r^6*x20^5) + ((-105 + 298*r^2)*s^2*(x1 - x10)^4*(x2 - x20)^2)/(11520*r^6*x20^5) + ((35 - 228*r^2)*sqrt(1 - r^2)*s*(x1 - x10)^3*(x2 - x20)^3)/(2880*r^6*x20^5) - ((105 - 1368*r^2 + 2008*r^4)*(x1 - x10)^2*(x2 - x20)^4)/(11520*r^6*x20^5) + (sqrt(1 - r^2)*(21 - 428*r^2 + 1152*r^4)*(x1 - x10)*(x2 - x20)^5)/(5760*r^6*s*x20^5) - ((7 - 214*r^2 + 1152*r^4)*(x2 - x20)^6)/(11520*r^6*s^2*x20^5) - (s^2*(x1 - x10)^4*(x2 - x20))/(64*r^4*x20^4) + (sqrt(1 - r^2)*s*(x1 - x10)^3*(x2 - x20)^2)/(16*r^4*x20^4) + (3*(-1 + 2*r^2)*(x1 - x10)^2*(x2 - x20)^3)/(32*r^4*x20^4) + ((1 - 4*r^2)*sqrt(1 - r^2)*(x1 - x10)*(x2 - x20)^4)/(16*r^4*s*x20^4) + ((-1 + 8*r^2)*(x2 - x20)^5)/(64*r^4*s^2*x20^4) + (s^2*(x1 - x10)^4)/(96*r^4*x20^3) - (sqrt(1 - r^2)*s*(x1 - x10)^3*(x2 - x20))/(24*r^4*x20^3) + ((3 - 10*r^2)*(x1 - x10)^2*(x2 - x20)^2)/(48*r^4*x20^3) + (sqrt(1 - r^2)*(-1 + 8*r^2)*(x1 - x10)*(x2 - x20)^3)/(24*r^4*s*x20^3) + ((1 - 16*r^2)*(x2 - x20)^4)/(96*r^4*s^2*x20^3) + ((x1 - x10)^2*(x2 - x20))/(4*r^2*x20^2) - (sqrt(1 - r^2)*(x1 - x10)*(x2 - x20)^2)/(2*r^2*s*x20^2) + (x2 - x20)^3/(4*r^2*s^2*x20^2) - (s^2*(x1 - x10)^2 - 2*sqrt(1 - r^2)*s*(x1 - x10)*(x2 - x20) + (x2 - x20)^2)/(2*r^2*s^2*x20)
# c0 <- (7*s^4*(x1-x10)^4)/(1920*r^4*x20^4)-(1/(480*r^4*x20^4))*(s*(30*a*sqrt(1-r^2)+s*(-30*m+7*sqrt(1-r^2)*s))*(x1-x10)^3*(x2-x20))+(1/(2880*r^4*x20^4))*((-540*a*(-1+r^2)+s*(-540*m*sqrt(1-r^2)+(63-160*r^2)*s))*(x1-x10)^2*(x2-x20)^2)+(1/(1440*r^4*s*x20^4))*((270*a*sqrt(1-r^2)*(-1+2*r^2)+s*(-270*m*(-1+2*r^2)+sqrt(1-r^2)*(-21+118*r^2)*s))*(x1-x10)*(x2-x20)^3)+(1/(5760*r^4*s^2*x20^4))*((-360*a*(-1+5*r^2)+s*(360*m*sqrt(1-r^2)*(-1+4*r^2)+(21-236*r^2)*s))*(x2-x20)^4)+(s*(a*sqrt(1-r^2)-m*s)*(x1-x10)^3)/(24*r^4*x20^3)+((3*a*(-1+r^2)+s*(3*m*sqrt(1-r^2)+r^2*s))*(x1-x10)^2*(x2-x20))/(24*r^4*x20^3)+(1/(24*r^4*s*x20^3))*((a*(3-10*r^2)*sqrt(1-r^2)+s*(m*(-3+10*r^2)-2*r^2*sqrt(1-r^2)*s))*(x1-x10)*(x2-x20)^2)+((a*(-1+9*r^2)+s*(m*(1-8*r^2)*sqrt(1-r^2)+r^2*s))*(x2-x20)^3)/(24*r^4*s^2*x20^3)-(s^2*(x1-x10)^2)/(24*r^2*x20^2)+((6*a*sqrt(1-r^2)-6*m*s+sqrt(1-r^2)*s^2)*(x1-x10)*(x2-x20))/(12*r^2*s*x20^2)-((12*a-12*m*sqrt(1-r^2)*s+s^2)*(x2-x20)^2)/(24*r^2*s^2*x20^2)-((x1-x10)*(2*a*sqrt(1-r^2)-2*m*s-2*b*sqrt(1-r^2)*x20+s*x20))/(2*r^2*s*x20)+((x2-x20)*(2*a-2*m*sqrt(1-r^2)*s-2*b*x20+sqrt(1-r^2)*s*x20))/(2*r^2*s^2*x20)
# c1 <- (s*((-a)*sqrt(1-r^2)+m*s)*(x1-x10))/(12*r^2*x20^2)+(1/(2880*r^4*x20^3))*((x1-x10)^2*(-60*a^2*(-3+2*r^2)+180*m^2*s^2+2*r^2*s^4-60*a*s*(6*m*sqrt(1-r^2)+r^2*s)-60*b^2*x20^2+60*b*sqrt(1-r^2)*s*x20^2-15*s^2*x20^2))+(1/(48*r^2*s^2*x20^2))*((x2-x20)*(12*a^2+12*m^2*s^2-4*m*sqrt(1-r^2)*s^3+2*r^2*s^4-4*a*s*(6*m*sqrt(1-r^2)+(-1+4*r^2)*s)-12*b^2*x20^2+12*b*sqrt(1-r^2)*s*x20^2-3*s^2*x20^2))+(1/(2880*r^4*s^2*x20^3))*((x2-x20)^2*(-180*a^2*(-1+4*r^2)-60*m^2*(-3+10*r^2)*s^2+240*m*r^2*sqrt(1-r^2)*s^3+2*r^2*s^4-96*r^4*s^4+60*a*s*(2*m*sqrt(1-r^2)*(-3+10*r^2)+r^2*(-5+14*r^2)*s)-60*b^2*x20^2+120*b^2*r^2*x20^2+60*b*sqrt(1-r^2)*s*x20^2-120*b*r^2*sqrt(1-r^2)*s*x20^2-15*s^2*x20^2+30*r^2*s^2*x20^2))+(1/(1440*r^4*s*x20^3))*((x1-x10)*(x2-x20)*(60*a^2*sqrt(1-r^2)*(-3+2*r^2)-180*m^2*sqrt(1-r^2)*s^2-120*m*r^2*s^3-2*r^2*sqrt(1-r^2)*s^4+180*a*s*(-2*m*(-1+r^2)+r^2*sqrt(1-r^2)*s)+60*b^2*sqrt(1-r^2)*x20^2-60*b*s*x20^2+60*b*r^2*s*x20^2+15*sqrt(1-r^2)*s^2*x20^2))-(1/(24*r^2*s^2*x20))*(12*a^2+12*m^2*s^2+2*r^2*s^4-12*m*s*(-2*b*sqrt(1-r^2)+s)*x20+12*b^2*x20^2-12*b*sqrt(1-r^2)*s*x20^2+3*s^2*x20^2-12*a*(2*m*sqrt(1-r^2)*s+r^2*s^2+2*b*x20-sqrt(1-r^2)*s*x20))
# c2 <- -((60*a^2*(1+r^2)+60*m^2*s^2+23*r^2*s^4-120*a*s*(m*sqrt(1-r^2)+r^2*s))/(720*r^2*x20^2))
Dv <- (1/2)*log(1 - r^2) + log(s) + log(x2)
cm1 <- (7*s^4*(x1 - x10)^6)/(11520*(-1 + r^2)^3*x20^5) - (7*r*s^3*(x1 - x10)^5*(x2 - x20))/(1920*(-1 + r^2)^3*x20^5) + ((-193 + 298*r^2)*s^2*(x1 - x10)^4*(x2 - x20)^2)/(11520*(-1 + r^2)^3*x20^5) + (r*(193 - 228*r^2)*s*(x1 - x10)^3*(x2 - x20)^3)/(2880*(-1 + r^2)^3*x20^5) + ((745 - 2648*r^2 + 2008*r^4)*(x1 - x10)^2*(x2 - x20)^4)/(11520*(-1 + r^2)^3*x20^5) + ((-745*r + 1876*r^3 - 1152*r^5)*(x1 - x10)*(x2 - x20)^5)/(5760*(-1 + r^2)^3*s*x20^5) + ((945 - 2090*r^2 + 1152*r^4)*(x2 - x20)^6)/(11520*(-1 + r^2)^3*s^2*x20^5) - (s^2*(x1 - x10)^4*(x2 - x20))/(64*(-1 + r^2)^2*x20^4) + (r*s*(x1 - x10)^3*(x2 - x20)^2)/(16*(-1 + r^2)^2*x20^4) + ((3 - 6*r^2)*(x1 - x10)^2*(x2 - x20)^3)/(32*(-1 + r^2)^2*x20^4) + (r*(-3 + 4*r^2)*(x1 - x10)*(x2 - x20)^4)/(16*(-1 + r^2)^2*s*x20^4) + ((7 - 8*r^2)*(x2 - x20)^5)/(64*(-1 + r^2)^2*s^2*x20^4) + (s^2*(x1 - x10)^4)/(96*(-1 + r^2)^2*x20^3) - (r*s*(x1 - x10)^3*(x2 - x20))/ (24*(-1 + r^2)^2*x20^3) + ((-7 + 10*r^2)*(x1 - x10)^2*(x2 - x20)^2)/ (48*(-1 + r^2)^2*x20^3) + ((7*r - 8*r^3)*(x1 - x10)*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s*x20^3) + ((-15 + 16*r^2)*(x2 - x20)^4)/ (96*(-1 + r^2)^2*s^2*x20^3) - ((x1 - x10)^2*(x2 - x20))/(4*(-1 + r^2)*x20^2) + (r*(x1 - x10)*(x2 - x20)^2)/(2*(-1 + r^2)*s*x20^2) - (x2 - x20)^3/(4*(-1 + r^2)*s^2*x20^2) + (s^2*(x1 - x10)^2 - 2*r*s*(x1 - x10)*(x2 - x20) + (x2 - x20)^2)/ (2*(-1 + r^2)*s^2*x20)
c0 <- (7*s^4*(x1 - x10)^4)/(1920*(-1 + r^2)^2*x20^4) - (s*(30*a*r + s*(-30*m + 7*r*s))*(x1 - x10)^3*(x2 - x20))/ (480*(-1 + r^2)^2*x20^4) + ((540*a*r^2 + s*(-540*m*r + (-97 + 160*r^2)*s))* (x1 - x10)^2*(x2 - x20)^2)/(2880*(-1 + r^2)^2*x20^4) + ((-270*a*r*(-1 + 2*r^2) + s*(270*m*(-1 + 2*r^2) + r*(97 - 118*r^2)*s))* (x1 - x10)*(x2 - x20)^3)/(1440*(-1 + r^2)^2*s*x20^4) + ((360*a*(-4 + 5*r^2) + s*(-360*m*r*(-3 + 4*r^2) + (-215 + 236*r^2)*s))* (x2 - x20)^4)/(5760*(-1 + r^2)^2*s^2*x20^4) + (s*(a*r - m*s)*(x1 - x10)^3)/ (24*(-1 + r^2)^2*x20^3) + ((-3*a*r^2 + s*(3*m*r + s - r^2*s))* (x1 - x10)^2*(x2 - x20))/(24*(-1 + r^2)^2*x20^3) + ((a*r*(-7 + 10*r^2) + s*(m*(7 - 10*r^2) + 2*r*(-1 + r^2)*s))*(x1 - x10)* (x2 - x20)^2)/(24*(-1 + r^2)^2*s*x20^3) + ((a*(8 - 9*r^2) + s*(m*r*(-7 + 8*r^2) + s - r^2*s))*(x2 - x20)^3)/ (24*(-1 + r^2)^2*s^2*x20^3) + (s^2*(x1 - x10)^2)/(24*(-1 + r^2)*x20^2) + ((12*a + s*(-12*m*r + s))*(x2 - x20)^2)/(24*(-1 + r^2)*s^2*x20^2) - ((x1 - x10)*(2*a*r - 2*m*s - 2*b*r*x20 + s*x20))/ (2*s*x20 - 2*r^2*s*x20) + ((x2 - x20)*(2*a - 2*m*r*s - 2*b*x20 + r*s*x20))/(2*s^2*x20 - 2*r^2*s^2*x20) + ((6*a*r - 6*m*s + r*s^2)*(x1 - x10)*(x2 - x20))/ (12*s*x20^2 - 12*r^2*s*x20^2)
c1 <- (s*(a*r - m*s)*(x1 - x10))/(12*(-1 + r^2)*x20^2) + ((x1 - x10)^2*(60*a^2*(1 + 2*r^2) + 180*m^2*s^2 + 2*s^4 - 2*r^2*s^4 - 60*a*s*(6*m*r + s - r^2*s) - 60*b^2*x20^2 + 60*b*r*s*x20^2 - 15*s^2*x20^2))/(2880*(-1 + r^2)^2*x20^3) + ((x2 - x20)*(-12*a^2 - 12*m^2*s^2 + 4*m*r*s^3 - 2*s^4 + 2*r^2*s^4 + 4*a*s*(6*m*r + (3 - 4*r^2)*s) + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2))/(48*(-1 + r^2)*s^2*x20^2) + (1/(2880*(-1 + r^2)^2*s^2*x20^3))* ((x2 - x20)^2*(180*a^2*(-3 + 4*r^2) + 60*m^2*(-7 + 10*r^2)*s^2 - 240*m*r*(-1 + r^2)*s^3 - 94*s^4 + 190*r^2*s^4 - 96*r^4*s^4 + 60*a*s*(m*(14*r - 20*r^3) + (9 - 23*r^2 + 14*r^4)*s) + 60*b^2*x20^2 - 120*b^2*r^2*x20^2 - 60*b*r*s*x20^2 + 120*b*r^3*s*x20^2 + 15*s^2*x20^2 - 30*r^2*s^2*x20^2)) + (1/(24*(-1 + r^2)*s^2*x20))* (12*a^2 + 12*m^2*s^2 + 2*s^4 - 2*r^2*s^4 + 12*m*(2*b*r - s)*s* x20 + 12*b^2*x20^2 - 12*b*r*s*x20^2 + 3*s^2*x20^2 - 12*a*(2*m*r*s + s^2 - r^2*s^2 + 2*b*x20 - r*s*x20)) + (1/(1440*(-1 + r^2)^2*s*x20^3))*((x1 - x10)*(x2 - x20)*(-60*a^2*(r + 2*r^3) - 180*m^2*r*s^2 + 120*m*(-1 + r^2)*s^3 + 180*a*r*s* (2*m*r + s - r^2*s) + r*(2*(-1 + r^2)*s^4 + 60*b^2*x20^2 - 60*b*r*s*x20^2 + 15*s^2*x20^2)))
c2 <- -((60*a^2*(-2 + r^2) - 60*m^2*s^2 + 23*(-1 + r^2)*s^4 + 120*a*s*(m*r + s - r^2*s))/(720*(-1 + r^2)*x20^2))
output <- list()
output$llk <- -log(2*del*pi) - Dv + cm1/del + c0 + del*c1 + (1/2)*del^2*c2
# print(paste('r',r,sep='='))
# print(paste('s',s,sep='='))
# print(paste('x2',x2,sep='='))
# a1 <- param[1]
# b1 <- -0.5
# a2 <- param[2]
# b2 <- param[3]
# s <- param[4]
# r <- param[5]
#
# Dv <- 0.5*log((1-r^2)*s*x2)
# c1 <- (r*s*a2*b1+r*s*a1*b2-s^2*a1*b1-a2*b2)/((1-r^2)*s^2) + (2*r*s*b1*b2-s^2*b1^2-b2^2)*x20/(2*(1-r^2)*s^2) - (s^4-r^2*s^4+6*s^2*a1^2-6*s^2*a2+6*r^2*s^2*a2 - 12*r*s*a1*a2 +6*a2^2)/(12*(1-r^2)*s^2*x20)
# c0 <- (x2-x20)*(x20*b2-r*s*x20*b1 -r*s*a1+a2)/((1-r^2)*s^2*x20) -(x1-x10)*(r*x20*b2-s*x20*b1-s*a1+r*a2)/((1-r^2)*s*x20) -s^2*(x1-x10)^2/(24*(1-r^2)*x20^2) - (x2-x20)^2*(s*(s-12*r*a1) + 12*a2)/(24*(1-r^2)*s^2*x20^2) + (x2-x20)*(x1-x10)*(r*s^2-6*s*a1 + 6*r*a2)/(12*s*x20^2*(1-r^2))
# cm1 <- -((x2-x20)^2 - 2*r*s*(x2-x20)*(x1-x10) + s^2*(x1-x10)^2)/(2*(1-r^2)*s^2*x20) + (x2-x20)^3/(4*(1-r^2)*s^2*x20^2) - r*(x2-x20)^2*(x1-x10)/(2*(1-r^2)*s*x20^2) + (x2-x20)*(x1-x10)^2/(4*(1-r^2)*x20^2) + (7*r-8*r^3)*(x2-x20)^3*(x1-x10)/(24*(1-r^2)^2*s*x20^3) - (7-10*r^2)*(x2-x20)^2*(x1-x10)^2/(48*(1-r^2)^2*x20^3) - r*s*(x2-x20)*(x1-x10)^3/(24*(1-r^2)^2*x20^3) + (s^2*(x1-x10)^4)/(96*(1-r^2)^2*x20^3) - (15-16*r^2)*(x2-x20)^4/(96*(1-r^2)^2*s^2*x20^3)
# output <- list()
# output$llk <- -log(2*del*pi) - Dv + cm1/del + c0 + del*c1
# print(paste('Dv',Dv,sep='='))
# print(paste('c1',c1,sep='='))
# print(paste('c0',c0,sep='='))
# print(paste('cm1',cm1,sep='='))
return(output)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.