A. Library, Data, dan Langkah Kerja

pacman::p_load(tseries, forecast, TSA, xts, dygraphs, tidyverse, lubridate, ggplot2, dplyr, formattable, DiagrammeR, psd, readr, knitr, kableExtra, purrr, webshot, FitAR, TSstudio, glue, plotly, ggplot2, ggbreak) # Package loudout for Timeseries

## package 'glue' successfully unpacked and MD5 sums checked
## package 'DiagrammeR' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Client\AppData\Local\Temp\RtmpkXzamm\downloaded_packages

library(hrbrthemes)
library(viridis)
library(gapminder)

# data yang digunakan
data <- read_csv("rawdata.csv", col_types = cols(date = col_date(format = "%m/%d/%Y")))

data %>% kbl(booktabs = T) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))

Untuk melakukan proses peramalan library yang dibutuhkan dan data yang dijadikan sebagai objek peramalan sebagai berikut. Kolom month menunjukkan bulan.

C. Pembagian Data

training <- data.frame(data[1:42,])
testing <- data.frame(data[43:52,])

ts_training <- ts(training$qty, freq = 365.25/7, start = c(2021,1))

D. Stasioneritas

# adf test
adf.test(ts_training)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_training
## Dickey-Fuller = -0.93136, Lag order = 3, p-value = 0.9368
## alternative hypothesis: stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -3.8007, Lag order = 3, p-value = 0.02963
## alternative hypothesis: stationary

E. Pembentukan Model

Secara keseluruhan akan coba dibentuk empat model peramalan. Keempat model yang dibangun tersebut adalah: (1) Arima, (2) Error, trend, seasonal model-ETS, (3) TBATS, dan (4) Exponential Smoothing. Dari hasi eksploratori data sebelumnya memang diketahui bahwa selain datanya tidak terlalu banyak, juga cenderung tidak memiliki efek seasonal, namun model (2) dan (3) tetap digunakan sebagai opsi. Keempat model tersebut kemudian akan dievaluasi dan dibandingkan keakurasiannya dalam meramal. Model dengan tingkat akurasi paling baik akan dipilih sebagai model terbaik dan digunakan untuk meramalkan nilai ke depan.

Untuk mengevaluasi beberapa model, dibentuk terlebih dahulu fungsi analisis sisaan yang bertujuan untuk melihat apakah sisaan yang dihasilkan model signifikan atau tidak.

analisis_sisaan <- function(model, std = TRUE){
  if (std == TRUE){
    res.model = rstandard(model)
  }
  else
    {
    res.model = residuals(model)
  }
  
  par(mfrow=c(2,2))
  plot(res.model,type='o',ylab='Standardised residuals',
       main="Time series plot of standardised residuals")
  abline(h=0)
  qqnorm(res.model,main="QQ plot of standardised residuals")
  qqline(res.model, col = 2)
  acf(res.model,main="ACF of standardised residuals")
  k=0
  LBQPlot(res.model, lag.max = length(model$residuals)-1 ,
          StartLag = k + 1, k = 0, SquaredQ = FALSE)
  par(mfrow=c(1,1))
}

# sources from
# https://www.titanwolf.org/Network/q/9119c26d-885e-4722-9037-3322eac273da/y

# model arima
ts_training %>% auto.arima(d = 1, stationary = F) 

## Series: . 
## ARIMA(1,1,2) 
## 
## Coefficients:
##           ar1     ma1      ma2
##       -0.9787  0.2116  -0.6226
## s.e.   0.0383  0.1467   0.1492
## 
## sigma^2 estimated as 17884:  log likelihood=-257.98
## AIC=523.97   AICc=525.08   BIC=530.82

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.978724   0.038295 -25.5577 < 2.2e-16 ***
## ma1  0.211581   0.146703   1.4422    0.1492    
## ma2 -0.622621   0.149207  -4.1729 3.008e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,2)
## Q* = 5.5067, df = 5, p-value = 0.3572
## 
## Model df: 3.   Total lags used: 8

options(warn = -2)

ts_training %>% Arima(order=c(2,1,2))

## Series: . 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1     ar2     ma1      ma2
##       -0.8793  0.0965  0.1517  -0.6578
## s.e.   0.2499  0.2399  0.1984   0.1610
## 
## sigma^2 estimated as 18272:  log likelihood=-257.9
## AIC=525.81   AICc=527.52   BIC=534.37

ts_training %>% Arima(order=c(1,1,3))

## Series: . 
## ARIMA(1,1,3) 
## 
## Coefficients:
##           ar1     ma1      ma2      ma3
##       -0.9775  0.2974  -0.6518  -0.1307
## s.e.   0.0360  0.2086   0.1480   0.2324
## 
## sigma^2 estimated as 18155:  log likelihood=-257.82
## AIC=525.63   AICc=527.35   BIC=534.2

