Description

This reports provides material forecasting demand using Arima, Exponential Smoothing, Croston, and Syntetos Boyland Approximation method. The dataset using in this report for modeling is real case data in one of the state owned company in Indonesia.

The report as structured as follows:
1. Data Extraction
2. Exploratory Data Analysis
3. Data Preparation
4. Modeling
5. Evaluation
6. Conclusion

1. Data Extraction

Import necessary libraries.

library(readxl)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.2     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.2
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(fpp3)
## -- Attaching packages -------------------------------------------- fpp3 0.4.0 --
## v lubridate   1.7.10     v feasts      0.2.2 
## v tsibble     1.0.1      v fable       0.3.1 
## v tsibbledata 0.3.0
## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date()    masks base::date()
## x dplyr::filter()      masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval()  masks lubridate::interval()
## x dplyr::lag()         masks stats::lag()
## x tsibble::setdiff()   masks base::setdiff()
## x tsibble::union()     masks base::union()
library(dplyr)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(flexclust)
## Loading required package: grid
## Loading required package: lattice
## Loading required package: modeltools
## Loading required package: stats4
## 
## Attaching package: 'modeltools'
## The following object is masked from 'package:fabletools':
## 
##     refit

Library readlxl: for import excel files into R.
Library tidyverse: for select unique/distinct rows from data frame.
Library fpp3: for visualization of correlation coefficient matrix
Library gridExtra: for plotting multiple graphs. Library caret: for One Hod Encoding

For the next step, we can load the data and selecting relevant columns. Then, see the dataframe’s structure.

dt2 <- read_excel("Data/Data Kerja Praktik.XLSX", sheet = "FINAL", skip = 1) %>%
  select("material" = `Short text`, "harga" = Harga, "tingkat_kekritisan" = `Tingkat kekritisan`, "rutin_nonrutin" = `Rutin/Non Rutin`, Jan:Des, `2014`:`2018`, `2019` = `2019...30`, "ROQ" = ROQ)
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in C1881 / R1881C3: got '(blank)'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in C1882 / R1882C3: got 'Grand Total'
## New names:
## * `2019` -> `2019...24`
## * `2019` -> `2019...30`
dt2 <- dt2[2:(nrow(dt2)-2),]

dt <- dt2 %>%
  group_by(material) %>% 
  distinct(material, .keep_all = T)

