# R/hessian.R In LSTS: Locally Stationary Time Series

#### Documented in hessian

#' @title Hessian Matrix
#' @description Numerical aproximation of the Hessian of a function.
#' @details
#' Computes the numerical approximation of the Hessian of \code{f}, evaluated at
#' \code{x0}.
#' Usually needs to pass additional parameters (e.g. data).  N.B. this uses no
#' numerical sophistication.
#' @param f (type: numeric) name of function that defines log likelihood
#' (or negative of it).
#' @param x0 (type: numeric) scalar or vector of parameters that give the point
#' at which you want the hessian estimated (usually will be the mle).
#' @param ... Additional arguments to be passed to the function.
#' @examples
#' # Variance of the maximum likelihood estimator for mu parameter in
#' # gaussian data
#' loglik <- function(series, x, sd = 1) {
#'   -sum(log(dnorm(series, mean = x, sd = sd)))
#' }
#' sqrt(c(var(malleco) / length(malleco), diag(solve(hessian(
#'   f = loglik, x = mean(malleco), series = malleco,
#'   sd = sd(malleco)
#' )))))
#' @return
#' An \eqn{n \times n} matrix of 2nd derivatives, where \eqn{n} is the length of
#' \code{x0}.
#' @export
hessian <- function(f, x0, ...) {
n <- length(x0)
mdelta <- matrix(0., nrow = n, ncol = n)
hess <- mdelta
delta <- 0.0001 * (sign(x0) + (x0 == 0.)) * pmax(abs(x0), 0.01)
diag(mdelta) <- delta
f0 <- f(x0, ...)
for (i in 1.:n) {
grad[i] <- f(x0 + mdelta[, i], ...)
}
for (i in 1.:n) {
for (j in i:n) {
hess[i, j] <- f(x0 + mdelta[, i] + mdelta[, j], ...) - grad[i] - grad[j]
hess[j, i] <- hess[i, j]
}
}
(hess + f0) / outer(delta, delta, "*")
}


## Try the LSTS package in your browser

Any scripts or data that you put into this service are public.

LSTS documentation built on July 29, 2021, 5:07 p.m.