CHAPTER 9 OUTPUT

This shows the outputs from Chapter 9 using R. The textbook is “A Modern Approach to Regression with R” by Simon J. Sheather (2008). The R code provided with the book has been updated.

#March 17, 2009

#Please change the file path in the command below to coincide with where you have stored the data files
setwd("C:/gokul/ModernApproach/RCodes")

library(readxl)
confood2 <- read_excel("MARData.xlsx", sheet  =  "confood2")

attach(confood2)

#Figure 9.1 on page 306
par(mfrow = c(1,1))
plot(log(Price),log(Sales),xlab = expression(log(Price[t])),ylab = expression(log(Sales[t])),
pch = Promotion+2,col = Promotion+1)
legend(-0.3, 8.5,legend = c("No","Yes"),pch = 2:3,col = 1:2,title = "Promotion")

#Figure 9.2 on page 307
plot(Week,log(Sales),type = 'o',ylab = expression(log(Sales[t])),xlab = "Week, t",
pch = Promotion+2,col = Promotion+1)
legend(0, 8.5,legend = c("No","Yes"),pch = 2:3,col = 1:2,title = "Promotion")

#Figure 9.3 on page 307
plot(log(SalesLag1),log(Sales),ylab = expression(log(Sales[t])),
xlab = expression(log(Sales[t-1])))

#Figure 9.4 on page 308
acf(log(Sales))

#Figure 9.5 on page 309
lsm1 <- lm(log(Sales)~log(Price)+Promotion+Week,data = confood2)
StanRes1 <- rstandard(lsm1)
par(mfrow = c(2,2))
plot(log(Price),StanRes1,ylab = "Standardized Residuals",xlab = expression(log(Price[t])))
plot(Week,StanRes1,ylab = "Standardized Residuals",xlab = "Week, t",type = 'o')
plot(Promotion,StanRes1,ylab = "Standardized Residuals")
plot(lsm1$fitted.values,StanRes1,ylab = "Standardized Residuals",xlab = "Fitted Values")

#Figure 9.6 on page 310
par(mfrow = c(1,1))
acf(StanRes1,main = "Series Standardized Residuals")

#R output on page 313
#install.packages("nlme")
library(nlme)

m1 <- gls(log(Sales)~log(Price)+Promotion+Week,correlation = corAR1(form = ~Week),data = confood2,method = "ML")
summary(m1)
## Generalized least squares fit by maximum likelihood
##   Model: log(Sales) ~ log(Price) + Promotion + Week 
##   Data: confood2 
##        AIC     BIC   logLik
##   6.537739 18.2452 2.731131
## 
## Correlation Structure: AR(1)
##  Formula: ~Week 
##  Parameter estimate(s):
##       Phi 
## 0.5503593 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept)  4.675667 0.2383703 19.615142   0.000
## log(Price)  -4.327391 0.5625564 -7.692368   0.000
## Promotion    0.584650 0.1671113  3.498565   0.001
## Week         0.012517 0.0046692  2.680813   0.010
## 
##  Correlation: 
##            (Intr) lg(Pr) Promtn
## log(Price)  0.807              
## Promotion   0.559  0.682       
## Week       -0.625 -0.157 -0.206
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.9473082 -0.6095076  0.1031472  0.5769989  2.9558179 
## 
## Residual standard error: 0.2740294 
## Degrees of freedom: 52 total; 48 residual
intervals(m1)
## Approximate 95% confidence intervals
## 
##  Coefficients:
##                    lower        est.       upper
## (Intercept)  4.196391300  4.67566686  5.15494243
## log(Price)  -5.458486702 -4.32739122 -3.19629575
## Promotion    0.248649971  0.58464986  0.92064974
## Week         0.003129195  0.01251724  0.02190529
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.2867307 0.5503593 0.7365028
## 
##  Residual standard error:
##     lower      est.     upper 
## 0.2113301 0.2740294 0.3553309
#Figure 9.7 on page 314 
acf(m1$residuals,main = "Series GLS Residuals")

