Hourly Demand Forecasting of Scotty Ride-Sharing Service

Github : Hourly Demand Forecasting of Scotty Ride-sharing Service

1 Intro

Scotty Technologies Inc. (“Scotty”) adalah sebuah perusahaan start-up teknologi yang didirikan pada tahun 2017 di Istanbul, Turkey. Salah satu layanan utama dari Scotty adalah motorcycle ride-sharing atau yang akrab kita ketahui Ojek Online. Menurut informasi yang dilansir dari markets.businessinsider.com, Scotty berencana untuk menjadi super-app pertama di Turkey dengan tekonologi dan model bisnis yang distruptif, serta berencana untuk melanjutkan pertumbuhan yang menjanjikan dan mencapai profitabilitas dalam waktu dekat. Berbicara mengenai pertumbuhan, tentu saja tidak luput dari demand. Oleh karena itu, projek ini ditujukan untuk membuat model analisa dan memprediksi demand transaksi per-jam pada layanan motorcycle ride-sharing di Scotty. Dalam proses pembuatan model ini, kita menggunakan data transaksi Scotty dari 2017-10-01 sampai 2017-12-02.

2 Problem Identification

Problem yang diselesaikan yaitu membuat model untuk memprediksi demand transaksi per-jam pada layanan motorcycle ride-sharing di Scotty. Namun yang hendak diprediksi dalam rentang waktu berapa lama? Mari kita cek dulu supaya dapat membantu proses pemodelan yang akan kita buat.

# Read Data Test
scotty_actual_test <- read_csv("data_input/scotty_ts/data-test.csv")
glimpse(scotty_actual_test)

#> Observations: 504
#> Variables: 3
#> $ src_sub_area <chr> "sxk97", "sxk97", "sxk97", "sxk97", "sxk97", "sxk97", ...
#> $ datetime     <dttm> 2017-12-03 00:00:00, 2017-12-03 01:00:00, 2017-12-03 ...
#> $ demand       <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...

Struktut data diatas belum sesuai. Mari sesuaikan struktur datanya dan meilhat rentang waktu yang hendak diprediksi dahulu:

scotty_actual_test <- scotty_actual_test %>% 
   mutate(
      # Transaksi Per-Jam
      datetime = floor_date(datetime, unit="hour"),
      demand = as.integer(demand)
   )
summary(scotty_actual_test)

#>  src_sub_area          datetime                       demand   
#>  Length:504         Min.   :2017-12-03 00:00:00   Min.   : NA  
#>  Class :character   1st Qu.:2017-12-04 17:45:00   1st Qu.: NA  
#>  Mode  :character   Median :2017-12-06 11:30:00   Median : NA  
#>                     Mean   :2017-12-06 11:30:00   Mean   :NaN  
#>                     3rd Qu.:2017-12-08 05:15:00   3rd Qu.: NA  
#>                     Max.   :2017-12-09 23:00:00   Max.   : NA  
#>                                                   NA's   :504

Data diatas merupakan summary dari data yang hendak diprediksi. Berdasarkan informasi tersebut, dapat disimpullkan bawah masalah yang perlu kita selesaikan yaitu melakukan prediksi demand per-jam di sub area sxk97, sxk9e dan sxk9 dalam rentang waktu dari 2017-12-03 00:00 sampai 2017-12-09 23:00 atau 24 jam * 1 minggu. Oke, mari kita mulai proses nya.

3 Data Preparation

3.1 Read Data

Berikut ini adalah data yang kita punya untuk melakukan analisa dan prediksi.

scotty_input <- read_csv("data_input/scotty_ts/data-train.csv")

5 Top Line Data

head(scotty_input)

5 Bottom Line Data

tail(scotty_input,5)

---

3.2 Variable Description

Variable	Description
id	Transaction id
trip_id	Trip id
driver_id	Driver id
rider_id	Rider id
start_time	Transaction Time
src_lat	Request source latitude
src_lon	Request source longitude
src_area	Request source area
src_sub_area	Request source sub-area
dest_lat	Requested destination latitude
dest_lon	Requested destination longitude
dest_area	Requested destination area
dest_sub_area	Requested destination sub-area
distance	Trip distance (in KM)
status	Trip status (all status considered as a demand)
confirmed_time_sec	Time different from request to confirmed (in seconds)

---

3.3 Missing & Duplicate Value

Missing Value
Jika dilihat terdapat missing value pada variabel trip_id dan driver_id, namun tidak menjadi masalah karena pada case ini kita butuh data waktu, sub area dan demand.

colSums(is.na(scotty_input))

#>                 id            trip_id          driver_id           rider_id 
#>                  0               4676               4676                  0 
#>         start_time            src_lat            src_lon           src_area 
#>                  0                  0                  0                  0 
#>       src_sub_area           dest_lat           dest_lon          dest_area 
#>                  0                  0                  0                  0 
#>      dest_sub_area           distance             status confirmed_time_sec 
#>                  0                  0                  0                  0

Duplicate Value
Dataset ini tidak memiliki data transaksi yang duplikat, sehingga bisa kita lanjutkan ke tahap pre-processing.

data.frame(jumlah.observasi.awal=length(scotty_input$id),
           jumlah.observasi.unik=length(unique(scotty_input$id)))

---

4 Data Pre-processing

4.1 Data Structure

glimpse(scotty_input)

#> Observations: 90,113
#> Variables: 16
#> $ id                 <chr> "59d005e1ffcfa261708ce9cd", "59d0066affcfa261708...
#> $ trip_id            <chr> "59d005e9cb564761a8fe5d3e", NA, "59d006c131e39c6...
#> $ driver_id          <chr> "59a892c5568be44b2734f276", NA, "599dc0dfa5b4fd5...
#> $ rider_id           <chr> "59ad2d6efba75a581666b506", "59cd704bcf482f6ce2f...
#> $ start_time         <dttm> 2017-10-01 00:00:17, 2017-10-01 00:02:34, 2017-...
#> $ src_lat            <dbl> 41.07047, 41.07487, 41.04995, 41.05287, 41.06760...
#> $ src_lon            <dbl> 29.01945, 28.99528, 29.03107, 28.99522, 28.98827...
#> $ src_area           <chr> "sxk9", "sxk9", "sxk9", "sxk9", "sxk9", "sxk9", ...
#> $ src_sub_area       <chr> "sxk9s", "sxk9e", "sxk9s", "sxk9e", "sxk9e", "sx...
#> $ dest_lat           <dbl> 41.11716, 41.08351, 41.04495, 41.08140, 41.02125...
#> $ dest_lon           <dbl> 29.03650, 29.00228, 28.98192, 28.98197, 29.11316...
#> $ dest_area          <chr> "sxk9", "sxk9", "sxk9", "sxk9", "sxk9", "sxk9", ...
#> $ dest_sub_area      <chr> "sxk9u", "sxk9e", "sxk9e", "sxk9e", "sxk9q", "sx...
#> $ distance           <dbl> 5.379250, 1.126098, 4.169492, 3.358296, 11.69357...
#> $ status             <chr> "confirmed", "nodrivers", "confirmed", "confirme...
#> $ confirmed_time_sec <dbl> 8, 0, 32, 65, 0, 27, 32, 30, 0, 24, 9, 0, 118, 1...

