R: FORECASTING AND TIME SERIES # 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")