Bachelorarbeit

1 Theorieteil

1.1 R-Pakete

#R-Paket für REDWINE-Datensatz
library(tsapp) 
#R-Paket für Funktionen zur Analyse von Zeitreihen
library(forecast)

1.2 Visualisierung der Zeitreihe

#Auswahl Datensatz
data("REDWINE")
#Visualisierung der Zeitreihe für Abbildung 1
plot(REDWINE, xlab = "Jahr", ylab = "Liter (in 1000)")

#multiplikative Zeitreihenzerlegung für Abbildung 2
plot(decompose(REDWINE, type = "multiplicative"), xlab = "Jahr")

1.3 Korrelogramme und weißes Rauschen

#Korrelogramm für Abbildung 3
ggAcf(REDWINE)

#Erstellung der Zeitreihe weißes Rauschen
#set.seed() für Reproduzierbarkeit von Abbildung 4
#rnorm(50) für Normalverteilung und Zeitperiode bis 50
set.seed(30)
weißes_Rauschen <- ts(rnorm(50))
autoplot(weißes_Rauschen)

#Korrelogramm: weißes Rauschen
ggAcf(weißes_Rauschen)

1.4 Naive Prognose

#Anzeige der Zeitreihenwerte für Tabelle 1 und Tabelle 2
show(REDWINE)

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1980  464  675  703  887 1139 1077 1318 1260 1120  963  996  960
## 1981  530  883  894 1045 1199 1287 1565 1577 1076  918 1008 1063
## 1982  544  635  804  980 1018 1064 1404 1286 1104  999  996 1015
## 1983  615  722  832  977 1270 1437 1520 1708 1151  934 1159 1209
## 1984  699  830  996 1124 1458 1270 1753 2258 1208 1241 1265 1828
## 1985  809  997 1164 1205 1538 1513 1378 2083 1357 1536 1526 1376
## 1986  779 1005 1193 1522 1539 1546 2116 2326 1596 1356 1553 1613
## 1987  814 1150 1225 1691 1759 1754 2100 2062 2012 1897 1964 2186
## 1988  966 1549 1538 1612 2078 2137 2907 2249 1883 1739 1828 1868
## 1989 1138 1430 1809 1763 2200 2067 2503 2141 2103 1972 2181 2344
## 1990  970 1199 1718 1683 2025 2051 2439 2353 2230 1852 2147 2286
## 1991 1007 1665 1642 1518 1831 2207 2822 2393 2306 1785 2047 2171
## 1992 1212 1335 2011 1860 1954 2152 2835 2224 2182 1992 2389 2724
## 1993  891 1247 2017 2257 2255 2255 3057 3330 1896 2096 2374 2535
## 1994 1041 1728 2201 2455 2204 2660 3670 2665 2639 2226 2586 2684
## 1995 1185 1749 2459 2618 2585 3310 3923

1.6 ARMA-Prozesse

1.5 Prognose: Neuronale Netze

#Prognose REDWINE mithilfe neuronaler Netze (Abbildung 8)
fit <- nnetar(REDWINE, lambda=0)
autoplot(forecast(fit, h=30), xlab = "Jahr", ylab = "Liter in 1000")

2 Praxisteil

2.1 R-Pakete

library("ggplot2") #für Visualisierung
library("prettydoc") #Designs für Rmarkdown
library("zoo") #für monatliche Klassifizierung der Daten
library("fpp2") #für "seasonplots"
library("caret") #Datenaufteilung und Zeitreihenanalysen
library("dplyr")
library("data.table") #Erstellung Tabelle
library("stargazer") #Summary
library("prophet") #Prophet-Modell
library("tseries")
library("ggthemes")
library("xts")

2.2 Deskriptive Statistik

# Einlesen und Öffnen des Datensatzes
test <- read.csv("~/Desktop/Uni/7. Semester/David BA/demand-forecasting-kernels-only/test.csv")
train <- read.csv("~/Desktop/Uni/7. Semester/David BA/demand-forecasting-kernels-only/train.csv")

