1 Description of the Dataset

This dataset contains the sale of Perrin Freres Champagne in millions of dollars every month from January, 1964 to September, 1972.

champagne = na.omit(read.csv("C:/Users/qinfa/Desktop/school/STA 321/champagne.csv"))

2 Create the Time Series Object

We will proceed by holding back the last 12 data points for testing our final model and creating the time series object.

train.champagne <- champagne[,2][1:93]
test.champagne <- champagne[,2][94:105]

champagne.ts = ts(champagne[,2][1:93], start = c(1964, 1), frequency = 12)
par(mar=c(2,2,2,2))
plot(champagne.ts, main="Champagne Sales from January, 1964 to September, 1972 (Millions of Dollars)", ylab="Monthly Rate", xlab="")

3 Fit Three Types of Smoothing Models

We continue by creating three types of smoothing models and fitting them.

fit1 = ses(champagne.ts, h=12)           
fit2 = holt(champagne.ts,damped=TRUE, h=12 )                   
fit3 = hw(champagne.ts, h=12, seasonal="additive",damped=TRUE)
fit4 = hw(champagne.ts, h=12, seasonal="multiplicative",damped=TRUE)

accuracy.table = round(rbind(accuracy(fit1), accuracy(fit2), accuracy(fit3), accuracy(fit4)),4)
row.names(accuracy.table)=c("SES","Holt Add. Damped", "HW Add. Damp", "HW Exp. Damp")
kable(accuracy.table, caption = "The accuracy measures of various exponential smoothing models 
      based on the training data")
The accuracy measures of various exponential smoothing models based on the training data
ME RMSE MAE MPE MAPE MASE ACF1
SES 275.3871 2454.6724 1666.6698 -15.8673 38.9449 2.5287 0.4216
Holt Add. Damped -1.8720 2530.9036 1730.2236 -17.8534 44.4789 2.6251 0.0692
HW Add. Damp 7.6758 786.7418 557.9193 -3.0582 13.7124 0.8465 0.2598
HW Exp. Damp 25.2781 569.8484 416.6006 -1.3267 10.3711 0.6321 -0.0469

From this table, we see that the Holts-Winter model with multiplicative seasonality and damping gave us the best result. This is consistent with our plot of the data and its patterns, as the peaks seem to grow larger with each year.

We continue and plot these models.

par(mfrow=c(1,1), mar=c(3,4,3,1))
###### plot the original data
pred.id = 94:105
plot(1:93, train.champagne, lwd=2,type="o", ylab="Champagne Sales", xlab="", 
     xlim=c(1, 140), ylim=c(1950, 14000), cex=0.5,
     main="Non-seasonal Smoothing Models")
lines(pred.id, fit1$mean, col="red")
lines(pred.id, fit2$mean, col="blue")
lines(pred.id, fit3$mean, col="purple")
lines(pred.id, fit4$mean, col="navy")
##
points(pred.id, fit1$mean, pch=16, col="red", cex = 0.5)
points(pred.id, fit2$mean, pch=17, col="blue", cex = 0.5)
points(pred.id, fit3$mean, pch=19, col="purple", cex = 0.5)
points(pred.id, fit4$mean, pch=21, col="navy", cex = 0.5)
#points(fit0, col="black", pch=1)
legend("bottomright", lty=1, col=c("red","blue","purple", "navy"),pch=c(16,17,19,21),
   c("SES","Holt Linea Dampedr","HW Additive Damped", "HW Multiplicative Damped"), 
   cex = 0.7, bty="n")

The plot shows that the simple exponential model and Holtz additive model with dampening did relatively poorly, while the Holtz-Winter Additive Damped and Multiplicative Damped models did significantly better. This matches with the accuracy results we derived from the table.

We continue by creating a table to see the accuracy of the different exponential smoothing methods based on the test dataset.

acc.fun = function(test.data, mod.obj){
  PE=100*(test.data-mod.obj$mean)/mod.obj$mean
  MAPE = mean(abs(PE))
  ###
  E=test.data-mod.obj$mean
  MSE=mean(E^2)
  ###
  accuracy.metric=c(MSE=MSE, MAPE=MAPE)
  accuracy.metric
}

pred.accuracy = rbind(SES =acc.fun(test.data=test.champagne, mod.obj=fit1),
                      Holt.Add.Damp =acc.fun(test.data=test.champagne, mod.obj=fit2),
                      HW.Add.Damp =acc.fun(test.data=test.champagne, mod.obj=fit3),
                      HW.Exp.Damp =acc.fun(test.data=test.champagne, mod.obj=fit4))
kable(pred.accuracy, caption="The accuracy measures of various exponential smoothing models 
      based on the testing data")
The accuracy measures of various exponential smoothing models based on the testing data
MSE MAPE
SES 8519347.8 38.164400
Holt.Add.Damp 8442531.2 37.844518
HW.Add.Damp 128425.8 6.382174
HW.Exp.Damp 115041.5 6.621129

This table shows that the damped model with multiplicative and additive seasonality did very similarly, with the MSE of the multiplicative seasonal model being slightly smaller. Looking at the plot, we can see that while the peaks seemed to grow larger for the first part of the plot, they then decreased later on, as is visible for the last year of the data. This might explain the similar results.

We finish by updating the model based on the entire dataset.

ts.champagne=ts(champagne[,2][1:105], start=c(1964, 1), frequency = 12)
final.model = hw(ts.champagne,h=12, seasonal="multiplicative", damped=TRUE) 
smoothing.parameter = final.model$model$par[1:3]
kable(smoothing.parameter, caption="Estimated values of the smoothing parameters in
      Holt-Winters damped trend with multiplicative seasonality")
Estimated values of the smoothing parameters in Holt-Winters damped trend with multiplicative seasonality
x
alpha 0.3825703
beta 0.0001020
gamma 0.0001006