str(dt)
## grouped_df [1,861 x 23] (S3: grouped_df/tbl_df/tbl/data.frame)
##  $ material          : chr [1:1861] "CABLE PWR;NYAF;1X4mm2;0.6/1kV;OH" "INVERTER;;380/500V;3P;0.75KW" "INVERTER;;380/480VAC;6P;15KW" "INVERTER;;380/480VAC;6P;5.5KW" ...
##  $ harga             : num [1:1861] 380000 3153700 47000000 22000000 34500000 ...
##  $ tingkat_kekritisan: chr [1:1861] "C" "C" "C" "C" ...
##  $ rutin_nonrutin    : chr [1:1861] "NON RUTIN" "NON RUTIN" "NON RUTIN" "NON RUTIN" ...
##  $ Jan               : num [1:1861] 0 0 0 0 1 0 0 0 0 0 ...
##  $ Peb               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Mar               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Apr               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ May               : num [1:1861] 0 0 0 0 7 0 0 0 0 0 ...
##  $ Jun               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Jul               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Aug               : num [1:1861] 0 0 0 0 3 0 0 0 0 0 ...
##  $ Sep               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Okt               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Nop               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Des               : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ 2014              : num [1:1861] 0 4 0 0 0 0 0 0 0 0 ...
##  $ 2015              : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ 2016              : num [1:1861] 3 0 2 2 0 0 0 0 3 0 ...
##  $ 2017              : num [1:1861] 0 0 0 0 0 0 0 0 0 3 ...
##  $ 2018              : num [1:1861] 0 0 0 0 0 0 0 0 0 0 ...
##  $ 2019              : num [1:1861] 0 0 0 0 1 0 0 0 0 0 ...
##  $ ROQ               : num [1:1861] 0 0 0 0 6.08 ...
##  - attr(*, "groups")= tibble [1,861 x 2] (S3: tbl_df/tbl/data.frame)
##   ..$ material: chr [1:1861] "ACTUATOR ACC;POSITIONER SIPART PS2" "ACTUATOR;AUM;SA075AM011;GS1003" "ACTUATOR;AUM;SA075AM011GF1003;" "ACTUATOR;CCI;FFWCV;85104503" ...
##   ..$ .rows   : list<int> [1:1861] 
##   .. ..$ : int 602
##   .. ..$ : int 237
##   .. ..$ : int 249
##   .. ..$ : int 149
##   .. ..$ : int 142
##   .. ..$ : int 136
##   .. ..$ : int 239
##   .. ..$ : int 947
##   .. ..$ : int 432
##   .. ..$ : int 358
##   .. ..$ : int 433
##   .. ..$ : int 434
##   .. ..$ : int 584
##   .. ..$ : int 565
##   .. ..$ : int 482
##   .. ..$ : int 483
##   .. ..$ : int 484
##   .. ..$ : int 485
##   .. ..$ : int 486
##   .. ..$ : int 487
##   .. ..$ : int 586
##   .. ..$ : int 572
##   .. ..$ : int 243
##   .. ..$ : int 244
##   .. ..$ : int 245
##   .. ..$ : int 246
##   .. ..$ : int 1613
##   .. ..$ : int 1758
##   .. ..$ : int 1757
##   .. ..$ : int 1756
##   .. ..$ : int 1750
##   .. ..$ : int 1754
##   .. ..$ : int 1462
##   .. ..$ : int 1250
##   .. ..$ : int 463
##   .. ..$ : int 464
##   .. ..$ : int 465
##   .. ..$ : int 466
##   .. ..$ : int 467
##   .. ..$ : int 468
##   .. ..$ : int 461
##   .. ..$ : int 456
##   .. ..$ : int 454
##   .. ..$ : int 455
##   .. ..$ : int 452
##   .. ..$ : int 453
##   .. ..$ : int 462
##   .. ..$ : int 98
##   .. ..$ : int 99
##   .. ..$ : int 100
##   .. ..$ : int 101
##   .. ..$ : int 102
##   .. ..$ : int 103
##   .. ..$ : int 104
##   .. ..$ : int 691
##   .. ..$ : int 442
##   .. ..$ : int 441
##   .. ..$ : int 435
##   .. ..$ : int 436
##   .. ..$ : int 437
##   .. ..$ : int 438
##   .. ..$ : int 447
##   .. ..$ : int 439
##   .. ..$ : int 440
##   .. ..$ : int 444
##   .. ..$ : int 448
##   .. ..$ : int 449
##   .. ..$ : int 450
##   .. ..$ : int 451
##   .. ..$ : int 457
##   .. ..$ : int 458
##   .. ..$ : int 459
##   .. ..$ : int 573
##   .. ..$ : int 596
##   .. ..$ : int 597
##   .. ..$ : int 590
##   .. ..$ : int 551
##   .. ..$ : int 591
##   .. ..$ : int 592
##   .. ..$ : int 593
##   .. ..$ : int 598
##   .. ..$ : int 709
##   .. ..$ : int 552
##   .. ..$ : int 548
##   .. ..$ : int 556
##   .. ..$ : int 555
##   .. ..$ : int 579
##   .. ..$ : int 566
##   .. ..$ : int 576
##   .. ..$ : int 577
##   .. ..$ : int 488
##   .. ..$ : int 489
##   .. ..$ : int 490
##   .. ..$ : int 491
##   .. ..$ : int 492
##   .. ..$ : int 493
##   .. ..$ : int 543
##   .. ..$ : int 494
##   .. ..$ : int 547
##   .. .. [list output truncated]
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE

The dataset has 1861 observations and 23 variables.

2. Exploratory Data Analysis

After in the previous stage we can find out the conclusions from each column, then we can further analyze it using exploratory data analysis (EDA) techniques. The purpose of this EDA stage is to be able to find out information more clearly and efficiently through graph plots of the data to be predicted.

dt2$tingkat_kekritisan <- as.factor(dt2$tingkat_kekritisan)
dt2$rutin_nonrutin <- as.factor(dt2$rutin_nonrutin)
dt2$tingkat_kekritisan[dt2$tingkat_kekritisan == "c"]
##  [1] <NA> <NA> c    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [16] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [31] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: A B c C
dt2 <- dt2 %>%
  mutate(tingkat_kekritisan =
           replace(tingkat_kekritisan,
                   tingkat_kekritisan == "c",
                   "C"))