#R output on page 318
g <- lm(log(Sales)~log(Price)+Promotion+Week,data = confood2)
rho <- 0.5504
x <- model.matrix(g)
Sigma <- with(confood2, diag(length(Week)))
Sigma <- rho^abs(row(Sigma)-col(Sigma))
sm <- chol(Sigma)
smi <- solve(t(sm))
xstar <- smi %*% x
ystar <- with(confood2, smi %*% log(Sales))
m1tls <- lm(ystar ~ xstar-1) 
summary(m1tls)
## 
## Call:
## lm(formula = ystar ~ xstar - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99263 -0.12518  0.01946  0.13055  0.67210 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## xstar(Intercept)  4.67566    0.23838  19.614  < 2e-16 ***
## xstarlog(Price)  -4.32741    0.56256  -7.692 6.44e-10 ***
## xstarPromotion    0.58464    0.16711   3.499  0.00102 ** 
## xstarWeek         0.01252    0.00467   2.681  0.01004 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2852 on 48 degrees of freedom
## Multiple R-squared:  0.9943, Adjusted R-squared:  0.9939 
## F-statistic:  2106 on 4 and 48 DF,  p-value: < 2.2e-16
#Figure 9.8 on page 318
par(mfrow = c(2,2))
plot(xstar[,1],ystar,xlab = "Intercept*",ylab = "log(Sales)*")
identify(xstar[,1],ystar,Week)
## integer(0)
plot(xstar[,2],ystar,xlab = "log(Price)*",ylab = "log(Sales)*")
identify(xstar[,2],ystar,confood2$Week)
## integer(0)
plot(xstar[,3],ystar,xlab = "Promotion*",ylab = "log(Sales)*")
identify(xstar[,3],ystar,confood2$Week)
## integer(0)
plot(xstar[,4],ystar,xlab = "Week*",ylab = "log(Sales)*")
identify(xstar[,4],ystar,Week)

## integer(0)
#Figure 9.9 on page 319
StanRes1 <- rstandard(m1tls)
par(mfrow = c(1,1))
acf(StanRes1,main = "Series Standardized LS Residuals")

#Figure 9.10 on page 320
par(mfrow = c(2,2))
plot(xstar[,2],StanRes1,ylab = "Standardized LS Residuals",xlab = "log(Price)*")
plot(xstar[,4],StanRes1,ylab = "Standardized LS Residuals",xlab = "Week*",type = 'o')
identify(xstar[,4],StanRes1,labels = Week,cex = 0.75)
## integer(0)
plot(xstar[,3],StanRes1,ylab = "Standardized LS Residuals",xlab = "Promotion*")
plot(m1tls$fitted.values,StanRes1,ylab = "Standardized LS Residuals",xlab = "Fitted Values*")

#Figure 9.11 on page 320
par(mfrow = c(2,2))
plot(m1tls)
abline(v = 2*4/length(Week),lty = 2)

detach(confood2)
BayArea <- read_excel("MARData.xlsx", sheet  =  "BayArea")
attach(BayArea)

#Figure 9.12 on page 321
pairs(InterestRate~LoansClosed+VacancyIndex)

#Figure 9.13 on page 322
m1 <- lm(InterestRate~LoansClosed+VacancyIndex)
summary(m1)
## 
## Call:
## lm(formula = InterestRate ~ LoansClosed + VacancyIndex)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22079 -0.08295 -0.01989  0.13077  0.19638 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.580687   0.312308   24.27 4.74e-14 ***
## LoansClosed  -0.007756   0.001249   -6.21 1.25e-05 ***
## VacancyIndex -0.176577   0.132733   -1.33    0.202    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1432 on 16 degrees of freedom
## Multiple R-squared:  0.8344, Adjusted R-squared:  0.8137 
## F-statistic: 40.31 on 2 and 16 DF,  p-value: 5.651e-07
StanRes1 <- rstandard(m1)
mres1 <- lm(StanRes1~LoansClosed+I(LoansClosed^2))
a1 <- mres1$coeff[1]
a2 <- mres1$coeff[2]
a3 <- mres1$coeff[3]
mres2 <- lm(StanRes1~m1$fitted.values+I(m1$fitted.values^2))
b1 <- mres2$coeff[1]
b2 <- mres2$coeff[2]
b3 <- mres2$coeff[3]
par(mfrow = c(2,2))
plot(LoansClosed,StanRes1,ylab = "Standardized Residuals")
curve(a1 + a2*x + a3*x^2, add  =  TRUE, col  =  "blue",lty = 2)
plot(VacancyIndex,StanRes1,ylab = "Standardized Residuals")
plot(m1$fitted.values,StanRes1,ylab = "Standardized Residuals",xlab = "Fitted Values")
curve(b1 + b2*x + b3*x^2, add  =  TRUE, col  =  "blue",lty = 2)
acf(m1$residuals,main = "Standardized LS Residuals")