Dataset yang dimiliki terdiri dari 90.133 observasi dan 16 variabel. Jika dilihat dari struktur datanya, beberapa variabel memiliki tipe data yang belum sesuai. Namun, karena projek ini bertujuan untuk melakukan time-series forecasting terhadap demand per-jam pada layanan motorcycle ride-sharing di Scotty, maka data yang kita perlukan adalah waktu transaksi per-jam (start_time) dan demand (count(id)) pada setiap Sub-Ara (src_sub_area). Maka dari Berikut proses dan datanya:

scotty_input <- scotty_input %>% 
   mutate(
      # Transaksi Per-Jam
      datetime = floor_date(start_time, unit="hour")
   ) %>% 
   group_by(src_sub_area, datetime) %>% 
   summarise(demand = length(status)) %>% 
   ungroup()

scotty_input

4.2 Missing Sequential Datetime

Dalam melakukan analisa dan prediksi time series ada beberapa hal yang harus terpenuhi antara lain:
1. Data tidak boleh ada yang missing
2. Data harus terurut berdasarkan periode waktunya
3. Tidak boleh ada waktu atau periode yang terlewat/bolong

Terkait case ini, Jam opearasional Scotty yaitu 24jam/hari, sehingga kita harus memastikan seluruh data jam sudah lengkap secara sekuensial dari rentang tanggal minimum hingga tanggal maximum. Jika terdapat deret waktu yang kosong, maka dapat dilakukan imputasi data dengan asumsi bahwa pada waktu tersebut memang tidak terdapat transaksi, sehingga nilai demand = 0.

scotty_input <- scotty_input %>% 
  group_by(src_sub_area) %>% 
  # PAD tanggal sekuensial
  padr::pad(start_val = min(scotty_input$datetime), end_val = max(scotty_input$datetime)) %>% 
  ungroup() %>% 
  distinct() %>% 
  mutate(
    # replace NA Value menjadi 0
    demand = ifelse(is.na(demand),0,demand)
  ) %>% 
  arrange(datetime)

Berikut data hari yang tidak memiliki transaksi:

scotty_input %>% filter(demand==0)

Data diatas menunjukan data waktu tanpa transaksi dari sub Area sxk97, sxk9e dan sxk9s. Total waktu tanpa transaksi yaitu 311 Jam dimana setiap baris merepresentasikan 1 jam.

4.3 Data Summary

summary(scotty_input)

#>  src_sub_area          datetime                       demand      
#>  Length:4536        Min.   :2017-10-01 00:00:00   Min.   :  0.00  
#>  Class :character   1st Qu.:2017-10-16 17:45:00   1st Qu.:  5.00  
#>  Mode  :character   Median :2017-11-01 11:30:00   Median : 15.00  
#>                     Mean   :2017-11-01 11:30:00   Mean   : 19.87  
#>                     3rd Qu.:2017-11-17 05:15:00   3rd Qu.: 27.00  
#>                     Max.   :2017-12-02 23:00:00   Max.   :217.00

Informasi di atas merupakan summary keseluruhan data yang kita butuhkan untuk melakukan melakukan forecasting. Dapat diketahui bahwa dataset ini merupakan data transaksi dari 3 sub area yaitu sxk97, sxk9e dan sxk9s yang terjadi dari 2017-10-01 00:00 sampai 2017-12-02 23:00 dengan jumlah minimum 0 demand dan jumlah maximum 217 demand. Demand 0 diasumsikan bahwa pada waktu tertentu memang tidak terdapat transaksi.

---

5 Exploratory Data Analysis

5.1 Demand Analysis by Sub-Area

Berikut ini visualisasi demand berdasarkan sub-area dari dari 2017-10-01 00:00 sampai 2017-12-02 23:00. Titik/Poin berwarna merah, hijau dan biru menandakan data deret waktu tersebut kosong atau tidak terdapat transaksi pada waktu tersebut.

scotty_demand <- scotty_input 

scotty_demand <- scotty_demand %>% 
    mutate(
      weekdays = weekdays(datetime),
      popup = glue(
      "Category: {src_sub_area}
      Transaction Time: {datetime}
      Day: {weekdays}
      Demand: {demand}"
    )) 


plot_scotty_demand <- ggplot(scotty_demand,aes(x=datetime, y=demand),show.legend = FALSE)+
   #geom_point()+
  geom_line(show.legend = FALSE,size=0.3)+
  geom_point(aes(text=popup),size=0.05)+
  geom_point(data = scotty_demand %>% filter(demand==0),
             aes(x=datetime, y=demand, color=src_sub_area, text=popup),size=0.5, show.legend = FALSE)+
  labs(
    x = NULL, 
    y = NULL,
    title = "Demand by Sub-Area"
  ) +
  scale_x_datetime(date_breaks = "7 day", labels = date_format("%b %d"))+
  facet_wrap(~ src_sub_area, scale = "free", ncol = 1)+
  theme_minimal()+
  theme(
    legend.position = "none", 
    plot.title = element_text(hjust = 0.5),
    panel.spacing = unit(1, "lines"))

ggplotly(p = plot_scotty_demand, tooltip="text") %>% 
  layout(showlegend=FALSE)

Pada grafik diatas dapat dilihat bahwa jumlah demand dari ketiga sub-area relatif sama dan juga terdapat pola seasonal meskipun belum terlihat jelas detailnya. Selain itu, juga dapat dilihat demand pada Jumat,03 November 2017 pukul 18.00 dan Jumat, 24 November 2017 pukul 18.00 jauh lebih tinggi dari hari lainnya. Sayangnya data yang kita miliki belum dapat menjawab hal tersebut, namun mari kita fokus ke demand. Berikut adalah demand per-harinya:

# Daily Demand
polar_day <- scotty_demand %>% 
  mutate(day=weekdays(datetime)) %>% 
  group_by(src_sub_area,day) %>% 
  summarise(demand_per_day=sum(demand)) %>% 
  ungroup() %>% 
  mutate(
    day = ordered(day, levels=c("Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"))
  )

serror <- function(x) sqrt(var(x)/length(x)) 

ggplot(polar_day %>% filter(src_sub_area=="sxk97"), aes(x=day, y=demand_per_day, fill = day)) +
  geom_bar(width = 1, stat = "identity", color = "white", show.legend = FALSE) +
  geom_errorbar(aes(ymin = demand_per_day - serror(demand_per_day), 
                    ymax = demand_per_day + serror(demand_per_day), 
                    color = day), 
                    width = .2) + 
  scale_y_continuous(breaks = 0:nlevels(day)) +
  labs(
    title = "Daily Demand in Sub-Area sxk97"
  )+
 theme_minimal() +
  theme(axis.title = element_blank(),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5, size=12, face="bold"),
        axis.text.y = element_blank(),
        axis.text.x=element_text(face="bold"))+
  coord_polar() -> day_polar1