p1a <- ggplot(data=dt2, aes(x = tingkat_kekritisan)) +
  geom_bar(fill = rainbow(4),
           color = "blue")+
  stat_count(geom = "text", color = "black", size = 5,
             aes(label = ..count..), position=position_stack(vjust = 0.5)) +
  labs(title = "Plot Tingkat kekritisan")
p2a <- ggplot(data=dt2, aes(x = rutin_nonrutin)) +
  geom_bar(fill = rainbow(2),
           color = "blue")+
  stat_count(geom = "text", color = "black", size = 5,
             aes(label = ..count..), position=position_stack(vjust = 0.5)) +
  labs(title = "Plot Sifat Material")

plot_sifat <- grid.arrange(p1a, p2a, ncol = 2, 
                           top=textGrob("Plot Material",gp=gpar(fontsize=30,font=3)))

Plot visual data on the classification and properties of materials based on material list data. Seen in the classification of materials, materials with a criticality level of C have a greater number, of which there are 1552 materials and based on their nature, non-routine materials have more quantities than non-routine ones. There are 1789 non-routine materials and 88 routine materials. In 2020, there are 1737 non-routine materials and 88 routine materials.

3. Data Preparation

After going through the process of cleaning and merging data, the next preprocessing stage is to transform the data. This process is used to convert the previous data format into the desired format in the model.

yearly <- dt %>%
  select(material, `2014`:`2019`) %>%
  pivot_longer(`2014`:`2019`, names_to = "tahun", values_to = "n") %>%
  mutate(tahun = as.integer(tahun)) %>%
  arrange(tahun, material) %>%
  distinct(tahun, material, .keep_all = T) %>%
  as_tsibble(index = tahun, key = material)

The previous data dimensions totaling 1862 observations and 23 variables divided into 2 data, namely training data and test data. Columns in 2014-2018 were changed to 1 column thereby increasing the total number of observations. The variable n is the value of the rate of use in each year which later this value will be included in the forecasting test. Based on these results, it was found that there were 11166 observations and 3 data variables. The variable consists of the name of the material, the year, and n (the value of the rate of use).

4. Modelling

The modeling stage will discuss the model that will be used for forecasting. The transformed material data will be divided first into training data and test data. Training data is data that we build to train the algorithm to make predictions that later the program can read patterns based on the data provided. Test data is data given to test the results of predicted values to measure the accuracy of values based on actual data.

4.1 Generate Test Design

Training data is the usage rate taken in the 2014-2018 period. As for the test data, the total usage rate will be taken in 2019. After the data is distributed, the total training data consists of 9305 observations and the test data has 1861 observation data.

# Split training = 6, test = 1 ----
train_ts <- yearly %>%
  group_by_key(material) %>%
  slice(1:5)
test_ts <- yearly %>%
  group_by_key(material) %>%
  slice(6)

train_ts
## # A tsibble: 9,305 x 3 [1Y]
## # Key:       material [1,861]
## # Groups:    material [1,861]
##    material                           tahun     n
##    <chr>                              <int> <dbl>
##  1 ACTUATOR ACC;POSITIONER SIPART PS2  2014    12
##  2 ACTUATOR ACC;POSITIONER SIPART PS2  2015     1
##  3 ACTUATOR ACC;POSITIONER SIPART PS2  2016     0
##  4 ACTUATOR ACC;POSITIONER SIPART PS2  2017     1
##  5 ACTUATOR ACC;POSITIONER SIPART PS2  2018     0
##  6 ACTUATOR;AUM;SA075AM011;GS1003      2014     0
##  7 ACTUATOR;AUM;SA075AM011;GS1003      2015     0
##  8 ACTUATOR;AUM;SA075AM011;GS1003      2016     2
##  9 ACTUATOR;AUM;SA075AM011;GS1003      2017     0
## 10 ACTUATOR;AUM;SA075AM011;GS1003      2018     0
## # ... with 9,295 more rows
test_ts
## # A tsibble: 1,861 x 3 [1Y]
## # Key:       material [1,861]
## # Groups:    material [1,861]
##    material                           tahun     n
##    <chr>                              <int> <dbl>
##  1 ACTUATOR ACC;POSITIONER SIPART PS2  2019     7
##  2 ACTUATOR;AUM;SA075AM011;GS1003      2019     0
##  3 ACTUATOR;AUM;SA075AM011GF1003;      2019     0
##  4 ACTUATOR;CCI;FFWCV;85104503         2019     0
##  5 ACTUATOR;CCI;MFWCV;30202002SOFT     2019     0
##  6 ACTUATOR;CCI;TBPVALE;500RMSOFT      2019     0
##  7 ACTUATOR;CDC;63000604D3;PSLG15      2019     0
##  8 ACTUATOR;EOM;63000604D3;ED0100      2019    12
##  9 ACTUATOR;MOTOR;;213050413           2019     0
## 10 ACTUATOR;MOTOR;A0000012A;48RPM      2019     0
## # ... with 1,851 more rows

