############################
# Author: Joshua M. Tebbs #
# Date: 20 Dec 2009 #
# Update: 25 Jul 2011 #
# Purpose: STAT 520 R code #
# CHAPTER 3 #
############################
# Example 3.4.
# Global temperature data from 1900
# Page 53
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
plot(globaltemps.1900,ylab="Global temperature deviations (since 1900)",xlab="Year",type="o")
# Example 3.4.
# Global temperature data from 1900
# Straight line model fit with output
# Page 54
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
summary(fit)
plot(globaltemps.1900,ylab="Global temperature deviations (since 1900)",xlab="Year",type="o")
abline(fit)
# Example 3.4.
# Global temperature data from 1900
# Residuals from straight line model fit
# Page 55
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
plot(resid(fit),ylab="Residuals",xlab="Year",type="o")
# Example 3.4.
# Global temperature data from 1900
# Plot of first data differences
# Page 57
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
plot(diff(globaltemps.1900),ylab="Global temperature deviation differences",xlab="Year",type="o")
# Example 3.5.
# Gold price data
# Page 59
data(gold)
plot(gold,ylab="Price",xlab="Time",type="o")
# Example 3.5.
# Gold price data
# Quadratic model fit with output
# Page 59-60
data(gold)
t <- time(gold)
t.sq <- t^2
fit <- lm(gold ~ t + t.sq)
summary(fit)
plot(gold,ylab="Price",xlab="Time",type="o")
curve(expr = fit$coefficients[1] +
fit$coefficients[2]*x +
fit$coefficients[3]*x^2, col = "black",
lty = "solid", lwd = 1, add = TRUE)
# Example 3.5.
# Gold price data
# Residuals from quadratic model fit
# Page 61
data(gold)
t <- time(gold)
t.sq <- t^2
fit <- lm(gold ~ t + t.sq)
plot(resid(fit),ylab="Residuals",xlab="Time",type="o")
# Example 3.6.
# Beer sales data
# Monthly plotting symbols added
# Regression output included
# Page 62-63
data(beersales)
b.sales<-window(beersales,start=1980)
plot(b.sales,ylab="Sales",xlab="Year",type='l')
points(y=b.sales,x=time(b.sales),pch=as.vector(season(b.sales)))
month.<-season(b.sales)
fit<-lm(b.sales ~ month.-1)
summary(fit)
# Example 3.6.
# Beer sales data
# Residuals from seasonal means model fit
# Page 64
data(beersales)
b.sales<-window(beersales,start=1980)
fit<-lm(b.sales ~ month.-1)
plot(resid(fit),ylab="Residuals",xlab="Time",type="o")
# Example 3.6.
# Beer sales data
# Cosine trend model fit and output
# Page 66-67
data(beersales)
b.sales<-window(beersales,start=1980)
har. <- harmonic(b.sales,1)
fit <- lm(b.sales~har.)
summary(fit)
plot(ts(fitted(fit),freq=12,start=c(1980,1)),ylab="Sales",type='l',ylim=range(c(fitted(fit),b.sales)))
points(b.sales)
# Example 3.6.
# Beer sales data
# Residuals from cosine trend model fit
data(beersales)
# Page 68
b.sales<-window(beersales,start=1980)
har. <- harmonic(b.sales,1)
fit <- lm(b.sales~har.)
plot(resid(fit),ylab="Residuals",xlab="Time",type="o")
# Example 3.4.
# Global temperature data from 1900
# Standardised residuals from straight line model fit
# Histogram and qq plot
# Figures were constructed separately
# Page 72
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
hist(rstudent(fit),main="Histogram of standardized residuals",xlab="Standardized residuals")
qqnorm(rstudent(fit),main="QQ plot of standardized residuals")
# Example 3.4.
# Global temperature data from 1900
# Shapiro-Wilk test for normality
# Page 73
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
shapiro.test(rstudent(fit))
# Example 3.4.
# Global temperature data from 1900
# Standardised residuals from straight line model fit
# Horizontal line added at 0
# Page 74
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
plot(rstudent(fit),ylab="Residuals",xlab="Year",type="o")
abline(h=0)
# Example 3.4.
# Global temperature data from 1900
# Runs test for independence on standardised residuals
# Page 75
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
runs(rstudent(fit))
# Example 3.4.
# Global temperature data from 1900
# Calculate sample ACF for standardised residuals
# Page 78
globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)
globaltemps.1900 <- window(globaltemps,start=1900)
fit <- lm(globaltemps.1900~time(globaltemps.1900))
acf(rstudent(fit),main="Sample ACF for standardized residuals")
# Simulation exercise
# Simulate white noise processes and plot ACFs
# Page 79
w.n.1 <- rnorm(100,0,1)
w.n.2 <- rnorm(100,0,1)
par(mfrow=c(2,2))
plot(w.n.1,ylab="White noise process.1",xlab="Time",type="o")
acf(w.n.1,main="Sample ACF")
plot(w.n.2,ylab="White noise process.2",xlab="Time",type="o")
acf(w.n.2,main="Sample ACF")