library(forecast)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: timeDate
## This is forecast 7.3
Read in Data and summarize
Appliance.data <- read.csv("ApplianceShipmentsReordered.csv")
shipments.ts <- ts(Appliance.data$Shipments, start = c(1985,1), end = c(1989, 4), freq = 4)
head(Appliance.data)
## Quarter Shipments
## 1 Q1 - 1985 4009
## 2 Q2 - 1985 4321
## 3 Q3 - 1985 4224
## 4 Q4 - 1985 3944
## 5 Q1 - 1986 4123
## 6 Q2 - 1986 4522
summary(Appliance.data)
## Quarter Shipments
## Q1 - 1985: 1 Min. :3944
## Q1 - 1986: 1 1st Qu.:4240
## Q1 - 1987: 1 Median :4489
## Q1 - 1988: 1 Mean :4425
## Q1 - 1989: 1 3rd Qu.:4588
## Q2 - 1985: 1 Max. :4900
## (Other) :14
tail(Appliance.data)
## Quarter Shipments
## 15 Q3 - 1988 4417
## 16 Q4 - 1988 4258
## 17 Q1 - 1989 4245
## 18 Q2 - 1989 4900
## 19 Q3 - 1989 4585
## 20 Q4 - 1989 4533
dim(Appliance.data)
## [1] 20 2
Plot
plot(shipments.ts, xlab = "Time", ylab = "shipments", ylim = c(3900, 5000), bty = "l")
Plot with trend and plot with zoom
shipments.lm <- tslm(shipments.ts ~ trend + I(trend^2))
par(mfrow = c(2, 1))
plot(shipments.ts, xlab = "Time", ylab = "shipments", ylim = c(3900, 5000), bty = "l")
lines(shipments.lm$fitted, lwd = 2)
shipments.ts.zoom <- window(shipments.ts, start = c(1986, 1), end = c(1988, 4))
plot(shipments.ts.zoom, xlab = "Time", ylab = "shipments", ylim = c(3900, 5000), bty = "l")
Forecasts in the validation period from a quadratic trend model estimated from the training period
#quadratic model training and validation lengths
nValid <- 4
shipments.ts
## Qtr1 Qtr2 Qtr3 Qtr4
## 1985 4009 4321 4224 3944
## 1986 4123 4522 4657 4030
## 1987 4493 4806 4551 4485
## 1988 4595 4799 4417 4258
## 1989 4245 4900 4585 4533
nTrain <- length(shipments.ts) - nValid
#quad model training and valid in windows
train.ts <- window(shipments.ts, start = c(1985, 1), end = c(1985, nTrain))
train.ts
## Qtr1 Qtr2 Qtr3 Qtr4
## 1985 4009 4321 4224 3944
## 1986 4123 4522 4657 4030
## 1987 4493 4806 4551 4485
## 1988 4595 4799 4417 4258
valid.ts <- window(shipments.ts, start = c(1985, nTrain + 1), end = c(1985, nTrain + nValid))
valid.ts
## Qtr1 Qtr2 Qtr3 Qtr4
## 1989 4245 4900 4585 4533
#quadratic trend
shipments.lm <- tslm(train.ts ~ trend + I(trend^2))
#forecast based on trend
shipments.lm.pred <- forecast(shipments.lm, h = nValid, level = 0)
#quadratic trend forecast
shipments.lm.pred
## Point Forecast Lo 0 Hi 0
## 1989 Q1 4409.468 4409.468 4409.468
## 1989 Q2 4355.424 4355.424 4355.424
## 1989 Q3 4291.984 4291.984 4291.984
## 1989 Q4 4219.148 4219.148 4219.148
#quadratic error
valid.ts - shipments.lm.pred$mean
## Qtr1 Qtr2 Qtr3 Qtr4
## 1989 -164.4679 544.5758 293.0159 313.8523
plot(shipments.lm.pred, ylim = c(3900, 5000), ylab = "shipments", xlab = "Time", bty = "l", xaxt = "n", xlim = c(1985,1990.25), main ="", flty = 3)
axis(1, at = seq(1985, 1990, 1), labels = format(seq(1985, 1990, 1)))
lines(shipments.lm$fitted, lwd = 1)
lines(valid.ts, lwd = 3)
Accuracy measures for validation period forecast
shipments.lm.pred
## Point Forecast Lo 0 Hi 0
## 1989 Q1 4409.468 4409.468 4409.468
## 1989 Q2 4355.424 4355.424 4355.424
## 1989 Q3 4291.984 4291.984 4291.984
## 1989 Q4 4219.148 4219.148 4219.148
accuracy(shipments.lm.pred$mean, valid.ts)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 246.744 356.3615 328.978 5.138469 7.075664 -0.3992743 0.9384219
Names
names(shipments.lm.pred)
## [1] "model" "mean" "lower" "upper" "level"
## [6] "x" "method" "newdata" "residuals" "fitted"
shipments.lm.pred$residuals
## Qtr1 Qtr2 Qtr3 Qtr4
## 1985 12.73775 228.43971 44.53803 -312.96730
## 1986 -202.07626 138.21113 223.89489 -443.02500
## 1987 -10.54853 281.32430 14.59349 -53.74097
## 1988 63.32094 283.77920 -72.36618 -196.11520
Quadratic Error
valid.ts - shipments.lm.pred$mean
## Qtr1 Qtr2 Qtr3 Qtr4
## 1989 -164.4679 544.5758 293.0159 313.8523
Histogram of Errors
hist(shipments.lm.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="")
#Create a time series out of the data.
shipments.ets.AAN <- ets(shipments.ts, model = "AAN")
AAN <- ets(valid.ts, model = "AAN")
ANNfc <- forecast(AAN, h = nValid, level = 0)
#Fit Model 1 to the time series.
shipments.ets.MMN <- ets(shipments.ts, model = "MMN", damped = FALSE)
shipments.ets.MMN
## ETS(M,M,N)
##
## Call:
## ets(y = shipments.ts, model = "MMN", damped = FALSE)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
##
## Initial states:
## l = 4023.4958
## b = 1.0083
##
## sigma: 0.0546
##
## AIC AICc BIC
## 289.0256 293.3113 294.0043
MMN <- ets(valid.ts, model = "MMN")
MMNfc <- forecast(MMN, h = nValid, level = 0)
#Fit Model 2.
shipments.ets.MMdN <- ets(shipments.ts, model = "MMN", damped = TRUE)
shipments.ets.MMdN
## ETS(M,Md,N)
##
## Call:
## ets(y = shipments.ts, model = "MMN", damped = TRUE)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
## phi = 0.8728
##
## Initial states:
## l = 4020.7947
## b = 1.0208
##
## sigma: 0.0488
##
## AIC AICc BIC
## 286.9602 293.4218 292.9346
MMdN <- ets(valid.ts, model = "MMN")
MMdNfc <- forecast(MMdN, h = nValid, level = 0)
#Fit Model 3.
shipments.ets.AAN.pred <- forecast(shipments.ets.AAN, h = nValid, level = c(0.2, 0.4, 0.6, 0.8))
#Accuracy AAN
accuracy(shipments.ets.AAN.pred, nValid)
## ME RMSE MAE MPE MAPE
## Training set 9.291158 219.3889 183.8221 -3.318428e-02 4.168075e+00
## Test set -4622.039089 4622.0391 4622.0391 -1.155510e+05 1.155510e+05
## MASE ACF1
## Training set 0.6962957 0.002451674
## Test set 17.5077238 NA
shipments.ets.MMN.pred <- forecast(shipments.ets.MMN, h = nValid, level = c(0.2, 0.4, 0.6, 0.8))
#Accuracy MMN
accuracy(shipments.ets.MMN.pred, nValid)
## ME RMSE MAE MPE MAPE
## Training set 32.87174 240.2824 211.2008 5.292109e-01 4.746086e+00
## Test set -4779.84338 4779.8434 4779.8434 -1.194961e+05 1.194961e+05
## MASE ACF1
## Training set 0.8000029 0.1258353
## Test set 18.1054674 NA
shipments.ets.MMdN.pred <- forecast(shipments.ets.MMdN, h = nValid, level = c(0.2, 0.4, 0.6, 0.8))
shipments.ets.MMdN.pred
## Point Forecast Lo 20 Hi 20 Lo 40 Hi 40 Lo 60
## 1990 Q1 4593.628 4534.654 4646.666 4478.961 4708.926 4407.844
## 1990 Q2 4599.058 4542.643 4652.771 4480.317 4709.938 4409.610
## 1990 Q3 4603.801 4547.746 4663.214 4488.329 4724.089 4417.389
## 1990 Q4 4607.946 4551.108 4658.680 4490.349 4720.409 4416.139
## Hi 60 Lo 80 Hi 80
## 1990 Q1 4779.928 4313.730 4884.779
## 1990 Q2 4775.458 4316.073 4876.773
## 1990 Q3 4796.393 4313.576 4896.166
## 1990 Q4 4793.088 4325.578 4889.928
#Accuracy MMdN
accuracy(shipments.ets.MMdN.pred)
## ME RMSE MAE MPE MAPE MASE
## Training set -3.6349 217.0341 177.1341 -0.3320249 4.032877 0.8303971
## ACF1
## Training set -0.01185158
par(mfrow = c(1, 3))
#This command sets the plot window to show 1 row of 3 plots.
plot(shipments.ets.AAN.pred, xlab = "Quarter", ylab = "shipments (in millions)", ylim = c(3900, 5000))
plot(shipments.ets.MMN.pred, xlab = "Quarter", ylab="shipments (in millions)", ylim = c(3900, 5000))
plot(shipments.ets.MMdN.pred, xlab = "Quarter", ylab="shipments (in millions)", ylim = c(3900, 5000))
fixed.nValid <- 4
fixed.nTrain <- length(shipments.ts) - fixed.nValid
stepsAhead <- 1
error <- rep(0, fixed.nValid - stepsAhead + 1)
percent.error <- rep(0, fixed.nValid - stepsAhead + 1)
for(j in fixed.nTrain:(fixed.nTrain + fixed.nValid - stepsAhead)) {
train.ts <- window(shipments.ts, start = c(1985, 1), end = c(1985, j))
valid.ts <- window(shipments.ts, start = c(1985, j + stepsAhead), end = c(1985, j + stepsAhead))
naive.pred <- naive(train.ts, h = stepsAhead)
error[j - fixed.nTrain + 1] <- valid.ts - naive.pred$mean[stepsAhead]
percent.error[j - fixed.nTrain + 1] <- error[j - fixed.nTrain + 1]/ valid.ts
}
mean(abs(error))
## [1] 258.75
sqrt(mean(error^2))
## [1] 364.3909
mean(abs(percent.error))
## [1] 0.0542274
fixed.nValid <- 4
fixed.nTrain <- length(shipments.ts) - fixed.nValid
train.ts <- window(shipments.ts, start = c(1985, 1), end = c(1985, fixed.nTrain))
valid.ts <- window(shipments.ts, start = c(1985, fixed.nTrain + 1), end = c(1985, fixed.nTrain + fixed.nValid))
naive.pred <- naive(train.ts, h = fixed.nValid)
snaive.pred <- snaive(train.ts, h = fixed.nValid)
accuracy(naive.pred, valid.ts)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 16.60 304.8922 265.40 0.1556967 6.077216 1.264311 -0.1617179
## Test set 307.75 385.6446 314.25 6.4985932 6.651715 1.497023 -0.4693744
## Theil's U
## Training set NA
## Test set 1.049607
accuracy(snaive.pred, valid.ts)
## ME RMSE MAE MPE MAPE MASE
## Training set 130.9167 252.3287 209.9167 2.8592405 4.665890 1.000000
## Test set 48.5000 243.1820 223.5000 0.8867438 5.009241 1.064708
## ACF1 Theil's U
## Training set 0.08095064 NA
## Test set 0.05467639 0.436126
Quad Error
Quad <- hist(shipments.lm.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="", ylim = c(0,8))
# Use stored hist object to help set up density curve
multiplier <- Quad$counts / Quad$density
# Need to ignore NA from 1985
Quadmydensity <- density(shipments.lm.pred$residuals, na.rm=TRUE)
Quadmydensity$y <- Quadmydensity$y * multiplier[1]
# Add the density curve
lines(Quadmydensity)
AAN Error
AANhist <- hist(shipments.ets.AAN.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="")
# Use stored hist object to help set up density curve
multiplier <- AANhist$counts / AANhist$density
# Need to ignore NA from 1985
AANmydensity <- density(shipments.ets.AAN.pred$residuals, na.rm=TRUE)
AANmydensity$y <- AANmydensity$y * multiplier[1]
# Add the density curve
lines(AANmydensity)
MMN Error
MMNhist <- hist(shipments.ets.MMN.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="")
# Use stored hist object to help set up density curve
multiplier <- MMNhist$counts / MMNhist$density
# Need to ignore NA from 1985
MMNmydensity <- density(shipments.ets.MMN.pred$residuals, na.rm=TRUE)
MMNmydensity$y <- MMNmydensity$y * multiplier[1]
# Add the density curve
lines(MMNmydensity)
MMdN Error
MMdNhist <- hist(shipments.ets.MMdN.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="")
# Use stored hist object to help set up density curve
multiplier <- MMdNhist$counts / MMdNhist$density
# Need to ignore NA from 1985
MMdNmydensity <- density(shipments.ets.MMdN.pred$residuals, na.rm=TRUE)
MMdNmydensity$y <- AANmydensity$y * multiplier[1]
# Add the density curve
lines(MMdNmydensity)
Naive Error
Naivehist <- hist(naive.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="")
multiplier <- Naivehist$counts / Naivehist$density
# Need to ignore NA from 1985
Naivemydensity <- density(naive.pred$residuals, na.rm=TRUE)
Naivemydensity$y <- Naivemydensity$y * multiplier[1]
# Add the density curve
lines(Naivemydensity)
SNaive Error
SNaivehist <- hist(snaive.pred$residuals, ylab = "Frequency", xlab = "Forecast Error", bty = "l", main ="")
multiplier <- SNaivehist$counts / SNaivehist$density
SNaivemydensity <- density(snaive.pred$residuals, na.rm=TRUE)
SNaivemydensity$y <- SNaivemydensity$y * multiplier[1]
lines(SNaivemydensity)
Qu1 <- mean(4245, 4258, 4595, 4409.468)
Qu2 <- mean(4900, 4258, 4799, 4355.434)
Qu3 <- mean(4585, 4258, 4417, 4291.984)
Qu4 <- mean(4533,4258,4258,4219.148)
Qu1
## [1] 4245
Qu2
## [1] 4900
Qu3
## [1] 4585
Qu4
## [1] 4533
Combined <- c(4245,4990, 4585, 4533)
Combined.ts <- ts(Combined, start = c(1989,1), end = c(1989, 4), freq = 4)
plot(valid.ts, bty="l", xaxt="n", xlab="The Year 1989", yaxt="n", ylab="Shipments", ylim = c(4000, 8000))
axis(1, at=seq(1989,1989.75,0.25), labels=c("Quarter 1", "Quarter 2", "Quarter 3", "Quarter 4"))
axis(2, las=2)
# Now add the forecasts and make the line red and dashed
lines(naive.pred$mean, col=2, lty=1)
# Add a legend
lines(snaive.pred$mean, col=3, lty=1)
lines(shipments.lm.pred$mean, col = 4, lty = 1)
lines(Combined.ts, col = 6, lty = 2)
# Add a legend
legend(1989,8000, c("Actual","Naive", "SNaive", "Quad Trend", "Combined"), col=1:6, lty=1:5)
Four Forecasts
head(valid.ts)
## [1] 4245 4900 4585 4533
head(naive.pred$mean)
## [1] 4258 4258 4258 4258
head(snaive.pred$mean)
## [1] 4595 4799 4417 4258
head(shipments.lm.pred$mean)
## 1 2 3 4
## 4409.468 4355.424 4291.984 4219.148