The transformation of the data entered into the dataframe (yearly) is then divided to determine the training data and test data used for testing. In the training data there are 5 points of request time, while the test data has 1 point of request time.

4.2 Build Model

The methods used in this forecasting are Croston method, Syntetos Boyland Approximation, Arima and Exponential Smoothing. These four models will predict the value of training which will automatically learn the pattern of data movement for each material.

# Fitting Model ----
model_fit <- train_ts %>%
  model(
    croston = CROSTON(n, type = "croston"),
    sba = CROSTON(n, type = "sba"),
    arima = ARIMA(n),
    ets = ETS(n)
  )
## Warning: 1315 errors (1 unique) encountered for croston
## [1315] At least two non-zero values are required to use Croston's method.
## Warning: 1315 errors (1 unique) encountered for sba
## [1315] At least two non-zero values are required to use Croston's method.
model_fit
## # A mable: 1,861 x 5
## # Key:     material [1,861]
##    material                     croston          sba          arima          ets
##    <chr>                        <model>      <model>        <model>      <model>
##  1 ACTUATOR ACC;POSITIONE~    <croston>    <croston> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  2 ACTUATOR;AUM;SA075AM01~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  3 ACTUATOR;AUM;SA075AM01~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  4 ACTUATOR;CCI;FFWCV;851~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  5 ACTUATOR;CCI;MFWCV;302~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  6 ACTUATOR;CCI;TBPVALE;5~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  7 ACTUATOR;CDC;63000604D~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  8 ACTUATOR;EOM;63000604D~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
##  9 ACTUATOR;MOTOR;;213050~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
## 10 ACTUATOR;MOTOR;A000001~ <NULL model> <NULL model> <ARIMA(0,0,0)> <ETS(A,N,N)>
## # ... with 1,851 more rows

It is known that the arima model produces ARIMA(0,0,0) which means that the data has a random walk nature where the variable value generated from time to time does not show a predictive pattern at all. Exponential smoothing produces an ETS(A,N,N) model which means that the data does not have a trend and seasonal level in the data. The model built on Croston and SBA, there are some unpredictable materials that result in (NULL model). This happens because the Croston and SBA methods cannot predict when the use of materials in the 2014-2018 time frame is used for less than 2 years, so the material that cannot be predicted will be defined as (NULL model).

# Filter material dengan nilai null pada semua model ----
not_null_model <- model_fit %>%
  filter(!is_null_model(sba) | !is_null_model(croston) | 
           !is_null_model(arima) | is_null_model(ets))

# Menghitung Akurasi ----
akurasi <- forecast(not_null_model, h = 2) %>%
  accuracy(test_ts) %>%
  select(.model, material, RMSE)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
## 1 observation is missing at 2020
akurasi_summary <- akurasi %>%
  na.omit() %>%
  group_by(.model) %>%
  summarise(RMSE_mean = mean(RMSE))

dt_prep <- dt %>% 
  select(1, "harga" = 2, "tingkat_kekritisan" = 3, "rutin_nonrutin" = 4, "2019" = 22, "ROQ" = 23)