ggplot(polar_day %>% filter(src_sub_area=="sxk9e"), aes(x=day, y=demand_per_day, fill = day)) +
  geom_bar(width = 1, stat = "identity", color = "white", show.legend = FALSE) +
  geom_errorbar(aes(ymin = demand_per_day - serror(demand_per_day), 
                    ymax = demand_per_day + serror(demand_per_day), 
                    color = day), 
                    width = .2) + 
  scale_y_continuous(breaks = 0:nlevels(day)) +
  labs(
    title = "Daily Demand in Sub-Area sxk9e"
  )+
  theme_minimal() +
  theme(axis.title = element_blank(),
        legend.position = "none",
         plot.title = element_text(hjust = 0.5, size=12, face="bold"),
        axis.text.y = element_blank(),
        axis.text.x=element_text(face="bold"))+
  coord_polar() ->day_polar2


ggplot(polar_day %>% filter(src_sub_area=="sxk9s"), aes(x=day, y=demand_per_day, fill = day)) +
  geom_bar(width = 1, stat = "identity", color = "white", show.legend = FALSE) +
  geom_errorbar(aes(ymin = demand_per_day - serror(demand_per_day), 
                    ymax = demand_per_day + serror(demand_per_day), 
                    color = day), 
                    width = .2) + 
  scale_y_continuous(breaks = 0:nlevels(day)) +
  labs(
    title = "Daily Demand in Sub-Area sxk9s"
  )+
 theme_minimal() +
  theme(axis.title = element_blank(),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5,size=12, face="bold"),
        axis.text.y = element_blank(),
        axis.text.x=element_text(face="bold"))+
  coord_polar() ->day_polar3


grid.arrange(day_polar1,day_polar2,day_polar3, ncol = 3)

Grafik diatas memberikan informasi total demmand per-hari dari masing-masing Sub-Area. Secara keseluruhan, ketiga Sub-Area menunjukan pola demand per-hari yang sama, dimana total demand paling besar yaitu pada hari Jumat dan total demand paling kecil yaitu pada hari Minggu. Mari kita lihat data demand per-jamnya berikut:

# Hourly Demand
polar_hourly <- scotty_demand %>% 
  mutate(hour=hour(datetime)) %>% 
  group_by(src_sub_area,hour) %>% 
  summarise(demand_per_hour=sum(demand)) %>% 
  ungroup() %>% 
  mutate(
    hour = ordered(hour,levels=c("0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17",
                              "18","19","20","21","22","23"))
  )

ggplot(polar_hourly %>% filter(src_sub_area=="sxk97"), aes(x=hour, y=demand_per_hour, fill = hour)) +
  geom_bar(width = 1, stat = "identity", color = "white", show.legend = FALSE) +
  geom_errorbar(aes(ymin = demand_per_hour - serror(demand_per_hour), 
                    ymax = demand_per_hour + serror(demand_per_hour), 
                    color = hour), 
                    width = .2) + 
  scale_y_continuous(breaks = 0:nlevels(hour)) +
  labs(
    title = "Hourly Demand in Sub-Area sxk97",
    subtitle = "24 Hour Format"
  )+
 theme_minimal() +
  theme(axis.title = element_blank(),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5,size=12, face="bold"),
        plot.subtitle = element_text(hjust = 0.5,size=11),
        axis.text.y = element_blank(),
        axis.text.x=element_text(size=11, face="bold"))+
  coord_polar() -> hour_polar1


ggplot(polar_hourly %>% filter(src_sub_area=="sxk9e"), aes(x=hour, y=demand_per_hour, fill = hour)) +
  geom_bar(width = 1, stat = "identity", color = "white", show.legend = FALSE) +
  geom_errorbar(aes(ymin = demand_per_hour - serror(demand_per_hour), 
                    ymax = demand_per_hour + serror(demand_per_hour), 
                    color = hour), 
                    width = .2) + 
  scale_y_continuous(breaks = 0:nlevels(hour)) +
  labs(
    title = "Hourly Demand in Sub-Area sxk9e",
    subtitle = "24 Hour Format"
  )+
 theme_minimal() +
  theme(axis.title = element_blank(),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5,size=12, face="bold"),
        plot.subtitle = element_text(hjust = 0.5,size=11),
        axis.text.y = element_blank(),
        axis.text.x=element_text(size=11, face="bold"))+
  coord_polar() -> hour_polar2


ggplot(polar_hourly %>% filter(src_sub_area=="sxk9s"), aes(x=hour, y=demand_per_hour, fill = hour)) +
  geom_bar(width = 1, stat = "identity", color = "white", show.legend = FALSE) +
  geom_errorbar(aes(ymin = demand_per_hour - serror(demand_per_hour), 
                    ymax = demand_per_hour + serror(demand_per_hour), 
                    color = hour), 
                    width = .2) + 
  scale_y_continuous(breaks = 0:nlevels(hour)) +
  labs(
    title = "Hourly Demand in Sub-Area sxk9s",
    subtitle = "24 Hour Format"
  )+
 theme_minimal() +
  theme(axis.title = element_blank(),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5,size=12, face="bold"),
        plot.subtitle = element_text(hjust = 0.5,size=11),
        axis.text.y = element_blank(),
        axis.text.x=element_text(size=11, face="bold"))+
  coord_polar() -> hour_polar3


grid.arrange(hour_polar1,hour_polar2,hour_polar3, ncol = 3)

Grafik diatas memberikan informasi total demmand per-jam dari masing-masing Sub-Area. Secara keseluruhan, ketiga Sub-Area menunjukan pola demand per-jam yang sama, dimana total demand relatif tinggi dari pukul 17:00 hingga pukul 19:00 dan total demmand relatif rendah dari pukul 01:00 hingga pukul 05:00. Namun, yang berbeda adalah total demand sub-area sxk9e juga cukup tinggi pada pukul 08:00, hal ini berbeda dengan sub-area lainnya.

---

5.2 Seasonal Analysis

Sebelum melakukan forecasting kita perlu mengetahui informasi tren dan seasonal yang ada pada data. Salah satu cara untuk mengetahuinya yaitu melakukan Decompose. Decompose merupakan tahapan dalam analisis time series yang digunakan untuk menguraikan beberapa komponen yang berupa:

1. Data : data objek time-series yang diamati.
   
2. Trend : pola data secara general, cenderung untuk naik atau turun.
   
3. Seasonal : pola musiman yang membentuk pola berulang pada periode waktu yang tetap.
   
4. Residual/Error : pola yang tidak dapat ditangkap dalam trend dan seasonal.
   

Jika pada hasil decompose, trend masih membentuk sebuah pola maka dapat dicurigai masih ada seasonality yang belum ditangkap oleh frekuensi data. Seharusnya trend cenderung naik atau cendurung turun atau membentuk garis yang smooth.

