Library Call
We only need 4 libraries to do any Holt Winters and ARIMA modelling on the NY birth data.
library(tidyverse) #untuk mengolah data
library(forecast) #karena error saat buat ggplot: "Objects of type ts not supported by autoplot."
library(MLmetrics) #MAPE()
library(tseries) #adf.testing for objective differencing decisionsData Reading
birth <- read.csv("data_input/nybirth.csv")
head(birth)#> date births
#> 1 1946-01-01 26.663
#> 2 1946-02-01 23.598
#> 3 1946-03-01 26.931
#> 4 1946-04-01 24.740
#> 5 1946-05-01 25.806
#> 6 1946-06-01 24.364Data Pre-processing
Converting Dataframe to TS Class
birth_ts <- ts(data = birth$births , start = 1946, frequency = 12)
birth_ts#> Jan Feb Mar Apr May Jun Jul Aug Sep Oct
#> 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
#> 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
#> 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
#> 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
#> 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
#> 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
#> 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
#> 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
#> 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
#> 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
#> 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
#> 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
#> 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
#> 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
#> Nov Dec
#> 1946 21.672 21.870
#> 1947 21.759 22.073
#> 1948 21.059 21.573
#> 1949 21.519 22.025
#> 1950 22.084 22.991
#> 1951 22.964 23.981
#> 1952 23.162 24.707
#> 1953 25.246 25.180
#> 1954 24.712 25.688
#> 1955 25.693 26.881
#> 1956 26.291 26.987
#> 1957 26.634 27.735
#> 1958 25.912 26.619
#> 1959 26.992 27.897Cross-Validation / Train-Test Splitting
Separating the last 2 years (or 12*2 months) of data into testing data.
# train
birth_train <- head(birth_ts, -12*2)
# test
birth_test <- tail(birth_ts, 12*2)Feature Selection
Holt-Winters
Visualizing TS Data
birth_ts %>% autoplot()
Looking at the data seemingly not having any significant change in the magnitude or local maximum or minimum in each seasonal, it seems that our Holt-Winters model is best using additive, but let’s see until we try the model if this is true.
Decompose Plots
Checking both of the decompose plots of any helpful, distinctive signs of using additive of multiplicative.
birth_train %>%
decompose(type = "additive") %>%
autoplot()
birth_train %>%
decompose(type = "multiplicative") %>%
autoplot()
Holt-Winters using Additive Seasonality
# auto additive
birth_hw_auto <- HoltWinters(x = birth_train)
# manual additive
birth_hw_manual <- HoltWinters(x = birth_train, alpha = 0.9, beta = 0.1, gamma = 0.5)Holt-Winters using Multiplicative Seasonality
# auto multiplicative
birth_hw_auto_m <- HoltWinters(x = birth_train, seasonal = "multiplicative")
# manual multiplicative
birth_hw_manual_m <- HoltWinters(x = birth_train, seasonal = "multiplicative", alpha = 0.9, beta = 0.1, gamma = 0.5)SARIMA with Automatic Orders
ARIMA and Seasonal ARIMA (SARIMA) has parameters so that we can tune our model to forecast the best outcome, have the least error, and/or have the least information loss.
The parameters of ARIMA are: (p,d,q)(P,D,Q)[F] p = how many data are used to smooth the auto-regression d = how many differencing of the original data is needed to make it stationary q = how many data are used to smooth the residuals using moving average
P = same with p, but along the seasonal frequency (F) D = same with d, but for seasonal differencing, before the overall differencing is done (using seasonal frequency (F) as parameter) Q = same with q, but along the seasonal frequency (F)
F = seasonal frequency or how many times the data seems to repeat some patterns.
birth_auto_arima <- auto.arima(birth_train)
birth_auto_arima#> Series: birth_train
#> ARIMA(0,1,2)(0,1,2)[12]
#>
#> Coefficients:
#> ma1 ma2 sma1 sma2
#> -0.1054 -0.2263 -1.1221 0.3499
#> s.e. 0.0883 0.0822 0.1095 0.1127
#>
#> sigma^2 estimated as 0.4161: log likelihood=-135.08
#> AIC=280.15 AICc=280.63 BIC=294.53Pre-test model evaluation using accuracy:
accuracy(birth_auto_arima)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.1329977 0.6057793 0.450476 0.5482743 1.825522 0.4734038
#> ACF1
#> Training set -0.01184833MAPE of 1.85%. Not bad!
SARIMA with Manual Orders
Original Data
birth_train %>%
adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -5.6971, Lag order = 5, p-value = 0.01
#> alternative hypothesis: stationaryIt appears that doing no differencing at all is still acceptable by the ADF hypothesis test, using alpha 0.05, it shows that p-value is < alpha, therefore the data is “stationary”.
birth_train %>%
autoplot()
But since we can see some trend and seasonal inside the data, we need to try differencing the data again. We will still consider the 0 differencing in both seasonal and overall. (D = 0, d = 0)
Seasonal Frequency
We are using 12 as the frequency when converting the birth dataframe into TS class, therefore, F=12.
Differencing
D = 1
birth_train %>%
diff(lag = 12) %>%
autoplot()
Seems like the seasonal is gone. Checking with ADF.Test:
birth_train %>%
diff(lag = 12) %>%
adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -4.9415, Lag order = 5, p-value = 0.01
#> alternative hypothesis: stationaryIt’s still acceptable.
D = 1, d = 1
birth_train %>%
diff(lag = 12) %>%
diff(lag = 1) %>%
autoplot()
It looks much more centralized on the mean or X-axis, therefore should be more stationary and lacking trend and seasonality.
birth_train %>%
diff(lag = 12) %>%
diff(lag = 1) %>%
adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -6.5978, Lag order = 5, p-value = 0.01
#> alternative hypothesis: stationaryStill an acceptable result.
Therefore, we will consider D=0/1, and d=0/1
PACF and ACF Plots
Since the data is visually more stationary using one seasonal differencing and one overall differencing, therefore we will use it as the basis to look at PACF and ACF plots, to determine the rest of ARIMA parameters.
birth_train %>%
diff(lag = 12) %>%
diff(lag = 1) %>%
tsdisplay(lag.max = 36)
#### PACF Plot: Determining P and p Since we have the seasonal frequency of 12, let’s look at each lag multiples of 12 in PACF plot. 12 passes the blue line, 24 also, while 36 is barely below the line. This seems fitting for a ‘dies down’ category, therefore we will be choosing P = 2.
For p, we should look in only the first few lags. It seems more fitting as a ‘cut off’ category, therefore we will try p = 2 or p = 3 as both are passing the blue line.
ACF Plot: Determining Q and q
Again, looking for seasonal lag multiples of 12 in ACF. Seems that only 12 passes the blue line, therefore Q = 1. Similar to PACF, the first few lags in ACF that passes the blue line is only 2 and 3, so, we will try q = 2 or q = 3.
Possibilities for Manual Models
Concluding our observations, here are the possibilities of (P,D,Q)(p,d,q): P = 2 D = 0/1 Q = 1
p = 2/3 d = 0/1 q = 2/3
There are 4 different parameters with each 2 different possibilities. Trying all combinations would result in 16 models.
# birth_arima202201 <- arima(x = birth_train, order = c(2,0,2), seasonal = c(2,0,1)) #this results in an error "non-finite finite-difference value"
birth_arima203201 <- arima(x = birth_train, order = c(2,0,3), seasonal = c(2,0,1))
# birth_arima212201 <- arima(x = birth_train, order = c(2,1,2), seasonal = c(2,0,1)) #this results in an error "non-stationary seasonal AR part from CSS"
birth_arima213201 <- arima(x = birth_train, order = c(2,1,3), seasonal = c(2,0,1))
birth_arima302201 <- arima(x = birth_train, order = c(3,0,2), seasonal = c(2,0,1))
birth_arima303201 <- arima(x = birth_train, order = c(3,0,3), seasonal = c(2,0,1))
birth_arima312201 <- arima(x = birth_train, order = c(3,1,2), seasonal = c(2,0,1))
birth_arima313201 <- arima(x = birth_train, order = c(3,1,3), seasonal = c(2,0,1))
birth_arima202211 <- arima(x = birth_train, order = c(2,0,2), seasonal = c(2,1,1))
birth_arima203211 <- arima(x = birth_train, order = c(2,0,3), seasonal = c(2,1,1))
birth_arima212211 <- arima(x = birth_train, order = c(2,1,2), seasonal = c(2,1,1))
birth_arima213211 <- arima(x = birth_train, order = c(2,1,3), seasonal = c(2,1,1))
birth_arima302211 <- arima(x = birth_train, order = c(3,0,2), seasonal = c(2,1,1))
birth_arima303211 <- arima(x = birth_train, order = c(3,0,3), seasonal = c(2,1,1))
birth_arima312211 <- arima(x = birth_train, order = c(3,1,2), seasonal = c(2,1,1))
birth_arima313211 <- arima(x = birth_train, order = c(3,1,3), seasonal = c(2,1,1))# accuracy(birth_arima202201)
accuracy(birth_arima203201)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.04672472 0.636629 0.4989282 0.1356397 2.029055 0.4229358
#> ACF1
#> Training set 0.008008026# accuracy(birth_arima212201)
accuracy(birth_arima213201)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.1210151 0.6314041 0.4899203 0.4678223 1.980338 0.4152999
#> ACF1
#> Training set -0.0347919accuracy(birth_arima302201)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.03076363 0.6217371 0.5012895 0.0484277 2.036431 0.4249375
#> ACF1
#> Training set -0.005484707accuracy(birth_arima303201)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.0421287 0.6207129 0.4989447 0.1122587 2.02679 0.4229498
#> ACF1
#> Training set -0.04979023accuracy(birth_arima312201)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.1182959 0.631228 0.4919142 0.4567391 1.987358 0.4169901
#> ACF1
#> Training set -0.0341704accuracy(birth_arima313201)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.1262301 0.6103531 0.4739159 0.4861465 1.912912 0.4017331
#> ACF1
#> Training set -0.0402518accuracy(birth_arima202211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.06024549 0.6001871 0.4512522 0.1877161 1.838544 0.3825214
#> ACF1
#> Training set -0.04604288accuracy(birth_arima203211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.07070011 0.5910344 0.4460204 0.2495491 1.818078 0.3780864
#> ACF1
#> Training set -0.1041361accuracy(birth_arima212211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.09770979 0.5786945 0.4387213 0.4022654 1.781933 0.3718991
#> ACF1
#> Training set -0.0948942accuracy(birth_arima213211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.1168173 0.5847007 0.437083 0.4778213 1.775556 0.3705104
#> ACF1
#> Training set -0.01359231accuracy(birth_arima302211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.04899447 0.5774071 0.4458807 0.1625226 1.814789 0.3779681
#> ACF1
#> Training set -0.1277111accuracy(birth_arima303211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.06857755 0.586897 0.4443339 0.2410027 1.811111 0.3766569
#> ACF1
#> Training set -0.09678015accuracy(birth_arima312211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.09476621 0.5693816 0.4233738 0.3806519 1.722433 0.3588892
#> ACF1
#> Training set -0.02140837accuracy(birth_arima313211)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.123292 0.5692868 0.4284582 0.5079665 1.738566 0.3631992
#> ACF1
#> Training set -0.006509432Only 4 of the first models have such higher errors of above 4% (or error when fitting the models) much more than the rest, which is 1.7 - 2.0% MAPE error.
Forecasting
Holt-Winters Models
birth_fc_auto <- forecast(object = birth_hw_auto, h = 2*12)
birth_fc_manual <- forecast(object = birth_hw_manual, h = 2*12)
birth_fc_auto_m <- forecast(object = birth_hw_auto_m, h = 2*12)
birth_fc_manual_m <- forecast(object = birth_hw_manual_m, h = 2*12)ARIMA Models
birth_fc_arima_auto <- forecast(birth_auto_arima, h = 2*12)
# birth_fc_arima202201 <- forecast(birth_arima202201, h = 2*12)
birth_fc_arima203201 <- forecast(birth_arima203201, h = 2*12)
# birth_fc_arima212201 <- forecast(birth_arima212201, h = 2*12)
birth_fc_arima213201 <- forecast(birth_arima213201, h = 2*12)
birth_fc_arima302201 <- forecast(birth_arima302201, h = 2*12)
birth_fc_arima303201 <- forecast(birth_arima303201, h = 2*12)
birth_fc_arima312201 <- forecast(birth_arima312201, h = 2*12)
birth_fc_arima313201 <- forecast(birth_arima313201, h = 2*12)
birth_fc_arima202211 <- forecast(birth_arima202211, h = 2*12)
birth_fc_arima203211 <- forecast(birth_arima203211, h = 2*12)
birth_fc_arima212211 <- forecast(birth_arima212211, h = 2*12)
birth_fc_arima213211 <- forecast(birth_arima213211, h = 2*12)
birth_fc_arima302211 <- forecast(birth_arima302211, h = 2*12)
birth_fc_arima303211 <- forecast(birth_arima303211, h = 2*12)
birth_fc_arima312211 <- forecast(birth_arima312211, h = 2*12)
birth_fc_arima313211 <- forecast(birth_arima313211, h = 2*12)Model Evaluation
Holt-Winters Models
MAPE(y_pred = birth_fc_auto$mean, y_true = birth_test)*100#> [1] 4.674416MAPE(y_pred = birth_fc_manual$mean, y_true = birth_test)*100#> [1] 4.479165MAPE(y_pred = birth_fc_auto_m$mean, y_true = birth_test)*100#> [1] 5.330195MAPE(y_pred = birth_fc_manual_m$mean, y_true = birth_test)*100#> [1] 4.922051The best MAPE was predicted at 4.47% by the manual parameters of Holt-Winters (alpha = 0.9, beta = 0.1, gamma = 0.5).
ARIMA Models
MAPE(y_pred = birth_fc_arima_auto$mean, y_true = birth_test)*100#> [1] 4.403916# MAPE(y_pred = birth_fc_arima202201$mean, y_true = birth_test)*100
MAPE(y_pred = birth_fc_arima203201$mean, y_true = birth_test)*100#> [1] 3.260017# MAPE(y_pred = birth_fc_arima212201$mean, y_true = birth_test)*100
MAPE(y_pred = birth_fc_arima213201$mean, y_true = birth_test)*100#> [1] 3.764099MAPE(y_pred = birth_fc_arima302201$mean, y_true = birth_test)*100#> [1] 2.35525MAPE(y_pred = birth_fc_arima303201$mean, y_true = birth_test)*100#> [1] 3.070731MAPE(y_pred = birth_fc_arima312201$mean, y_true = birth_test)*100#> [1] 3.839012MAPE(y_pred = birth_fc_arima313201$mean, y_true = birth_test)*100#> [1] 3.222905MAPE(y_pred = birth_fc_arima202211$mean, y_true = birth_test)*100#> [1] 2.42024MAPE(y_pred = birth_fc_arima203211$mean, y_true = birth_test)*100#> [1] 3.476886MAPE(y_pred = birth_fc_arima212211$mean, y_true = birth_test)*100#> [1] 5.121851MAPE(y_pred = birth_fc_arima213211$mean, y_true = birth_test)*100#> [1] 5.322314MAPE(y_pred = birth_fc_arima302211$mean, y_true = birth_test)*100#> [1] 2.811651MAPE(y_pred = birth_fc_arima303211$mean, y_true = birth_test)*100#> [1] 3.579452MAPE(y_pred = birth_fc_arima312211$mean, y_true = birth_test)*100#> [1] 8.472663MAPE(y_pred = birth_fc_arima313211$mean, y_true = birth_test)*100#> [1] 4.995938Rather than the auto-ARIMA, most of our manual parametered models seems to perform better. The Auto-ARIMA performs at 4.40% MAPE error, while birth_arima302201 or ARIMA(3,0,2)(2,0,1)[12] performs better at 2.36% error, and close by birth_arima202211 or ARIMA(2,0,2)(2,1,1)[12] with 2.42% error.
Overall
Best MAPE of Holt-Winters is model birth_fc_manual or additive HoltWinters(alpha = 0.9, beta = 0.1, gamma = 0.5) with 4.47% MAPE error. Meanwhile, our best performing (S)ARIMA model is model birth_arima302201 or ARIMA(3,0,2)(2,0,1)[12] with 2.36% MAPE error.
It seems like automatic parametered model is not better than our manually parametered ones, be it Holt-Winters or (S)ARIMA. Also, ARIMA seems to perform better to forecast the next 2 years data with this dataset, having less error than Holt-Winters.
Data Visualization of Forecast Comparing
Without Forecasted Data
Let’s compare the model’s fitted values on the training data.
birth_ts %>%
head(-12*2) %>%
autoplot()+
autolayer(object = birth_fc_manual$fitted, lwd = 1, series = "HW Add Auto")+
autolayer(object = birth_fc_arima302201$fitted, lwd = 1, series = "ARIMA (3,0,2)(2,0,1)")
Both are very close, though ARIMA seems most fitting at most times.
Only Forecasted Data
We compare only the last 2 years of data, the most recent, and should have been unseen from any of the models’ training process.
birth_ts %>%
window(start = c(1958,1), end = c(1959,12)) %>%
autoplot()+
autolayer(object = birth_fc_manual$mean, lwd = 1, series = "HW Add Auto")+
autolayer(object = birth_fc_arima302201$mean, lwd = 1, series = "ARIMA (3,0,2)(2,0,1)")
Again, ARIMA is surely most fitting model compared to our Holt-Winters. But other finely tuned models may perform better than the one we found here.