ts_training %>% Arima(order=c(2,1,2)) %>% lmtest::coeftest()

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.87930    0.24990 -3.5187 0.0004337 ***
## ar2  0.09651    0.23992  0.4023 0.6874912    
## ma1  0.15167    0.19843  0.7643 0.4446663    
## ma2 -0.65781    0.16096 -4.0869 4.372e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ts_training %>% Arima(order=c(1,1,3)) %>% lmtest::coeftest()

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.977479   0.036038 -27.1235 < 2.2e-16 ***
## ma1  0.297424   0.208585   1.4259    0.1539    
## ma2 -0.651813   0.148023  -4.4035 1.065e-05 ***
## ma3 -0.130719   0.232395  -0.5625    0.5738    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ts_training %>% Arima(order=c(2,1,2)) %>% analisis_sisaan()

ts_training %>% Arima(order=c(1,1,3)) %>% analisis_sisaan()

ts_training %>% Arima(order=c(2,1,2)) %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)
## Q* = 6.0517, df = 4, p-value = 0.1953
## 
## Model df: 4.   Total lags used: 8

ts_training %>% Arima(order=c(1,1,3)) %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,3)
## Q* = 6.7122, df = 4, p-value = 0.1519
## 
## Model df: 4.   Total lags used: 8

model_arima <- ts_training %>% auto.arima(d = 1, stationary = F)
model_arima$aic

## [1] 523.9669

# ets
model_ets <- ts_training %>% ets()
model_ets$aic

## [1] 575.2265

# tbats
model_tbats <- ts_training %>% tbats(use.trend = F, use.damped.trend = F, use.box.cox = F, use.arma.errors = T)
model_tbats$AIC

## [1] 546.071

# Holt-Winters
model_exsmot <- ts_training %>% HoltWinters(gamma = F, beta = 0.1, alpha = 0.3)
model_exsmot

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = ., alpha = 0.3, beta = 0.1, gamma = F)
## 
## Smoothing parameters:
##  alpha: 0.3
##  beta : 0.1
##  gamma: FALSE
## 
## Coefficients:
##        [,1]
## a 427.88254
## b  19.98501

F. Peramal (Plot) dan Evaluasi Model

# model arima

f_arima <- forecast(model_arima, h = 10) # h = testing data
plot_forecast(f_arima, title = "Model Arima")

# saveRDS(f_arima, file = "SFK.rds")

# model ets
f_ets <- forecast(model_ets, h = 10)
plot_forecast(f_ets, title = "Model ETS")

# model tbats
f_tbats <- forecast(model_tbats, h = 10)
plot_forecast(f_tbats, title = "Model TBATS")

# model exponential smoothing

f_exsmot <- forecast(model_exsmot, h = 10)
plot_forecast(f_exsmot, title = "Model Exponential Smoothing")

df_arima <- f_arima$mean %>% as.data.frame() %>%
  mutate(across(where(is.numeric), ~ round(., 0)))
colnames(df_arima) <- "arima"

df_ets <- f_ets$mean %>% as.data.frame() %>%
  mutate(across(where(is.numeric), ~ round(., 0)))
colnames(df_ets) <- "ets"

df_tbats <- f_tbats$mean %>% as.data.frame() %>%
  mutate(across(where(is.numeric), ~ round(., 0)))
colnames(df_tbats) <- "tbats"

df_exsmot <- f_exsmot$mean %>% as.data.frame() %>%
  mutate(across(where(is.numeric), ~ round(., 0)))
colnames(df_exsmot) <- "holtwinters"

kompilasi <- cbind(testing, df_arima, df_ets, df_exsmot, df_tbats)

kompilasi %>% select(-3)  %>% kbl(booktabs = T) %>%
  kable_styling(latex_options = c("striped", "HOLD_position"))

ak_arima <- accuracy(kompilasi$qty, kompilasi$arima, d = 1, D = 0)
ak_ets <- accuracy(kompilasi$qty, kompilasi$ets, d = 1, D = 0)
ak_holtwinters <- accuracy(kompilasi$qty, kompilasi$holtwinters, d = 1, D = 0)
ak_tbats <- accuracy(kompilasi$qty, kompilasi$tbats, d = 1, D = 0)

ak_compile <- rbind(ak_arima, ak_ets, ak_holtwinters, ak_tbats)
rownames(ak_compile) <- c("arima","ets", "holtwinters", "tbats")

ak_compile %>% as.data.frame() %>%  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  kbl(booktabs = T) %>% kable_styling(latex_options = c("striped", "HOLD_position"))

G. Peramalan

Misalnya kita ingin melakukan peramalan selama periode 3 bulan ke depan (tahun 2022) atau 12 pekan. Hal tersebut dapat dilakukan sebagai berikut:

tanggal <- c("21-10-25", "21-11-01", "21-11-08", "21-11-15", "21-11-22", "21-11-29", "21-12-06", "21-12-13", "21-12-20", "21-12-27", "22-01-03", "22-01-10", "22-01-17", "22-01-24", "22-01-31", "22-02-07", "22-02-14", "22-02-21", "22-02-28", "22-03-07", "22-03-14", "22-03-21")

ymd(tanggal)

