Time Series and Forecasting For Predict Sales
Ruli Erinton
29 June 2022
Intro
on this occasion we will do LBB to predict the sale of souvenirs for the next two years, this is done so that we prevent the occurrence of over stock or Out of Stock in our retail.
Data Preperation
Load required library.
library(dplyr) # data wrangling
library(lubridate) # date manipulation
library(forecast) # time series library
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
library(tseries) # adf.test
library(fpp) # usconsumtion
library(TSstudio) # mempercantik visualisasi timeseries
library(ggplot2)
library(tidyr)
library(padr) # complete data frame
library(zoo) # Missing value imputationLoad data to perform timeseries and forecasting, We will use the sales souvenir data in Indonesia which is in the fancy.dat file. From this data, we will use the 2 columns we need, namely the month column to show the time and sales as the Y value that we observe to create the ts object.
souvenir <- scan(file="fancy.dat")Inspect Data
after we have successfully imported our data, we will do a data inspection to find out contents our data, actually we can use the view() function to view the contents of the data but it will take time to see the whole data so we use a function that sees the head() only.
head(souvenir)#> [1] 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74colSums(is.na(data.frame(souvenir)))#> souvenir
#> 0The results of the check state that the data does not have NA data, and nothing is missing, the next step is to change the data type to ts data type.
Creating Time Series Object
to process a data so that it can be predicted we must convert the
data into a TS object so that the machine can understand and be useful
for selecting the type of model in the Forecasting that we
make.
souvenir_ts <- ts(data = souvenir, start = c(1987, 1),frequency = 12)
souvenir_ts#> Jan Feb Mar Apr May Jun Jul
#> 1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61
#> 1988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12
#> 1989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62
#> 1990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22
#> 1991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55
#> 1992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78
#> 1993 10243.24 11266.88 21826.84 17357.33 15997.79 18601.53 26155.15
#> Aug Sep Oct Nov Dec
#> 1987 3566.34 5021.82 6423.48 7600.60 19756.21
#> 1988 4752.15 5496.43 5835.10 12600.08 28541.72
#> 1989 8176.62 8573.17 9690.50 15151.84 34061.01
#> 1990 7979.25 8093.06 8476.70 17914.66 30114.41
#> 1991 12552.22 11637.39 13606.89 21822.11 45060.69
#> 1992 19888.61 23933.38 25391.35 36024.80 80721.71
#> 1993 28586.52 30505.41 30821.33 46634.38 104660.67Decomposition
Decomposition is the stage of decomposing time series data into its components. The components in the time series:
Trend (Tt): general data pattern, tends to go up or down. If there is a trend there is still a pattern, it means that there is a pattern that has not been decomposed properly.
Seasonal (St): a seasonal pattern that forms a repeating pattern over a fixed period of time
Error/Remainder/Random (Et): patterns that cannot be caught in trends and seasonality
Let’s do our data visualization by highlighting the CO2 emission level.
souvenir_ts %>% autoplot()
Insights:
- there is seasonality in the data
- The data is multiplicative as seen from the increasing seasonal pattern
- Uptrend
From the time series data, it can be seen that it has a multiplecative pattern because the variance changes -> the ‘valley/hill’ deviation increases/decreases
souvenir_decomposition <- decompose(souvenir_ts, type = "multiplicative")
souvenir_decomposition %>% autoplot()
From the plot results above, it can be said that the seasonal data pattern is multiplecative, so I can assume that the data has a seasonal pattern, while the trend data pattern also does not have an uptrend.
Splitting Data
We will split the data to create a learning model and test the data with unseen data on souvenir sales.
# test using `tail()`
souvenir_test <- tail(souvenir_ts, n=12)
# train using `head()`
souvenir_train <- head(souvenir_ts, n = -12)
souvenir_train#> Jan Feb Mar Apr May Jun Jul Aug
#> 1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61 3566.34
#> 1988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12 4752.15
#> 1989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62 8176.62
#> 1990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22 7979.25
#> 1991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55 12552.22
#> 1992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78 19888.61
#> Sep Oct Nov Dec
#> 1987 5021.82 6423.48 7600.60 19756.21
#> 1988 5496.43 5835.10 12600.08 28541.72
#> 1989 8573.17 9690.50 15151.84 34061.01
#> 1990 8093.06 8476.70 17914.66 30114.41
#> 1991 11637.39 13606.89 21822.11 45060.69
#> 1992 23933.38 25391.35 36024.80 80721.71Model Fitting
##Simple Moving Average
A method that uses a moving average for forecasting. Because it uses an average, the weights used are the same for each observation in the past. This method is often used for data that does not contain a trend and is seasonal (the data moves around the average).
We use two models for comparison, namely Simple Moving Average and Exponential Smoothing, then we make a model with Simple Moving Average.
model_sma <- SMA(souvenir_train, n = 10)
model_sma#> Jan Feb Mar Apr May Jun Jul
#> 1987 NA NA NA NA NA NA NA
#> 1988 6023.286 6188.381 6535.599 6644.728 6799.855 6938.355 7054.085
#> 1989 8377.987 8467.647 8873.317 8908.661 8939.992 9127.857 9313.176
#> 1990 10735.187 10786.167 11379.049 11352.946 11378.896 11283.709 11238.514
#> 1991 10672.958 10683.004 10885.969 11045.611 11105.726 11328.749 11647.098
#> 1992 15132.390 15235.242 15818.845 15956.630 15762.231 15815.218 16324.757
#> Aug Sep Oct Nov Dec
#> 1987 NA NA 3727.929 4321.508 6057.376
#> 1988 6886.952 6676.535 5284.424 6294.451 8628.799
#> 1989 9547.328 9144.637 7259.515 8302.997 11138.835
#> 1990 11067.389 10361.511 7803.080 9002.436 11432.419
#> 1991 12054.650 11426.923 9776.171 11475.718 15334.764
#> 1992 16952.929 17164.056 15197.122 18038.099 25125.301##Exponential Smoothing
The SMA method only considers n observations in the past to forecast both trend and seasonal patterns, which tend not to be caught, so another method is needed that takes into account all past data, namely the exponential method. There are three types of exponential smoothing models:
Simple Exponential Smoothing: TS data without trend and seasonality
Holt’s Exponential (Double Exponential Smoothing): TS data without seasonality
Holt’s Winters Exponential (Triple Exponential Smoothing): TS data containing errors, trends, and seasonality
because from the Decomposition results we find that our data has seasonal patterns, trends and errors, the appropriate model is Holt’s Winters Exponential (Triple Exponential Smoothing): TS data containing errors, trends, and seasonality.
#Model Holt Winter
model_hw <- HoltWinters(souvenir_train, seasonal = "multiplicative")
model_hw#> Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
#>
#> Call:
#> HoltWinters(x = souvenir_train, seasonal = "multiplicative")
#>
#> Smoothing parameters:
#> alpha: 0.3469842
#> beta : 0.07501578
#> gamma: 0.5711478
#>
#> Coefficients:
#> [,1]
#> a 27582.5247693
#> b 807.9173318
#> s1 0.4676809
#> s2 0.6030877
#> s3 0.9592971
#> s4 0.7056491
#> s5 0.6754072
#> s6 0.7911595
#> s7 0.8912119
#> s8 0.8937350
#> s9 0.8856460
#> s10 0.9134738
#> s11 1.3925717
#> s12 2.9151023Forecasting
After we create a model with Simple Moving Average and Exponential Smoothing, we will do forecasting and also visualize the forecasting results, we will compare the results.
#Forecast Model Holt Winters
forecast_souvenir_hw <- forecast(model_hw, h = 12)
forecast_souvenir_hw#> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
#> Jan 1993 13277.67 11237.89 15317.44 10158.10 16397.24
#> Feb 1993 17609.17 15156.03 20062.32 13857.41 21360.94
#> Mar 1993 28784.94 25268.64 32301.23 23407.23 34162.65
#> Apr 1993 21744.01 18610.83 24877.19 16952.23 26535.79
#> May 1993 21357.80 18034.65 24680.95 16275.49 26440.12
#> Jun 1993 25657.32 21616.65 29698.00 19477.65 31837.00
#> Jul 1993 29622.05 24836.99 34407.12 22303.93 36940.18
#> Aug 1993 30427.98 25296.39 35559.56 22579.90 38276.06
#> Sep 1993 30868.11 25438.75 36297.48 22564.62 39171.61
#> Oct 1993 32576.03 26643.50 38508.56 23503.01 41649.05
#> Nov 1993 50786.55 41491.30 60081.80 36570.70 65002.41
#> Dec 1993 108667.82 88629.46 128706.18 78021.79 139313.85#Forecast Model SMA
forecast_souvenir_sma <- forecast(model_sma, h = 12)
forecast_souvenir_sma#> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
#> Jan 1993 24552.27 23618.93 25485.61 23124.85 25979.70
#> Feb 1993 25095.38 23748.61 26442.14 23035.68 27155.07
#> Mar 1993 26363.00 24582.70 28143.31 23640.26 29085.75
#> Apr 1993 26555.84 24416.05 28695.62 23283.32 29828.35
#> May 1993 26792.04 24295.62 29288.45 22974.09 30609.98
#> Jun 1993 27078.61 24221.43 29935.80 22708.92 31448.30
#> Jul 1993 27522.34 24283.30 30761.37 22568.66 32476.01
#> Aug 1993 28101.90 24455.49 31748.30 22525.20 33678.59
#> Sep 1993 27358.89 23480.49 31237.30 21427.39 33290.40
#> Oct 1993 22588.93 19116.46 26061.39 17278.25 27899.60
#> Nov 1993 26361.12 21993.97 30728.27 19682.14 33040.10
#> Dec 1993 35601.15 29278.39 41923.90 25931.33 45270.97Visualization of forecast results using model_hw
souvenir_train %>%
autoplot() +
autolayer(souvenir_test, series = "Data test") + #Actual Data test
autolayer(model_hw$fitted[,1], series = "Fitted Values") + # Prediction on Data Train
autolayer(forecast_souvenir_hw$mean, series = "Prediksi Data test") # Prediction on Data Test
after visualizing the black grid data is the actual data from the tarin data and the green one is the prediction result from model_hw if we look briefly the pattern follows the actual data on the train data, it’s just that there are still some points that are still off, and for testing in the test data from the visualization at the beginning of the year that you want to predict is still a bit far away, but it has the same pattern and at the end of the year you want to predict it is not much different from the actual test data, I assume the model made with holts winter forecasts well visual.
Visualization of forecast results using model_sma
souvenir_train %>%
autoplot() +
autolayer(souvenir_test, series = "Data test") + #Actual Data test
autolayer(model_sma, series = "Fitted Values") + # Prediction on Data Train
autolayer(forecast_souvenir_sma$mean, series = "Prediksi Data test") # Prediction on Data Test
From the visualization, we can see that the prediction of the model with SMA is very far compared to using Holt Winters, this is because SMA only performs average calculations over a certain period of time, in contrast to Exponential Smoothing which must take into account seasonal, trend, and error patterns so that the resulting model is far away. better visually.
Model Evalution
We have got the best model using Holt Winters, but how to see the
error there is how big the error is from our model, we will use the
accuracy function to see the MAPE in our model, whether it
is overfitting or not.
# model evaluation on train data
accuracy(model_hw$fitted[,1], souvenir_train)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set 189.1763 2200.411 1512.406 -0.1299785 13.13863 0.1829319 0.3931926# model evaluation on test data
accuracy(forecast_souvenir_hw$mean, souvenir_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -4060.199 4541.688 4060.199 -21.15116 21.15116 0.5148249 0.6278038from the results we get that the MAPE value in the test on the data train is quite good, it’s just that the test on the MAPE test data has increased which I can assume the model that we make is overfitting in this case.
What if we try to use the SARIMA model
#do a stationary check
adf.test(souvenir_train)#>
#> Augmented Dickey-Fuller Test
#>
#> data: souvenir_train
#> Dickey-Fuller = -1.2891, Lag order = 4, p-value = 0.8648
#> alternative hypothesis: stationarybecause p-value > 0.05, it can be concluded that the data is not stationary, for that we need to do differencing
diff(souvenir_train, lag = 12) %>% diff() %>% diff() %>% adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -5.9747, Lag order = 3, p-value = 0.01
#> alternative hypothesis: stationaryAfter doing the differencing 2 times, we get a p-value < 0.05, so it can be said that the data is stationary. lalu kita akan membuat model SARIMA dengan cara auto.
model_auto_sarima <- auto.arima(souvenir_train, seasonal = T)
model_auto_sarima#> Series: souvenir_train
#> ARIMA(0,1,2)(0,1,1)[12]
#>
#> Coefficients:
#> ma1 ma2 sma1
#> -0.7050 0.4720 0.6439
#> s.e. 0.1416 0.1617 0.2139
#>
#> sigma^2 = 13525127: log likelihood = -570.12
#> AIC=1148.23 AICc=1148.97 BIC=1156.54Then we will do a forecast and see the results of the visualization
forecast_sarima <- forecast(model_auto_sarima, h = 12)
forecast_sarima#> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
#> Jan 1993 26076.18 21349.82 30802.55 18847.83 33304.53
#> Feb 1993 34451.07 29523.78 39378.36 26915.43 41986.71
#> Mar 1993 41905.51 35791.16 48019.87 32554.41 51256.62
#> Apr 1993 35208.21 28102.40 42314.02 24340.81 46075.61
#> May 1993 32173.89 24198.94 40148.83 19977.26 44370.52
#> Jun 1993 36267.49 27509.24 45025.74 22872.90 49662.08
#> Jul 1993 41656.64 32179.61 51133.67 27162.77 56150.51
#> Aug 1993 44835.05 34690.04 54980.07 29319.59 60350.52
#> Sep 1993 52399.53 41627.87 63171.19 35925.70 68873.36
#> Oct 1993 52525.49 41161.70 63889.28 35146.07 69904.91
#> Nov 1993 66741.26 54814.69 78667.82 48501.15 84981.37
#> Dec 1993 118785.10 106320.57 131249.63 99722.25 137847.96souvenir_train %>%
autoplot() +
autolayer(souvenir_test, series = "Data test") + #Actual Data test
autolayer(forecast_sarima$mean, series = "Prediksi Data test") # Prediction on Data Test
from the results of the model it still doesn’t look so good compared to using Holt Winters but if the pattern is a little according to my assumptions, let’s see the MAPE value.
accuracy(forecast_sarima$mean, souvenir_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -18364.03 18579.26 18364.03 -88.8476 88.8476 0.09059368 2.157673If from MAPE it’s even worse when compared to using the Holt Winters model, my assumption is that this happens because SARIMA should be used for random data or has no pattern at all while the data we have has seasonal multiplecative patterns, increasing trends, and some errors so that Holt winters model is suitable for data types that have these 3 things.
Conclusion
From the model that we have created using Holt Winters, it is the best because visually it is in accordance with the plot in the test data and when viewed from the MAPE value, it has the lowest value from other models. and in the case of the SARIMA model, visually it follows the test data pattern but the MAPE value is still bad, for that we cannot assume the model is only visual, we also have to look at MAPE on the model.