#Ergründung der Struktur der Trainingsdaten
str(train)

## 'data.frame':    913000 obs. of  4 variables:
##  $ date : chr  "2013-01-01" "2013-01-02" "2013-01-03" "2013-01-04" ...
##  $ store: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ item : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ sales: int  13 11 14 13 10 12 10 9 12 9 ...

#Ergründung der Struktur der Testdaten
str(test)

## 'data.frame':    45000 obs. of  4 variables:
##  $ id   : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ date : chr  "2018-01-01" "2018-01-02" "2018-01-03" "2018-01-04" ...
##  $ store: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ item : int  1 1 1 1 1 1 1 1 1 1 ...

#Summary-Tabelle
stargazer(train, type="text")

## 
## ===========================================================
## Statistic    N     Mean  St. Dev. Min Pctl(25) Pctl(75) Max
## -----------------------------------------------------------
## store     913,000 5.500   2.872    1     3        8     10 
## item      913,000 25.500  14.431   1     13       38    50 
## sales     913,000 52.250  28.801   0     30       70    231
## -----------------------------------------------------------

#Aggregation nach Store
agg = aggregate(train,
                by = list(train$store), 
                FUN = mean)

#Durchschnittliche Verkäufe pro Store pro Tag
ggplot(agg, aes(x = store, y = sales, fill = store)) +
  geom_bar(stat = "identity", fill = "steelblue")+
  geom_text(aes(label=round(sales)), color="white", vjust=1.5, size=3.5)+
  theme_classic() +
  labs(x = "Filiale",
       y = "durchschn. Verkäufe",
       title = "Durchschnittliche Verkäufe pro Tag") +
  scale_x_discrete(name = "Filiale(store)",
                   limits = c("1", "2", "3", "4", "5",
                              "6", "7", "8", "9", "10"))

#Aggregation nach item
agg2 = aggregate(train,
                by = list(train$item), 
                FUN = mean)

#Aufteilung in 4 Gruppen für übersichtliche Diagramme
agg2$group <- as.numeric(cut(agg2$item, 4))

#Erstellung Diagramm: durchschn. Verkäufe pro Item pro Tag
ggplot(agg2, aes(x = factor(item), y = sales)) + 
  geom_bar(stat = "identity", fill = "steelblue") + 
  geom_text(aes(label=round(sales)),color="white", vjust=1.5, size=2.5) +
  theme(axis.text.x = element_text( vjust = 0.5)) +
  labs(x = "item",
       y = "durchschn. Verkäufe",
       title = "Durchschnittliche Verkäufe pro item pro Tag") +
  facet_wrap(~group, ncol = 2, scales = "free_x")

#Histogramm von Sales (Häufigkeitsverteilung der sales)
ggplot(train, aes(x=sales))+
  geom_histogram(fill="steelblue", alpha=.7)+
  labs(x="sales", y="Häufigkeit", title="Histogramm von sales")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

2.3 Datenaufteilung

#Ähnlichkeiten der verschiedenen Zeitreihen

#Erstellung Zeitreihe store 2 & item 15
store2_item15 <- setDT(train)[store == 2 & item == 15]

ts_store2_item15 <- ts(store2_item15[,c("sales")],     
           start = c(2013,1),
           frequency = 365)

#Erstellung Zeitreihe store 10 & item 6
store10_item6 <- setDT(train)[store == 10 & item == 6]

ts_store10_item6 <- ts(store10_item6[,c("sales")],     
           start = c(2013,1),
           frequency = 365)

#Erstellung Zeitreihe store 4 & item 1
store4_item1 <- setDT(train)[store == 4 & item == 1]

ts_store4_item1 <- ts(store4_item1[,c("sales")],     
           start = c(2013,1),
           frequency = 365)

#Erstellung Zeitreihe store 7 & item 12
store7_item12 <- setDT(train)[store == 7 & item == 12]