#R output on page 323
library(nlme)
m1 <- gls(InterestRate~LoansClosed+VacancyIndex,correlation = corAR1(form = ~Month),data = BayArea,method = "ML")
acf(m1$residuals,main = "GLS Residuals")
summary(m1)
## Generalized least squares fit by maximum likelihood
##   Model: InterestRate ~ LoansClosed + VacancyIndex 
##   Data: BayArea 
##         AIC       BIC   logLik
##   -35.30833 -30.58613 22.65416
## 
## Correlation Structure: AR(1)
##  Formula: ~Month 
##  Parameter estimate(s):
##       Phi 
## 0.9572093 
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept)   7.122990 0.4182065 17.032232  0.0000
## LoansClosed  -0.003432 0.0011940 -2.874452  0.0110
## VacancyIndex -0.076340 0.1307842 -0.583710  0.5676
## 
##  Correlation: 
##              (Intr) LnsCls
## LoansClosed  -0.316       
## VacancyIndex -0.822  0.117
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.3439744 -1.1533069 -0.8047658  0.4467079  1.1783687 
## 
## Residual standard error: 0.2377426 
## Degrees of freedom: 19 total; 16 residual
intervals(m1)
## Approximate 95% confidence intervals
## 
##  Coefficients:
##                     lower         est.         upper
## (Intercept)   6.236431638  7.122989795  8.0095479516
## LoansClosed  -0.005963412 -0.003432182 -0.0009009516
## VacancyIndex -0.353590009 -0.076339971  0.2009100658
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.5284235 0.9572093 0.9969063
## 
##  Residual standard error:
##      lower       est.      upper 
## 0.06868792 0.23774259 0.82287451
#R output on page 323
g <- lm(InterestRate~LoansClosed+VacancyIndex,data = BayArea)
rho <- 0.9572093
x <- model.matrix(g)
Sigma <- diag(length(InterestRate))
Sigma <- rho^abs(row(Sigma)-col(Sigma))
sm <- chol(Sigma)
smi <- solve(t(sm))
xstar <- smi %*% x
ystar <- smi %*% InterestRate
m1tls <- lm(ystar ~ xstar-1) 
summary(m1tls)
## 
## Call:
## lm(formula = ystar ~ xstar - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27439 -0.12946  0.00934  0.25772  0.52132 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## xstar(Intercept)   7.122990   0.418207  17.032 1.12e-11 ***
## xstarLoansClosed  -0.003432   0.001194  -2.874    0.011 *  
## xstarVacancyIndex -0.076340   0.130784  -0.584    0.568    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2591 on 16 degrees of freedom
## Multiple R-squared:  0.9831, Adjusted R-squared:  0.9799 
## F-statistic: 309.8 on 3 and 16 DF,  p-value: 2.227e-14
#Figure 9.14 on page 324
par(mfrow = c(2,2))

plot(xstar[,1],ystar,xlab = "Intercept*",ylab = "InterestRate*")
identify(xstar[,1],ystar,Month)
## integer(0)
plot(xstar[,2],ystar,xlab = "LoansClosed*",ylab = "InterestRate*")
identify(xstar[,2],ystar,Month)
## integer(0)
plot(xstar[,3],ystar,xlab = "VacancyIndex*",ylab = "InterestRate*")
identify(xstar[,3],ystar,Month)
## integer(0)
plot(xstar[,2],xstar[,3],xlab = "LoansClosed*",ylab = "VacancyIndex*")
identify(xstar[,2],xstar[,3],Month)

