esreg: Joint Quantile and Expected Shortfall Regression

Description

Estimates a joint linear regression model for the pair (VaR, ES): $$Q_\alpha(Y | X) = X'\beta_q$$ $$ES_\alpha(Y | X) = X'\beta_e$$

Usage

esreg(formula, data, alpha, g1 = 2L, g2 = 1L, target = "rho",
  shift_data = TRUE, method = "ils", control = list(terminate_after = 10,
  max.time = 10, box = 10))

Arguments

formula

Formula object, e.g.: y ~ x1 + x2 + ...

data

data.frame that holds the variables. If missing the data is extracted from the environment.

alpha

Probability level

1, 2 (see G1_fun, G1_prime_fun). Default is 1.

1, 2, 3, 4, 5 (see G2_curly_fun, G2_fun, G2_prime_fun). Default is 2.

target

The functions to be optimized: either the loss (rho) or the moment function (psi). Optimization of the loss function is strongly recommended.

shift_data

If g2 is 1, 2 or 3, we can either estimate the model without or with shifting of the Y variable. We either risk biased estimates (no shifting) or slightly different estimates due to the changed loss function (shifting). Defaults to shifting to avoid biased estimates.

method

Iterated local search (ils) or Simulated annealing (sa). Defaults to ils for reasons of computation speed.

control

A list with control parameters passed to either the ils or sa:

terminate_after - Stop the iterated local search if there is no improvement within max_step consecutive steps. Default is 10.
max.time - Maximum running time of the sa optimizer.
box - Box around the parameters for the sa optimizer.

Value

An esreg object

References

A Joint Quantile and Expected Shortfall Regression Framework

Examples

Run this code

# NOT RUN {
# Simulate data (DGP-(2) in the linked paper)
set.seed(0)
x <- rchisq(2000, df=1)
y <- -x + (1 + 0.5 * x) * rnorm(1000)

# True quantile and expected shortfall regression parameters (for alpha=0.025)
alpha=0.025
true_pars <- c(-1.959964, -1.979982, -2.337803, -2.168901)

# Estimate the model using the standard settings
fit <- esreg(y ~ x, alpha=alpha)

# Compare the different variance-covariance estimators
cov1 <- vcov(object=fit, sparsity="iid", cond_var="ind")
cov2 <- vcov(object=fit, sparsity="nid", cond_var="scl_N")
cov3 <- vcov(object=fit, sparsity="nid", cond_var="scl_sp")

print("Comparison of the variance-covariance estimators")
print(rbind(Truth=true_pars,
            Estimate=coef(fit),
            SE_iid_ind=sqrt(diag(cov1)),
            SE_ind_N=sqrt(diag(cov2)),
            SE_ind_sp=sqrt(diag(cov3))))

# Compares estimates using different G2 functions
fit1 <- esreg(y ~ x, alpha=alpha, g2=1)
fit2 <- esreg(y ~ x, alpha=alpha, g2=2)
fit3 <- esreg(y ~ x, alpha=alpha, g2=3)
fit4 <- esreg(y ~ x, alpha=alpha, g2=4)
fit5 <- esreg(y ~ x, alpha=alpha, g2=5)
fits <- sapply(list(fit1, fit2, fit3, fit4, fit5), coef)
colnames(fits) <- sapply(1:5, function(i) esreg:::.G_function_names(1, i)[2])
print("Comparison of the five G2 functions")
print(rbind(Truth=true_pars, t(fits)))

# Compare the M- and Z-estimator
fit_m <- esreg(y ~ x, alpha=alpha, target="rho")
fit_z <- esreg(y ~ x, alpha=alpha, target="psi")
print("Comparison of the M- and Z-estimator")
print(t(cbind(Truth=true_pars, M=coef(fit_m), Z=coef(fit_z))))
# }