ts_store7_item12 <- ts(store7_item12[,c("sales")],     
           start = c(2013,1),
           frequency = 365)

par(mfrow=c(2,2))
plot(ts_store2_item15, xlab="Jahr", ylab="durchschn. Verkäufe", main="store2_item15")
plot(ts_store10_item6, xlab="Jahr", ylab="durchschn. Verkäufe", main="store10_item6")
plot(ts_store4_item1, xlab="Jahr", ylab="durchschn. Verkäufe", main="store4_item1")
plot(ts_store7_item12, xlab="Jahr", ylab="durchschn. Verkäufe", main="store7_item12")

#Datenaufteilung 80/20
train_window =window(ts_store2_item15, start=2013, end=c(2016,365))
test_window = window(ts_store2_item15, start=2017, end=c(2017,365))
Trainingsdaten_80 <- ts(train_window, frequency = 365, start = c(2013,1))
plot(Trainingsdaten_80)

Testdaten_20 <- ts(test_window, frequency = 365, start = c(2017,1))
plot(Testdaten_20)

autoplot(Trainingsdaten_80, xlab="Jahr", ylab="durchschn. Verkäufe")+
  autolayer(Testdaten_20)+
  autolayer(Trainingsdaten_80)+
  guides(colour=guide_legend(title="Daten"))

#Korrelogramm der Trainingsdaten
#Zeitreihenzerlung der Trainingsdaten
Korrelogramm_train <- ggAcf(ts_store2_item15)

plot(Korrelogramm_train)

plot(decompose(ts_store2_item15))

2.3 Anwendung der Prognosemodelle

2.3.1 Naive Prognose

#für kürzere Prognose h-Wert anpassen 3Monate: h=90
#Naive Prognose, h= Anzahl an Prognose-Perioden
naive_prognose <- naive(Trainingsdaten_80, h = 90, c(80, 95))

autoplot(naive_prognose)+
  xlab("Jahr")+
  ylab("durchschn. Verkäufe")

accuracy(naive_prognose, Testdaten_20)

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.01987663 24.13404 17.97464 -2.7204683 18.01048 0.8307543
## Test set     4.46666667 19.96775 16.28889  0.8087245 16.31557 0.7528420
##                    ACF1 Theil's U
## Training set -0.3884003        NA
## Test set      0.3594551 0.8731591

#Anzeige Prognosewerte 
#dunkle Markierung: Zeigt an, wo der Wert in 80% der Fälle ist
#helle Markierung: Zeigt an, wo der Wert in 95% der Fälle ist
#show(naive_prognose) <- wegen Übersicht Befehl nicht durchgeführt

2.3.2 Naive Prognose (saisonal)

#Naive Prognose (saisonal), h= Anzahl an Prognose-Perioden
#Prognosewert = Wert der gleichen Saison aus Vorperiode
#für kürzere Prognose h-Wert anpassen 3Monate: h=90
naive_prognose_saisonal <- snaive(Trainingsdaten_80, h = 365, level=c(80,95))

autoplot(naive_prognose_saisonal, 
  xlab="Jahr", 
  ylab="durchschn. Verkäufe")

accuracy(naive_prognose_saisonal, Testdaten_20)

##                    ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 9.181735 26.37250 21.63653 5.546645 19.79364 1.0000000 -0.1512906
## Test set     4.778082 26.02591 20.19178 1.693608 16.83065 0.9332264 -0.2428726
##              Theil's U
## Training set        NA
## Test set     0.9787598

#Anzeige Prognosewerte 
#dunkle Markierung: Zeigt an, wo der Wert in 80% der Fälle ist
#helle Markierung: Zeigt an, wo der Wert in 95% der Fälle ist
#show(naive_prognose_saisonal) <- wegen Übersicht Befehl nicht durchgeführt

2.3.3 einfache exponentielle Glättung

