Statistics 5601 (Geyer, Fall 2013) The Sandwich Estimator

General Instructions

To do each example, just click the Submit button. You do not have to type in any R instructions or specify a dataset. That's already done for you.

Theory

This page goes with notes, the PDF form of which is found on the course notes page.

We just repeat the R computations in that note in Rweb form here.

The Cauchy Location Problem

The Cauchy Location Problem with Sandwich

R statements logl <- expression(- log(1 + (x - theta)^2)) loglfun <- deriv3(logl, "theta", c("theta", "x")) loglfun theta.start <- median(x) ##### does work vectorwise ##### loglfun(theta.start, x = x) fred <- function(theta, x) { sally <- loglfun(theta, x) result <- sum(- sally) attr(result, "gradient") <- sum(- attr(sally, "gradient")) return(result) } out <- nlm(fred, median(x), x = x) print(out$estimate) print(out$gradient) vfun <- function(theta, x) { sally <- loglfun(theta, x) sum(attr(sally, "gradient")^2) } jfun <- function(theta, x) { sally <- loglfun(theta, x) sum(- attr(sally, "hessian")) } theta.hat <- out$estimate Vhat <- vfun(theta.hat, x) Jhat <- jfun(theta.hat, x) print(Vhat) print(Jhat) conf.level <- 0.95 crit <- qnorm((1 + conf.level) / 2) theta.hat + crit * c(-1, 1) * sqrt(Vhat / Jhat^2)

Dataset URL

For security reasons, Rweb forbids the use of the eval function. This makes the code in the form above much messier than the code in the notes, but both calculate the same things.

The Cauchy Location-Scale Problem with Sandwich

R statements logl <- expression(- log(1 + ((x - mu) / sigma)^2) - log(sigma)) loglfun <- deriv3(logl, c("mu", "sigma"), c("mu", "sigma", "x")) loglfun mu.start <- median(x) sigma.start <- IQR(x) / 2 ##### does work vectorwise ##### loglfun(mu.start, sigma.start, x = x) fred <- function(theta, x) { mu <- theta[1] sigma <- theta[2] sally <- loglfun(mu, sigma, x) result <- sum(- sally) attr(result, "gradient") <- apply(- attr(sally, "gradient"), 2, sum) return(result) } theta.start <- c(mu.start, sigma.start) out <- nlm(fred, theta.start, x = x) print(out$estimate) print(out$gradient) vfun <- function(theta, x) { mu <- theta[1] sigma <- theta[2] sally <- loglfun(mu, sigma, x) herman <- attr(sally, "gradient") t(herman) %*% herman } jfun <- function(theta, x) { mu <- theta[1] sigma <- theta[2] sally <- loglfun(mu, sigma, x) apply(- attr(sally, "hessian"), c(2, 3), sum) } theta.hat <- out$estimate Vhat <- vfun(theta.hat, x) Jhat <- jfun(theta.hat, x) Mhat <- solve(Jhat) %*% Vhat %*% solve(Jhat) print(Vhat) print(Jhat) print(Mhat) conf.level <- 0.95 crit <- qnorm((1 + conf.level) / 2) theta.hat[1] + crit * c(-1, 1) * sqrt(Mhat[1, 1]) theta.hat[2] + crit * c(-1, 1) * sqrt(Mhat[2, 2])