hessian: Hessian Matrix

Description Usage Arguments Value Author(s) Examples

View source: R/hessian.R

Description

Computes numerical approximation to Hessian of f, evaluated at x0. Usually need to pass additional parameters (e.g. data). N.B. this uses no numerical sophistication.

Usage

1
hessian(f, x0, ...)

Arguments

f

name of function that defines log likelihood (or negative of it).

x0

scalar or vector of parameters that give the point at which you want the hessian estimated (usually will be the mle).

...

additional arguments passed to f.

Value

An n x n matrix of 2nd derivatives, where n is the length of x0.

Author(s)

Wilfredo Palma <[email protected]>

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
## Ex 1:  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)))
}
n = 500
series = rnorm(500, mean = 10, sd = 2)
sqrt(c(var(series)/n,diag(solve(hessian(f = loglik, x = mean(series), series = series, 
                          sd = sd(series))))))


## Ex 2:  Variance of the maximum likelihood estimator for phi parameter AR(1)
## in gaussian data
loglik = function(series, x, sd = 1){
n = length(series)
-(sum(log(dnorm(series[2:n], mean = x * series[1:(n-1)], sd = sd))) + 
log(dnorm(series[1], mean = 0, sd = sqrt(sd^2/(1-x^2)))))
}
n = 1000
series = arima.sim(n, model = list(c(1,0,0), ar = 0.7))
fit = arima(series, c(1,0,0), include.mean = FALSE)
round(c(fit$var.coef,diag(solve(hessian(f = loglik, x = fit$coef, series = series,
      sd = sqrt(fit$sigma2))))),6)


## Ex 3:  Variance of the whittle maximum likelihood estimator for phi parameter
## AR(1) in gaussian  data
loglik = function(series, x, sd = 1){
n = length(series)
aux = periodogram(series, plot = FALSE)
lambda = aux$lambda
I = aux$periodogram
f = fdensity(ar = x, sd = sd, lambda = lambda)
lik = sum(log(f)+I/f)
lik = lik/n
}
n = 500
series = arima.sim(n, model = list(c(1,0,0), ar = 0.7))
fit = arima(series, c(1,0,0), include.mean = FALSE)
round(c(fit$var.coef,diag(solve(hessian(f = loglik, x = fit$coef, series = series,
      sd = sqrt(fit$sigma2))))/n),4)

LSTS documentation built on May 29, 2017, 6:36 p.m.