#einfache exponentielle Glättung 
#für kürzere Prognose h-Wert anpassen 3 Monate: h=90
einf_exp_glättung_ts <- ses(Trainingsdaten_80, h=365, alpha=0.4, level=c(80,95))
autoplot(einf_exp_glättung_ts, xlab="Jahr", ylab="durchschn. Verkäufe", 
         title="Naive Prognose (saisonal) Store 2 Item 15")

#show(einf_exp_glättung_ts) <- Werte anzeigen lassen

accuracy(einf_exp_glättung_ts, Testdaten_20)

##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  0.03897858 20.62337 16.04304 -3.018734 15.76091 0.7414792
## Test set     38.17736637 48.60660 40.50330 25.859420 28.98618 1.8719871
##                      ACF1 Theil's U
## Training set -0.007125855        NA
## Test set      0.608289250  1.713129

2.3.4 Holt-Winters-Verfahren

#gesamte Sicht start=2013 -> Fenster (2013-2017)
#Holt-Winters-Verfahren benötigt wöchentliche Sicht
#Window(), für bessere Skalierung der Grafik
windowtsDataWeekly <- window(Trainingsdaten_80, start= 2013, end = c(2016,365))
tsDataWeekly <- ts(windowtsDataWeekly, frequency = 7)
#für kürzere Prognose h-Wert anpassen 3 Monate: h=90
fit1 <- hw(tsDataWeekly, h=365, seasonal = "additive")
fit2 <- hw(tsDataWeekly, h=365, seasonal = "multiplicative")
autoplot(tsDataWeekly) +
  autolayer(fit1, series="HW additive Prognosen", PI=FALSE) +
  autolayer(fit2, series="HW multiplicative Prognosen", PI=FALSE) +
  xlab("Woche") +
  ylab("durchschn. Verkäufe") +
  ggtitle("HW Prognose") +
  guides(colour=guide_legend(title="Prognose"))

#für hw-Verfahren werden Testdaten wöchentlich benötigt
#start=209, da 4 Jahre * 52 Wochen = 208 (Trainingsdaten)
testDataWeekly <- ts(Testdaten_20, frequency = 7, start=209)
accuracy(fit1, testDataWeekly) #additiv HW

##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set -0.1104765 12.32162  9.686705 -1.229303  9.31917 0.7475453
## Test set     33.6494184 49.97770 39.569646 21.379560 28.81701 3.0536805
##                    ACF1 Theil's U
## Training set 0.05114907        NA
## Test set     0.51342671  1.775544

accuracy(fit2, testDataWeekly) #multiplikativ HW

##                       ME     RMSE       MAE      MPE      MAPE      MASE
## Training set -0.01517778 11.98593  9.393619 -1.11151  8.980088 0.7249272
## Test set     36.54109705 50.37117 40.834554 23.92737 29.471798 3.1512963
##                    ACF1 Theil's U
## Training set 0.02684915        NA
## Test set     0.54487498  1.774124

#HoltWinters(windowtsDataWeekly)$alpha
#HoltWinters(windowtsDataWeekly)$beta
#HoltWinters(windowtsDataWeekly)$gamma
#tail(HoltWinters(windowtsDataWeekly)$fitted)
#show(fit1)
#show(fit2)

#Detailsicht start=2016 -> kleineres Fenster (2016-2017)
#Holt-Winters-Verfahren benötigt wöchentliche Sicht
#Window(), für bessere Skalierung der Grafik
windowtsDataWeekly <- window(Trainingsdaten_80, start= 2016, end= c(2016, 365))
tsDataWeekly <- ts(windowtsDataWeekly, frequency = 7)
#für kürzere Prognose h-Wert anpassen 3 Monate: h=90
fit1 <- hw(tsDataWeekly, h=365, seasonal = "additive")
fit2 <- hw(tsDataWeekly, h=365, seasonal = "multiplicative")
autoplot(tsDataWeekly) +
  autolayer(fit1, series="HW additive Prognosen", PI=FALSE) +
  autolayer(fit2, series="HW multiplicative Prognosen", PI=FALSE) +
  xlab("Woche") +
  ylab("durchschn. Verkäufe") +
  ggtitle("HW Prognose") +
  guides(colour=guide_legend(title="Prognose"))