## integer(0)
#Figure 9.15 on page 325
StanRes1 <- rstandard(m1tls)
par(mfrow = c(2,2))
acf(StanRes1,main = "Standardized LSResiduals")
plot(xstar[,2],StanRes1,ylab = "Standardized LS Residuals",xlab = "LoansClosed*")
identify(xstar[,2],StanRes1,Month)
## integer(0)
plot(xstar[,3],StanRes1,ylab = "Standardized LS Residuals",xlab = "VacancyIndex*")
identify(xstar[,3],StanRes1,Month)
## integer(0)
plot(m1tls$fitted.values,StanRes1,ylab = "Standardized Residuals",xlab = "Fitted Values*")
identify(m1tls$fitted.values,StanRes1,Month)

## integer(0)
detach(BayArea)
#################EXERCISES

#Exercise 9.4.1
boxoffice <- read_excel("MARData.xlsx", sheet  =  "boxoffice")
attach(boxoffice)

#Figure 9.16 on page 326
YearsS1975 <- year - 1975
lsm1 <- lm(GrossBoxOffice~YearsS1975,data = boxoffice)
StanRes1 <- rstandard(lsm1)
par(mfrow = c(2,2))
plot(YearsS1975,GrossBoxOffice,ylab = "Gross Box Office ($M)",xlab = "Years since 1975")
abline(lsm1,lty = 2)
plot(YearsS1975,StanRes1,ylab = "Standardized Residuals",xlab = "Years since 1975")
acf(StanRes1,main = "Series Standardized Residuals")

#R output on page 327
library(nlme)
m1 <- gls(GrossBoxOffice~YearsS1975,correlation = corAR1(form = ~YearsS1975),data = boxoffice,method = "ML")
summary(m1)
## Generalized least squares fit by maximum likelihood
##   Model: GrossBoxOffice ~ YearsS1975 
##   Data: boxoffice 
##        AIC      BIC    logLik
##   330.3893 336.2522 -161.1947
## 
## Correlation Structure: AR(1)
##  Formula: ~YearsS1975 
##  Parameter estimate(s):
##       Phi 
## 0.8782065 
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept)  4.514082  72.74393 0.062054  0.9509
## YearsS1975  27.075395   3.44766 7.853259  0.0000
## 
##  Correlation: 
##            (Intr)
## YearsS1975 -0.782
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.93420837 -1.38592020  0.01822281  0.33202594  1.54269773 
## 
## Residual standard error: 76.16492 
## Degrees of freedom: 32 total; 30 residual
#R output on page 327
g <- lm(GrossBoxOffice~YearsS1975,data = boxoffice)
rho <- 0.8782065
x <- model.matrix(g)
Sigma <- diag(length(YearsS1975))
Sigma <- rho^abs(row(Sigma)-col(Sigma))
sm <- chol(Sigma)
smi <- solve(t(sm))
xstar <- smi %*% x
ystar <- smi %*% GrossBoxOffice
m1tls <- lm(ystar ~ xstar-1) 
summary(m1tls)
## 
## Call:
## lm(formula = ystar ~ xstar - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -214.235  -42.370    0.902   33.011  202.415 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## xstar(Intercept)    4.514     72.744   0.062    0.951    
## xstarYearsS1975    27.075      3.448   7.853 9.17e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 78.66 on 30 degrees of freedom
## Multiple R-squared:  0.8427, Adjusted R-squared:  0.8322 
## F-statistic: 80.37 on 2 and 30 DF,  p-value: 8.919e-13
#Figure 9.17 on page 328
StanRes1 <- rstandard(m1tls)
mres2 <- lm(StanRes1~m1tls$fitted.values+I(m1tls$fitted.values^2)+I(m1tls$fitted.values^3))
b1 <- mres2$coeff[1]
b2 <- mres2$coeff[2]
b3 <- mres2$coeff[3]
b4 <- mres2$coeff[4]
mres3 <- lm(StanRes1~m1tls$fitted.values+I(m1tls$fitted.values^2)+I(m1tls$fitted.values^3)+I(m1tls$fitted.values^4)+I(m1tls$fitted.values^5))
par(mfrow = c(1,2))

plot(m1tls$fitted.values,StanRes1,ylab = "Standardized LS Residuals",xlab = "Fitted Values*")
curve(b1 + b2*x + b3*x^2 + + b4*x^3, add  =  TRUE,lty = 2)
acf(StanRes1,main = "Stand LS Residuals")

detach(boxoffice)