Computation of standard errors

Objects of maxlogL class (outputs frommaxlogL 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}}) = I^{-1}(\hat{\boldsymbol{\theta}}) = C(\hat{\boldsymbol{\theta}}), \end{equation}\]

where \(I(\hat{\boldsymbol{\theta}})\) is the 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:

  • Firstly, it is estimated through optim, with option hessian = TRUE in maxlogL function.
  • If the previous implementation fails, 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))
#>   0:     43502.610:  1.00000
#>   1:     32479.822:  2.00000
#>   2:     32247.813:  1.91870
#>   3:     32106.589:  1.76062
#>   4:     32097.456:  1.79477
#>   5:     32097.298:  1.79091
#>   6:     32097.298:  1.79079
#>   7:     32097.298:  1.79079
summary(theta_1)
#> _______________________________________________________________
#> Optimization routine: nlminb 
#> Standard Error calculation: Hessian from optim 
#> _______________________________________________________________
#>       AIC     BIC
#>   64194.6 64194.6
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   5.99419    0.04239   141.4   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators 
#> ---

## Hessian
print(theta_1$fit$hessian)
#>          [,1]
#> [1,] 556.6325

## Standard errors
print(theta_1$fit$StdE)
#> [1] 0.04238534
print(theta_1$outputs$StdE_Method)
#> [1] "Hessian from optim"

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))
#>   0:     43502.610:  1.00000
#>   1:     32479.822:  2.00000
#>   2:     32247.813:  1.91870
#>   3:     32106.589:  1.76062
#>   4:     32097.456:  1.79477
#>   5:     32097.298:  1.79091
#>   6:     32097.298:  1.79079
#>   7:     32097.298:  1.79079
bootstrap_maxlogL(theta_2, R = 200)
#> 
#> ...Bootstrap computation of Standard Error. Please, wait a few minutes...
#> 
#> 
#>  --> Done <---
summary(theta_2)
#> _______________________________________________________________
#> Optimization routine: nlminb 
#> Standard Error calculation: Bootstrap 
#> _______________________________________________________________
#>       AIC     BIC
#>   64194.6 64194.6
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   5.99419    0.04309   139.1   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators 
#> ---

## Hessian
print(theta_2$fit$hessian)
#>          [,1]
#> [1,] 556.6325

## Standard errors
print(theta_2$fit$StdE)
#> [1] 0.04309105
print(theta_2$outputs$StdE_Method)
#> [1] "Bootstrap"

Notice that Standard Errors calculated with optim (\(0.042385\)) and those calculated with bootstrap implementation (\(0.043091\)) are approximately equals, but no identical.