It can be concluded that the method produced by ARIMA, Croston, Exponential Smoothing, and SBA has a greater error value than the method used by the company. The method produced by arima has an average error of 10.6 and Exponential Smoothing produces an average error of 22.5, while the average error value for the croston and SBA methods cannot be generated because there are some materials that cannot produce a value. prediction. The method produced by the company itself has a smaller error of 2.07. In this study, the forecasting results from the ARIMA, Croston, Exponential Smoothing, SBA and existing methods will be combined by choosing the method that can produce the smallest error value.

5. Evaluation

This evaluation consists of an analysis of the calculation of the error value, the value of inventory level, service level and warehouse costs based on the existing policy method, analysis method Arima, Exponential Smoothing, Croston, and Syntetos Boyland Approximation.

# Calculate forecast value ----
peramalan <- forecast(not_null_model, h = 2)

# Tabel report ----
t1 <- peramalan %>%
  as_tibble() %>%
  filter(tahun == 2020) %>%
  mutate(prediksi_2020 = paste0(.model, "_", as.character(tahun))) %>%
  select(material, .mean, prediksi_2020) %>%
  pivot_wider(names_from = prediksi_2020, values_from = .mean)

t2 <- akurasi %>%
  select(.model, material, RMSE) %>%
  pivot_longer(RMSE, names_to = "matriks", values_to = "nilai_akurasi") %>%
  mutate(model_akurasi = paste0(.model, "_", matriks)) %>%
  select(material, nilai_akurasi, model_akurasi) %>%
  pivot_wider(names_from = model_akurasi, values_from = nilai_akurasi)

TO get the forecast value, we can enter the results of the model that has been made for each material. in table t1, will display the results of forecasting values for each method used. while for table t2 will display the error value resulting from each method.

Next we will assume that the model that returns NULL values is 0, because the model cannot generate forecasted values. After that we can also combine the four models used by choosing the smallest error value from each forecasting model and it will be entered into the optimal method table.

# Memasukkan 4 model dan memilih model yang terbaik
## Memasukkan data peramalan yang sudah dibersihkan di Excel
dt_prep <- read_excel("Data/dt_prep.xlsx")

akurasi2 <- dt_prep %>%
  select(material, ROQ.e) %>%
  na.omit() %>%
  pivot_longer(ROQ.e, names_to = ".model", values_to = "RMSE") %>%
  bind_rows(akurasi) %>%
  group_by(.model)

akurasi3 <- akurasi2 %>%
  group_by(material) %>%
  slice_min(order_by = RMSE, n = 1)

akurasi2 %>%
  group_by(material) %>%
  slice_min(order_by = RMSE, n = 1) %>%
  ungroup() %>%
  distinct(material, .keep_all = TRUE) %>%
  summarise(mean(RMSE))
## # A tibble: 1 x 1
##   `mean(RMSE)`
##          <dbl>
## 1         5.03
metode_optimal <- akurasi2 %>%
  group_by(material) %>%
  slice_min(order_by = RMSE, n = 1) %>%
  ungroup() %>%
  distinct(material, .keep_all = TRUE)

The report table contains the optimal model for forecasting materials. The optimal model is determined based on the smallest error value in each forecast tested based on demand data in 2019. After we combine these 5 methods, we get an average error value of 5.03 points. from the value obtained, it is proven that by combining the 5 models used, it can reduce the error value in the company’s material forecast which was previously worth 7,32.

nilai_3_peramalan <- dt_prep %>%
  select(material, ROQ) %>%
  left_join(t1)
## Joining, by = "material"
nilai_3_peramalan$ROQ <- as.numeric(nilai_3_peramalan$ROQ)
## Warning: NAs introduced by coercion
nilai_3_peramalan %>%
  pivot_longer(cols = 2:6, names_to = ".model", values_to = "nilai_prediksi") %>%
  right_join(metode_optimal, by = c("material", ".model")) %>%
  ungroup() %>%
  arrange(.model)