Sesuai dengan masalah yang hendak diselesaikan yaitu melakukan prediksi demand per-jam pada masing-masing sub-area selama 1 minggu, maka mari kita coba decompose secara harian dan mingguan terhadap salah satu sub-area.

sxk9s <- scotty_input %>% filter(src_sub_area == "sxk9s") %>% select(-src_sub_area)

sxk9s_ts_weekly <- sxk9s %>% .$demand %>% ts(frequency = 24*7)
  
plot_sxk9s_ts_weekly <- sxk9s_ts_weekly %>% 
  tail(24*7*4) %>% 
  decompose() %>% 
   autoplot()+ 
    theme_minimal()+
    theme(
          title = element_text(size=11),
          axis.title=element_text(size=9, face="bold"),
          axis.text.x=element_text(size=9),
          axis.text.y=element_text(size=9),
          legend.position = "none",
          panel.spacing = unit(4, "pt")
          )


ggplotly(plot_sxk9s_ts_weekly) %>% 
  layout(title = list(text = paste0('Decomposition on Weekly in Sub-Area sxk9s',
                                    '<br>',
                                    '<sup>',
                                    'Deompose single seasonal time-series object using frequency = 24*7 (Weekly)',
                                    '</sup>')))

Hasil decompose objek single sesonal time-series diatas menggunakan frekuensi Mingguan pada data 1 bulan terkahir. Beberapa informasi yang bisa kita dapat yaitu:

1. Panel data merupakan data observasi time-series scotty pada data 1 bulan terkahir. Data ini sebenarnya sama seperti yang ditampilkan pada poin 5.1 Demand Analysis, namun karena pola kurang terlihat pada poin 5.1 tersebut maka kita perkecil analisanya menggunakan data 1 bulan terakhir.
   
2. Panel seasonal merupakan komponen seasonal. Jika dilihat terdapat 4 bagian grafik yang memuncak paling tinggi dan yang lainnya rendah. Berhubung data yang kita time-series yang kita buat menggunakan frekuensi mingguan, maka hal ini dapat mengindikasikan bahwa pola seasonal yang terjadi sebenarnya secara harian dan mingguan.
   
3. Panel trend merupakan komponen trend. Pola trend yang ditampilkan masih belum smooth dan terkesan masih membentuk pola pada beberapa bagian, mungkin saja terdapat pola seasonal yang belum ditangkap. Hal ini masuk akal dengan yang ditampilkan pada panel seasonal, sehingga mengindikasikan bahwa data kita merupakan multiple seasonal time series object.
   
4. Panel Remainder/Error merupakan komponen error yang menunjukan pola yang tidak dapat ditangkap dalam trend dan seasonal.
   

Hasil decompose diatas mengindikasikan bahwa data time-series yang kita observasi memiliki pola Multi-Seasonal. Berdasarkan pola seasonalnya menunjukan bahwa data ini memiliki seasonal harian dan mingguan, sehingga mari kita mebuat data multiple seasonal time series dan melakukan decompose kembali.

sxk9s_msts_weekly <- sxk9s %>% .$demand %>%  
  msts(seasonal.periods = c(24, 24*7))

plot_sxk9s_msts_weekly <- sxk9s_msts_weekly %>% 
  tail(24*7*4) %>% 
   mstl() %>% 
   autoplot()+ 
   labs(
      title = "Decomposition on Daily and Weekly in Sub-Area sxk9s"
    )+
    theme_minimal()+
    theme(
          title = element_text(size=12),
          axis.title=element_text(size=9, face="bold"),
          axis.text.x=element_text(size=9),
          axis.text.y=element_text(size=9),
          legend.position = "none",
          panel.spacing = unit(4, "pt")
          ) 

ggplotly(plot_sxk9s_msts_weekly) %>% 
    layout(title = list(text = paste0('Decomposition on Daily & Weekly in Sub-Area sxk9s',
                                    '<br>',
                                    '<sup>',
                                    'Deompose multiple seasonal time-series object using seasonal.periods = c(24, 24*7) (Daily & Weekly)',
                                    '</sup>')))

Pada hasil decompose objek multiple seasonal time-series diatas, dapat dilihat pola tren terlihat lebih smooth dan tidak membentuk pola. Hal ini mengindikasikan bahwa kemungkinan besar keseluruhan seasonal sudah ditangkap. Untuk lebih jelasnya, berikut ini visualisasi seasonal perjam-nya dalam harian dan mingguan pada data observasi kita:

sxk9s_msts_weekly_dc <- mstl(sxk9s_msts_weekly %>% tail(24*7*4))

plot_seasonal_msts <- as.data.frame(sxk9s_msts_weekly_dc) %>% 
  mutate(
    datetime = sxk9s %>% tail(24*7*4) %>% .$datetime
  ) %>% 
  mutate(
    day = wday(datetime, label = TRUE, abbr = FALSE),
    hour = as.factor(hour(datetime))
  ) %>% 
  group_by(day, hour) %>% 
  summarise(
    seasonal = sum(Seasonal24 + Seasonal168)
  ) %>% 
  ggplot(mapping = aes(x = hour, y = seasonal)) +
  geom_col(aes(fill = day)) +
  scale_fill_viridis_d(option = "VIRIDIS") +
  labs(
    title = "Seasonal Analysis based Daily and Weekly in Sub-Area sxk9s",
    fill = "Day of Week"
  ) + 
  theme_minimal() +
  theme(
        title = element_text(size=11),
        axis.title=element_text(size=9, face="bold"),
        axis.text.x=element_text(size=9),
        axis.text.y=element_text(size=9),
        legend.position = "right"
        ) 

ggplotly(plot_seasonal_msts)

Pada chart diatas dapat dilihat bahwa high season terjadi pada setiap pukul 17:00 sampai 19:00 dan setiap hari Jumat.

---

6 Cross Validation

Dalam case ini kita akan mencoba membuat beberapa model, tentunya kita akan membagi seluruh data yang kita punya menjadi data train dan data test. Problem yang diselesaikan yaitu memprediksi demand perjam dalam 1 minggu sehingga saya akan membagi data test menggunakan data 1 minggu terakhir dan sisanya sebagai data train. Berikut interval waktu di data train dan data test:

scotty_data <- scotty_input

#data test 1 minggu
test_size <- 24*7

test_end <- max(scotty_data$datetime)
test_start <- test_end - hours(test_size) + hours(1)

train_end <- test_start - hours(1)
train_start <- min(scotty_data$datetime)

# get the interval of each samples
intrain <- interval(train_start, train_end)
intest <- interval(test_start, test_end)

data.frame("interval_data_train"=intrain,
              "interval_data_test"=intest)

#head(scotty_data)

Berikut ini visualisasi pembagian data train dan data test dalam deret waktu:

viz_scotty_data <- scotty_data %>% 
  mutate(sample = case_when(
        datetime %within% intrain ~ "train",
        datetime %within% intest ~ "test")) %>% 
  drop_na() %>% 
  mutate(sample = factor(sample, levels = c("train", "test")))

