Using the Winter-Holts method,model the Quarterly beer sales data and predict for the next 2 years.

Load the data set

beer.sales <- read.csv("D:/PGP BA-BI Course Materials/TIME SERIES FORECASTING/GA/beer.csv")

The beer sales data is read into the R and saved in the object ‘beer.sales’.The next step is to store the data in a time series object using ts() function.

Plotting of Periodogram - to check for frequency

library(TSA)
## Warning: package 'TSA' was built under R version 3.4.3
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 3.4.3
## Loading required package: locfit
## Warning: package 'locfit' was built under R version 3.4.3
## locfit 1.5-9.1    2013-03-22
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-17. For overview type 'help("mgcv-package")'.
## Loading required package: tseries
## Warning: package 'tseries' was built under R version 3.4.3
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
periodogram(beer.sales)

The above periodogram plot is to identify the dominant periods(or frequncies) of a time series. From this plot, the time period is identified as 0.25 and the frequency is calculated as 1/0.25 = 4. Therefore, we can conclude that the behaviour in the series is quarterly

Making as Timeseries data

beersales.ts <- ts(beer.sales, frequency = 4, start = c(2000,1))
beersales.ts
##       Qtr1  Qtr2  Qtr3  Qtr4
## 2000 284.4 212.8 226.9 308.4
## 2001 262.0 227.9 236.1 320.4
## 2002 271.9 232.8 237.0 313.4
## 2003 261.4 226.8 249.9 314.3
## 2004 286.1 226.5 260.4 311.4
## 2005 294.7 232.6 257.2 339.2
## 2006 279.1 249.8 269.8 345.7
## 2007 293.8 254.7 277.5 363.4
## 2008 313.4 272.8 300.1 369.5
## 2009 330.8 287.8 305.9 386.1
## 2010 335.2 288.0 308.3 402.3
## 2011 352.8 316.1 324.9 404.8
## 2012 393.0 318.9 327.0 442.3
## 2013 383.1 331.6 361.4 445.9
## 2014 386.6 357.2 373.6 466.2
## 2015 409.6 369.8 378.6 487.0
## 2016 419.2 376.7 392.8 506.1
## 2017 458.4 387.4 426.9 525.0

Assuming the data is from the year 2000, the start year is set as 2000 with 1 as increasing year.

Plotting Time Series

plot.ts(beersales.ts,main = "Timeseries of Beer Sales", col = "blue")
abline(reg = lm(beersales.ts~time(beersales.ts)))

We can see from the plot that there is a linear trend and the magnitude of seasonal pattern is almost constant as the time increases, hence the timeseries can be described using Additive Holt-Winter’s method and also forecast future sales.

Test for Stationarity

kpss.test(beersales.ts)
## Warning in kpss.test(beersales.ts): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  beersales.ts
## KPSS Level = 3.0456, Truncation lag parameter = 1, p-value = 0.01

The above test validates that the given data is not stationary,hence the time series has to be decomposed i.e trend and seasonal components has to be separated to forecast the beer sales.

Applying Holt- winter’s additive method

library(forecast)
## Warning: package 'forecast' was built under R version 3.4.3
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:nlme':
## 
##     getResponse
fit <- hw(beersales.ts, seasonal = "additive")
fitted(fit)
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2000 264.5711 222.4177 235.8788 314.6716
## 2001 269.1815 219.2323 232.9907 312.6972
## 2002 268.0774 222.4238 236.2082 317.5723
## 2003 272.1912 227.0247 238.0433 318.7928
## 2004 272.0273 230.6423 244.4301 323.1572
## 2005 280.1063 235.9443 254.1045 327.2834
## 2006 291.0935 242.5409 262.7210 338.2571
## 2007 297.8437 253.8770 274.2161 349.9664
## 2008 308.1938 266.3948 288.5314 367.9393
## 2009 325.1499 284.5868 308.4243 385.7316
## 2010 343.9186 302.0324 323.1390 399.0237
## 2011 354.2837 310.6274 332.0563 411.6082
## 2012 365.1922 324.7481 343.7037 422.8812
## 2013 384.8971 336.8464 353.6417 441.2810
## 2014 398.9884 349.6714 369.6807 456.7315
## 2015 411.8143 367.3324 387.0117 474.6557
## 2016 427.8514 383.8941 400.7029 491.7647
## 2017 440.7737 398.7928 415.8050 512.9939

The data is smoothed by applying Holt-Winter’s additive method.Above is the smoothed or predicted values of the given data.

Plotting the smoothed data

plot.ts(beersales.ts,main = "Smoothed Timeseries of Beer Sales", col = "blue")
lines(fitted(fit),col = "red")

The original data and data smoothed with Holt-Winter’s method is plotted.

Estimates of model parameters

fit$model
## Holt-Winters' additive method 
## 
## Call:
##  hw(y = beersales.ts, seasonal = "additive") 
## 
##   Smoothing parameters:
##     alpha = 0.0395 
##     beta  = 0.0395 
##     gamma = 0.1854 
## 
##   Initial states:
##     l = 258.2137 
##     b = 0.0239 
##     s=55.077 -23.9983 -37.4121 6.3335
## 
##   sigma:  9.738
## 
##      AIC     AICc      BIC 
## 653.6687 656.5719 674.1587

From the above results, we are able to find the overall smoothing parmater (alpha), trend smoothing parameter(beta) and seasonal smoothing parameter(gamma). The initial values of level, trend and seasonality is interpreted.

Beer Sales : Forecast for next two years

forecast(fit,8)
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2018 Q1       463.6812 451.2014 476.1609 444.5951 482.7672
## 2018 Q2       416.0852 403.5665 428.6039 396.9395 435.2309
## 2018 Q3       438.6314 426.0255 451.2373 419.3523 457.9105
## 2018 Q4       536.0253 523.2657 548.7849 516.5112 555.5394
## 2019 Q1       484.0115 470.3853 497.6376 463.1721 504.8509
## 2019 Q2       436.4155 422.4713 450.3597 415.0897 457.7413
## 2019 Q3       458.9617 444.5959 473.3275 436.9911 480.9323
## 2019 Q4       556.3556 541.4571 571.2541 533.5704 579.1409

The forecasted values of beer sales for the next two years i.e assuming it is for 2018,2019 is displayaed. The values are calculated at 80% and 95% confidence interval.

Plot of forecasted beer sales

plot(forecast(fit,8))

Decomposing the additive time series data

states <- fit$model$states[,1:3]
colnames(states) <- cbind('Level','Trend','Seasonality')
plot(states,col = "blue", main = "Decompostion of time series")

The above shows the plot of estimates of level,trend and seasonal component of the time series data.

Measuring Forecast Accuracy

Residual Plot

plot(residuals(fit))

Thus, we can see that the above residual plot doesn’t show any pattern in it. Hence, we can conclude that the forecasted values are correct.