############################
# 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 = "http://people.stat.sc.edu/hansont/stat520/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 = "globaltemps = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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
# Page 68
data(beersales)
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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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")