ggplotly(
 viz_scotty_data %>% 
  ggplot(aes(x = datetime, y = demand, colour = sample)) +
    geom_line() +
    labs(
      x = NULL, 
      y = NULL,
      title = "Demand by Sub-Area on Data Train Vs Data Test"
    ) +
    scale_x_datetime(date_breaks = "7 day", labels = date_format("%b %d"))+
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_minimal()+
    theme(
      legend.position = "bottom", 
      plot.title = element_text(hjust = 0.5),
      panel.spacing = unit(1, "lines"))+
 my_viridis_color()
  ) %>% 
  layout(
    xaxis = list(showticklabels = FALSE),
    legend=list(orientation = "h",
                   y = -0.1, x = 0.4)
  )

Dapat dilihat pada visualisasi di atas, dimana jumlah demand di 3 Sub-Area berbeda-beda dan pada waktu-waktu tertentu memiliki selisih jauh berbeda. Hal ini dapat mengindikasikan outlier, sehingga perlu dilakukan scaling data train pada masing-masing Sub-Area sebelum memulai pemodelan supaya tidak sensitif pada data outlier. Berikut prosesnya:

1. Convert data to wide format

# Spread data menjadi 3 Sub-Area 
scotty_data <- scotty_data %>%
  spread(src_sub_area, demand)

head(scotty_data)

2. Scalling data
Sebelumya saya sudah mencoba scalling menggunakan transformasi Log. Memang outlier-nya lebih minimum, namun hasil evaluasi errornya lebih besar dibandingkan jika scalling menggunakan transformasi square root. Sehingga, saya memutuskan scalling menggunakan transformasi square root. Berikut sampel hasil scalling:

# mean(scotty_data$sxk97)
# mean(scotty_data$sxk9e)
# mean(scotty_data$sxk9s)
# 
# sd(scotty_data$sxk97)
# sd(scotty_data$sxk9e)
# sd(scotty_data$sxk9s)
# 
# var(scotty_data$sxk97)
# var(scotty_data$sxk9e)
# var(scotty_data$sxk9s)
#
# boxplot(scotty_data$sxk97)
# boxplot(scotty_data$sxk9e)
# boxplot(scotty_data$sxk9s)

# recipes: square root transformation- sqrt, center, scale
scalling <- recipe(~ ., filter(scotty_data, datetime %within% intrain)) %>%
  step_sqrt(all_numeric()) %>%
  step_center(all_numeric()) %>%
  step_scale(all_numeric()) %>%
  prep()

# preview the bake results
scotty_data <- bake(scalling, scotty_data)

# log transform
#scotty_data$sxk97 <- log(scotty_data$sxk97+1)
#scotty_data$sxk9e <- log(scotty_data$sxk9e+1)
#scotty_data$sxk9s <- log(scotty_data$sxk9s+1)


head(scotty_data)

Data yang sudah di-scalling tentunya nanti perlu dikembalikan lagi menjadi nilai aslinya, maka dibuat fungsi revert_back:

revert_back <- function(vector, recipe, varname) {
  # store recipe values
  rec_center <- recipe$steps[[2]]$means[varname]
  rec_scale <- recipe$steps[[3]]$sds[varname]
  
  # convert back based on the recipe
  results <- (vector * rec_scale + rec_center) ^ 2
  
  # add additional adjustment if necessary
  results <- round(results)
  
  # return the results
  results
}


# Exp dari Log +1
# revert_back_log <- function(vector) {
#   # store recipe values
#   #rec_center <- recipe$steps[[2]]$means[varname]
#   #rec_scale <- recipe$steps[[3]]$sds[varname]
#
#   # convert back based on the recipe
#   #results <- (vector * rec_scale + rec_center) ^ 2

#   results <- exp(vector)-1
#   # add additional adjustment if necessary
#   results <- round(results)
#
#   # return the results
#   results
# }

3. Convert data to long format
Berikut ini merupakan struktur data kita setelah dilakukan cross validation train-test dan scalling:

variabel datetime merupakan data waktu transaksi. variabel sample merupakan data tipe factor yang terbagi train dan test, sehingga kita dapat mengidentifikasi data yang hendak digunakan dalam pemodelan dan evaluasi model. variabel src_sub_area merupakan data nama dari Sub-Area. variabel demand merupakan data demand.

scotty_data <- scotty_data %>% 
  gather(src_sub_area, demand, -datetime)

head(scotty_data)

7 Nested Modelling

Pada pemodelan ini kita akan menggunakan fungsi purr dan untuk menggunakan fungsi purr, maka kita perlu merubah dataset kita yang sebelumnya berbentuk table menjadi nested table, berikut lebih jelasnya:

# Split Train-Test Sets
scotty_data <- scotty_data %>% 
  mutate(sample = case_when(
        datetime %within% intrain ~ "train",
        datetime %within% intest ~ "test")) %>% 
  drop_na() %>% 
  mutate(sample = factor(sample, levels = c("train", "test")))

# nested table
scotty_data <- scotty_data %>%
  group_by(src_sub_area, sample) %>%
  nest(.key = "data") %>%
  spread(sample, data)

scotty_data

Data diatas merupakan format nested table, variabel src_sub_area merupakan nama sub-area, variabel train merupakan bentuk nested table yang isinya merupakan data datetime dan demand untuk data train, dan variabel test juga merupakan bentuk nested table yang isinya merupakan data datetime dan demand untuk data test.

7.1 Data Time-Series Object List

Berdasarkan proses decompose sebelumnya, kita sudah melihat bahwa data kita memiliki Mutiple Seasonal sehingga data kita akan lebih cocok apabila dikonversi menjadi objek msts dengan seasonal harian dan mingguan, namun dengan menggunakan konsep nested table kita juga dapat mengecek dengan mudah apabila data kita dikonversi ke objek ts. Hal tersebut akan sangat membantu kita untuk membdandingkan hasilnya. Sehingga, kita akan mencoba melakukan pemodelan menggunakan single seasonal time-series dengan frequency harian dan multiple seasonal time-series dengan seasonal harian dan mingguan pada masing-masing Sub-Area.

# Membuat list objek data time-series
# replace_na(lag(x = as.vector(x$demand), n = 7),0)
data_model <- list(
  # single seasonal time-series dengan frequency harian
  ts = function(x) ts(x$demand, frequency = 24),
  #ts = function(x) ts(replace_na(lag(x = as.vector(x$demand), n = 7),0), frequency = 24),
  
  # multiple seasonal time-series dengan seasonal harian dan mingguan
  msts = function(x) msts(x$demand, seasonal.periods = c(24, 24*7))
  #msts = function(x) msts(replace_na(lag(x = as.vector(x$demand), n = 7),0), seasonal.periods = c(24, 24*7))
 
)

# convert to nested table
data_model <- data_model %>%
  rep(
    length(unique(scotty_data$src_sub_area))
  ) %>%
  enframe("data_fun_name", "data_fun") %>%
  mutate(src_sub_area =
    sort(
      rep(unique(scotty_data$src_sub_area), length(unique(.$data_fun_name)))
         )
  )

