## R Packages Utilized
library(readr)
library(ggplot2)
library(forecast)
library(fGarch)
library(tseries)
library(bayesGARCH)
library(rugarch)
## Importing the data set
## Source: https://archive.ics.uci.edu/ml/datasets/Daily+Demand+Forecasting+Orders
## Abstract: The dataset was collected during 60 days, this is a real database of a brazilian logistics company.
## Added column called DayTotal to capture total number of days [=(Week-1)*7+Day]
library(readr)
bank <- read_csv("C:/Users/bryce_anderson/Desktop/Boston College/Predictive Analytics and Forecasting/Week 3 (PA&F)/Daily_Demand_Forecasting_Orders3.csv")
## Creating a time series of the data set
bankTS <- ts(bank)
library(ggplot2)
attach(bank)
plot(Target~DayTotal, data=bank, xlab="Total Number of Days", ylab="Target (Total Orders)", main="Target (Total Orders) over Number of Days", col="DodgerBlue")
abline(lm(Target~DayTotal), col="orange")
## Creating target forecast
bankFTarget <- (bankTS[,"Target"])
bankFTarget
## Time Series:
## Start = 1
## End = 60
## Frequency = 1
## [1] 539.577 224.675 129.412 317.120 210.517 207.364 263.043 248.958
## [9] 344.291 248.428 281.420 243.568 308.178 363.402 336.872 246.992
## [17] 308.880 233.126 404.380 298.560 229.249 236.304 297.174 409.401
## [25] 231.035 238.826 235.598 242.112 490.790 289.657 298.459 323.603
## [33] 616.453 346.035 307.645 253.847 530.944 333.359 306.356 416.830
## [41] 415.187 268.002 234.503 234.724 230.064 357.394 259.246 244.235
## [49] 402.607 255.061 342.606 268.640 188.601 202.022 213.509 316.849
## [57] 286.412 303.447 304.950 331.900
## Set a lag of 5 to represent a week of orders
autoplot(bankFTarget) + ylab("Daily Target Orders") + xlab("Day")
bankFTarget %>% diff(lag=5) %>% ggtsdisplay()
bankFTarget %>% diff(lag=5) %>% diff() %>% ggtsdisplay()
bankFTarget %>%
Arima(order=c(0,1,1), seasonal=c(0,1,1)) %>%
residuals() %>% ggtsdisplay()
arima1 <- Arima(bankFTarget, order=c(0,1,3), seasonal=c(0,1,1))
checkresiduals(arima1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,3)
## Q* = 7.6437, df = 7, p-value = 0.3651
##
## Model df: 3. Total lags used: 10
arima1 %>% forecast(h=12) %>% autoplot()
auto.arima(bankFTarget)
## Series: bankFTarget
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 300.8733
## s.e. 11.4708
##
## sigma^2 estimated as 8029: log likelihood=-354.35
## AIC=712.71 AICc=712.92 BIC=716.9
# Generate forecasts from an ETS model
fit.ets <- ets(bankFTarget)
fit.ets
## ETS(A,N,N)
##
## Call:
## ets(y = bankFTarget)
##
## Smoothing parameters:
## alpha = 1e-04
##
## Initial states:
## l = 300.8709
##
## sigma: 90.3757
##
## AIC AICc BIC
## 790.1036 790.5322 796.3866
checkresiduals(fit.ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 7.467, df = 8, p-value = 0.4872
##
## Model df: 2. Total lags used: 10
ets1 <- bankFTarget %>% ets() %>% forecast(h=12) %>% autoplot()
ets1
gBank <- garch(bankFTarget)
##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 7.225673e+03 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 5.072e+02
## 1 2 3.703e+02 2.70e-01 2.82e+00 6.8e-05 1.4e+03 1.0e+00 2.02e+03
## 2 4 3.701e+02 3.97e-04 5.43e-04 3.5e-06 1.9e+00 5.0e-02 2.05e-02
## 3 5 3.700e+02 2.67e-04 1.57e-04 2.2e-06 0.0e+00 3.3e-02 1.57e-04
## 4 6 3.699e+02 3.30e-04 3.86e-04 7.5e-06 2.1e+00 1.1e-01 5.95e-04
## 5 7 3.697e+02 5.47e-04 2.82e-04 4.5e-06 0.0e+00 8.4e-02 2.82e-04
## 6 9 3.691e+02 1.77e-03 1.61e-03 1.9e-05 0.0e+00 3.3e-01 2.84e-03
## 7 11 3.688e+02 8.34e-04 9.28e-04 7.3e-06 1.8e+00 1.2e-01 6.41e-03
## 8 14 3.679e+02 2.45e-03 2.50e-03 2.6e-05 3.2e-01 4.9e-01 7.31e-03
## 9 17 3.678e+02 1.41e-04 2.54e-04 4.7e-07 4.6e+00 8.8e-03 8.28e-04
## 10 18 3.678e+02 2.35e-05 2.20e-05 5.2e-07 2.0e+00 8.8e-03 9.06e-04
## 11 21 3.677e+02 2.78e-04 4.42e-04 7.5e-06 9.3e-01 1.5e-01 9.62e-04
## 12 24 3.677e+02 1.58e-05 5.20e-05 1.1e-07 3.5e+00 2.3e-03 1.37e-04
## 13 25 3.677e+02 7.91e-06 7.41e-06 1.4e-07 1.9e+00 2.3e-03 3.52e-05
## 14 27 3.677e+02 5.39e-06 5.20e-06 2.3e-07 9.8e-01 4.6e-03 4.84e-05
## 15 29 3.677e+02 1.12e-06 1.08e-06 4.6e-08 2.0e+00 9.2e-04 1.42e-03
## 16 31 3.677e+02 2.30e-06 2.23e-06 9.3e-08 2.0e+00 1.8e-03 1.26e-02
## 17 33 3.677e+02 4.76e-07 4.64e-07 1.9e-08 2.0e+00 3.7e-04 4.47e-01
## 18 35 3.677e+02 9.59e-07 9.34e-07 3.9e-08 2.0e+00 7.4e-04 1.79e+00
## 19 37 3.677e+02 1.86e-07 1.85e-07 8.3e-09 1.7e+01 1.5e-04 2.34e+00
## 20 39 3.677e+02 3.49e-07 3.43e-07 1.7e-08 2.1e+00 3.0e-04 1.90e+00
## 21 41 3.677e+02 7.23e-07 7.09e-07 3.5e-08 2.0e+00 5.9e-04 2.13e+00
## 22 44 3.677e+02 2.26e-08 2.26e-08 8.2e-10 1.4e+01 1.2e-05 2.38e+00
## 23 46 3.677e+02 4.33e-09 4.33e-09 1.6e-10 7.3e+01 2.4e-06 2.45e+00
## 24 48 3.677e+02 8.59e-10 8.59e-10 3.3e-11 3.7e+02 4.7e-07 2.45e+00
## 25 50 3.677e+02 1.71e-10 1.71e-10 6.5e-12 1.8e+03 9.5e-08 2.45e+00
## 26 53 3.677e+02 1.37e-09 1.37e-09 5.2e-11 5.9e+01 7.6e-07 2.45e+00
## 27 56 3.677e+02 2.73e-11 2.73e-11 1.0e-12 2.0e+00 1.5e-08 -1.30e-03
## 28 58 3.677e+02 5.46e-11 5.46e-11 2.1e-12 2.0e+00 3.0e-08 -1.30e-03
## 29 60 3.677e+02 1.09e-11 1.09e-11 4.2e-13 2.0e+00 6.1e-09 -1.30e-03
## 30 62 3.677e+02 2.18e-11 2.18e-11 8.4e-13 2.0e+00 1.2e-08 -1.30e-03
## 31 64 3.677e+02 4.37e-12 4.37e-12 1.7e-13 2.0e+00 2.4e-09 -1.30e-03
## 32 66 3.677e+02 8.74e-12 8.74e-12 3.3e-13 2.0e+00 4.8e-09 -1.30e-03
## 33 68 3.677e+02 1.75e-12 1.75e-12 6.7e-14 2.0e+00 9.7e-10 -1.30e-03
## 34 70 3.677e+02 3.50e-13 3.49e-13 1.3e-14 2.0e+00 1.9e-10 -1.30e-03
## 35 72 3.677e+02 6.98e-13 6.99e-13 2.7e-14 2.0e+00 3.9e-10 -1.30e-03
## 36 74 3.677e+02 1.40e-12 1.40e-12 5.4e-14 2.0e+00 7.8e-10 -1.30e-03
## 37 76 3.677e+02 -2.72e+07 2.80e-13 1.1e-14 2.0e+00 1.6e-10 -1.30e-03
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION 3.676807e+02 RELDX 1.072e-14
## FUNC. EVALS 76 GRAD. EVALS 37
## PRELDF 2.795e-13 NPRELDF -1.303e-03
##
## I FINAL X(I) D(I) G(I)
##
## 1 7.225675e+03 1.000e+00 1.298e-05
## 2 2.257120e-11 1.000e+00 6.624e-01
## 3 9.274881e-01 1.000e+00 -2.606e-02
gBank
##
## Call:
## garch(x = bankFTarget)
##
## Coefficient(s):
## a0 a1 b1
## 7.226e+03 2.257e-11 9.275e-01
plot(gBank)
summary(gBank)
##
## Call:
## garch(x = bankFTarget)
##
## Model:
## GARCH(1,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## 0.4100 0.7526 0.9073 1.0537 1.9528
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 7.226e+03 NA NA NA
## a1 2.257e-11 NA NA NA
## b1 9.275e-01 NA NA NA
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 38.971, df = 2, p-value = 3.448e-09
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.32515, df = 1, p-value = 0.5685
gBankFit <- garchFit(formula=~garch(1,1), bankFTarget)
##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 1)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 1
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: norm
## h.start: 2
## llh.start: 1
## Length of Series: 60
## Recursion Init: mci
## Series Scale: 89.60204
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -33.57884644 33.57885 3.357885 TRUE
## omega 0.00000100 100.00000 0.100000 TRUE
## alpha1 0.00000001 1.00000 0.100000 TRUE
## gamma1 -0.99999999 1.00000 0.100000 FALSE
## beta1 0.00000001 1.00000 0.800000 TRUE
## delta 0.00000000 2.00000 2.000000 FALSE
## skew 0.10000000 10.00000 1.000000 FALSE
## shape 1.00000000 10.00000 4.000000 FALSE
## Index List of Parameters to be Optimized:
## mu omega alpha1 beta1
## 1 2 3 5
## Persistence: 0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 85.333452: 3.35788 0.100000 0.100000 0.800000
## 1: 85.243691: 3.35514 0.100365 0.0913358 0.795148
## 2: 85.114007: 3.35081 0.119374 0.0808993 0.799501
## 3: 84.781066: 3.35152 0.136544 0.0410139 0.790682
## 4: 84.553231: 3.36789 0.165345 0.0119589 0.806992
## 5: 84.527582: 3.36267 0.162669 0.00680226 0.804288
## 6: 84.514483: 3.35094 0.165856 0.00529963 0.808631
## 7: 84.503716: 3.35793 0.164672 1.00000e-08 0.810563
## 8: 84.499437: 3.35818 0.166334 0.000386139 0.811999
## 9: 84.499398: 3.36305 0.166872 1.00000e-08 0.810989
## 10: 84.498678: 3.36077 0.167368 1.00000e-08 0.811415
## 11: 84.498563: 3.35780 0.167147 1.00000e-08 0.811172
## 12: 84.498493: 3.35815 0.167924 1.00000e-08 0.810647
## 13: 84.498448: 3.35837 0.170780 1.00000e-08 0.807387
## 14: 84.498446: 3.35840 0.171371 1.00000e-08 0.806771
## 15: 84.498446: 3.35840 0.171390 1.00000e-08 0.806752
##
## Final Estimate of the Negative LLH:
## LLH: 354.2211 norm LLH: 5.903686
## mu omega alpha1 beta1
## 3.009195e+02 1.376011e+03 1.000000e-08 8.067518e-01
##
## R-optimhess Difference Approximated Hessian Matrix:
## mu omega alpha1 beta1
## mu -7.671016e-03 -3.171017e-06 -0.2202708 -0.01131133
## omega -3.171017e-06 -1.368050e-05 -0.1039463 -0.09771529
## alpha1 -2.202708e-01 -1.039463e-01 -979.1185298 -764.89878455
## beta1 -1.131133e-02 -9.771529e-02 -764.8987845 -703.81379434
## attr(,"time")
## Time difference of 0.01504111 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
## Time difference of 0.034091 secs
gBankFit
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = bankFTarget)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x0000000017ead1c0>
## [data = bankFTarget]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## 3.0092e+02 1.3760e+03 1.0000e-08 8.0675e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 3.009e+02 1.187e+01 25.344 <2e-16 ***
## omega 1.376e+03 3.611e+03 0.381 0.703
## alpha1 1.000e-08 1.021e-01 0.000 1.000
## beta1 8.068e-01 5.731e-01 1.408 0.159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -354.2211 normalized: -5.903686
##
## Description:
## Tue Nov 26 22:28:30 2019 by user: bryce_anderson
##plot(gBankFit) ## This function does not run in R Markdown but allows you to plot different selections
## Make a plot selection (or 0 to exit):
## 1: Time Series 2: Conditional SD
## 3: Series with 2 Conditional SD Superimposed 4: ACF of Observations
## 5: ACF of Squared Observations 6: Cross Correlation
## 7: Residuals 8: Conditional SDs
## 9: Standardized Residuals 10: ACF of Standardized Residuals
## 11: ACF of Squared Standardized Residuals 12: Cross Correlation between r^2 and r
## 13: QQ-Plot of Standardized Residuals
## Bayes Garch Model
The results from the ETS/ARIMA models and the GARCH model are similar but show some small differences in both. The residual time series appears to be very similar while the histogram of the residuals differs some. We also see some similarities between the ACF graphs in each, but the GARCH model shows different lag peaks. The GARCH models produces a different output in it’s summary and provides some interesting tests including the Box-L Jung Test with squared residual values.