##  [1] "2021-10-25" "2021-11-01" "2021-11-08" "2021-11-15" "2021-11-22"
##  [6] "2021-11-29" "2021-12-06" "2021-12-13" "2021-12-20" "2021-12-27"
## [11] "2022-01-03" "2022-01-10" "2022-01-17" "2022-01-24" "2022-01-31"
## [16] "2022-02-07" "2022-02-14" "2022-02-21" "2022-02-28" "2022-03-07"
## [21] "2022-03-14" "2022-03-21"

df_tanggal <- tanggal %>% as.data.frame.Date() 
df_tanggal <- df_tanggal[-c(1:10),] %>% as.data.frame.Date() 
colnames(df_tanggal) <- "date_yymmdd"

data_forecast <- f_arima_2022$mean %>% as.data.frame()
data_forecast <- data_forecast[-c(1:10),] %>% as.data.frame()
colnames(data_forecast) <- "forecast"

Maka peramalan untuk periode 12 pekan ke depan adalah:

data_forecast_2022 <- cbind(df_tanggal, data_forecast)

data_forecast_2022 %>% mutate(across(where(is.numeric), ~ round(., 1))) %>%
  kbl(booktabs = T) %>% kable_styling(latex_options = c("striped", "HOLD_position"))

library(readxl)
ul <- read_excel("detail_jualan.xls", skip = 7)
head(ul)

## # A tibble: 6 x 25
##      NO FAKTUR       TANGGAL             KODE_CUS  CUSTOMER  DIVISI SALES GUDANG
##   <dbl> <chr>        <dttm>              <chr>     <chr>     <chr>  <chr> <chr> 
## 1     1 0001/F/01/21 2021-01-04 00:00:00 "\"10006~ SINAR TI~ D-RMF  SAL-M GUD R~
## 2     2 0001/F/02/21 2021-02-01 00:00:00 "\"10000~ HARYANTO~ D-RMF  SAL-Y GUD R~
## 3     3 0001/F/03/21 2021-03-02 00:00:00 "\"10024~ AGUNG SE~ D-RMF  SAL-Y GUD R~
## 4     4 0001/F/04/21 2021-04-01 00:00:00 "\"10003~ CITRA MA~ D-RMF  SAL-Y GUD R~
## 5     5 0001/F/05/21 2021-05-03 00:00:00 "\"17007~ PHANG EN~ D-RMF  SAL-6 GUD R~
## 6     6 0001/F/06/21 2021-06-01 00:00:00 "\"18000~ WHITE HO~ D-RMF  SAL-6 GUD R~
## # ... with 17 more variables: JTH_TEMPO <dttm>, KODE_BRG <chr>, NAMA_BRG <chr>,
## #   BANYAK <dbl>, HARGA <dbl>, KURS <dbl>, DISC <dbl>, DISCOUNT <dbl>,
## #   TOTAL <dbl>, TOTAL FAKTUR <dbl>, DPP <dbl>, PPN <dbl>, FAKGDG <lgl>,
## #   KETERANGAN <chr>, KET ITEM <lgl>, PO <lgl>, TGLPO <dttm>

ul_clean_discard <-
  ul %>% select(
    -c(
      NO,
      FAKTUR,
      KODE_CUS,
      DIVISI,
      GUDANG,
      KODE_BRG,
      DISCOUNT,
      FAKGDG,
      KURS,
      PO,
      KETERANGAN,
      KODE_BRG,
      TGLPO,
      JTH_TEMPO,
      DISC,
      `KET ITEM`
    )
  )  

ul_clean <-
  ul_clean_discard %>% mutate_at(c("SALES","NAMA_BRG", "CUSTOMER"), as.factor) %>%
  filter_if(is.numeric, all_vars((.) != 0))
head(ul_clean)

## # A tibble: 6 x 10
##   TANGGAL             CUSTOMER        SALES  NAMA_BRG       BANYAK  HARGA  TOTAL
##   <dttm>              <fct>           <fct>  <fct>           <dbl>  <dbl>  <dbl>
## 1 2021-01-04 00:00:00 SINAR TIMUR - ~ SAL-M  GRACE F 55050~      1 168300 1.14e5
## 2 2021-02-01 00:00:00 HARYANTO - JAK~ SAL-Y  R.GRANIT GT91~      1 490050 3.03e5
## 3 2021-03-02 00:00:00 AGUNG SENTRA K~ SAL-Y  R.GRANIT GT61~      1 314050 1.97e5
## 4 2021-04-01 00:00:00 CITRA MANDIRI   SAL-Y  INTERLOK GL63~     22 188100 2.82e6
## 5 2021-06-01 00:00:00 WHITE HOUSE PD~ SAL-6  R.GRESS F 227~      2 133650 1.85e5
## 6 2021-07-02 00:00:00 MAKMUR KERAMIK~ SAL-B6 R.GRANIT GT63~      1 348150 2.14e5
## # ... with 3 more variables: TOTAL FAKTUR <dbl>, DPP <dbl>, PPN <dbl>

str(ul_clean)