data_model

Data diatas menunjukan data time-series pada masing-masing sub-area. mari kita satukan dengan nested table sebelumnya untuk mempermudah pemodelan.

# join time-series object to data
scotty_data <- scotty_data %>%
  left_join(data_model)

scotty_data

Data diatas merupakan nested table antara data train dan data test pada masing-masing sub-area yang digabungkan dengan objek time series yang sudah dibuat. Sampai dengan proses ini, berarti kita sudah memiliki objek single seasonal time-series dengan frequency harian dan multiple seasonal time-series dengan seasonal harian dan mingguan pada masing-masing Sub-Area. Mari kita siapkan model-model untuk masing-masing objek.

7.2 Time Series Model List

objek time-series yang kita buat memiliki efek seasonal, sehingga seharusnya hanya metode yang dapat mengatasi efek seasonal yang kita gunakan, namun mari kita coba beberapa model dan diterapkan pada setiap objek time-series yang sudah dibuat. Metode yang akan digunakan yaitu:

1. ETS
   
2. Holt Winter Exponential Smoothing
   
3. Seasonal ARIMA
   
4. STLM method=“ets”
   
5. STLM method=“arima”
   
6. TBATS
   

seluruh model diatas akan kita join/merge dengan nested table dari objek time-series yang sudah kita buat sebelumnya supaya dapat di generate dengan mudah. Berikut proses dan hasilnya:

#models list
models <- list(
  ets = function(x) ets(x),
  #Model Holt’s Winter Exponential Smoothing
  holt_winters = function(x) HoltWinters(x),
  # Seasonal ARIMA (SARIMA)
  auto_arima = function(x) {
    while (adf.test(x)$p.value > 0.05) {
       x = diff(x)
     }
    auto.arima(x, seasonal = TRUE)
  },
  #STLM ets 
  stlm_ets = function(x) stlm(x, method="ets"),
  #STLM arima
  stlm_arima = function(x) stlm(x, method="arima"),
  
  #TBATS
  tbats = function(x) tbats(x, use.box.cox = FALSE, 
                  use.trend = TRUE, 
                  use.damped.trend = NULL,
                  use.arma.errors = TRUE,
                  use.parallel = TRUE,
                  num.cores = 6)
)

#Nested Model by sub-area
models <- models %>%
  rep(length(unique(scotty_data$src_sub_area))) %>%
  enframe("model_name", "model") %>%
  mutate(src_sub_area =
      sort(rep(unique(scotty_data$src_sub_area), length(unique(.$model_name))))
  )


#join time-series object and model
 scotty_models <- scotty_data %>%
   left_join(models)
#    filter(
#      !(model_name == "ets" & data_fun_name == "msts"),
#      !(model_name == "auto_arima" & data_fun_name == "msts"),
#      !(model_name == "tbats" & data_fun_name == "ts"),
#      !(model_name == "tbats_nat" & data_fun_name == "ts")
#      )
#   #   !(model_name == "stlm_arima" & data_fun_name == "msts")
#   # )

scotty_models

7.3 Execute Nested Model Fitting

Chunk di bawah ini digunakan untuk menggenerate seluruh model yang dibuat. Chunk ini di-set eval=FALSE karena cukup memakan waktu. Hasil generate model disimpan dengan format file .rds.

scotty_ts_model <- scotty_models %>%
  mutate(
    params = map(train, ~ list(x = .x)),
    data = invoke_map(data_fun, params),
    params = map(data, ~ list(x = .x)),
    fitted = invoke_map(model, params)
  ) %>%
  select(-data, -params)

# 
# tes <- msts(scotty_models$train[[1]]$demand, seasonal.periods = c(24,24*7))
# length(tes)
# 
# msts(replace_na(lag(x = as.vector(scotty_models$train[[1]]$demand), n = 7),0),seasonal.periods = c(24,24*7))
# 
# data.frame(
#   x = c(tes),
#   lag1 = lag(x = as.vector(tes), n = 1),
#   lag2 = lag(x = as.vector(tes), n = 2),
#   lag3 = lag(x = as.vector(tes), n = 3),
#   lag4 = lag(x = as.vector(tes), n = 4),
#   lag5 = lag(x = as.vector(tes), n = 5),
#   lag6 = lag(x = as.vector(tes), n = 6),
#   lag7 = lag(x = as.vector(tes), n = 7),
#   lag8 = lag(x = as.vector(tes), n = 8),
#   lag9 = lag(x = as.vector(tes), n = 9)
# ) %>%  na.omit() %>%
#   GGally::ggcorr(label = T)
# 
# replace_na(lag(x = as.vector(scotty_models$train[[1]]$demand), n = 7),0)

models <- scotty_ts_model
wd <-  as.character(getwd())
saveRDS(object=models, file=paste(paste(wd,"/scotty_ts_model/",sep = ""),"scotty_ts_model.rds",sep=""))

Mari kita breakdown dulu. Nested table yang sebelumnya kita buat terdiri dari 3 Sub-Area, dan setiap sub-area dibuat menjadi 2 objek time series yaitu single seasonal dan multiple seasonal, kemudian kita membuat 6 model yang akan digunakan berdasarkan sub-area dan jenis objek time-series nya. Sehingga, totalnya kita sudah membuat 36 model time-series dimana 18 model untuk single seasonal object dan 18 model untuk multiple seasonal object. Berikut ini hasil generate seluruh model yang sudah dibuat dan siap untuk digunakan.

scotty_ts_model <- readRDS("scotty_ts_model/scotty_ts_model.rds")
head(scotty_ts_model)

7.4 Model Evaluation

Berikut ini adalah hasil evaluasi seluruh model berdasarkan nilai Mean Absolute Error (MAE). Jika dilihat secara manual, nilai MAE paling kecil pada setiap sub-area dihasilkan oleh model TBATS menggunakan objek Multiple Seasonal Time Series. Berikut ini data MAE jika kita melakukan forecasting terhadap data train dan juga data test yang sudah kita siapkan.

scotty_model_eval <- scotty_ts_model %>% 
  mutate(
    #evaluasi model terhadap data train
    mae_train = map(fitted, ~ forecast(.x, h = test_size)) %>%
                      map2_dbl(train, ~ mae_vec(truth = revert_back(.y$demand,scalling,src_sub_area), 
                                                estimate = revert_back(.x$fitted,scalling,src_sub_area))),
    #evaluasi model terhadap data test
     mae_test = map(fitted, ~ forecast(.x, h = test_size)) %>%
                      map2_dbl(test, ~ mae_vec(truth = revert_back(.y$demand,scalling,src_sub_area), 
                                               estimate = revert_back(.x$mean,scalling,src_sub_area)))
    ) %>% 
  arrange(src_sub_area, mae_test, mae_train) 


scotty_model_eval %>% select(src_sub_area,data_fun_name, model_name,fitted,mae_train,mae_test)

7.5 Forecasting on Data Test

Berikut ini adalah proses dan hasil forecasting dari seluruh model terhadap data test:

