t.ols <- beta2hat/se.ols
t.ols
ifelse(t.ols > t.crit, "Reject", "Fail to reject")
t.white <- beta2hat/se.white
t.white
ifelse(t.white > t.crit, "Reject", "Fail to reject")
t.jackknife <- beta2hat/se.jackknife
t.jackknife
ifelse(t.jackknife > t.crit, "Reject", "Fail to reject")
t.bootstrap <- beta2hat/se.bootstrap
t.bootstrap
ifelse(t.bootstrap > t.crit, "Reject", "Fail to reject")
## ----brksbss,keep.source=TRUE--------------------------------------------
require(bootstrap)
set.seed(42)
n <- 10
M <- 1000
bias.jack.ml <- numeric(M)
bias.jack.ls <- numeric(M)
for(m in 1:M) {
x <- rnorm(n)
## Jackknife the ML and LS estimators of the sample variance.
var.jack.ml <- function(mydat) {
sum((mydat-mean(mydat))^2)/length(mydat)
}
var.jack.ls <- function(mydat) {
sum((mydat-mean(mydat))^2)/(length(mydat)-1)
}
output.jack.ml <- jackknife(x,var.jack.ml)
bias.jack.ml[m] <- output.jack.ml$jack.bias
output.jack.ls <- jackknife(x,var.jack.ls)
bias.jack.ls[m] <- output.jack.ls$jack.bias
}
## The true bias of the ML estimator is -1/n...
true.bias <- -1/n
true.bias
mean(bias.jack.ml)
## The true bias of the OLS estimator is 0...
mean(bias.jack.ls)
# ## (the first meq constraints are treated as equality constraints,
## ## all further as inequality constraints).
## b <- solve.QP(Dmat=D,dvec=d,Amat=A,bvec=b0,meq=1)$solution
## ----qpexam--------------------------------------------------------------
require(quadprog)
## Assume we want to minimize: -(0 5 0) %*% b + 1/2 b^T b
## under the constraints:      A^T b >= b0
## with b0 = (-8,2,0)^T
## and      (-4  2  0)
##      A = (-3  1 -2)
##          ( 0  0  1)
## we can use solve.QP as follows:
##
Dmat       <- matrix(0,3,3)
diag(Dmat) <- 1
dvec       <- c(0,5,0)
Amat       <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
bvec       <- c(-8,2,0)
solve.QP(Dmat,dvec,Amat,bvec)$solution
## # Problem Set {-}
## ----qandqibmclose,fig.asp=0.5-------------------------------------------
par(mfrow=c(1,3),cex=1/0.83)
require(fpp2)
data(ibmclose)
plot(ibmclose)
acf(ibmclose)
ndiffs(ibmclose)
## ----rceas---------------------------------------------------------------
set.seed(42)
y <- arima.sim(n=100,list(ar=0.6,sd=1))
plot(y,ylab="$y_t$")
## ----gtwas---------------------------------------------------------------
set.seed(42)
y1 <- arima.sim(n=100,list(ar=0.6,ma=0.6,sd=1))
y2 <- arima.sim(n=100,list(ar=c(-0.5,0.3),sd=1))
plot(y1,ylab="$y_t$")
lines(y2,col=1,lty=2)
## ----fpp2cond------------------------------------------------------------
require(fpp2)
data(condmilk)
condmilk
frequency(condmilk)
plot(condmilk)
## ----sdsp----------------------------------------------------------------
tsdisplay(diff(condmilk,differences=1,lag=12))
## ----acld----------------------------------------------------------------
tsdisplay(diff(condmilk,differences=2,lag=12))
## ----sugar---------------------------------------------------------------
fit1 = Arima(condmilk, c(0,1,0), c(0,1,1))
fit2 = Arima(condmilk, c(0,1,0), c(2,1,0))
fit3 = Arima(condmilk, c(5,1,0), c(0,1,1))
fit4 = Arima(condmilk, c(5,1,0), c(2,1,0))
fit5 = Arima(condmilk, c(0,1,5), c(0,1,1))
fit6 = Arima(condmilk, c(0,1,5), c(2,1,0))
## ----fitai---------------------------------------------------------------
fit.aic <- c(AIC(fit1), AIC(fit2), AIC(fit3), AIC(fit4), AIC(fit5), AIC(fit6))
which.min(fit.aic)
## ----ratsd---------------------------------------------------------------
hist(residuals(fit1))
tsdisplay(residuals(fit1))
Box.test(residuals(fit1), lag = 24, type = "Ljung-Box", fitdf = 1)
## ----fcp1----------------------------------------------------------------
fcast1 = forecast(fit1, h = 24)
plot(fcast1)
set.seed(42)
n <- 1000
x <- sort(rnorm(n))
dgp <- 1+x^2
y <- dgp + rnorm(n)
model <- lm(y~x)
summary(model)
## ----quaddgp,echo=FALSE,fig.cap="Data, $g(x)$, and Linear Model Fit (Top), $dg(x)/dx$ and $\\hat\\beta_1$ (Bottom).",fig.asp=1.5----
par(mfrow=c(2,1),cex=1/0.83)
plot(x,y,cex=0.25,xlab="$x$",ylab="$y$")
lines(x,dgp,col="grey")
abline(model,col=1,lty=2)
legend("top",c("Data","$g(x)$","$\\hat\\beta_0+\\hat\\beta_1x$"),
col=c(1,"grey",1),
lty=c(NA,1,2),
pch=c(1,NA,NA),
pt.cex=0.25,
bty="n")
plot(x,2*x,col="grey",type="l",xlab="$x$",ylab="$dg(x)/dx$")
abline(h=coef(model)[2],col=1,lty=2)
legend("topleft",c("$dg(x)/dx$",paste("$\\hat\\beta_1=",format(coef(model)[2],digits=2),"$",sep="")),
col=c(1,1),
lty=1:2,
bty="n")
## ----mcdgp,echo=FALSE,fig.cap="Monte Carlo DGP."-------------------------
n <- 1000
x <- seq(-3.5,3.5,length=n)
dgp <- x^6
dgp <- 1+dgp/sd(dgp)
y <- dgp + rnorm(n,sd=0.5)
plot(x,y,cex=0.25,col="grey",xlab="$x$",ylab="$y$")
lines(x,dgp)
## ----smcmma,echo=FALSE---------------------------------------------------
## A simple Monte Carlo that compares Mallows model averaging (MMA), model
## selection (Cp, Mallows's criterion), and model assertion. We need to
## load the quadprog package to compute the model averaging weights.
require(quadprog)
## Set the seed, number of observations (n), number of candidate models (M)
## and number of Monte Carlo replications (num.reps).
set.seed(42)
n <- 100
M <- 10
num.reps <- 1000
## Create some matrices and vectors for storing results
fitted.mat <- matrix(NA,n,M)
residual.mat <- matrix(NA,n,M)
mse <- matrix(NA,num.reps,M+2)
cp <- numeric(M)
cp.vec <- numeric(M)
## Conduct num.reps Monte Carlo replications
for(m in 1:num.reps) {
## For each replication simulate the DGP
x <- rnorm(n)
dgp <- x^6
dgp <- 1+dgp/sd(dgp)
y <- dgp + rnorm(n,sd=0.5)
## X has M columns which are used to generate the
## candidate models, e.g., column 1 for model 1, columns
## 1 and 2 for model 2 etc.
X <- poly(x,degree=M)
## For each replication fit the M candidate models
for(i in 1:M) {
model <- lm(y~X[,1:i])
fitted.mat[,i] <- fitted(model)
residual.mat[,i] <- residuals(model)
## Number of parameters in model i (intercept + i columns)
k <- i+1
sum.sqres <- sum(residuals(model)^2)
sigsq <- sum.sqres/(n-k)
## Mallows's model selection criterion
cp[i] <- sum.sqres+2*sigsq*k
}
## Construct the model averaging estimator
## The largest model contains all M columns of X
sigsq <- (t(residual.mat[,M]) %*% residual.mat[,M])/(n-M)
## Set up the quadratic programming problem
Amat <- cbind(rep(1,M),diag(1,M,M))
bvec <- c(1,rep(0,M))
Dmat <- t(residual.mat)%*%residual.mat
if(qr(Dmat)$rank<M) Dmat <- Dmat + diag(1e-10,M,M)
K <- 2:(M+1)
dvec <- -2*sigsq*K
## Solve the quadratic program for the model averaging weights
w.hat.mma <- solve.QP(Dmat,dvec,Amat,bvec,1)$solution
## Compute the model averaging estimator
yhat.mma <- fitted.mat%*%w.hat.mma
## Compute MSEs for each approach (averaging, selection, assertion)
mse.mma <- mean((dgp-yhat.mma)^2)
mse.Cp <- mean((dgp-fitted.mat[,which.min(cp)])^2)
mse.assertion <- colMeans(sweep(fitted.mat,1,dgp)^2)
mse[m,] <- c(mse.mma,mse.Cp,mse.assertion)
## Keep track of which candidate model was selected by Cp
cp.vec[m] <- which.min(cp)
}
## ----msefig,echo=FALSE,fig.cap="Model Selection, Averaging, and Assertion with 10 Candidate Models (orthogonal polynomials of order 1--10)."----
mse <- data.frame(mse)
names(mse) <- c("Averaging","Selection",paste("k =",1:M))
boxplot(mse,outline=FALSE,notch=TRUE,ylab="MSE",las=2)
mean.mse <- apply(mse,2,mean)
median.mse <- apply(mse,2,median)
## ----msetab,echo=FALSE---------------------------------------------------
options(scipen=9)
foo <- data.frame(formatC(mean.mse,format="f",digits=4),
formatC(mean.mse/mean.mse[rank(mean.mse)==2],format="f",digits=2),
rank(mean.mse))
knitr::kable(foo,
caption="Mean MSE and Ranking of MSE Performance (k = 6 is the oracle model).",
align=c("r","r","r"),
col.names=c("Mean MSE","Relative Efficiency","Rank"))
## ----mseproptab,echo=FALSE-----------------------------------------------
foo <- prop.table(table(cp.vec))
names(foo) <- paste("k =",names(foo))
knitr::kable(data.frame(names(foo),format(as.numeric(foo),digits=2)),
align=c("l","c"),
caption="Model Selection Proportion Among the Candidate Models (k = 6 is the oracle model).",
col.names=c("Model","Selection Proportion"))
## ----maslrp--------------------------------------------------------------
require(MASS)
set.seed(42)
n <- 1000
x <- rnorm(n)
## Compute the largest desired model using raw polynomials
X <- data.frame(poly(x,10,raw=TRUE))
## Generate the DGP (column is polynomial order in x,
## so dgp <- X[,3] means x^3)
dgp <- X[,3]
y <- dgp/sd(dgp) + rnorm(n,sd=0.5)
## First, estimate the full model and examine its AIC value
model <- lm(y~.,data=X)
formula(model)
## AIC
AIC(model)
## Next run stepwise AIC model selection
model.aic <- stepAIC(model,k=2,trace=0,direction="both",
scope=list(upper=~.,lower=~1))
## Check the model that was selected and its AIC value
formula(model.aic)
AIC(model.aic)
## Repeat using the BIC criterion
BIC(model)
## Stepwise BIC model selection
model.bic <- stepAIC(model,k=log(n),trace=0,direction="both",
scope=list(upper=~.,lower=~1))
formula(model.bic)
BIC(model.bic)
## ----cvdoouems-----------------------------------------------------------
set.seed(42)
n <- 1000
x <- rnorm(n)
## Compute the largest desired model using raw polynomials
X <- poly(x,10,raw=TRUE)
## Generate the DGP (column is polynomial order in x,
## so dgp <- X[,3] means x^3)
dgp <- X[,3]
y <- dgp/sd(dgp) + rnorm(n,sd=0.5)
cv.vec <- numeric()
model <- list()
## For simplicity we consider models with one predictor, x^j,
## j=1,...,10
for(i in 1:ncol(X)) {
model[[i]]<-lm(y~X[,i,drop=FALSE])
cv.vec[i]<-mean(residuals(model[[i]])^2/(1-hatvalues(model[[i]]))^2)
}
which.min(cv.vec)
BIC(model[[which.min(cv.vec)]])
summary(model[[which.min(cv.vec)]])
## ----erbmapex------------------------------------------------------------
## DGP
set.seed(42)
n <- 100
x1 <- sort(runif(n,-1,1))
x2 <- runif(n,-1,1)
dgp <- 1 -2*x1**3 + x2
dgp.x2.med <- 1 -2*x1**3 + median(x2)
dgp.deriv.x1 <- -6*x1**2
y <- dgp + rnorm(n,sd=0.5)
## Models do not include the true model (lm(y~I(x1^3)+x2)) nor a model
## with highest degree that for the true model, they are all
## simple (raw) misspecified polynomials
model.1 <- lm(y~x1+x2)
model.2 <- lm(y~poly(x1,2,raw=TRUE)+poly(x2,2,raw=TRUE))
model.3 <- lm(y~poly(x1,4,raw=TRUE)+poly(x2,4,raw=TRUE))
model.4 <- lm(y~poly(x1,5,raw=TRUE)+poly(x2,5,raw=TRUE))
## Derivatives withrespect to x1
model.1.deriv.x1 <- rep(coef(model.1)[2],n)
model.2.deriv.x1 <- coef(model.2)[2]+2*coef(model.2)[3]*x1
model.3.deriv.x1 <- coef(model.3)[2]+2*coef(model.3)[3]*x1+
3*coef(model.3)[4]*x1**2+
4*coef(model.3)[5]*x1**3
model.4.deriv.x1 <- coef(model.4)[2]+2*coef(model.4)[3]*x1+
3*coef(model.4)[4]*x1**2+
4*coef(model.4)[5]*x1**3+
5*coef(model.4)[6]*x1**4
## Create the residual matrix (residual.mat), rank vector (K), and
## variance of the largest dimensional model
residual.mat <- cbind(residuals(model.1),
residuals(model.2),
residuals(model.3),
residuals(model.4))
K <- c(model.1$rank,model.2$rank,model.3$rank,model.4$rank)
sigsq <- summary(model.4)$sigma^2
## Solve the quadratic program for the Mallows model average weights
require(quadprog)
M <- ncol(residual.mat)
D <- t(residual.mat)%*%residual.mat
if(qr(D)$rank<M) D <- D + diag(1e-10,M,M)
A <- cbind(rep(1,M),diag(1,M,M))
b0 <- c(1,rep(0,M))
d <- -sigsq*K
b <- solve.QP(Dmat=D,dvec=d,Amat=A,bvec=b0,meq=1)$solution
## Compute the MMA derivative estimate and MMA fitted values
## For these models there are no interaction terms, so no need to
## control for the value of x2 for the derivative
ma.deriv.x1 <- cbind(model.1.deriv.x1,
model.2.deriv.x1,
model.3.deriv.x1,
model.4.deriv.x1)%*%b
## Surface depends on value of x2, set to median to see partial plot
ma.fit<-cbind(predict(model.1,newdata=data.frame(x1,x2=median(x2))),
predict(model.2,newdata=data.frame(x1,x2=median(x2))),
predict(model.3,newdata=data.frame(x1,x2=median(x2))),
predict(model.4,newdata=data.frame(x1,x2=median(x2))))%*%b
## ----erbmapexfig,echo=FALSE,fig.cap="Model Averaging Illustration (true DGP is $y_i=1-2x_{i1}^3+x_{i2}+\\epsilon_i$, correctly specified model is not among the set of candidate models).",fig.asp=1.5----
par(mfrow=c(2,1),cex=1/0.83)
plot(x1,ma.fit,xlab="$x_1$",ylab="$y$",type="l")
lines(x1,dgp.x2.med,col=1,lty=2)
legend("topleft",c("DGP","MMA estimate $\\hat y_m$"),col=c(1,1),lty=1:2,bty="n")
plot(x1,dgp.deriv.x1,xlab="$x_1$",ylab="$d\\hat y_m/d\\,x_1$",type="l")
lines(x1,ma.deriv.x1,col=1,lty=2)
legend("topleft",c("DGP marginal effect","MMA marginal effect"),col=c(1,1),lty=1:2,bty="n")
## ----exrraeumma,eval=FALSE-----------------------------------------------
## ## Linear example
## dgp <- 1 -2*x1 + x2
## dgp.x2.med <- 1 -2*x1 + median(x2)
## dgp.deriv.x1 <- rep(-2,n)
## ## Order 5 polynomial
## dgp <- 1 -2*x1**5 + x2
## dgp.x2.med <- 1 -2*x1**5 + median(x2)
## dgp.deriv.x1 <- -10*x1**4
## ----fBactqbcCc,fig.cap="A Quadratic Bézier Curve.",echo=FALSE-----------
B <- function(x,b0,b1,b2) { (1-x)^2*b0+(1-x)*x*b1+x^2*b2 }
x <- seq(0,1,length=1000)
b0 <- 1
b1 <- -1
b2 <- 2
plot(x,B(x,b0,b1,b2),ylab="B(x)",type="l",lwd=2,cex.lab=1.25)
## ----pbubcfTBc,fig.cap="The Quadratic Bézier Curve Bases.", echo=FALSE----
Bernstein <- function(n,i,x) { factorial(n)/(factorial(n-i)*factorial(i))*(1-x)^{n-i}*x^i }
x <- seq(0,1,length=100)
degree <- 2
plot(x,Bernstein(degree,0,x),type="l",lwd=2,ylab="B(x)",col=1)
for(i in 1:degree) lines(x,Bernstein(degree,i,x),lty=i+1,lwd=2,col=1)
## ----tpmmfssfrb,eval=TRUE,echo=TRUE--------------------------------------
## Take 4 predictors and their polynomial bases
x1 <- poly(rnorm(100),2)
x2 <- poly(rnorm(100),3)
x3 <- poly(rnorm(100),4)
x4 <- poly(rnorm(100),5)
## Additive dimension
dim(cbind(x1,x2,x3,x4))
## Generalized Taylor dimension
dim(ma:::taylor.model.matrix(list(x1,x2,x3,x4)))
## # Model Uncertainty Answers {-}
## # Problem Set {-}
## ----ansfcckammups,echo=TRUE---------------------------------------------
require(ma)
require(robustbase)
set.seed(42)
M <- 1000
N <- seq(500,1500,by=100)
nobs <- numeric()
mse <- numeric()
i <- 1
for(m in 1:M) {
for(n in N) {
x1 <- runif(n,-1,1)
x2 <- runif(n,-1,1)
dgp <- x1 + x2**2 + x1*x2
y <- dgp + rnorm(n,sd=0.5*sd(dgp))
mse[i] <- mean((fitted(lm.ma(y~x1+x2,verbose=FALSE))-
dgp)^2)
nobs[i] <- n
i <- i + 1
}
require(ma)
require(robustbase)
robustbase
install.packages("ma")
## ----ansfcckammups,echo=TRUE---------------------------------------------
require(ma)
require(robustbase)
set.seed(42)
M <- 1000
N <- seq(500,1500,by=100)
nobs <- numeric()
mse <- numeric()
i <- 1
for(m in 1:M) {
for(n in N) {
x1 <- runif(n,-1,1)
x2 <- runif(n,-1,1)
dgp <- x1 + x2**2 + x1*x2
y <- dgp + rnorm(n,sd=0.5*sd(dgp))
mse[i] <- mean((fitted(lm.ma(y~x1+x2,verbose=FALSE))-
dgp)^2)
nobs[i] <- n
i <- i + 1
}
install.packages("ansfcckammups")
## ----ansfcckammups,echo=TRUE---------------------------------------------
require(ma)
require(robustbase)
set.seed(42)
M <- 1000
N <- seq(500,1500,by=100)
nobs <- numeric()
mse <- numeric()
i <- 1
for(m in 1:M) {
for(n in N) {
x1 <- runif(n,-1,1)
x2 <- runif(n,-1,1)
dgp <- x1 + x2**2 + x1*x2
y <- dgp + rnorm(n,sd=0.5*sd(dgp))
mse[i] <- mean((fitted(lm.ma(y~x1+x2,verbose=FALSE))-
dgp)^2)
nobs[i] <- n
i <- i + 1
}
## # Problem Set {-}
library(lm.ma)
for(m in 1:M) {
for(n in N) {
x1 <- runif(n,-1,1)
x2 <- runif(n,-1,1)
dgp <- x1 + x2**2 + x1*x2
y <- dgp + rnorm(n,sd=0.5*sd(dgp))
mse[i] <- mean((fitted(lm.ma(y~x1+x2,verbose=FALSE))-
dgp)^2)
nobs[i] <- n
i <- i + 1
}
## ----mmaximevsl,eval=FALSE-----------------------------------------------
## ## R code to compute the Mallows Model Averaging (MMA) weight vector.
## ## Presumed inputs are residual.mat (E), rank vector (K), and variance
## ## of the residuals from the largest model (sigsq).
## require(quadprog)
## M <- ncol(residual.mat)
## ## The D matrix in b'D b (M x M).
## D <- t(residual.mat)%*%residual.mat
## ## Ensure that D is well-conditioned.
## if(qr(D)$rank<M) D <- D + diag(1e-10,M,M)
## ## Create the constraint matrix. The first column of A has the adding
## ## up constraint (weights sum to 1), the following M columns the
## ## non-negativity constraints (each weight >= 0).
## A <- cbind(rep(1,M),diag(1,M,M))
## ## The constraint vector b0.
## b0 <- c(1,rep(0,M))
## ## The d vector in d'b (1 x M).
## d <- -sigsq*K
## ## Solve the quadratic program for the MMA weight vector telling
## ## the solver that the first constraint is an equality constraint
## ## (the first meq constraints are treated as equality constraints,
## ## all further as inequality constraints).
## b <- solve.QP(Dmat=D,dvec=d,Amat=A,bvec=b0,meq=1)$solution
## ----qpexam--------------------------------------------------------------
require(quadprog)
## Assume we want to minimize: -(0 5 0) %*% b + 1/2 b^T b
## under the constraints:      A^T b >= b0
## with b0 = (-8,2,0)^T
## and      (-4  2  0)
##      A = (-3  1 -2)
##          ( 0  0  1)
## we can use solve.QP as follows:
##
Dmat       <- matrix(0,3,3)
diag(Dmat) <- 1
dvec       <- c(0,5,0)
Amat       <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
bvec       <- c(-8,2,0)
solve.QP(Dmat,dvec,Amat,bvec)$solution
## ----fBactqbcC,fig.cap="A Quadratic Bézier Curve.",echo=FALSE------------
B <- function(x,b0,b1,b2) { (1-x)^2*b0+(1-x)*x*b1+x^2*b2 }
x <- seq(0,1,length=1000)
b0 <- 1
b1 <- -1
b2 <- 2
plot(x,B(x,b0,b1,b2),ylab="B(x)",type="l",lwd=2,cex.lab=1.25)
## ----pbubcfTB,fig.cap="Quadratic Bézier curve basis functions.",echo=FALSE----
Bernstein <- function(n,i,x) { factorial(n)/(factorial(n-i)*factorial(i))*(1-x)^{n-i}*x^i }
x <- seq(0,1,length=100)
degree <- 2
plot(x,Bernstein(degree,0,x),type="l",lwd=2,ylab="B(x)",col=1)
for(i in 1:degree) lines(x,Bernstein(degree,i,x),lty=i+1,lwd=2,col=1)
## ----cTFrmsfptB,fig.asp=1.25,fig.cap="Third degree B-spline with three interior knots along with its first derivative.",echo=FALSE----
require(crs)
par(mfrow=c(2,1),cex=1/0.83)
degree <- 3
m <- degree+1
## nbreak is the total number of knots (2 indicates 2 endpoints, 0 interior)
nbreak <- 2+3
## N is the number of interior knots
N <- nbreak-2
x <- seq(0,1,length=1000)
B <- gsl.bs(x,degree=degree,nbreak=nbreak,intercept=TRUE)
matplot(x,B,type="l",lwd=2)
deriv <- 1
B.deriv <- gsl.bs(x,degree=degree,nbreak=nbreak,deriv=deriv,intercept=TRUE)
matplot(x,B.deriv,type="l",lwd=2)