## tibble [44,642 x 10] (S3: tbl_df/tbl/data.frame)
##  $ TANGGAL     : POSIXct[1:44642], format: "2021-01-04" "2021-02-01" ...
##  $ CUSTOMER    : Factor w/ 684 levels "ACC - CIBUBUR*",..: 591 187 6 98 676 288 126 207 24 437 ...
##  $ SALES       : Factor w/ 19 levels "SAL-2","SAL-20",..: 18 19 19 19 6 10 11 11 4 11 ...
##  $ NAMA_BRG    : Factor w/ 1332 levels "ASIA CLIFFORD GREY 25x50 KW 1",..: 25 656 509 209 894 585 403 192 327 23 ...
##  $ BANYAK      : num [1:44642] 1 1 1 22 2 1 1 14 5 12 ...
##  $ HARGA       : num [1:44642] 168300 490050 314050 188100 133650 ...
##  $ TOTAL       : num [1:44642] 113603 302606 196972 2824322 185426 ...
##  $ TOTAL FAKTUR: num [1:44642] 113603 302606 196972 2824322 185426 ...
##  $ DPP         : num [1:44642] 103275 275096 179066 2567565 168569 ...
##  $ PPN         : num [1:44642] 10328 27510 17907 256757 16857 ...

dim(ul_clean)

## [1] 44642    10

ul_clean <- 
ul_clean[order(ul_clean$TANGGAL),]
ul_clean

## # A tibble: 44,642 x 10
##    TANGGAL             CUSTOMER        SALES  NAMA_BRG      BANYAK  HARGA  TOTAL
##    <dttm>              <fct>           <fct>  <fct>          <dbl>  <dbl>  <dbl>
##  1 2021-01-04 00:00:00 SINAR TIMUR - ~ SAL-M  GRACE F 5505~      1 168300 1.14e5
##  2 2021-01-04 00:00:00 ANEKA BANGUNAN~ SAL-M  INTERLOK GL6~      1 187000 1.28e5
##  3 2021-01-04 00:00:00 PHANG ENG LIE   SAL-6  ROMAN F 2631~      4 147400 4.48e5
##  4 2021-01-04 00:00:00 TERATAI KERAMIK SAL-B6 R.GRES F 229~      4 153450 4.26e5
##  5 2021-01-04 00:00:00 TERATAI KERAMIK SAL-B6 PRECUT WA 33~     25 239250 4.26e6
##  6 2021-01-04 00:00:00 PT.SINAR BARU ~ SAL-B6 R.GRANIT GTA~     13 275000 2.24e6
##  7 2021-01-04 00:00:00 MMS - CINERE    SAL-B6 INTERLOK GL6~     12 229900 1.88e6
##  8 2021-01-04 00:00:00 SINAR MAJU - K~ SAL-B6 RECTIFIED W ~     15 201300 2.15e6
##  9 2021-01-04 00:00:00 SINAR MAJU - K~ SAL-B6 GRACE F 3320~     10 177650 1.27e6
## 10 2021-01-04 00:00:00 PELITA HARAPAN~ SAL-8  STEPTILE BGT~      7 193600 1.36e6
## # ... with 44,632 more rows, and 3 more variables: TOTAL FAKTUR <dbl>,
## #   DPP <dbl>, PPN <dbl>

plot_penjualan_persale <- ul_clean %>%
  group_by(SALES) %>%
  summarise(total = mean (`TOTAL FAKTUR`)) %>%
  arrange(desc(total)) %>%
  mutate(label = glue("Rata-rata penjualan = {total}")) %>%
  ggplot(plot, mapping = aes(
    x = reorder(SALES, desc(total)),
    y = total,
    fill = SALES,
    text = label
  )) + geom_bar(stat = "identity") + xlab("Sales Group") + ylab("Total Rata-rata Penjualan Perbulan") + theme_minimal() + theme(legend.position = "none") + guides(x = guide_axis(n.dodge = 2), y.sec = guide_axis())

ggplotly(plot_penjualan_persale, tooltip = "label")

weekly_quantity <-
  ul_clean %>% mutate(date = week(TANGGAL)) %>%  # converted become weekly
  group_by(date, SALES) %>% # group by berdasarkan category tiap team sales dan perminggu
  summarise(quantity = sum(BANYAK),
            # transaksi yang sukses
            total = sum(`TOTAL FAKTUR`)) # total perminggu tiap team sales

penjualan_tahunan <-
  weekly_quantity %>%
  mutate(
    label = glue(
      "Jumlah Transaksi = {quantity}
    Total Rupiah = {comma(total)}
    Weekly = {date}
    Team Sales = {SALES}"
    )
  ) %>%
filter() %>% 
  ggplot(mapping = aes(
    x = date,
    y = quantity,
    group = SALES,
    color = SALES,
    text = label
  )) + geom_line(lwd = 1) + geom_point()
ggtitle("Team sales penjualan 2020-2021") +
  theme_minimal()

## NULL

ggplotly(penjualan_tahunan, tooltip = "label")