scotty_forecast <- scotty_model_eval %>%
  mutate(
    forecast =
      map(fitted, ~ forecast(.x, h = test_size)) %>%
      map2(test, ~ tibble(
        datetime = .y$datetime,
        demand = as.vector(.x$mean)
      )),
    key = paste(data_fun_name, model_name, sep = "-")
  )

# Extract Value using unnest
scotty_forecast <- scotty_forecast %>% 
  select(src_sub_area, key, actual = test, forecast) %>%
  spread(key, forecast) %>%
  gather(key, value, -src_sub_area) %>% 
  unnest(value) %>% 
  mutate(demand = revert_back(demand,scalling,src_sub_area)) %>% 
  mutate(
    popup = glue(
      "Data Model: {key}
      Sub Area: {src_sub_area}
      Datetime: {datetime}
      Demand: {demand}"
    )
  ) 
  
scotty_forecast %>%
 select(-popup) %>% 
 spread(key,demand) 

Untuk lebih jelasnya, berikut visualisasi terhadap hasil forecasting diatas:

plot_scotty_forecast <- scotty_forecast %>% 
  ggplot(aes(x = datetime, y = demand)) +
    #geom_line(data = scotty_input,aes(y = demand),alpha = 0.2,size = 0.8)+
    geom_line(data = scotty_forecast %>% filter(key == "actual"),aes(y = demand),alpha = 0.2,size = 0.8)+
    geom_point(data = scotty_forecast %>% filter(key == "actual"),aes(y = demand, text=popup),alpha = 0.1,size=0.1)+
    geom_line(data = scotty_forecast %>% filter(key != "actual"),aes(frame = key,col = key)) +
    geom_point(data = scotty_forecast %>% filter(key != "actual"),
               aes(frame = key,col = key, text=popup),alpha = 0.1, size=0.1)+
    geom_point(data = scotty_forecast %>% filter(key != "actual",demand==0),
               aes(frame = key,col = "black", text=popup),alpa=0.5, size=1)+
    labs(
      x = "Date", 
      y = "Demand",
      title = "Comparison of Time Series Forecasting Models Result on Data Test", 
      frame = "Models"
    ) +
    scale_x_datetime(date_breaks = "1 day", labels = date_format("%b %d"))+
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    # scale_color_manual(values = c(color,color,color,color,color,color,color,color,color)) +
    theme_minimal()+
    theme(legend.position = "none", 
      plot.title = element_text(hjust = 0.5),
      panel.spacing = unit(1, "lines"))+
    scale_color_viridis_d()

ggplotly(plot_scotty_forecast, tooltip = "text", height = 560) %>% 
  animation_opts(
    2000, eas = "elastic", redraw = TRUE
  ) %>% 
  animation_slider(
    currentvalue = list(prefix = "Model: ", font = list(color=mycolor_hex("choc")))
  )

Pada chart diatas, garis abu-abu menunjukan data test yang hendak di-forecast, garis berwarna lainnya merupakan representasi hasil forecast dari model yang sudah kita buat dan titik/poin merepresentasikan hasil forecast dengan demand = 0 (waktu tanpa transaksi). Jika dilihat secara Visual, model TBATS (msts-tbats) menggunakan objek Multiple Seasonal Time Series terlihat cenderung paling bisa mengikuti pola garis dari data test pada setiap sub-area. Namun mari kita cek model terbaik dari setiap sub-area.

7.6 Automated Best Model Selection

Berikut ini adalah model terbaik yang dapat diterapkan pada ketiga sub-area:

# Minimum error for each sub-area
scotty_best_model <- scotty_model_eval %>%
  select(-fitted) %>% # remove unused
  group_by(src_sub_area) %>%
  filter(mae_test == min(mae_test)) %>%
  ungroup()


scotty_best_model %>% select(src_sub_area,data_fun_name, model_name,mae_train,mae_test)

Informasi diatas menunjukan bahwa error paling kecil pada ketiga sub-area dihasilkan oleh model TBATS menggunakan objek Multiple Seasonal Time Series dengan nilai MAE < 10, sehingga model inilah yang paling tepat untuk kita terapkan. Namun, perlu dingat sebelumnya kita membuat model ini menggunakan data train dan melakukan splitting menjadi tranning set dan test set untuk menentukan model terbaik. Untuk melakukan forecast terhadap data aktual maka kita harus menggabungkannya kembali seluruh datanya supaya deret waktu lengkap dan baru bisa menerapkan ketiga model terbaik diatas. Maka, berikut prosesnya:

# Menggabungkan data test dan train
scotty_final_model <- scotty_best_model %>%
  mutate(
    fulldata = map2(train, test, ~ bind_rows(.x, .y))) %>%
  select(src_sub_area, fulldata, everything(), -train, -test)

Chunk di bawah ini digunakan untuk menggenerate final model yang dibuat. Chunk ini di-set eval=FALSE karena cukup memakan waktu. Hasil generate model disimpan dengan format file .rds.

#Execute Nested Final Model
scotty_final_model <- scotty_final_model %>%
  mutate(
    params = map(fulldata, ~ list(x = .x)),
    data = invoke_map(data_fun, params),
    params = map(data, ~ list(x = .x)),
    fitted = invoke_map(model, params)
  ) %>%
  select(-data, -params)

final_models <- scotty_final_model
wd <-  as.character(getwd())
saveRDS(object=final_models, file=paste(paste(wd,"/scotty_ts_model/",sep = ""),"scotty_final_model.rds",sep=""))

---

8 Forecasting Result

Berdasarkan problem yang ingin diselesaikan maka kita perlu akan melakukan prediksi demand per-jam di sub area sxk97, sxk9e dan sxk9 dalam rentang waktu 1 minggu kedepan atau dari 2017-12-03 00:00 sampai 2017-12-09 23:00. Berikut ini hasil forecastingnya:

# Read Model
scotty_final_model <- readRDS("scotty_ts_model/scotty_final_model.rds")

# Forecasting 
time_size <- 24*7

next_week_prediction <- scotty_final_model %>% 
  mutate(forecast =
    map(fitted, ~ forecast(.x, h = time_size)) %>%
    map2(fulldata, ~ tibble(
      datetime = timetk::tk_make_future_timeseries(.y$datetime, time_size),
      demand = as.vector(.x$mean)
    ))
  )

# Extract Forecasting Result using unnest
next_week_prediction_result <- next_week_prediction %>% 
  select(src_sub_area, actual = fulldata, forecast) %>%
  gather(key, value, -src_sub_area) %>% 
  unnest(value) %>% 
  mutate(demand = revert_back(demand,scalling,src_sub_area)) 

# Display Forecasting Result  
next_week_prediction_result %>%
 filter(key=="forecast") %>% 
 spread(src_sub_area,demand) 

# Submit result
# result_submit <- next_week_prediction_result %>%
#     filter(key == "forecast") %>%
#     select(-key)
# 
# wd <-  as.character(getwd())
# filename <- paste0("gehm_scotty_",as.integer(Sys.time()),".csv")
# write_csv(result_submit,paste(paste(wd,"/data_result/",sep = ""),filename,sep=""))

