R : Copyright 2001, The R Development Core Team Version 1.3.1 (2001-08-31) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type `license()' or `licence()' for distribution details. R is a collaborative project with many contributors. Type `contributors()' for more information. Type `demo()' for some demos, `help()' for on-line help, or `help.start()' for a HTML browser interface to help. Type `q()' to quit R. > # 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) Call: lm(formula = dist ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -29.069 -9.525 -2.272 9.215 43.201 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.5791 6.7584 -2.601 0.0123 * speed 3.9324 0.4155 9.464 1.49e-12 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 15.38 on 48 degrees of freedom Multiple R-Squared: 0.6511, Adjusted R-squared: 0.6438 F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12 > 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) Analysis of Variance Table Model 1: dist ~ speed Model 2: dist ~ speed + I(speed^2) Res.Df RSS Df Sum of Sq F Pr(>F) 1 48 11353.5 2 47 10824.7 1 528.8 2.296 0.1364 > library(lmtest) > bptest(carfit) Breusch-Pagan-Test data: form1 BP = 4.3257, df = 2, p-value = 0.115 > dwtest(carfit) Durbin-Watson-Test data: DW = 1.6762, p-value = NA > carpred <- predict(carfit,se.fit=TRUE,interval="prediction") > carpred$fit[1,] fit lwr upr -1.84946 -34.49984 30.80092 > x <- rnorm(100) > y <- 1 + x + rnorm(100) > y <- as.integer(y<0) > probitfit <- glm(y ~ x, family = binomial(link = probit)) > summary(probitfit) Call: glm(formula = y ~ x, family = binomial(link = probit)) Deviance Residuals: Min 1Q Median 3Q Max -1.7926 -0.5329 -0.2885 -0.1010 2.4735 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.2010 0.2025 -5.931 3.01e-09 *** x -0.8413 0.2038 -4.128 3.66e-05 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 91.177 on 99 degrees of freedom Residual deviance: 67.376 on 98 degrees of freedom AIC: 71.376 Number of Fisher Scoring iterations: 5 > 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) Formula: dist ~ a + speed^b Parameters: Estimate Std. Error t value Pr(>|t|) a -3.35701 4.44957 -0.754 0.454 b 1.39114 0.02971 46.827 <2e-16 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 15.11 on 48 degrees of freedom Correlation of Parameter Estimates: a b -0.8771 > 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) Call: arima0(x = UKDriverDeaths, order = c(1, 0, 1), seasonal = list(order = c(0, 1, 1), period = 12), xreg = seatbelt) Coefficients: ar1 ma1 sma1 xreg1 0.9360 -0.6042 -0.8898 -346.9095 Approx standard errors: ar1 ma1 sma1 xreg1 0.0550 0.1328 0.0117 86.0320 sigma^2 estimated as 16767: log likelihood = -1139.96, aic = 2287.93 > 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)) ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** I INITIAL X(I) D(I) 1 0.569929E-04 0.100E+01 2 0.500000E-01 0.100E+01 3 0.500000E-01 0.100E+01 IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF 0 1 -0.806E+04 1 8 -0.806E+04 0.99E-05 0.26E-04 0.1E-04 0.2E+12 0.1E-05 0.27E+07 2 9 -0.806E+04 0.63E-06 0.64E-06 0.1E-04 0.2E+01 0.1E-05 0.10E+01 3 18 -0.807E+04 0.49E-03 0.84E-03 0.3E+00 0.2E+01 0.5E-01 0.10E+01 4 21 -0.808E+04 0.12E-02 0.10E-02 0.6E+00 0.2E+01 0.2E+00 0.11E+00 5 22 -0.809E+04 0.13E-02 0.25E-02 0.4E+00 0.2E+01 0.4E+00 0.92E+01 6 34 -0.809E+04 0.10E-02 0.23E-02 0.2E-05 0.4E+01 0.2E-05 0.42E-02 7 35 -0.809E+04 0.22E-04 0.18E-04 0.2E-05 0.2E+01 0.2E-05 0.89E-03 8 36 -0.809E+04 0.18E-05 0.19E-05 0.2E-05 0.2E+01 0.2E-05 0.12E-02 9 45 -0.810E+04 0.25E-03 0.27E-03 0.3E-01 0.7E+00 0.4E-01 0.12E-02 10 47 -0.810E+04 0.47E-03 0.47E-03 0.2E-01 0.1E+01 0.4E-01 0.44E-02 11 49 -0.811E+04 0.10E-02 0.11E-02 0.4E-01 0.1E+01 0.8E-01 0.16E-01 12 51 -0.811E+04 0.19E-03 0.28E-03 0.8E-02 0.2E+01 0.2E-01 0.42E-02 13 54 -0.812E+04 0.80E-03 0.83E-03 0.2E-01 0.8E+00 0.5E-01 0.44E-02 14 55 -0.812E+04 0.63E-03 0.93E-03 0.2E-01 0.1E+01 0.5E-01 0.30E-02 15 57 -0.813E+04 0.47E-03 0.68E-03 0.8E-02 0.1E+01 0.2E-01 0.91E-03 16 58 -0.813E+04 0.35E-04 0.44E-04 0.7E-02 0.6E+00 0.2E-01 0.55E-04 17 59 -0.813E+04 0.18E-06 0.19E-04 0.4E-02 0.0E+00 0.9E-02 0.19E-04 18 61 -0.813E+04 0.64E-05 0.17E-04 0.5E-03 0.1E+01 0.1E-02 0.25E-04 19 62 -0.813E+04 0.59E-05 0.70E-05 0.5E-03 0.2E+01 0.1E-02 0.15E-04 20 64 -0.813E+04 0.17E-05 0.21E-05 0.1E-02 0.4E+00 0.2E-02 0.23E-05 21 65 -0.813E+04 0.64E-07 0.62E-07 0.4E-04 0.0E+00 0.9E-04 0.62E-07 22 66 -0.813E+04 0.12E-08 0.97E-09 0.3E-04 0.0E+00 0.7E-04 0.97E-09 23 67 -0.813E+04 0.87E-11 0.11E-11 0.1E-05 0.0E+00 0.3E-05 0.11E-11 24 68 -0.813E+04 0.74E-13 0.89E-16 0.8E-08 0.0E+00 0.2E-07 0.89E-16 25 69 -0.813E+04 0.39E-14 0.42E-18 0.4E-09 0.0E+00 0.8E-09 0.42E-18 26 70 -0.813E+04 0.00E+00 0.43E-23 0.1E-11 0.0E+00 0.3E-11 0.43E-23 ***** X- AND RELATIVE FUNCTION CONVERGENCE ***** FUNCTION -0.812580E+04 RELDX 0.125E-11 FUNC. EVALS 70 GRAD. EVALS 26 PRELDF 0.433E-23 NPRELDF 0.433E-23 I FINAL X(I) D(I) G(I) 1 0.872214E-06 0.100E+01 0.195E-02 2 0.453209E-01 0.100E+01 0.114E-06 3 0.941866E+00 0.100E+01 0.113E-06 > summary(ftse.garch) Call: garch(x = ftse, order = c(1, 1)) Model: GARCH(1,1) Residuals: Min 1Q Median 3Q Max -4.57065 -0.57878 0.01228 0.69132 6.67351 Coefficient(s): Estimate Std. Error t value Pr(>|t|) a0 8.722e-07 3.083e-07 2.829 0.00467 ** a1 4.532e-02 6.791e-03 6.674 2.49e-11 *** b1 9.419e-01 1.018e-02 92.478 < 2e-16 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Diagnostic Tests: Jarque Bera Test data: Residuals X-squared = 215.3627, df = 2, p-value = < 2.2e-16 Box-Ljung test data: Squared.Residuals X-squared = 0.1001, df = 1, p-value = 0.7517 > library(MASS) Attaching package `MASS': The following object(s) are masked _by_ .GlobalEnv : vcov > 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) lower upper [1,] 2.1449125 4.0137381 [2,] -0.2896639 -0.1842742 > 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) >