In Jaimemosg/EstimationTools: Maximum Likelihood Estimation for Probability Functions from Data Sets

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Computation of standard errors

Objects of maxlogL class (outputs from maxlogL and maxlogLreg) stores the estimated parameters of probability density/mass functions by Maximum Likelihood. The variance-covariance matrix is computed from Fisher information matrix, which is obtained by means of the Inverse Hessian matrix of estimators:

\begin{equation} Var(\hat{\boldsymbol{\theta}}) = \mathcal{J}^{-1}(\hat{\boldsymbol{\theta}}) = C(\hat{\boldsymbol{\theta}}), \end{equation}

where $\mathcal{J}(\hat{\boldsymbol{\theta}})$ is the observed Fisher Information Matrix. Hence, the standard errors can be calculated as the square root of the diagonal elements of matrix $C$, as follows:

\begin{equation} SE(\hat{\boldsymbol{\theta}}) = \sqrt{C_{jj}(\hat{\boldsymbol{\theta}})}, \end{equation}

To install the package, type the following commands:

if (!require('devtools')) install.packages('devtools')
devtools::install_github('Jaimemosg/EstimationTools', force = TRUE)

In EstimationTools Hessian matrix is computed in the following way:

• If StdE_Method = optim, it is estimated through the optim function (with option hessian = TRUE under the hood in maxlogL or maxlogLreg function).
• If the previous implementation fails or if the user chooses StdE_Method = numDeriv, it is calculated with hessian function from numDeriv package.
• If both of the previous methods fail, then standard errors are computed by bootstrapping with the function bootstrap_maxlogL.

Additionally, EstimationTools allows implementation of bootstrap for standard error estimation, even if the Hessian computation does not fail.

Standard Error with maxlogL function

Lets fit the following distribution:

$$ \begin{aligned} X &\sim N(\mu, \:\sigma^2) \ \mu &= 160 \quad (\verb|mean|) \ \sigma &= 6 \quad (\verb|sd|) \end{aligned} $$

The following chunk illustrates the fitting with Hessian computation via optim:

library(EstimationTools)

x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                   link = list(over = "sd", fun = "log_link"),
                   fixed = list(mean = 160))
summary(theta_1)

## Hessian
print(theta_1$fit$hessian)

## Standard errors
print(theta_1$fit$StdE)
print(theta_1$outputs$StdE_Method)

a <- theta_1$fit$StdE

Note that Hessian was computed with no issues. Now, lets check the aforementioned feature in maxlogL: the user can implement bootstrap algorithm available in bootstrap_maxlogL function. To illustrate this, we are going to create another object theta_2:

# Bootstrap
theta_2 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                   link = list(over = "sd", fun = "log_link"),
                   fixed = list(mean = 160))
bootstrap_maxlogL(theta_2, R = 200)
summary(theta_2)

## Hessian
print(theta_2$fit$hessian)

## Standard errors
print(theta_2$fit$StdE)
print(theta_2$outputs$StdE_Method)

b <- theta_2$fit$StdE

Notice that Standard Errors calculated with optim ($r round(a, 6)$) and those calculated with bootstrap implementation ($r round(b, 6)$) are approximately equals, but no identical.

Jaimemosg/EstimationTools documentation built on Oct. 23, 2023, 10 a.m.