pembeli_terbanyak_12ribu <-
  ul_clean %>% mutate(bulan = month(TANGGAL)) %>%
  group_by(CUSTOMER) %>% # mengroupkan bulan nya
  summarise(quantity = sum(BANYAK),
            PPN = sum(PPN),
            DPP = sum(DPP)) %>%
  filter(quantity > 12000) %>%
  mutate(label = glue("Pembeli = {CUSTOMER}
                                 Total Transaksi = {quantity}")) %>% 
ggplot(
  mapping = aes(
    x = reorder(CUSTOMER, desc(quantity)),
    y = quantity,
    fill = quantity,
    text = label
  )
) + geom_bar(stat = "identity") + theme(legend.position = "none") + labs(x = "Nama Perusahaan", y = "Jumlah transaksi") + guides(x = guide_axis(n.dodge = 2), y.sec = guide_axis()) + scale_fill_continuous()

ggplotly(pembeli_terbanyak_12ribu, tooltip = "label")

weekly_buyers <- 
ul_clean %>% transmute(TANGGAL,BANYAK) %>% 
  mutate(week = week(TANGGAL)) %>% 
  group_by(week) %>% 
  summarise(stock = sum(BANYAK))
weekly_buyers

## # A tibble: 52 x 2
##     week stock
##    <dbl> <dbl>
##  1     1  8762
##  2     2 16530
##  3     3 12663
##  4     4 16440
##  5     5 15091
##  6     6 15691
##  7     7  9937
##  8     8 14848
##  9     9 16140
## 10    10 11907
## # ... with 42 more rows

training1 <- weekly_buyers[1:42,]
testing1<- weekly_buyers[43:52,]

ts_1 <- ts(training1$stock, frequency = 365/7, start = c(2021,1)) # perminggu dalam 1 tahun
adf.test(ts_1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_1
## Dickey-Fuller = -1.9808, Lag order = 3, p-value = 0.5812
## alternative hypothesis: stationary

tail(ul_clean)

## # A tibble: 6 x 10
##   TANGGAL             CUSTOMER         SALES NAMA_BRG      BANYAK  HARGA   TOTAL
##   <dttm>              <fct>            <fct> <fct>          <dbl>  <dbl>   <dbl>
## 1 2021-12-31 00:00:00 AGUS SUDJONO     SAL-~ R.GRANIT GT8~      1 827376  5.11e5
## 2 2021-12-31 00:00:00 SEDERHANA        SAL-~ RECTIFIED W ~     28 187550  3.66e6
## 3 2021-12-31 00:00:00 SEDERHANA        SAL-~ RECTIFIED F ~      4 198550  5.66e5
## 4 2021-12-31 00:00:00 CV. ANDALAN CAH~ SAL-~ STEPNOSING R~   3500   9482  3.32e7
## 5 2021-12-31 00:00:00 CAHAYA TERANG -~ SAL-~ RECTIFIED W ~      1 187550  1.34e5
## 6 2021-12-31 00:00:00 NUSANTARA BANGU~ SAL-M R.GRANIT GT6~      5 353650  1.11e6
## # ... with 3 more variables: TOTAL FAKTUR <dbl>, DPP <dbl>, PPN <dbl>

ts_1 %>% diff() %>% tsdisplay()

ts_1 %>% diff() %>% adf.test()

## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -4.5006, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

analisis_sisaan <- function(model, std = TRUE){
  if (std == TRUE){
    res.model = rstandard(model)
  }
  else
    {
    res.model = residuals(model)
  }
  
  par(mfrow=c(2,2))
  plot(res.model,type='o',ylab='Standardised residuals',
       main="Time series plot of standardised residuals")
  abline(h=0)
  qqnorm(res.model,main="QQ plot of standardised residuals")
  qqline(res.model, col = 2)
  acf(res.model,main="ACF of standardised residuals")
  k=0
  LBQPlot(res.model, lag.max = length(model$residuals)-1 ,
          StartLag = k + 1, k = 0, SquaredQ = FALSE)
  par(mfrow=c(1,1))
}

ts_1 %>% auto.arima(d = 1, stationary = F) 

## Series: . 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.1666  -0.7922
## s.e.  0.2195   0.1415
## 
## sigma^2 estimated as 15175555:  log likelihood=-396.49
## AIC=798.98   AICc=799.63   BIC=804.12

options(warn = -2)

ts_1 %>% auto.arima(d = 1, stationary = F) %>% lmtest::coeftest()

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.16656    0.21945  0.7590    0.4479    
## ma1 -0.79225    0.14150 -5.5988 2.158e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ts_1 %>% auto.arima(d = 1, stationary = F) %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 5.6409, df = 6, p-value = 0.4646
## 
## Model df: 2.   Total lags used: 8

ts_1 %>% auto.arima(d = 1, stationary = F) %>% analisis_sisaan()

# model bandingan overfitting

ACF1 <- ts_1 %>% diff() %>% Acf(plot = F)
PACF1 <- ts_1 %>% diff() %>% Pacf(plot = F)

par(mar=c(3,2,3,0))
plot(ACF1, main = "Stock penjualan perbulan XYZ Time Series ACF")

plot(PACF1, main = "Stock penjualan perbulan XYZ Time Series PACF")

options(warn = -2)

ts_1 %>% Arima(order=c(1,1,1))

## Series: . 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.1666  -0.7922
## s.e.  0.2195   0.1415
## 
## sigma^2 estimated as 15175555:  log likelihood=-396.49
## AIC=798.98   AICc=799.63   BIC=804.12

ts_1 %>% Arima(order=c(3,1,2))

## Series: . 
## ARIMA(3,1,2) 
## 
## Coefficients:
##           ar1     ar2      ar3      ma1      ma2
##       -0.2592  0.0294  -0.2897  -0.3469  -0.2098
## s.e.   0.3331  0.2699   0.1928   0.3337   0.2644
## 
## sigma^2 estimated as 15220344:  log likelihood=-395.02
## AIC=802.03   AICc=804.5   BIC=812.31

ts_1 %>% Arima(order=c(1,1,1)) %>% lmtest::coeftest()

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.16656    0.21945  0.7590    0.4479    
## ma1 -0.79225    0.14150 -5.5988 2.158e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ts_1 %>% Arima(order=c(3,1,2)) %>% lmtest::coeftest()

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.25923    0.33306 -0.7783   0.4364
## ar2  0.02942    0.26992  0.1090   0.9132
## ar3 -0.28975    0.19281 -1.5028   0.1329
## ma1 -0.34688    0.33370 -1.0395   0.2986
## ma2 -0.20978    0.26443 -0.7933   0.4276

ts_1 %>% Arima(order=c(1,1,1)) %>% analisis_sisaan()

ts_1 %>% Arima(order=c(3,1,2)) %>% analisis_sisaan()

ts_1 %>% Arima(order=c(1,1,1)) %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 5.6409, df = 6, p-value = 0.4646
## 
## Model df: 2.   Total lags used: 8

ts_1 %>% Arima(order=c(3,1,2)) %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,2)
## Q* = 3.5757, df = 3, p-value = 0.3111
## 
## Model df: 5.   Total lags used: 8

model_arima3 <- ts_1 %>% arima(order = c(1,1,1),seasonal = list(order = c(1,0,0), period = 12), method = "ML")
model_arima2 <- ts_1 %>% auto.arima(d = 1, stationary = F)
model_ets1 <- ts_1 %>% ets()
model_tbats1 <- ts_1 %>% tbats(use.trend = F, use.damped.trend = F,
                                 use.box.cox = F, use.arma.errors = T)
model_exsmot1 <- ts_1 %>% HoltWinters(gamma = F, beta = 0.1, alpha = 0.3)

auto.arima(ts_1, trace=TRUE)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2)            with drift         : Inf
##  ARIMA(0,1,0)            with drift         : 793.2523
##  ARIMA(1,1,0)            with drift         : 785.7222
##  ARIMA(0,1,1)            with drift         : 784.178
##  ARIMA(0,1,0)                               : 791.1848
##  ARIMA(1,1,1)            with drift         : 783.6381
##  ARIMA(2,1,1)            with drift         : Inf
##  ARIMA(1,1,2)            with drift         : 784.4443
##  ARIMA(0,1,2)            with drift         : 785.9104
##  ARIMA(2,1,0)            with drift         : 788.6544
##  ARIMA(1,1,1)                               : 781.3218
##  ARIMA(0,1,1)                               : 782.5907
##  ARIMA(1,1,0)                               : 783.4559
##  ARIMA(2,1,1)                               : 783.6508
##  ARIMA(1,1,2)                               : 782.3807
##  ARIMA(0,1,2)                               : 784.4344
##  ARIMA(2,1,0)                               : 786.3107
##  ARIMA(2,1,2)                               : 784.3377
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(1,1,1)                               : 799.6273
## 
##  Best model: ARIMA(1,1,1)

## Series: ts_1 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.1666  -0.7922
## s.e.  0.2195   0.1415
## 
## sigma^2 estimated as 15175555:  log likelihood=-396.49
## AIC=798.98   AICc=799.63   BIC=804.12

forecast_arima <- 
forecast(model_arima2, h = 10)
plot_forecast(forecast_arima,title =  "Auto Arima Model",color = "Purple", Xtitle = "Weekly",
              Ytitle = "Stock paid")

forecast_ets <- forecast(model_ets1, h = 10)
plot_forecast(forecast_ets, title =  "ETS Model",color = "RED", Xtitle = "Weekly",
              Ytitle = "Stock paid") 

ul_clean %>% group_by(CUSTOMER) %>% 
  summarise(Jumlah_Transaksi = n(),
            total_pengeluaran = sum(`TOTAL FAKTUR`),
            total_pengeluaran_in_M = round(total_pengeluaran/1000000000, digits = 2)) %>% 
  arrange(desc(total_pengeluaran))

## # A tibble: 674 x 4
##    CUSTOMER                   Jumlah_Transaksi total_pengeluar~ total_pengeluar~
##    <fct>                                 <int>            <dbl>            <dbl>
##  1 CV.KEMENANGAN JAYA TRADING             1414     27463956595.            27.5 
##  2 PT. SUMBER JAYA LANGGENG               1325      9457786110.             9.46
##  3 PT. SEHATI BANGUNAN ABADI               604      7140263517.             7.14
##  4 CV. HAKA ANUGRAH JAYA                   756      4560503585.             4.56
##  5 DEPOK KERAMIK - DEPOK*                  714      4189033101.             4.19
##  6 INTI GADING KERAMIK                     470      3856634187.             3.86
##  7 BERGUNA JAYA - P.POLIM*                 709      3796426639.             3.8 
##  8 CAHAYA INDAH K. - P.AREN*               715      3744965404.             3.74
##  9 BINTANG KERAMIK - JAKBAR*              1119      3599621570.             3.6 
## 10 WILLY HANSEN                            687      3410855618.             3.41
## # ... with 664 more rows

ul_clean %>% group_by(NAMA_BRG) %>% 
  summarise(Jumlah = sum(BANYAK),
            per_transaksi = n()) %>% 
  arrange(desc(Jumlah)) 

## # A tibble: 1,316 x 3
##    NAMA_BRG                 Jumlah per_transaksi
##    <fct>                     <dbl>         <int>
##  1 RECTIFIED W 63801R KW 1   12936           668
##  2 STEPNOSING SN 13606 KW 1  11187            13
##  3 R.GRANIT GT609860FR KW 1   9947           182
##  4 ROMAN W 40106 KW 1         9730           335
##  5 R.GRANIT GT609866FR KW 1   9343           222
##  6 RECTIFIED W 63364R KW 1    8963           297
##  7 R.GRES F 447106 KW1        8452           132
##  8 R.GRES F 447442 KW1        8114           112
##  9 RECTIFIED W 63765R KW 1    7919           129
## 10 GRACE F 557379 KW 1        7256            44
## # ... with 1,306 more rows

df1 <- 
ul %>% mutate(week = week(TANGGAL)) %>% 
  group_by(NAMA_BRG, week) %>% 
  summarise(quantity_barang = sum(BANYAK), # banyaknya quantity terjual 
            per_transaksi = n(),
            total_in_rp = sum(`TOTAL FAKTUR`)) %>% filter(NAMA_BRG == "RECTIFIED W 63801R KW 1") %>% mutate(NAMA_BRG = as.factor(NAMA_BRG))
head(df1)

## # A tibble: 6 x 5
## # Groups:   NAMA_BRG [1]
##   NAMA_BRG                 week quantity_barang per_transaksi total_in_rp
##   <fct>                   <dbl>           <dbl>         <int>       <dbl>
## 1 RECTIFIED W 63801R KW 1     1              93             6   67632896.
## 2 RECTIFIED W 63801R KW 1     2             349            15   73045701.
## 3 RECTIFIED W 63801R KW 1     3              52             9   16088016.
## 4 RECTIFIED W 63801R KW 1     4             294            16   51035688.
## 5 RECTIFIED W 63801R KW 1     5             225            11   60508205.
## 6 RECTIFIED W 63801R KW 1     6             132            12   37976807.

dataset_barang <- 
ul_clean %>% mutate(week = week(TANGGAL)) %>% 
  group_by(NAMA_BRG, week) %>% summarise(quantity_barang = sum(BANYAK), # banyaknya quantity terjual 
            per_transaksi = n(),
            total_in_rp = sum(`TOTAL FAKTUR`)) %>% mutate(NAMA_BRG = as.factor(NAMA_BRG)) %>% 
  arrange(desc(week))
head(dataset_barang)

## # A tibble: 6 x 5
## # Groups:   NAMA_BRG [6]
##   NAMA_BRG              week quantity_barang per_transaksi total_in_rp
##   <fct>                <dbl>           <dbl>         <int>       <dbl>
## 1 GRACE F 332132P KW 1    53               3             1     390308.
## 2 GRACE F 33554P KW 1     53               3             1    3469874.
## 3 GRACE F 550005 KW 1     53              32             1    3706560 
## 4 GRACE F 550012 KW 1     53               3             1     334008.
## 5 GRACE F 550500 KW 1     53              45             1    5212350 
## 6 GRACE F 550501 KW 1     53               4             1     463320

ul_clean %>% group_by(SALES, CUSTOMER) %>% 
  summarise(pembelian_tiap_sale_customer = n(),
            stock_terbeli = sum(BANYAK),
            total = sum(`TOTAL FAKTUR`),
            total_M = round(total/1000000000, digits = 2)) %>% 
  arrange(desc(total_M)) 

## # A tibble: 849 x 6
## # Groups:   SALES [18]
##    SALES  CUSTOMER         pembelian_tiap_sale_~ stock_terbeli     total total_M
##    <fct>  <fct>                            <int>         <dbl>     <dbl>   <dbl>
##  1 SAL-Y  CV.KEMENANGAN J~                   888         29368   1.63e10   16.3 
##  2 SAL-41 CV.KEMENANGAN J~                   388         20823   7.44e 9    7.44
##  3 SAL-Y  PT. SUMBER JAYA~                   743         20925   6.39e 9    6.39
##  4 SAL-L  CV. HAKA ANUGRA~                   755         10875   4.55e 9    4.55
##  5 SAL-B8 DEPOK KERAMIK -~                   714         14085   4.19e 9    4.19
##  6 SAL-B8 BERGUNA JAYA - ~                   709         12740   3.80e 9    3.8 
##  7 SAL-8  CV.KEMENANGAN J~                   130          5344   3.72e 9    3.72
##  8 SAL-B6 BINTANG KERAMIK~                  1115         12691   3.59e 9    3.59
##  9 SAL-8  WILLY HANSEN                       687         11146   3.41e 9    3.41
## 10 SAL-B8 CV. TIMUR BARU                     666         11501   3.30e 9    3.3 
## # ... with 839 more rows

dataset_viz <-
  dataset_barang %>% mutate(label = glue(
    "Weekly = {week}
  Total Rupiah = {comma(total_in_rp)}
  Product = {NAMA_BRG}"
  )) %>%
  filter(
    NAMA_BRG == "STEPNOSING SN 13606 KW 1" |
      NAMA_BRG == "RECTIFIED W 63801R KW 1" |
      NAMA_BRG == "R.GRANIT GT609860FR KW 1" |
      NAMA_BRG == "ROMAN W 40106 KW 1" |
      NAMA_BRG == "R.GRANIT GT609866FR KW 1"
  ) %>%
  ggplot(
    mapping = aes(
      x = week,
      y = quantity_barang,
      group = NAMA_BRG,
      fill = NAMA_BRG,
      color = NAMA_BRG,
      text = label
    )
  ) + geom_point() + geom_smooth(method = lm, se = F) + labs(title = "Regression Linear Product")

ggplotly(dataset_viz, tooltip = "label")

dataset <-
  ul_clean %>% filter(NAMA_BRG == "RECTIFIED W 63801R KW 1") %>%
  mutate(TANGGAL = TANGGAL %>% floor_date(unit = "week")) %>%
  group_by(TANGGAL) %>%
  summarise(var = sum(BANYAK)) %>%
  mutate(month = month(TANGGAL)) # mengubah jenis data pada kolom `month` menjadi factor
head(dataset, 10)

## # A tibble: 10 x 3
##    TANGGAL               var month
##    <dttm>              <dbl> <dbl>
##  1 2021-01-03 00:00:00   200     1
##  2 2021-01-10 00:00:00   265     1
##  3 2021-01-17 00:00:00    55     1
##  4 2021-01-24 00:00:00   447     1
##  5 2021-01-31 00:00:00   106     1
##  6 2021-02-07 00:00:00    72     2
##  7 2021-02-14 00:00:00   171     2
##  8 2021-02-21 00:00:00   360     2
##  9 2021-02-28 00:00:00   111     2
## 10 2021-03-07 00:00:00   393     3

dataset_sales <-
  xts(x = dataset$var, order.by = dataset$TANGGAL)
plot(dataset_sales, main = "Penjualan Sampel Produk Tahun 2021") # visualisasi data secara interaktif

plotly1 <- 
ul_clean %>% 
  group_by(TANGGAL, SALES, CUSTOMER) %>% # group by berdasarkan category tiap team sales dan perminggu
  summarise(quantity = sum(BANYAK),
            total_faktur = sum(`TOTAL FAKTUR`),
            product = (NAMA_BRG)) %>% mutate(label = glue (
              "Total barang terjual = {quantity}
    Total Rupiah = {comma(total_faktur)}
    Date = {TANGGAL}
    Team Sales = {SALES}
    Customer = {CUSTOMER}"
            )) %>% filter(SALES == "SAL-L" | SALES == "SAL-K" | SALES == "SAL-2") %>% 
  
  ggplot(mapping = aes(x = TANGGAL,
                       y = quantity,
                       color = SALES,
                       size = quantity,
                       text = label)) + geom_point(alpha=0.7) +
    scale_size(range = c(1.4, 10)) +scale_color_viridis(discrete=TRUE, guide=FALSE)
ggtitle("Team sales penjualan 2020-2021") +
  theme_minimal()

## NULL

ggplotly(plotly1, tooltip = "label")

ul_clean %>% group_by(SALES) %>% summarise(var = n()) %>% 
  arrange(desc(var))

## # A tibble: 18 x 2
##    SALES    var
##    <fct>  <int>
##  1 SAL-B8  9723
##  2 SAL-B6  7761
##  3 SAL-8   7176
##  4 SAL-Y   6731
##  5 SAL-M   5687
##  6 SAL-41  3257
##  7 SAL-6   1616
##  8 SAL-L   1251
##  9 SAL-2    607
## 10 SAL-C4   450
## 11 SAL-K    204
## 12 SAL-74    57
## 13 SAL-20    31
## 14 SAL-49    23
## 15 SAL-67    22
## 16 SAL-G     21
## 17 SAL-H     18
## 18 SAL-26     7