# The following file contains all examples found in # "Using R to Teach Econometrics" by J. Racine and # R. Hyndman. The output file `examples.out' was generated # by running all commands contained in this file using # R version 1.3.1 (11/2/01, J. Racine) data(cars) carfit <- lm(dist ~ speed, data=cars) summary(carfit) plot(speed ~ dist, data=cars) abline(carfit) vcov <- summary(carfit)$sigma^2 * summary(carfit)$cov.unscaled carfitunres <- lm(dist ~ speed + I(speed^2), data=cars) carfitres <- lm(dist ~ speed, data=cars) anova(carfitres,carfitunres) library(lmtest) bptest(carfit) dwtest(carfit) carpred <- predict(carfit,se.fit=TRUE,interval="prediction") carpred$fit[1,] x <- rnorm(100) y <- 1 + x + rnorm(100) y <- as.integer(y<0) probitfit <- glm(y ~ x, family = binomial(link = probit)) summary(probitfit) index <- probitfit$coef[1] + probitfit$coef[2]*x invmills <- dnorm(index)/pnorm(index) library(nls) carfitnls <- nls(dist ~ a + speed^b, start=list(a=1,b=1), data=cars) summary(carfitnls) library(ts) data(UKDriverDeaths) plot(UKDriverDeaths) diff.1 <- diff(UKDriverDeaths,12) plot(diff.1) par(mfrow=c(1,2)) acf(diff.1,lag.max=40) pacf(diff.1,lag.max=40) seatbelt <- as.matrix(c(rep(0,169),rep(1,23))) fit <- arima0(UKDriverDeaths,order=c(1,0,1),seasonal=list(order=c(0,1,1),period=12), xreg=seatbelt) print(fit) fcast <- predict(fit,n.ahead=12,newxreg=as.matrix(rep(1,12))) plot(UKDriverDeaths,xlim=c(1969,1986),ylim=c(900,2600)) lines(fcast$pred,col=2) lines(fcast$pred + 2*fcast$se,col=3) lines(fcast$pred - 2*fcast$se,col=3) library(tseries) data(EuStockMarkets) ftse <- diff(log(EuStockMarkets))[,"FTSE"] ftse.garch <- garch(ftse, order=c(1,1)) summary(ftse.garch) library(MASS) data(Boston) crim <- as.integer(Boston$crim > quantile(Boston$crim,0.8)) med.value <- Boston$medv logitreg <- function(beta,x,y) { eta <- exp(beta[1] + beta[2]*x) p <- eta/(1+eta) return(-sum(ifelse(y,log(p),log(1-p)))) } mle <- optim(c(0,0),logitreg,x=med.value,y=crim,hessian=TRUE) se <- sqrt(diag(solve(mle$hessian))) cbind(lower=mle$par - 1.96*se, upper=mle$par + 1.96*se) a <- seq(-10,15,l=50) b <- seq(-1,0.2,l=50) L <- matrix(NA,50,50) for(i in 1:50) for(j in 1:50) L[i,j] <- -logitreg(c(a[i],b[j]), med.value, crim) par(mfrow=c(1,2)) contour(a,b,L) persp(a,b,L)