if(!requireNamespace("lattice") || !requireNamespace("boot")) q()
###################################################
### chunk number 1: setup
###################################################
options(prompt = "R> ", continue = "+ ", width = 64,
digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE)
options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))},
twofig = function() {par(mfrow = c(1,2))},
threefig = function() {par(mfrow = c(1,3))},
fourfig = function() {par(mfrow = c(2,2))},
sixfig = function() {par(mfrow = c(3,2))}))
library("AER")
suppressWarnings(RNGversion("3.5.0"))
set.seed(1071)
###################################################
### chunk number 2: DGP
###################################################
dgp <- function(nobs = 15, model = c("trend", "dynamic"),
corr = 0, coef = c(0.25, -0.75), sd = 1)
{
model <- match.arg(model)
coef <- rep(coef, length.out = 2)
err <- as.vector(filter(rnorm(nobs, sd = sd), corr,
method = "recursive"))
if(model == "trend") {
x <- 1:nobs
y <- coef[1] + coef[2] * x + err
} else {
y <- rep(NA, nobs)
y[1] <- coef[1] + err[1]
for(i in 2:nobs)
y[i] <- coef[1] + coef[2] * y[i-1] + err[i]
x <- c(0, y[1:(nobs-1)])
}
return(data.frame(y = y, x = x))
}
###################################################
### chunk number 3: simpower
###################################################
simpower <- function(nrep = 100, size = 0.05, ...)
{
pval <- matrix(rep(NA, 2 * nrep), ncol = 2)
colnames(pval) <- c("dwtest", "bgtest")
for(i in 1:nrep) {
dat <- dgp(...)
pval[i,1] <- dwtest(y ~ x, data = dat,
alternative = "two.sided")$p.value
pval[i,2] <- bgtest(y ~ x, data = dat)$p.value
}
return(colMeans(pval < size))
}
###################################################
### chunk number 4: simulation-function
###################################################
simulation <- function(corr = c(0, 0.2, 0.4, 0.6, 0.8,
0.9, 0.95, 0.99), nobs = c(15, 30, 50),
model = c("trend", "dynamic"), ...)
{
prs <- expand.grid(corr = corr, nobs = nobs, model = model)
nprs <- nrow(prs)
pow <- matrix(rep(NA, 2 * nprs), ncol = 2)
for(i in 1:nprs) pow[i,] <- simpower(corr = prs[i,1],
nobs = prs[i,2], model = as.character(prs[i,3]), ...)
rval <- rbind(prs, prs)
rval$test <- factor(rep(1:2, c(nprs, nprs)),
labels = c("dwtest", "bgtest"))
rval$power <- c(pow[,1], pow[,2])
rval$nobs <- factor(rval$nobs)
return(rval)
}
###################################################
### chunk number 5: simulation
###################################################
set.seed(123)
psim <- simulation()
###################################################
### chunk number 6: simulation-table
###################################################
tab <- xtabs(power ~ corr + test + model + nobs, data = psim)
ftable(tab, row.vars = c("model", "nobs", "test"),
col.vars = "corr")
###################################################
### chunk number 7: simulation-visualization
###################################################
library("lattice")
xyplot(power ~ corr | model + nobs, groups = ~ test,
data = psim, type = "b")
###################################################
### chunk number 8: simulation-visualization1
###################################################
library("lattice")
trellis.par.set(theme = canonical.theme(color = FALSE))
print(xyplot(power ~ corr | model + nobs, groups = ~ test, data = psim, type = "b"))
###################################################
### chunk number 9: journals-lm
###################################################
data("Journals")
journals <- Journals[, c("subs", "price")]
journals$citeprice <- Journals$price/Journals$citations
jour_lm <- lm(log(subs) ~ log(citeprice), data = journals)
###################################################
### chunk number 10: journals-residuals-based-resampling-unused eval=FALSE
###################################################
## refit <- function(data, i) {
## d <- data
## d$subs <- exp(d$fitted + d$res[i])
## coef(lm(log(subs) ~ log(citeprice), data = d))
## }
###################################################
### chunk number 11: journals-case-based-resampling
###################################################
refit <- function(data, i)
coef(lm(log(subs) ~ log(citeprice), data = data[i,]))
###################################################
### chunk number 12: journals-boot
###################################################
library("boot")
set.seed(123)
jour_boot <- boot(journals, refit, R = 999)
###################################################
### chunk number 13: journals-boot-print
###################################################
jour_boot
###################################################
### chunk number 14: journals-lm-coeftest
###################################################
coeftest(jour_lm)
###################################################
### chunk number 15: journals-boot-ci
###################################################
boot.ci(jour_boot, index = 2, type = "basic")
###################################################
### chunk number 16: journals-lm-ci
###################################################
confint(jour_lm, parm = 2)
###################################################
### chunk number 17: ml-loglik
###################################################
data("Equipment", package = "AER")
nlogL <- function(par) {
beta <- par[1:3]
theta <- par[4]
sigma2 <- par[5]
Y <- with(Equipment, valueadded/firms)
K <- with(Equipment, capital/firms)
L <- with(Equipment, labor/firms)
rhs <- beta[1] + beta[2] * log(K) + beta[3] * log(L)
lhs <- log(Y) + theta * Y
rval <- sum(log(1 + theta * Y) - log(Y) +
dnorm(lhs, mean = rhs, sd = sqrt(sigma2), log = TRUE))
return(-rval)
}
###################################################
### chunk number 18: ml-0
###################################################
fm0 <- lm(log(valueadded/firms) ~ log(capital/firms) +
log(labor/firms), data = Equipment)
###################################################
### chunk number 19: ml-0-coef
###################################################
par0 <- as.vector(c(coef(fm0), 0, mean(residuals(fm0)^2)))
###################################################
### chunk number 20: ml-optim
###################################################
opt <- optim(par0, nlogL, hessian = TRUE)
###################################################
### chunk number 21: ml-optim-output
###################################################
opt$par
sqrt(diag(solve(opt$hessian)))[1:4]
-opt$value
###################################################
### chunk number 22: Sweave eval=FALSE
###################################################
## Sweave("Sweave-journals.Rnw")
###################################################
### chunk number 23: Stangle eval=FALSE
###################################################
## Stangle("Sweave-journals.Rnw")
###################################################
### chunk number 24: texi2dvi eval=FALSE
###################################################
## texi2dvi("Sweave-journals.tex", pdf = TRUE)
###################################################
### chunk number 25: vignette eval=FALSE
###################################################
## vignette("Sweave-journals", package = "AER")
Generated by dwww version 1.15 on Fri Jun 21 07:46:44 CEST 2024.