# -------------------------------
# Bootstrap for the probit model
# -------------------------------
# \donttest{
# ---
# Step 1
# Simulation of data
# ---
# Load required package
library("mnorm")
# Set seed for reproducibility
set.seed(123)
# The number of observations
n <- 100
# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
# Random errors
u <- rnorm(n = n, mean = 0, sd = 1)
# Coefficients
gamma <- c(-1, 2)
# Linear index
li <- gamma[1] * w1 + gamma[2] * w2
# Latent variable
z_star <- li + u
# Cuts
cuts <- c(-1, 0.5, 2)
# Observable ordered outcome
z <- rep(0, n)
z[(z_star > cuts[1]) & (z_star <= cuts[2])] <- 1
z[(z_star > cuts[2]) & (z_star <= cuts[3])] <- 2
z[z_star > cuts[3]] <- 3
table(z)
# Data
data <- data.frame(w1 = w1, w2 = w2, z = z)
# ---
# Step 2
# Estimation of the parameters
# ---
# Estimation
model <- msel(formula = z ~ w1 + w2, data = data)
summary(model)
# ---
# Step 3
# Bootstrap
# ---
# Perform bootstrap
bootstrap <- bootstrap_msel(model)
# Test the hypothesis that H0: gamma[2] = -2gamma[1]
# by using the t-test and with bootstrap p-values
fn_test <- function(object)
{
gamma <- coef(object, eq = 1)
return(gamma[2] + 2 * gamma[1])
}
b <- test_msel(object = model,
fn = fn_test,
test = "t",
method = "bootstrap",
ci = "percentile",
se_type = "bootstrap",
bootstrap = bootstrap)
summary(b)
# Replicate the analysis with the additional 20 bootstrap iterations
bootstrap2 <- bootstrap_msel(model, iter = 20)
bootstrap_new <- bootstrap_combine_msel(bootstrap, bootstrap2)
b2 <- test_msel(object = model,
fn = fn_test,
test = "t",
method = "bootstrap",
ci = "percentile",
se_type = "bootstrap",
bootstrap = bootstrap)
summary(b2)
# }
Run the code above in your browser using DataLab