Third order Delta method.R

# Third order Delta method 

bias_adj <- function(mu, sigma, fun){
  first_deriv <- pracma::fderiv(fun, mu)
  second_deriv <- pracma::fderiv(fun, mu, n = 2)
  third_deriv <- pracma::fderiv(fun, mu, n = 3)
  moment1 <- fun(mu)+ # first order
    1/2*second_deriv*sigma^2 # second order
  # the third order part evaluates to zero due to symmetry
  moment2 <- sqrt((first_deriv*sigma)^2+ # first order
                    1/2*(second_deriv*sigma^2)^2+ # second order
                    5/12*third_deriv^2*sigma^6+first_deriv*third_deriv*sigma^4) # third order
  return(c("mean" = moment1, "sd" = moment2))
}