## # A tibble: 1,861 x 4
##    material                             .model nilai_prediksi  RMSE
##    <chr>                                <chr>           <dbl> <dbl>
##  1 ACTUATOR;AUM;SA075AM011GF1003;       arima              NA     0
##  2 ACTUATOR;CDC;63000604D3;PSLG15       arima              NA     0
##  3 ACTUATOR;MOTOR;A0000012A;48RPM       arima              NA     0
##  4 BOILER ACC;IVC ROD END LINKAGE       arima              NA     2
##  5 BOILER ACC;TUBE 64X6X7000MM SA210-A1 arima              NA     0
##  6 BOILER ACC;TUBE SUPERHEAT GB12       arima              NA    15
##  7 BOLT&NUT;FUJ;M64X;ST114051AP2P       arima              NA     0
##  8 BUSH;DMW;VPFCM;R233                  arima              NA     0
##  9 CHEM;ANTISCALANT RO;PC 191 T         arima              NA    90
## 10 CHEM;HYDROCHLORIC ACID;32%           arima              NA   960
## # ... with 1,851 more rows
result_forecast <- nilai_3_peramalan %>%
  pivot_longer(cols = 2:6, names_to = ".model", values_to = "nilai_prediksi") %>%
  mutate(
    .model = str_remove(.model, "_2020"), 
    .model = str_replace(.model, "ROQ", "ROQ.e")) %>%
  right_join(metode_optimal, by = c("material", ".model")) %>%
  ungroup() %>%
  arrange(.model)

colSums(is.na(result_forecast))
##       material         .model nilai_prediksi           RMSE 
##              0              0              0              0
nrow(nilai_3_peramalan)
## [1] 1861
nrow(metode_optimal)
## [1] 1861
colSums(is.na(result_forecast))
##       material         .model nilai_prediksi           RMSE 
##              0              0              0              0
result_forecast2 <- result_forecast %>%
  left_join(dt_prep) %>%
  select(material, harga, tingkat_kekritisan, ROQ, rutin_nonrutin, `2019`,
         starts_with(".model"), starts_with("nilai_prediksi"))
## Joining, by = "material"
result_forecast2 <- result_forecast2 %>%
  mutate(ROQ =
           replace(ROQ,
                   is.na(result_forecast2$ROQ),
                   "0"))
result_forecast2$ROQ <- as.numeric(result_forecast2$ROQ)
## Warning: NAs introduced by coercion
result_forecast2$harga <- as.numeric(result_forecast2$harga)
## Warning: NAs introduced by coercion
result_forecast2$nilai_prediksi <- ceiling(result_forecast2$nilai_prediksi)
result_forecast2$ROQ <- ceiling(result_forecast2$ROQ)

The results of the merger will be entered into the result forecast 2 column to find out the difference in terms of the inventory value that will be ordered by the company in the next period.

Inventory Level

Determining the value of the inventory level will affect the allocation of material storage and also the costs that will be invested in purchasing materials. The less the amount of material ordered, the smaller the cost allocation used. The company targets that the amount of material ordered can be as optimal as possible in order to reduce the value of waste in the warehouse.

forecast2_arima <- filter(result_forecast2, .model == "arima")
forecast2_croston <- filter(result_forecast2, .model == "croston")
forecast2_ets <- filter(result_forecast2, .model == "ets")
forecast2_sba <- filter(result_forecast2, .model == "sba")
forecast2_roq <- filter(result_forecast2, .model == "ROQ.e")

sum(result_forecast2$`2019`)
## [1] 10878
sum(result_forecast2$ROQ)
## [1] NA
sum(result_forecast2$nilai_prediksi)
## [1] 13346

The data showing the total material after rounding based on the company’s forecasting results and combining 5 methods. The number of materials predicted by the company totaled 17,195 units of material and the forecasting results from the combination of the four methods amounted to 13,346 units of material. In the combination of these 5 methods, there is a percentage decrease in the number of material orders by 22.4%. It also shows that the company can save the amount of material ordered by 3,849 units.

Material Cost

Determination of warehouse costs is obtained by multiplying the number of forecasted materials with the price of each material obtained based on data from suppliers. These results will be tested to determine whether the proposed method can minimize the amount of budget allocated to warehousing.

## Perhitungan Anggaran'

# Anggaran berdasarkan peramalan
anggaran <- result_forecast2 %>% 
  mutate(
    anggaran_by_model = nilai_prediksi * harga,
    anggaran_by_roq  = ROQ * harga,
    anggaran_by_2019 = `2019` * harga
    # anggaran_by_ROQ
  ) %>% 
  pivot_longer(
    starts_with("anggaran_by"), 
    names_to = "metode_ramalan", 
    values_to = "hasil_peramalan"
  ) %>% 
  group_by(metode_ramalan) %>%
  summarise(total_anggaran = sum(hasil_peramalan, na.rm = T))