Berikut visualisasi forecast dari 2017-12-03 00:00 sampai 2017-12-09 23:00

vis_forecast <- next_week_prediction_result %>% 
  filter(month(datetime)>=11) %>% 
  mutate(
    popup = glue(
      "Data Model: {key}
      Sub Area: {src_sub_area}
      Datetime: {datetime}
      Demand: {demand}"
    )
  ) 

plot_viz_forecast <- vis_forecast %>% 
  ggplot(aes(x = datetime, y = demand, colour = key)) +
    geom_line() +
   geom_point(aes(text=popup),size=0.05)+
   geom_point(data = vis_forecast %>% filter(demand==0),
             aes(x=datetime, y=demand, text=popup),size=0.5, alpha=0.5, show.legend = FALSE)+
    labs(
      title = "Forecasting Result in Sub-Area sxk97, sxk9e and sxk9"
    ) +
    scale_x_datetime(date_breaks = "7 day", labels = date_format("%b %d"))+
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_minimal()+
    theme(
      legend.position = "bottom", 
      plot.title = element_text(hjust = 0.5),
      panel.spacing = unit(1, "lines"))+
    my_viridis_color()


ggplotly(
 plot_viz_forecast, tooltip = "text", height=560
  ) %>% 
  layout(
    legend=list(orientation = "h",
                   y = -0.1, x = 0.4)
  )%>% 
  layout(title = list(text = paste0('Hourly Demand Forecasting Result in Sub-Area sxk97, sxk9e and sxk9',
                                    '<br>',
                                    '<sup>',
                                    'Forecasting result from December 3,2017 00:00 to December 09, 2017 23:00',
                                    '</sup>')))

Berdasarkan hasil forecasting untuk 3 Desember 2017 pukul 00:00 sampai 9 Desember 2017 pukul 23:00, demand paling tinggi pada ketiga sub-area terjadi pada 8 Desember 2017 pukul 18:00, dimana demand di sxk87 sebanyak 85, demand di sxk9e sebanyak 106 dan demand di sxk9s sebanyak 91.

---

9 Conclusion

9.1 Assumption Checking

Assumption Checking : Auto Correlation and Normality of Residual

Jika dilihat pada plot ACF residual ketiga sub-area, masih terdapat lag yang keluar dari garis biru. Hal ini dapat disimpulkan masih ada residual yang berkorelasi sehingga menunjukan bahwa masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast. Kemudian, jika dilihat pada histogram residual pada ketiga sub-area, menunjukan bahwa residual kurang berdisitribusi normal. Memang residual cukup berpusat ditengah, namun masih terdapat data residual pada bagian kiri dan kanan yang cukup banyak. Terkait hal ini kita bisa melakukan improvement misalkan memilih train-test set dengan periode waktu berbeda, differencing, imputasi null value menggunakan nilai rata-rata, melakukan teknik scalling berbeda atau mengeksplorasi teknik-teknik pemodelan time-series forecasting lainnya.

9.1.1 Sub-Area sxk97

Residual from TBATS method applied for Sub-Area sxk97

# checkresiduals(scotty_final_model$fitted[[1]],plot=TRUE, main="Residual from TBATS in Sub-Area sxk97") 
# shapiro.test(x = residuals(scotty_final_model$fitted[[1]]))#sxk97

#par(mfrow=c(1,2), cex=1.1)
ggplotly(ggAcf(residuals(next_week_prediction$fitted[[2]]), lag.max = 504)+
  ggtitle("sxk97 : ACF of residuals"))

ggplotly(gghistogram(residuals(next_week_prediction$fitted[[2]]))+
           labs(title="sxk97 : Histogram of residuals",x="residuals"))

#vis_forecast[c(10,24,168,336,504),]

9.1.2 Sub-Area sxk9e

Residual from TBATS method applied for Sub-Area sxk9e

# checkresiduals(scotty_final_model$fitted[[2]],plot=TRUE, main="Residual from TBATS in Sub-Area sxk9e") 
# shapiro.test(x = residuals(scotty_final_model$fitted[[2]]))#sxk97
#par(mfrow=c(1,2), cex=1.1)
ggplotly(ggAcf(residuals(next_week_prediction$fitted[[2]]), lag.max = 504)+
  ggtitle("sxk9e : ACF of residuals"))

ggplotly(gghistogram(residuals(next_week_prediction$fitted[[2]]))+
           labs(title="sxk9e : Histogram of residuals",x="residuals") 
         )

9.1.3 Sub-Area sxk9s

Residual from TBATS method applied for Sub-Area sxk9s

#checkresiduals(scotty_final_model$fitted[[3]],plot=TRUE,df=NULL, main="Residual from TBATS in Sub-Area sxk9s") 
#shapiro.test(x = residuals(scotty_final_model$fitted[[3]]))#sxk97
par(mfrow=c(1,2), cex=1.1)
ggplotly(ggAcf(residuals(next_week_prediction$fitted[[3]]), lag.max = 504)+
  ggtitle("sxk9s : ACF of residuals"))

ggplotly(gghistogram(residuals(next_week_prediction$fitted[[3]]))+
      labs(title="sxk9s : Histogram of residuals",x="residuals")
      )

9.2 Summary

Hasil analisa data time-series yang sudah dilakukan menunjukan layanan Ride-sharing service Scotty di sub-area sxk97, sxk9e dan sxk9s memiliki seasonal harian dan mingguan yang dimana high season terjadi pada pada setiap pukul 17:00 sampai 19:00 dan setiap hari Jumat. Selain itu, dari hasil pemodelan yang menggunakan metode ETS, HoltWinter, SARIMA, STLM dan TBATS yang sudah di-uji, model TBATS menggunakan objek Multiple Seasonal Time Series menghasilkan nilai error paling kecil dengan nilai MAE < 10. Jika dilakukan forecasting menggunakan model TBATS untuk 3 Desember 2017 pukul 00:00 sampai 9 Desember 2017 pukul 23:00, hasilnya menunjukan demand paling tinggi pada ketiga sub-area terjadi pada 8 Desember 2017 pukul 18:00, dimana demand di sxk87 sebanyak 85, demand di sxk9e sebanyak 106 dan demand di sxk9s sebanyak 91.

---

10 Reference

Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia.[Diakses 24-27 Maret 2020]

Herlambang, R. Dimas Baghas. 2019. PURRR-Operly Fitting Multiple Time Series [Diakses 24-27 Maret 2020]

Prabhakaran, S. 2016. Time Series Analysis[Diakses 25 Maret 2020]

Prastiwi, Ajeng. 2019. Forecasting Time Series with Multiple Seasonal[Diakses 26 Maret 2020]

Dancho, M. 2019. Forecasting Using Multiple Models[Diakses 26 Maret 2020]

Sievert, C. 2019. Interactive web-based data visualization with R, plotly, and shiny[Diakses 27 Maret 2020]