2.3.5 Arima-Modell

#für kürzere Prognose h-Wert anpassen 3 Monate: h=90
Prognose_arima1 <- auto.arima(Trainingsdaten_80)
autoplot(forecast(Prognose_arima1, h=365), xlab = "Jahr", ylab= "durchschn. Verkäufe")

print(Prognose_arima1)

## Series: Trainingsdaten_80 
## ARIMA(5,1,2)(0,1,0)[365] 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4     ar5      ma1      ma2
##       -0.9110  -0.1920  -0.2059  -0.2456  0.0753  -0.1768  -0.7867
## s.e.   0.0386   0.0418   0.0412   0.0416  0.0329   0.0248   0.0237
## 
## sigma^2 estimated as 512.9:  log likelihood=-4967.91
## AIC=9951.83   AICc=9951.96   BIC=9991.81

accuracy(forecast(Prognose_arima1, h=365), Testdaten_20)

##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  0.04294924 19.54138 12.90186 -2.031587 12.21522 0.5962997
## Test set     -4.49933055 25.95744 18.64379 -6.087430 16.22927 0.8616809
##                     ACF1 Theil's U
## Training set  0.02394253        NA
## Test set     -0.24256707 0.9260199

2.3.6 Prophet-Modell

#Erstellen Zeitreihe item15 store2
train_subset_store2_item15=subset(train,train$store==2 & train$item==15)
#Erstellung window mit Zeitraum 2013-2016
setDT(train_subset_store2_item15)
store_2_item15_fenster <- train_subset_store2_item15[date <= "2016-12-31"]

#Durchführung von Prophet-Modell 1Jahresprognose 2017-2018
stats=data.frame(y=store_2_item15_fenster$sales
                 ,ds=store_2_item15_fenster$date)
stats=aggregate(stats$y,by=list(stats$ds),FUN=sum)
#head(stats)
colnames(stats)<- c("ds","y")
#für kürzere Prognose periods anpassen 3 Monate: periods = 90
model_prophet = prophet(stats)

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

#summary(model_prophet)
future = make_future_dataframe(model_prophet, periods = 90)
forecast = predict(model_prophet, future)
plot(model_prophet, forecast, xlab="Jahr", ylab="durchschn. Verkäufe")

#mithilfe von Excel wurdeen die Prognosewerte extrahiert
#Laden und umbenennen der Excel
#Prophet_prognosewerte <- R_prophet_modell_forecast
#prophet_train_set<- ts(Prophet_prognosewerte, frequency = 365, start=c(2017,1))
#accuracy für 1 Jahr 
#accuracy(prophet_train_set, Testdaten_20)
#accuracy 3 Monate
#Prophet_prognosewerte_3monate <- prophet_3monate_prognose
#prophet_train_set_3monate <- ts(Prophet_prognosewerte_3monate, frequency = 365,start=c(2017,1))
#accuracy(prophet_train_set_3monate, Testdaten_20)

2.3.7 Neuronale Netze

#für kürzere Prognose h-Wert anpassen 3 Monate: h=90
neuronales_netz1 <- nnetar(Trainingsdaten_80)
autoplot(forecast(neuronales_netz1, h=365), xlab = "Jahr", ylab = "durchschn. Verkäufe")

accuracy(forecast(neuronales_netz1, h=365), Testdaten_20)

##                        ME      RMSE       MAE        MPE      MAPE      MASE
## Training set -0.006384656  4.531664  3.454613 -0.4259944  3.213893 0.1596658
## Test set     -9.345518952 18.929416 15.211271 -9.4582690 13.251083 0.7030365
##                    ACF1 Theil's U
## Training set 0.03008802        NA
## Test set     0.44101031 0.8335747