Using the Winter-Holts method,model the Quarterly beer sales data and predict for the next 2 years.
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.
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
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.
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.
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.
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.
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.
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.
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(forecast(fit,8))
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.
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.