/
Bootstrap_implementation.Rmd
108 lines (81 loc) · 3.66 KB
/
Bootstrap_implementation.Rmd
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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
---
title: "Bootstrap implementation"
author: "Jaime Mosquera"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Bootstrap implementation}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
editor_options:
chunk_output_type: console
---
```{r setup, include = FALSE}
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:
```{r install, eval=FALSE}
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`:
```{r HessianOptim, warning=FALSE, message=FALSE}
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)
```
```{r Hessian1, echo=FALSE}
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`:
```{r HessianBootstrap}
# 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)
```
```{r Hessian2, echo=FALSE}
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.