Monte Carlo

Generates R number of parameter estimates from a vector of thetahat \(\left( \boldsymbol{\hat{\theta} }\right)\) and sampling variances and covariances of thetahat \(\widehat{\mathrm{Cov}} \left( \boldsymbol{\hat{\theta}} \right)\). If length of thetahat is greater than 1, a multivariate normal distribution is assumed. If length of thetahat is equal 1, a univariate normal distribution is assumed.

mc(thetahat, vcovhat, R = 20000L)

Arguments

thetahat	Numeric vector. Parameter estimates \(\boldsymbol{\hat{\theta}}\).
vcovhat	Symmetric matrix. Estimated sampling variances and covariances of parameter estimates \(\widehat{\mathrm{Cov}} \left( \boldsymbol{\hat{\theta}} \right)\).
R	Integer. Number of Monte Carlo replications.

Examples

R <- 20000L
# length 1
thetahat <- 100
vcovhat <- 1.5
mc_star_length_1 <- mc(
  thetahat = thetahat,
  vcovhat = vcovhat,
  R = R
)
str(mc_star_length_1)
#>  num [1:20000] 102 101 102 101 100 ...
hist(
  mc_star_length_1,
  main = expression(
    paste(
      "Histogram of ",
      hat(theta),
      "*"
    )
  ),
  xlab = expression(
    paste(
      hat(theta),
      "*"
    )
  )
)
qqnorm(mc_star_length_1)
qqline(mc_star_length_1)
# length greater than 1
alphahat <- 0.3386
betahat <- 0.4510
alphahat_betahat <- alphahat * betahat
varhat_alphahat <- 0.1224^2
varhat_betahat <- 0.1460^2
thetahat <- c(
  alphahat,
  betahat
)
vcovhat <- matrix(
  data = c(
    varhat_alphahat,
    0.00,
    0.00,
    varhat_betahat
  ),
  ncol = 2
)
mc_star_length_2 <- mc(
  thetahat = thetahat,
  vcovhat = vcovhat,
  R = R
)
str(mc_star_length_2)
#>  num [1:20000, 1:2] 0.387 0.395 0.392 0.319 0.499 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : NULL
#>   ..$ : NULL
alphahat_betahat_star <- mc_star_length_2[, 1] * mc_star_length_2[, 2]
hist(
  alphahat_betahat_star,
  main = expression(
    paste(
      "Histogram of ",
      hat(alpha),
      hat(beta),
      "*"
    )
  ),
  xlab = expression(
    paste(
      hat(alpha),
      hat(beta),
      "*"
    )
  )
)
qqnorm(alphahat_betahat_star)
qqline(alphahat_betahat_star)
wald(
  thetahat = alphahat_betahat,
  sehat = sd(alphahat_betahat_star)
)
#>    statistic            p           se      ci_0.05       ci_0.5       ci_2.5 
#>  1.998726138  0.045637993  0.076402963 -0.098697394 -0.044092392  0.002961543 
#>      ci_97.5      ci_99.5     ci_99.95 
#>  0.302455657  0.349509592  0.404114594