Determination of warehouse costs is obtained by multiplying the number of forecasting materials with the price of each material obtained based on data from suppliers. It is known that the total budget generated by using a combination of 5 methods is Rp. 35,830,075,396.00 where the budget is smaller than the total budget generated using the existing method of Rp. 36,012,143,757.00. the total savings that can be made by a combination of 5 methods is Rp. 182,068,361.00.

Service Level

Assessment at the service level is measured to determine how effective the forecasting method used to predict the amount of material demand is. The higher the service level value obtained, the better for the company. At this stage, we will create a new column named service level to calculate the value obtained for each material and ensure that there is no NA value.

# SERVICE LEVEL
service_level <- read_excel("Data/Data-hasil2 (1).xlsx")
## New names:
## * `` -> ...1
service_level <- service_level %>%
  mutate(`service level` =
           replace(`service level`,
                   is.na(service_level$`service level`),
                   "1"))

service_level <- service_level %>%
  mutate(`service level eksisting` =
           replace(`service level eksisting`,
                   is.na(service_level$`service level eksisting`),
                   "1"))

service_level$`service level` <- as.numeric(service_level$`service level`)
service_level$`service level eksisting` <- as.numeric(service_level$`service level eksisting`)

# Price for "NA" ==> 0
service_level$harga[service_level$harga == "NA"]
##  [1] "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA"
## [16] "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA"
## [31] "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA" "NA"
## [46] "NA" "NA" "NA" "NA"
service_level <- service_level %>%
  mutate(harga =
           replace(harga,
                   harga == "NA",
                   0))
# Criticallity for "NA" ==> 0
service_level$tingkat_kekritisan[service_level$tingkat_kekritisan == "0"]
## character(0)
service_level <- service_level %>%
  mutate(tingkat_kekritisan =
           replace(tingkat_kekritisan,
                   tingkat_kekritisan == "0",
                   "Unavailable"))

colSums(is.na(service_level))
##                    ...1                material                   harga 
##                       0                       0                       0 
##      tingkat_kekritisan                     ROQ          rutin_nonrutin 
##                       0                       0                       0 
##                    2019                  .model          nilai_prediksi 
##                       0                       0                       0 
##           service level service level eksisting 
##                       0                       0

After make sure that there is no NA value in it, then we can calculate the service level value for each material. we will separate the condition of the material that is overstocked and not so that we can find out whether this new method can reduce the value of overstock or not. service level values that are overstocked will be entered into the sl_over table and those that are not included in sl_target.

sl_target <- filter(service_level,
                    `service level` <= 1)
sl_target_eks <- filter(service_level,
                    `service level eksisting` <= 1)


mean(sl_target$`service level`)
## [1] 0.9533147
mean(sl_target_eks$`service level eksisting`)
## [1] 0.9506557
service_level <- arrange(service_level, desc(`service level`))

in the result, it can be seen that the combination of 5 methods can produce a service level value of 95.33%. Meanwhile, the company’s existing model produces a smaller value, namely 95.06%.

## Service Level Overstock  
sl_over <- filter(service_level,
                 `service level` > 1)
sl_over_eks <- filter(service_level,
                  `service level eksisting` > 1)
mean(sl_over$`service level`)
## [1] 1.72423
mean(sl_over_eks$`service level eksisting`)
## [1] 2.469007

The results show the average material that has overstock through service level values that have a value of more than 100%. The combination of 5 methods has a smaller average of 172.4% while the company’s existing model has a larger overstock value of 246.9%. it can be concluded that there is a 75% decrease in the value of overstocked materials.

6. Conclusion

  1. The combination of the 5 methods used can produce a more accurate forecast value for the company. a combination of 5 methods can minimize the percentage error value of 31.3%.
  2. The existing method produces an inventory level value of 17,195 material, service level value is 95.06%, and warehouse cost is Rp.36,012,143,757.00. The combination of 5 methods can produce a value of inventory level is 13,346 materials, service level value is 95.33%, and warehouse costs are Rp.35,830,075,396.00. It is known that the combination of 5 methods can reduce the value of inventory level by 22.4%, increase the value of service level by 0.27%, and decrease warehouse costs by 0.51%.