Statistics 5601 (Geyer, Spring 2006) Examples: Sandwich Estimator
Contents
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.
This page goes with notes, the PDF form
of which is found on the theory page.
We just repeat the R computations in that note in Rweb form here.
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)
External Data Entry
Enter a 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.
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])
External Data Entry
Enter a dataset URL :
Again, 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.