Food and Beverage Timeseries
Rangga Gemilang
1. PROBLEM IDENTIFICATION
Disini kita diminta untuk membuat report mengenai prediksi jumlah pengunjung untuk tiap jam dalam 7 hari kedepan (Senin, 19-Feb-2018 hingga Minggu, 25-Feb-2018) dan juga penjelesan mengenai seasonality yang terjadi berdasarkan data historis yang kita punya. Namun sebelum itu, mari kita lihat lebih jelas mengenai apa yang harus kita prediksi di permasalahan ini untuk memudahkan proses kita kedepannya.
# Read Data Test
fnb_test <- read_csv("data-test.csv")
glimpse(fnb_test)## Rows: 91
## Columns: 2
## $ datetime <dttm> 2018-02-19 10:00:00, 2018-02-19 11:00:00, 2018-02-19 12:00:0…
## $ visitor <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…Untuk memastikan rentang waktu yang diberikan, mari kita ubah
terlebih dahulu data datetime ke bentuk yang benar.
fnb_test <- fnb_test %>%
mutate(
datetime = floor_date(datetime, unit = "hour"),
visitor = as.integer(visitor)
)
summary(fnb_test)## datetime visitor
## Min. :2018-02-19 10:00:00 Min. : NA
## 1st Qu.:2018-02-20 19:30:00 1st Qu.: NA
## Median :2018-02-22 16:00:00 Median : NA
## Mean :2018-02-22 16:00:00 Mean :NaN
## 3rd Qu.:2018-02-24 12:30:00 3rd Qu.: NA
## Max. :2018-02-25 22:00:00 Max. : NA
## NA's :91Berdasarkan hasil summary() diatas, sudah dapat
dipastikan bahwa masalah yang perlu kita selesaikan adalah perdiksi
visitor tiap jam dalam rentang waktu 2018-02-19
10:00:00 sampai 2018-02-25 22:00:00 atau
13 jam * 1 minggu.
2. DATA PRE-PROCESSING
2.1 Read Data
fnb_train <- read_csv("data-train.csv")
glimpse(fnb_train)## Rows: 137,748
## Columns: 10
## $ transaction_date <dttm> 2017-12-01 13:32:46, 2017-12-01 13:32:46, 2017-12-01…
## $ receipt_number <chr> "A0026694", "A0026694", "A0026695", "A0026695", "A002…
## $ item_id <chr> "I10100139", "I10500037", "I10500044", "I10400009", "…
## $ item_group <chr> "noodle_dish", "drinks", "drinks", "side_dish", "drin…
## $ item_major_group <chr> "food", "beverages", "beverages", "food", "beverages"…
## $ quantity <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1,…
## $ price_usd <dbl> 7.33, 4.12, 2.02, 5.60, 3.01, 4.86, 6.34, 7.58, 4.12,…
## $ total_usd <dbl> 7.33, 4.12, 2.02, 5.60, 3.01, 4.86, 6.34, 7.58, 4.12,…
## $ payment_type <chr> "cash", "cash", "cash", "cash", "cash", "cash", "cash…
## $ sales_type <chr> "dine_in", "dine_in", "dine_in", "dine_in", "dine_in"…Melihat glimpse() dari data kita diatas, variabel
transaction_date merupakan variabel yang sama dengan
variabel datetime di data test kita tadi. Nanti akan kita
telaah lebih lanjut untuk pemrosesannya.
2.1.2 Data Descriptions
transaction_date: The timestamp of a transaction.receipt_number: The ID of a transaction.item_id: The ID of an item in a transaction.item_group: The group ID of an item in a transaction.item_major_group: The major-group ID of an item in a transaction.quantity: The quantity of purchased item.price_usd: The price of purchased item.total_usd: The total price of purchased item.payment_type: The payment method.sales_type: The sales method.
2.2 Data Wrangling
Jika kita lihat dari glimpse() data training kita tadi,
dapat dilihat bahwa tipe data transction_date nya belum
benar. Mengapa belum benar? Dikarenakan data tersebut belum berada di
satuan jam yang tepat yang dimana satuan
jam nya masih dalam bentukan yang detail. Satuan
jam yang kita inginkan adalah satuan jam yang
dapat menangkap berapa jumlah visitor di tiap jamnya. Maka dari
itu, langkah pertama yang harus kita lakukan adalah merapihkan
satuan jam kita dengan floor_date().
fnb <- fnb_train %>%
mutate(datetime = floor_date(transaction_date,
unit = "hours")) %>%
select(transaction_date,
receipt_number,
datetime)
paged_table(head(fnb))Sekarang kita sudah mempunyai kolom datetime dengan
satuan jam yang rapih, dan kita juga mengambil variable
receipt_number serta transaction_date nya
karena kita membutuhkan data receipt_number yang lengkap
agar nantinya bisa kita gunakan tiap receipt_number
yang berbeda sebagai penanda tiap visitor yang
berbeda. Untuk medapatkan jumlah visitor, kita
akan melakukan aggregasi atau summarise pada kolom
receipt_number kita.
fnb <- fnb %>%
group_by(datetime) %>%
summarise(visitor = n_distinct(receipt_number)) %>%
arrange(datetime)
paged_table(head(fnb))2.3 Missing Values & Padding
Setelah bentuk data kita sudah mirip dengan bentuk data test, kita akan cek missing values nya terlebih dahulu.
colSums(is.na(fnb))## datetime visitor
## 0 0Namun, ini bisa juga dikarenakan kita belum melakukan
padding untuk interval data kita. Yang dimaksud dengan
padding adalah kita menambahkan periode
waktu(datetime) berdasarkan kolom datetime
kita. Dan nanti datetime yang terlewat tersebut akan
diisikan nilai ( NA ) atau missing. Meskipun data yang
dibutuhkan hanyalah data pada pukul 10 pagi hingga 10
malam, kita akan tetap lakukan padding agar
tidak terlewat sedikitpun untuk interval waktu yang kita butuhkan. Nanti
setelah itu, baru kita akan melakukan filter() untuk
interval waktu yang dibutuhkan kembali.
Maka dari itu, mari kita lakukan padding terlebih
dahulu menggunakan fungsi pad().
# Padding
fnb <- fnb %>%
pad()
# Cek kembali missing values
colSums(is.na(fnb))## datetime visitor
## 0 663Dapat dilihat sekarang kolom visitor kita memiliki 663
data missing, untuk pengisian missing values tersebut
dapat bergantung pada perspektif bisnis kita masing-masing. Namun dalam
kasus ini, jika missing values diakibatkan oleh tutupnya toko,
maka kita akan mengisinya dengan nilai 0. Mengingat bahwa pembukaan toko
hanya berlangsung pada pukul 10 pagi hingga 10 malam.
Lalu pengisian missing values akan dilakukan dengan fungsi
na.fill().
fnb <- fnb %>%
mutate(visitor = na.fill(visitor, 0))
#Cek kembali missing values
colSums(is.na(fnb))## datetime visitor
## 0 0Setelah missing values sudah terisi, kita akan lanjut untuk memfilter data berdasarkan interval waktu yang kita butuhkan yaitu pada pukul 10:00:00 hingga 22:00:00.
fnb <- fnb %>%
filter(hour(datetime) >= 10 & hour(datetime) <= 22)
summary(fnb)## datetime visitor
## Min. :2017-12-01 13:00:00.00 Min. : 0.00
## 1st Qu.:2017-12-21 12:00:00.00 1st Qu.:16.00
## Median :2018-01-10 11:00:00.00 Median :26.00
## Mean :2018-01-10 06:45:25.17 Mean :28.99
## 3rd Qu.:2018-01-30 10:00:00.00 3rd Qu.:42.00
## Max. :2018-02-18 22:00:00.00 Max. :83.00Dilihat dari hasil summary() diatas, keseluruhan
dataset final kita dimulai dari 1 Desember,
2017 pukul 13:00 hingga 18 Februari,
2018 pukul 22:00.
2.4 Final Dataset
paged_table(fnb)3. EXPLORATORY DATA ANALYSIS
3.1 Basic Analysis
Sebelum kita menuju pada Seasonal Analysis, kita bisa melakukan Basic Analysis yang dimana analisis ini dilakukan pada data yang masih berbentuk mentah dan belum dijadikan menjadi objek timeseries. Maka dari itu kita akan mengambil kembali data mentah kita yang belum di padding.
fnb_raw <- read_csv("data-train.csv")
fnb_raw <- fnb_raw %>%
mutate(item_group = as.factor(item_group),
item_major_group = as.factor(item_major_group),
payment_type = as.factor(payment_type),
sales_type = as.factor(sales_type))Berikut adalah Basic Analysis yang bisa dilihat mengenai beberapa behaviour dasar dari Visitors kita:
# Set theme echarts untuk semua plot nantinya
e_common(font_family = "Cardo", theme = "infographic")
fnb_item <- fnb_raw %>%
select(item_group, quantity) %>%
group_by(item_group) %>%
summarise(Item = sum(quantity)) %>%
ungroup()
fnb_sales <- fnb_raw %>%
select(sales_type, receipt_number) %>%
group_by(sales_type, receipt_number) %>%
summarise(count = n()) %>%
ungroup()
fnb_sales <- fnb_sales %>%
select(sales_type) %>%
group_by(sales_type) %>%
summarise(Consumption = n())
fnb_pay <- fnb_raw %>%
select(payment_type, receipt_number) %>%
group_by(payment_type, receipt_number) %>%
summarise(count = n()) %>%
ungroup()
fnb_pay <- fnb_pay %>%
select(payment_type) %>%
group_by(payment_type) %>%
summarise(Payment = n())
item <- fnb_item %>%
e_charts(item_group) %>%
e_pie(Item,
itemStyle = list(
borderRadius = 0),
radius = c("25%", "70%"),
legend= FALSE) %>%
e_title("What Visitors Like to Buy?",
right = "10%",
left = "25%",
top = "1%") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)")
sales <- fnb_sales %>%
e_charts(sales_type) %>%
e_pie(Consumption,
itemStyle = list(
borderRadius = 1),
radius = c("25%", "70%"),
legend= FALSE) %>%
e_title("What Visitors Prefer?",
right = "10%",
left = "25%",
top = "1%") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)")
pay <- fnb_pay %>%
e_charts(payment_type) %>%
e_pie(Payment,
itemStyle = list(
borderRadius = 1),
radius = c("25%", "70%"),
legend= FALSE) %>%
e_title("How Visitors Pay?",
right = "10%",
left = "25%",
top = "1%") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)")
e_arrange(sales, item, pay, rows = 3)Dari plot diatas dapat disimpulkan bahwa:
Mayoritas Visitors lebih memilih untuk makan di tempat saat datang ke toko. Secara bisnis, ini juga akan mempengaruhi jumlah Visitors yang datang tiap jamnya, dikarenakan akan ada peristiwa dimana keadaan toko penuh dan tidak bisa menampung jumlah Visitors yang lebih banyak lagi. Namun, untuk mengetahui hal tersebut kita membutuhkan data yang lebih detail lagi mengenai kapasitas toko.
Makanan yang paling banyak terjual adalah makanan dengan tipe
rice dish&noodle dish. Sedangkan untuk minuman, yang paling banyak terjual adalah jenis minuman secara general selainKopi,Susu, danJus.Tipe pembayaran yang paling populer digunakan oleh Visitors adalah tipe pemabayaran
cash.
Setelah Basic Analysis dilakukan, kita akan melanjutkan Analysis berdasarkan waktu di tahap Seasonal Analysis selanjutnya.
3.2 Seasonality Analysis
Dalam proses Seasonal Analysis ini, kita ingin mengetahui mengenai trend dan seasonal yang terdapat dalam data kita. Untuk mencari informasi mengenai hal-hal tersebut, kita akan melakukan Decomposition. Ide utama dari Decomposition adalah untuk menguraikan ketiga komponen dari objek time series yaitu:
Trend : pola data secara general, cenderung untuk naik atau turun. Jika ada trend masih terdapat pola artinya masih ada pola yang belum terurai dengan baik.
Seasonal : pola musiman yang membentuk pola berulang pada periode waktu yang tetap.
Error/Reminder/Random : pola yang tidak dapat ditangkap dalam trend dan seasonal.
Jika pada hasil decompose, trend masih membetuk sebuah pola maka dapat dicurigai masih ada seasonality yang belum ditangkap. Seharusnya trend cenderung naik atau cenderung turun secara smooth. Maka dari itu, mari kita melakukan Decomposition dengan membuat objek timeseries terlebih dahulu.
# Data Preprocessing
# Membuat objek timeseries
fnb_ts <- ts(data = fnb$visitor,
frequency = 13)Jika kita bertanya mengapa frekuensinya 13? Seperti yang dijelaskan pada proses Problem Identification diatas, range waktu yang digunakan hanyalah range waktu pada saat toko buka hingga toko tutup yaitu pada pukul 10:00:00 hingga 22:00:00. Yang dimana range interval dari kedua waktu tersebut adalah 13 jika dihitung dalam satuan jam.
Setelah objek timeseries dibuat, kita akan memvisualisasi terlebih dahulu untuk melihat apakah data timeseries kita bersifat Additive atau Multiplicative.
df1 <- data.frame(date = as.Date(as.yearmon(time(fnb_ts))),
Y = as.matrix(fnb_ts))
df1 <- df1 %>%
mutate(date = fnb$datetime) %>%
rename(Visitor = Y)
df1 %>%
e_charts(x = date) %>%
e_line(Visitor,
symbol = "none",
legend = FALSE) %>%
e_title("Plot Historis u/ Jumlah Visitors per Jamnya ") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(date, axisPointer = list(show = TRUE)) %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
)Berdasarkan plot diatas, objek timeseries kita
memiliki varians dari polanya tetap atau konstan. Sehingga bisa kita
sebutkan bahwa data yang kita punya bersifat Additive dan Stasioner,
lalu kita akan memasukannya ke dalam fungsi decompose()
kita.
# Decompose objek TS
fnb_decomts <- fnb_ts %>%
decompose(type = "additive")
df_decom <- data.frame(date = fnb$datetime,
Visitor = fnb_decomts$x,
Trend = fnb_decomts$trend,
Seasonal = fnb_decomts$seasonal,
Remainder = fnb_decomts$random)
base_chart <- df_decom %>%
e_charts(x = date, height = 200) %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(date, axisPointer = list(
show = TRUE,
label = list(show = FALSE)
)) %>%
e_y_axis(
nameLocation = "end",
nameTextStyle = list(
color = "#666666",
fontWeight = "bold"
)
) %>%
e_group("decomp") %>%
e_grid(top = 35, bottom = 35)
y <- base_chart %>%
e_line(Visitor, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#e87c25")) %>%
e_axis_labels(y = "History")
trend <- base_chart %>%
e_line(Trend, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#fcce10")) %>%
e_axis_labels(y = "Trend")
seasonal <- base_chart %>%
e_line(Seasonal, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#27727b")) %>%
e_axis_labels(y = "Seasonality")
remainder <- base_chart %>%
e_bar(Remainder, legend = FALSE) %>%
e_datazoom(show = TRUE, toolbox = FALSE) %>%
e_color(c("#c1232b")) %>%
e_axis_labels(y = "Remainder") %>%
e_connect_group("decomp")
e_arrange(y, trend, seasonal, remainder, nrow = 4)Pada hasil decompose kita, kita mendapatkan informasi bahwa trend yang kita miliki masih berpola atau memiliki pattern. Oleh karena itu, bisa kita indikasikan bahwa data timeseries yang kita observasi memiliki pola Multi-Seasonal yang berarti masih ada seasonality yang belum berhasil kita tangkap.
Karena kita masih memiliki objek timeseries dalam satuan jam kita akan melakukan visualisasi terlebih dahulu agar trend pada satuan jam nya dapat dilihat.
fnb_dts <- fnb %>%
mutate(Hour = as.factor(hour(datetime)),
Seasonal = fnb_decomts$seasonal) %>%
distinct(Hour,
Seasonal) %>%
arrange(Hour) %>%
group_by(Hour)
fnb_dts %>%
e_charts(Hour) %>%
e_bar(Seasonal,
stack = T,
legend = F) %>%
e_title("Analisis Seasonal berdasarkan Hourly") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)",
formatter = htmlwidgets::JS("
function(params)
{
return `<strong style='color:Black;'>Time:</strong> ${params.value[0]}
<br/><strong style='color:Black;'>Seasonal:</strong> ${params.value[1]}`
} "))Dari plot diatas dapat disimpulkan bahwa:
Peak Season terjadi di tiap pukul 20:00.
Low Season terjadi di tiap pukul 10:00.
Trend Seasonal keseluruhannya terus meningkat tiap jamnya hingga mulainya pukul 21:00 yang dimana trend nya mulai menurun.
3.3 Multi-Seasonality Analysis
Untuk langkah Multi-Seasonal Analysis kali ini akan mirip
dengan yang kita lakukan pada Seasonal Analysis sebelumnya. Yang
berbeda hanya pada pembuatan objek timeseries yang
dimana kita akan menggunakan fungsi msts().
# Data Preprocessing
fnb_msts <- msts(data = fnb$visitor,
seasonal.periods = c(13,
13*7))Lalu kita lanjutkan dengan Decomposition dengan
menggunakan mstl() yang diperuntukkan untuk objek
timeseries yang multi-seasonal.
fnb_decomsts <- fnb_msts %>%
mstl()
df_decomsts <- as.data.frame(fnb_decomsts)
df_decomsts <- df_decomsts %>%
mutate(Date = fnb$datetime) %>%
rename(Visitor = Data)
base_chart <- df_decomsts %>%
e_charts(x = Date, height = 200) %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(
show = TRUE,
label = list(show = FALSE)
)) %>%
e_y_axis(
nameLocation = "end",
nameTextStyle = list(
color = "#666666",
fontWeight = "bold"
)
) %>%
e_group("decomp") %>%
e_grid(top = 35, bottom = 35)
y <- base_chart %>%
e_line(Visitor, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#b5c334")) %>%
e_axis_labels(y = "History")
trend <- base_chart %>%
e_line(Trend, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#e87c25")) %>%
e_axis_labels(y = "Trend")
seasonal13 <- base_chart %>%
e_line(Seasonal13, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#fcce10")) %>%
e_axis_labels(y = "Seasonal - Hourly")
seasonal91 <- base_chart %>%
e_line(Seasonal91, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#27727b")) %>%
e_axis_labels(y = "Seasonal - Weekly")
remainder <- base_chart %>%
e_bar(Remainder, legend = FALSE) %>%
e_datazoom(show = TRUE, toolbox = FALSE) %>%
e_color(c("#c1232b")) %>%
e_axis_labels(y = "Remainder") %>%
e_connect_group("decomp")
e_arrange(y, trend, seasonal13, seasonal91, remainder, nrow = 5)Bisa kita lihat bahwa data trend yang kita miliki sudah lebih smooth yang dimana ini menandakan bahwa pola-pola seasonality yang lainnya sudah tertangkap. Agar lebih jelasnya, kita sekarang memiliki Seasonal Mingguan yang telah berhasil ditangkap.
Seperti pada cara sebelumnya, kita akan membuat visualisasi kembali agar Seasonal Jam dan Seasonal Mingguan dapat lebih mudah dimengerti.
fnb_msts_df <- data.frame(fnb_decomsts)
fnb_msts_df <- fnb_msts_df %>%
mutate(datetime = fnb$datetime) %>%
mutate(Day = wday(datetime, label = TRUE, abbr = FALSE),
Hour = (hour(datetime))) %>%
group_by(Day, Hour) %>%
summarise(Seasonal = sum(Seasonal13 + Seasonal91)) %>%
ungroup() %>%
mutate(Day= as.factor(Day),
Hour= as.factor(Hour)) %>%
group_by(Day)
fnb_msts_df %>%
e_charts(Hour) %>%
e_bar(Seasonal,
stack = T) %>%
e_title("Analisis Seasonal berdasarkan Hourly & Weekly") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)",
formatter = htmlwidgets::JS("
function(params)
{
return `<strong style='color:Black;'>Day:</strong> ${params.seriesName}
<br/><strong style='color:Black;'>Time:</strong> ${params.value[0]}
<br/><strong style='color:Black;'>Seasonal:</strong> ${params.value[1]}`
} ")) %>%
e_legend(
orient = "horizontal",
bottom = "2%"
)Dari plot diatas dapat disimpulkan bahwa:
Peak Season terjadi di tiap pukul 20:00 pada hari Sabtu.
Low Season terjadi di tiap pukul 10:00 pada hari Jumat.
Trend Seasonal mulai meningkat pesat pada pukul 19:00 hingga pukul 22:00 dimana trend tersebut mulai menurun.
4. CROSS VALIDATION
Setelah Analysis timeseries dilakukan, kita akan melakukan Cross Validation pada objek timeseries kita. Tujuan utama dari Cross Validation itu sendiri adalah untuk proses evaluasi performa model kita nanti dalam memprediksi data baru.
Di tahap ini, kita akan membagi data kita menjadi 2 bagian yaitu data training dan data validation. Perbedaan tahap Cross Validation untuk data timeseries dengan data lainnya adalah kita tidak menyebarkan datanya kedalam dataset training dan validation secara random. Karena sangat tidak masuk akal jika kita ingin melakukan forecasting jumlah visitor yang datang di masa lalu berdasarkan data jumlah visitor di masa kedepannya.
Kita juga bisa melakukan K-Fold Cross
Validation pada data timeseries namun untuk saat ini, saya akan
menggunakan metode yang sangat basic diajarkan dikelas yaitu dengan
metode head() untuk data training dan
tail() untuk data validation.
Berdasarkan Problem Identification yang kita miliki, problem yang ingin diselesaikan adalah prediksi jumlah visitor tiap jamnya dalam interval waktu 1 minggu. Maka dari itu untuk data validation, kita akan menggunakan data 1 minggu terakhir dan sisanya sebagai data train.
fnb_val <- tail(fnb_msts,
13*7)
fnb_train <- head(fnb_msts,
length(fnb_msts) - 13*7) 5. MODEL FITTING
Dikarenakan objek timeseries kita memiliki indikasi Multiple Seasonality, sehingga kita akan menggunakan beberapa model Multi-Seasonal. Model-model tersebut adalah:
STLM ARIMA menggunakan objek MSTS.
STLM Holt-Winter menggunakan objek MSTS.
BATS Model menggunakan objek MSTS.
TBATS Model menggunakan objek MSTS.
SNAIVE Model menggunakan objek MSTS
5.1 STLM ARIMA
Merujuk pada artikel
ini, kita bisa menggunakan model STL(Seasonal and Trend
decomposition using Loess) untuk objek timeseries yang
memiliki pola Multi-Seasonal. STL secara konsep akan
melakukan smoothing terhadap data tetangga setiap masing-masing
observasi dengan memberikan bobot yang lebih berat terhadap data yang
dekat dengan observed data. Untuk proses modellingnya, kita bisa
menggunakan fungsi stlm().
Dikarenakan pada dasarnya model ARIMA diperuntukkan untuk objek
timeseries yang sifatnya non-seasonal, maka kita akan
menggunakan stlm(method = "arima") untuk pemodelan objek
timeseries kita yang multi-seasonal. Selain itu, jika
kita kembali pada Sub-Bab 4.2 Seasonality Analysis, data kita
juga telah bersifat Stasioner yang berarti sudah cocok
untuk digunakan dengan model ini.
# Model Fitting
fnb_arima <- stlm(fnb_train,
s.window = 13*7,
method = "arima")Setelah modelling dilakukan, kita akan melakukan
Forecasting dengan fungsi forecast().
# Forecasting
f_arima <- forecast(fnb_arima,
h = 13*7)Dilanjutkan dengan visualisasinya:
# Plotting
lower_arima <- as.data.frame(f_arima$lower)
upper_arima <- as.data.frame(f_arima$upper)
arima_forecast <- data.frame(Date = tail(fnb$datetime,
13*7),
Forecast = f_arima$mean,
Lower80 = lower_arima$"80%",
Upper80 = upper_arima$"80%",
Lower95 = lower_arima$"95%",
Upper95 = upper_arima$"95%")
arima_history <- data.frame(Date = head(fnb$datetime,
length(fnb_msts) - 13*7),
History = as.data.frame(fnb_train)
)
combined_arima <- bind_rows(
arima_history,
arima_forecast) %>%
rename(History = x)
combined_arima %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Forecast with SLTM-ARIMA") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan visualisasi perbandingan dari data aktual dan data forecast kita.
df_arima <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(f_arima$mean))
df_arima %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with STLM-ARIMA") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))5.2 STLM ETS
Untuk model exponential smoothing pada objek timeseries yang memiliki Seasonal lebih dari satu(Multiple Seasonal), maka kita bisa menggunakan model DHSW atau Double Seasonal Holt-Winters. Namun, dikarenakan data training kita memiliki nilai 0, maka model tersebut tidak bisa dijalankan dalam case ini.
Oleh karena itu, untuk pemodelan exponential smoothing pada objek
timeseries kita, akan kita ganti dengan salah satu
method STLM yaitu method="ets".
# # Model Fitting
fnb_ets <- stlm(fnb_train,
s.window = 13*7,
method = "ets")Setelah modelling dilakukan, kita akan melakukan
Forecasting dengan fungsi forecast().
# Forecasting
f_ets <- forecast(fnb_ets,
h = 13*7)Dilanjutkan dengan visualisasinya:
# Plotting
lower_ets <- as.data.frame(f_ets$lower)
upper_ets <- as.data.frame(f_ets$upper)
ets_forecast <- data.frame(Date = tail(fnb$datetime,
13*7),
Forecast = f_ets$mean,
Lower80 = lower_ets$"80%",
Upper80 = upper_ets$"80%",
Lower95 = lower_ets$"95%",
Upper95 = upper_ets$"95%")
ets_history <- data.frame(Date = head(fnb$datetime,
length(fnb_msts) - 13*7),
History = as.data.frame(fnb_train)
)
combined_ets <- bind_rows(
ets_history,
ets_forecast) %>%
rename(History = x)
combined_ets %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Forecast with STLM-ETS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan visualisasi perbandingan dari data aktual dan data forecast kita.
df_ets <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(f_ets$mean))
df_ets %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with STLM-ETS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))5.3 BATS
Model BATS memperlakukan objek timeseries yang seasonal sama seperti model ETS. Namun ada beberapa parameter tambahan yang dapat dipertimbangkan dalam BATS:
Box-cox Transformation.
ARMA errors.
Trend Components.
Seasonal Components.
Kelemahan dari model BATS itu sendiri adalah model tersebut tidak bisa memodelkan Seasonality yang kompleks. Yang dimaksud Seasonality kompleks adalah Seasonality yang tidak bersifat integer.
- Misalnya weekly time series dengan annual seasonal pattern 365.25/7 ≈ 52.18.
Maka dari itu, dikarenakan seasonal pattern kita bersifat integer, maka kita bisa menggunakan model ini untuk masuk dalam perbandingan nantinya.
# Model Fitting
fnb_bats <- fnb_train %>%
bats()Setelah modelling dilakukan, kita akan melakukan
Forecasting dengan fungsi forecast().
# Forecasting
f_bats <- forecast(fnb_bats,
h = 13*7)Dilanjutkan dengan visualisasinya:
# Plotting
lower_bats <- as.data.frame(f_bats$lower)
upper_bats <- as.data.frame(f_bats$upper)
bats_forecast <- data.frame(Date = tail(fnb$datetime,
13*7),
Forecast = f_bats$mean,
Lower80 = lower_bats$"80%",
Upper80 = upper_bats$"80%",
Lower95 = lower_bats$"95%",
Upper95 = upper_bats$"95%")
bats_history <- data.frame(Date = head(fnb$datetime,
length(fnb_msts) - 13*7),
History = as.data.frame(fnb_train)
)
combined_bats <- bind_rows(
bats_history,
bats_forecast) %>%
rename(History = x)
combined_bats %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Forecast with BATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan visualisasi perbandingan dari data aktual dan data forecast kita.
df_bats <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(f_bats$mean))
df_bats %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with BATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))5.4 TBATS
Perbedaan TBATS dan BATS ada pada Trigonometric functions yang terdapat pada model TBATS, kelebihan dari fungsi tersebut adalah memberikan fleksibilitas kepada model untuk memodelkan objek timeseries dengan Seasonality yang kompleks. Meski seasonality kita masih bersifat integer, kita akan tetap menggunakan model TBATS untuk perbandingan nantinya.
# Model Fitting
fnb_tbats <- fnb_train %>%
tbats()Setelah modelling dilakukan, kita akan melakukan
Forecasting dengan fungsi forecast().
# Forecasting
f_tbats <- forecast(fnb_tbats,
h = 13*7)Dilanjutkan dengan visualisasinya:
# Plotting
lower_tbats <- as.data.frame(f_tbats$lower)
upper_tbats <- as.data.frame(f_tbats$upper)
tbats_forecast <- data.frame(Date = tail(fnb$datetime,
13*7),
Forecast = f_tbats$mean,
Lower80 = lower_tbats$"80%",
Upper80 = upper_tbats$"80%",
Lower95 = lower_tbats$"95%",
Upper95 = upper_tbats$"95%")
tbats_history <- data.frame(Date = head(fnb$datetime,
length(fnb_msts) - 13*7),
History = as.data.frame(fnb_train)
)
combined_tbats <- bind_rows(
tbats_history,
tbats_forecast) %>%
rename(History = x)
combined_tbats %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Forecast with TBATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan visualisasi perbandingan dari data aktual dan data forecast kita.
df_tbats <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(f_tbats$mean))
df_tbats %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with TBATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))5.5 SNaive
Merunut pada artikel, kita juga bisa melakukan Seasonal Naive forecasting untuk objek timeseries kita yang Seasonalnya lebih dari 1.
# Model Fitting
fnb_snaive <- snaive(fnb_train,
h = 13*7)Setelah modelling dilakukan, kita akan melakukan
Forecasting dengan fungsi forecast().
# Forecasting
f_snaive <- forecast(fnb_snaive,
h = 13*7)Dilanjutkan dengan visualisasinya:
# Plotting
lower_snaive <- as.data.frame(f_snaive$lower)
upper_snaive <- as.data.frame(f_snaive$upper)
snaive_forecast <- data.frame(Date = tail(fnb$datetime,
13*7),
Forecast = f_snaive$mean,
Lower80 = lower_snaive$"80%",
Upper80 = upper_snaive$"80%",
Lower95 = lower_snaive$"95%",
Upper95 = upper_snaive$"95%")
snaive_history <- data.frame(Date = head(fnb$datetime,
length(fnb_msts) - 13*7),
History = as.data.frame(fnb_train)
)
combined_snaive <- bind_rows(
snaive_history,
snaive_forecast) %>%
rename(History = x)
combined_snaive %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Forecast with Snaive") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan visualisasi perbandingan dari data aktual dan data forecast kita.
df_snaive <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(f_snaive$mean))
df_snaive %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with Snaive") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))5B. PUTTING IT ALL TOGETHER
df_all <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
STLM_Arima = df_arima$Forecast,
STLM_ETS = df_ets$Forecast,
BATS = df_bats$Forecast,
TBATS = df_tbats$Forecast,
Snaive = df_snaive$Forecast)
df_all %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = F,
symbol = "none") %>%
e_line(STLM_Arima,
legend = F,
symbol = "none") %>%
e_line(STLM_ETS,
legend = F,
symbol = "none") %>%
e_line(BATS,
legend = F,
symbol = "none") %>%
e_line(TBATS,
legend = F,
symbol = "none") %>%
e_line(Snaive,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Perbandingan Forecast antar Model vs Actual") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))6. EVALUATION
Setelah Modelling dan Comparing dilakukan, kita akan melihat nilai error yang dihasilkan dari perbandingan yang dilakukan pada antar model. Ada beberapa jenis nilai error yang dapat dihasilkan dalam evaluasi ini, namun yang akan kita gunakan adalah nilai Mean Absolute Error (MAE) yang menunjukkan rata-rata dari nilai absolut error. MAE bisa diinterpretasikan sebagai seberapa besar penyimpangan hasil prediksi terhadap nilai aktualnya. Kelebihan lain dari MAE adalah lebih mudah dijelaskan ke orang tanpa latar belakang statistik.
6.1 STLM ARIMA
Untuk melakukan evaluasi model, dapat menggunakan fungsi
accuracy() dari package forecast. Yang
dilakukan adalah perhitungan nilai error yang dihasilkan dari
perbandingan selisih antara nilai forecast yang dihasilkan
dengan nilai aktual yang terdapat di data validasi
kita.
forecast::accuracy(f_arima$mean,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.143655 7.250587 5.624846 -Inf Inf 0.4100639 0Hasil diatas menunjukkan bahwa nilai MAE yang dihasilkan oleh model STLM ARIMA adalah sebesar 5.624846.
6.2 STLM ETS
Untuk melakukan evaluasi model, dapat menggunakan fungsi accuracy()
dari package forecast. Yang dilakukan adalah perhitungan
nilai error yang dihasilkan dari perbandingan selisih antara
nilai forecast yang dihasilkan dengan nilai aktual
yang terdapat di data validasi kita.
forecast::accuracy(f_ets$mean,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -7.767669 10.43721 8.846777 -Inf Inf 0.4196039 0Hasil diatas menunjukkan bahwa nilai MAE yang dihasilkan oleh model STLM ETS adalah sebesar 8.846777.
6.3 BATS
Untuk melakukan evaluasi model, dapat menggunakan fungsi accuracy()
dari package forecast. Yang dilakukan adalah perhitungan
nilai error yang dihasilkan dari perbandingan selisih antara
nilai forecast yang dihasilkan dengan nilai aktual
yang terdapat di data validasi kita.
forecast::accuracy(f_bats$mean,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -5.414 8.941681 7.230378 -Inf Inf 0.3903061 0Hasil diatas menunjukkan bahwa nilai MAE yang dihasilkan oleh model BATS adalah sebesar 7.230378.
6.4 TBATS
Untuk melakukan evaluasi model, dapat menggunakan fungsi accuracy()
dari package forecast. Yang dilakukan adalah perhitungan
nilai error yang dihasilkan dari perbandingan selisih antara
nilai forecast yang dihasilkan dengan nilai aktual
yang terdapat di data validasi kita.
forecast::accuracy(f_tbats$mean,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -3.783676 7.953826 6.533778 -Inf Inf 0.4875056 0Hasil diatas menunjukkan bahwa nilai MAE yang dihasilkan oleh model TBATS adalah sebesar 6.533778.
6.5 SNaive
Untuk melakukan evaluasi model, dapat menggunakan fungsi accuracy()
dari package forecast. Yang dilakukan adalah perhitungan
nilai error yang dihasilkan dari perbandingan selisih antara
nilai forecast yang dihasilkan dengan nilai aktual
yang terdapat di data validasi kita.
forecast::accuracy(f_snaive$mean,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -1.659341 9.388443 7.395604 -23.47607 44.4866 0.4183726 0Hasil diatas menunjukkan bahwa nilai MAE yang dihasilkan oleh model Snaive adalah sebesar 7.395604.
6.6 Kesimpulan
Dengan membandingkan semua nilai error yang dihasilkan, dapat disimpulkan bahwa model Multi-Seasonal ARIMA memiliki performa terbaik dengan MAE terkecil (5.624846). Oleh karena itu, kita akan menggunakan model ini untuk melakukan prediksi dengan data test kita.
7. PERFORMANCE WITH TEST DATA
Untuk melakukan prediction performance, kita bisa untuk memasukan hasil prediksi kita ke dalam data test yang telah disediakan.
fnb_test <- read_csv("data-test.csv")
fnb_test$visitor <- f_arima$mean
# Lihat apakah bentuk data sudah sesuai atau belum.
paged_table(fnb_test)Setelah data prediksi sudah dimasukkan, save data kembali ke bentuk csv agar bisa disubmit ke Algoritma Capstone Leaderboard.
write.csv(fnb_test,
file = "rangga_submission1.csv")Setelah data csv disubmit, kita bisa melihat apakah nilai MAE dari data prediksi kita sudah memenuhi syarat atau belum.
Berdasarkan dari gambar diatas, data prediksi kita mendapatkan nilai MAE sebesar 5.45 pada data test yang disediakan. Ini berarti dengan model Multiple-Seasonal ARIMA, kita telah memenuhi syarat penilaian:
8. ASSUMPTION
Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Mengapa menggunakan residual data? Karena dengan menggunakan residual data, kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model. Metode forecasting yang baik menghasilkan nilai residual berikut ini:
Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
Residual berdistribusi normal.
8.1 Auto-Correlation
H0: residual tidak berkorelasi.
H1: residual berkorelasi.
Nilai yang diinginkan p-value > 0.05 atau gagal tolak H0 yang berarti residual tidak berkorelasi.
Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time
series bisa melakukan uji Ljung-box dengan menggunakan
fungsi Box.test(residual model, type = "Ljung-Box).
arima_acf <- ggAcf(residuals(f_arima),
lag.max = 504) +
geom_segment(lineend = "butt",
color = "#c1232b") +
geom_hline(yintercept = 0,
color = "#c1232b") +
geom_hline(yintercept = c(0.0635, -0.0635),
color = "#27727b",
linetype = "dashed") +
theme_minimal() +
theme(
plot.background = element_rect(fill = "White",
color = "White"),
panel.background = element_rect(fill = "White"),
plot.title = element_text(colour = "#27727b",
family = "Cardo",
size = 15),
axis.title.x = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 13),
axis.title.y = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 13),
legend.background = element_rect(fill = "White",
color = "White"),
legend.key = element_rect(fill = "White",
color = "White"),
legend.title = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 11),
legend.text = element_text(color = "Black",
face = "bold",
family = "Cardo")
)
ggplotly(arima_acf)# menggunakan Ljung-Box test
Box.test(f_arima$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: f_arima$residuals
## X-squared = 0.00060238, df = 1, p-value = 0.9804Berdasarkan hasil diatas, nilai p-value yang didapatkan adalah sebesar 0.9804 yang dimana > 0.05. Oleh karena itu, dapat disimpulkan bahwa data forecast kita telah lolos uji asumsi auto-korelasi.
8.2 Normality
H0: residual menyebar normal.
H1: residual tidak menyebar normal.
Nilai yang diinginkan p-value > 0.05 atau gagal tolak H0 yang berarti residual menyebar normal.
Untuk mengecek normality residual pada hasil forecasting time series
kita bisa melakukan uji normality (shapiro test) dengan menggunakan
fungsi shapiro.test(residual model)
resarima <- data.frame(residual = f_arima$residuals)
resarima %>%
e_charts() %>%
e_histogram(residual,
name = "Error",
legend = F) %>%
e_density(residual,
areaStyle = list(opacity = .4,
color = "#27727b"),
itemStyle = list(color = "#27727b",
borderType = "dashed"),
y_index = 1,
legend = F) %>%
e_title(text = "Normality of Residuals STLM-ARIMA") %>%
e_hide_grid_lines() %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(axisLabel = list(fontFamily = "Cardo")) # menggunakan shapiro test
shapiro.test(f_arima$residuals)##
## Shapiro-Wilk normality test
##
## data: f_arima$residuals
## W = 0.98873, p-value = 0.000001169Berdasarkan hasil diatas, nilai p-value yang didapatkan adalah sebesar 0.000001169 yang dimana < 0.05. Oleh karena itu, dapat disimpulkan bahwa data forecast kita tidak lolos uji asumsi normality.
9. FIRST CONCLUSION
Dikarenakan model kita tidak lolos uji asumsi normality, kita dapat melakukan penambahan jumlah data untuk membangun model sehingga residual dapat berdistribusi normal. Merunut pada diskusi yang dilakukan oleh Pak Rob Hyndman, residual yang tidak normal dapat terjadi dikarenakan data observasi kita yang juga tidak berdistribusi normal. Cara untuk mengatasinya, berdasarkan petuah dari Mas Aksakal, beliau menyebutkan bahwa kita bisa melakukan:
Melakukan Log Transform pada data kita, sebelum di re-fit kembali ke dalam model. Setelah itu baru kita bisa uji normalitas kembali.
Menggunakan model yang tidak memerlukan asumsi normalitas. Namun saya belum bisa menemukan contoh model seperti ini untuk kasus timeseries.
Untuk cara nomor 1, tidak bisa kita lakukan dikarenakan ada data visitor bernilai 0, nilai tersebut akan berubah menjadi -Inf jika kita melakukan Log-Transform terhadap data. Akibat dari nilai -Inf tersebut adalah kita tidak bisa melakukan modelling untuk data kita. Berikut adalah contoh datanya:
log_train <- log(fnb_train)
head(log_train,100)## Multi-Seasonal Time Series:
## Start: 1 1
## Seasonal Periods: 13 91
## Data:
## [1] 2.7725887 3.6375862 3.2958369 3.3672958 3.7841896 3.9120230 4.1896547
## [8] 4.2484952 4.1431347 4.1431347 2.3025851 2.8332133 2.8903718 3.4657359
## [15] 3.0445224 3.6888795 3.5835189 3.5835189 3.7135721 4.2195077 4.1108739
## [22] 4.1271344 3.9889840 1.9459101 2.5649494 2.9957323 3.5553481 3.1354942
## [29] 3.4011974 3.4339872 3.7841896 4.0073332 4.1896547 3.8501476 3.8918203
## [36] 3.9889840 1.6094379 2.1972246 2.3025851 2.5649494 2.4849066 2.7080502
## [43] 2.9444390 3.2958369 2.9957323 3.6635616 3.9889840 4.0430513 3.9512437
## [50] 1.9459101 2.6390573 2.5649494 2.1972246 2.4849066 2.7725887 3.0910425
## [57] 3.1354942 3.5553481 4.2626799 4.0430513 3.9889840 3.5553481 1.3862944
## [64] 2.6390573 3.3322045 3.0910425 3.3672958 3.2188758 3.2580965 3.4011974
## [71] 3.5553481 3.8501476 3.8712010 4.0604430 3.4965076 0.6931472 2.3025851
## [78] 2.1972246 2.6390573 2.7725887 3.0445224 2.9957323 3.5263605 3.6888795
## [85] 4.0943446 4.0430513 3.9702919 3.7135721 -Inf -Inf -Inf
## [92] 2.5649494 2.9444390 3.1780538 3.3672958 3.2580965 3.8066625 3.7612001
## [99] 3.8286414 3.9889840Namun masih ada 1 cara lagi yang kita bisa lakukan yaitu melakukan scaling dengan nilai Square Root dari data kita. Kita akan mencoba cara tersebut dan melihat hasilnya apakah memenuhi syarat yang diberikan atau tidak.
10. DATA TRANSFORMATION
Dikarenakan model kita tidak lolos Uji Asumsi secara keseluruhan dapat disimpulkan bahwa hasil modeling kita belum cukup baik untuk menggambarkan dan menangkap informasi pada data. Maka dari itu, kita akan mencoba untuk melakukan Transformasi Data.
Pertama-tama, kita akan scaling data kita menggunakan nilai
square root-nya. Hal ini dapat dilakukan dengan menggunakan
fungsi dari recipe package. (Referensi bisa lihat disini)
Perbedaan dari Cross Validation yang pertama adalah, saat ini, kita melakukannya di awal. Kita hanya perlu data training yang sudah ditransformasi saja nantinya untuk proses Analysis timeseries. Kita tidak memerlukan data di interval waktu transformed data validation kita.
Mengapa begitu? Karena nanti pada tahap Evaluation kita akan unscale data forecast kita dan membandingkan akurasinya dengan data validation awal yang telah kita buat sebelumnya. Jika kita membandingkan transformed data forecast terhadap transformed data validation, nantinya Error yang dihasilkan akan sangat-sangat kecil terutama nilai MAE. (Mengingat semua transformed data akan mendekati 0 atau satuan yang sangat kecil)
# CROSS VALIDATION
# Buat objek penempatan baru agar yang awal tidak terubah-ubah
fnb_val2 <- tail(fnb,
13*7)
fnb_train2 <- head(fnb,
length(fnb) - 13*7)
# Transformasi Data
rec <- recipe(visitor~., data = fnb_train2) %>%
step_sqrt(all_outcomes()) %>%
step_center(all_outcomes()) %>%
step_scale(all_outcomes()) %>%
prep()Ubah Data Train dan Data Validation kita.
# Ubah Data Train
rec_train <- juice(rec)
# Ubah Data Validation yang hanya akan kita pakai untuk pembuatan plot nantinya.
rec_val <- bake(rec, fnb_val2)
paged_table(head(rec_train))Dapat dilihat data kita sudah ter-scaling. Dan kita akan membuat sebuah fungsi untuk mengembalikan nilai forecast kita nanti agar cocok dengan nilai yang ada di data validation kita.
rev_func <- function(x, rec) {
means <- rec$steps[[2]]$means[["visitor"]]
sds <- rec$steps[[3]]$sds[["visitor"]]
x <- (x * sds + means)^2
x <- round(x)
x <- ifelse(x < 0, 0, x)
x}11. ANOTHER EXPLORATORY DATA ANALYSIS
Setelah data sudah ter-scaling, dan berdasarkan pada
EDA pertama kita menunjukkan bahwa objek
timeseries kita adalah Multiple-Seasonal, maka kita
akan membuat objek baru dengan fungsi msts().
fnb_scale <- msts(rec_train$visitor,
seasonal.periods = c(13,
13*7))Lalu kita lanjutkan dengan Decomposition menggunakan
mstl().
fnb_decomsc <- fnb_scale %>%
mstl()
df_decomsc <- as.data.frame(fnb_decomsc)
df_decomsc <- df_decomsc %>%
mutate(Date = rec_train$datetime) %>%
rename(Visitor = Data)
base_chart <- df_decomsc %>%
e_charts(x = Date, height = 200) %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(
show = TRUE,
label = list(show = FALSE)
)) %>%
e_y_axis(
nameLocation = "end",
nameTextStyle = list(
color = "#666666",
fontWeight = "bold"
)
) %>%
e_group("decomp") %>%
e_grid(top = 35, bottom = 35)
y <- base_chart %>%
e_line(Visitor, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#b5c334")) %>%
e_axis_labels(y = "History")
trend <- base_chart %>%
e_line(Trend, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#e87c25")) %>%
e_axis_labels(y = "Trend")
seasonal13 <- base_chart %>%
e_line(Seasonal13, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#fcce10")) %>%
e_axis_labels(y = "Seasonal - Hourly")
seasonal91 <- base_chart %>%
e_line(Seasonal91, symbol = "none", legend = FALSE) %>%
e_datazoom(show = FALSE, toolbox = FALSE) %>%
e_color(c("#27727b")) %>%
e_axis_labels(y = "Seasonal - Weekly")
remainder <- base_chart %>%
e_bar(Remainder, legend = FALSE) %>%
e_datazoom(show = TRUE, toolbox = FALSE) %>%
e_color(c("#c1232b")) %>%
e_axis_labels(y = "Remainder") %>%
e_connect_group("decomp")
e_arrange(y, trend, seasonal13, seasonal91, remainder, nrow = 5)Dapat dipastikan lagi bahwa trend sudah tidak berpola dan para seasonality sudah tertangkap.
12. ANOTHER MODEL FITTING + EVALUATION
Sama seperti sebelumnya, kita akan melakukan modelling dengan 6 model:
STLM ARIMA menggunakan objek MSTS.
STLM Holt-Winter menggunakan objek MSTS.
BATS Model menggunakan objek MSTS.
TBATS Model menggunakan objek MSTS.
SNAIVE Model menggunakan objek MSTS.
12.1 STLM ARIMA
# Modelling
rec_arima <- stlm(fnb_scale,
s.window = 13*7,
method = "arima")Setelah modelling dilakukan, kita akan melakukan Forecasting dengan fungsi forecast().
# Forecasting
recf_arima <- forecast(rec_arima,
h = 13*7)
# Kembalikan Valuenya lagi ke bentuk semula dengan fungsi yang telah kita buat
# Untuk keperluan Evaluasi nanti.
rev_arima <- rev_func(recf_arima$mean,
rec = rec)Dilanjutkan dengan visualisasinya:
# Plotting
lower_arima <- as.data.frame(recf_arima$lower)
upper_arima <- as.data.frame(recf_arima$upper)
arima_forecast <- data.frame(Date = rec_val$datetime,
Forecast = recf_arima$mean,
Lower80 = lower_arima$"80%",
Upper80 = upper_arima$"80%",
Lower95 = lower_arima$"95%",
Upper95 = upper_arima$"95%")
arima_history <- data.frame(Date = rec_train$datetime,
History = as.data.frame(rec_train)
)
combined_arima <- bind_rows(
arima_history,
arima_forecast) %>%
rename(History = History.visitor)
combined_arima %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Transformed Data Forecast with STLM-ARIMA") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Visualisasi perbandingan dari data aktual dan data forecast kita.
df_arima <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(rev_arima))
df_arima %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with STLM-ARIMA") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan evaluasi untuk model kita:
forecast::accuracy(rev_arima, #Gunakan data yang telah diunscale
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -1.696629 7.196285 5.539326 -25.05003 35.59347 0.4054203 0Dapat dilihat dari hasil diatas bahwa nilai MAE yang kita hasilkan sebesar 5.539326 yang dimana lebih kecil dari model STLM-ARIMA pertama kita.
12.2 STLM ETS
# Modelling
rec_ets <- stlm(fnb_scale,
s.window = 13*7,
method = "ets")Setelah modelling dilakukan, kita akan melakukan Forecasting dengan fungsi forecast().
# Forecasting
recf_ets <- forecast(rec_ets,
h = 13*7)
# Kembalikan Valuenya lagi ke bentuk semula dengan fungsi yang telah kita buat
# Untuk keperluan Evaluasi nanti.
rev_ets <- rev_func(recf_ets$mean,
rec = rec)Dilanjutkan dengan visualisasinya:
# Plotting
lower_ets <- as.data.frame(recf_ets$lower)
upper_ets <- as.data.frame(recf_ets$upper)
ets_forecast <- data.frame(Date = rec_val$datetime,
Forecast = recf_ets$mean,
Lower80 = lower_ets$"80%",
Upper80 = upper_ets$"80%",
Lower95 = lower_ets$"95%",
Upper95 = upper_ets$"95%")
ets_history <- data.frame(Date = rec_train$datetime,
History = as.data.frame(rec_train)
)
combined_ets <- bind_rows(
ets_history,
ets_forecast) %>%
rename(History = History.visitor)
combined_ets %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Transformed Data Forecast with STLM-ETS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Visualisasi perbandingan dari data aktual dan data forecast kita.
df_ets <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(rev_ets))
df_ets %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with STLM-ETS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan evaluasi untuk model kita:
forecast::accuracy(rev_ets,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -5.213483 8.765356 6.988764 -43.26207 48.00071 0.4190816 0Dapat dilihat dari hasil diatas bahwa nilai MAE yang kita hasilkan sebesar 6.988764 yang dimana lebih kecil dari model STLM-ETS pertama kita.
12.3 BATS
# Modelling
rec_bats <- fnb_scale %>%
bats(use.box.cox = T,
use.trend = T,
use.damped.trend = T)Setelah modelling dilakukan, kita akan melakukan Forecasting dengan fungsi forecast().
# Forecasting
recf_bats <- forecast(rec_bats,
h = 13*7)
# Kembalikan Valuenya lagi ke bentuk semula dengan fungsi yang telah kita buat
# Untuk keperluan Evaluasi nanti.
rev_bats <- rev_func(recf_bats$mean,
rec = rec)Dilanjutkan dengan visualisasinya:
# Plotting
lower_bats <- as.data.frame(recf_bats$lower)
upper_bats <- as.data.frame(recf_bats$upper)
bats_forecast <- data.frame(Date = rec_val$datetime,
Forecast = recf_bats$mean,
Lower80 = lower_bats$"80%",
Upper80 = upper_bats$"80%",
Lower95 = lower_bats$"95%",
Upper95 = upper_bats$"95%")
bats_history <- data.frame(Date = rec_train$datetime,
History = as.data.frame(rec_train)
)
combined_bats <- bind_rows(
bats_history,
bats_forecast) %>%
rename(History = History.visitor)
combined_bats %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Transformed Data Forecast with BATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Visualisasi perbandingan dari data aktual dan data forecast kita.
df_bats <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(rev_bats))
df_bats %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with BATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan evaluasi untuk model kita:
forecast::accuracy(rev_bats,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -4.269663 8.361227 6.47191 -36.58648 42.60024 0.387254 0Dapat dilihat dari hasil diatas bahwa nilai MAE yang kita hasilkan sebesar 6.47191 yang dimana lebih kecil dari model BATS pertama kita.
12.4 TBATS
# Modelling
rec_tbats <- fnb_scale %>%
tbats(use.box.cox = T,
use.trend = T,
use.damped.trend = T)Setelah modelling dilakukan, kita akan melakukan Forecasting dengan fungsi forecast().
# Forecasting
recf_tbats <- forecast(rec_tbats,
h = 13*7)
# Kembalikan Valuenya lagi ke bentuk semula dengan fungsi yang telah kita buat
# Untuk keperluan Evaluasi nanti.
rev_tbats <- rev_func(recf_tbats$mean,
rec = rec)Dilanjutkan dengan visualisasinya:
# Plotting
lower_tbats <- as.data.frame(recf_tbats$lower)
upper_tbats <- as.data.frame(recf_tbats$upper)
tbats_forecast <- data.frame(Date = rec_val$datetime,
Forecast = recf_tbats$mean,
Lower80 = lower_tbats$"80%",
Upper80 = upper_tbats$"80%",
Lower95 = lower_tbats$"95%",
Upper95 = upper_tbats$"95%")
tbats_history <- data.frame(Date = rec_train$datetime,
History = as.data.frame(rec_train)
)
combined_tbats <- bind_rows(
tbats_history,
tbats_forecast) %>%
rename(History = History.visitor)
combined_tbats %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Transformed Data Forecast with TBATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Visualisasi perbandingan dari data aktual dan data forecast kita.
df_tbats <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(rev_tbats))
df_tbats %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with TBATS") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan evaluasi untuk model kita:
forecast::accuracy(rev_tbats,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.325843 7.437379 5.988764 -Inf Inf 0.492946 0Dapat dilihat dari hasil diatas bahwa nilai MAE yang kita hasilkan sebesar 5.988764 yang dimana lebih kecil dari model TBATS pertama kita.
12.5 SNaive
# Modelling
rec_snaive <- snaive(fnb_scale,
h = 13*7)Setelah modelling dilakukan, kita akan melakukan Forecasting dengan fungsi forecast().
# Forecasting
recf_snaive <- forecast(rec_snaive,
h = 13*7)
rev_snaive <- rev_func(recf_snaive$mean,
rec = rec)Dilanjutkan dengan visualisasinya:
# Plotting
lower_snaive <- as.data.frame(recf_snaive$lower)
upper_snaive <- as.data.frame(recf_snaive$upper)
snaive_forecast <- data.frame(Date = rec_val$datetime,
Forecast = recf_snaive$mean,
Lower80 = lower_snaive$"80%",
Upper80 = upper_snaive$"80%",
Lower95 = lower_snaive$"95%",
Upper95 = upper_snaive$"95%")
snaive_history <- data.frame(Date = rec_train$datetime,
History = as.data.frame(rec_train)
)
combined_snaive <- bind_rows(
snaive_history,
snaive_forecast) %>%
rename(History = History.visitor)
combined_snaive %>%
e_charts(x = Date) %>%
e_line(History,
legend = F,
symbol = "none") %>%
e_line(Forecast,
legend = F,
symbol = "none") %>%
e_line(Lower80,
legend = F,
symbol = "none") %>%
e_line(Upper80,
legend = F,
symbol = "none") %>%
e_line(Lower95,
legend = F,
symbol = "none") %>%
e_line(Upper95,
legend = F,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Transformed Data Forecast with SNaive") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Visualisasi perbandingan dari data aktual dan data forecast kita.
df_snaive <- data.frame(Date = fnb$datetime %>%
tail(13*7),
Actual = as.vector(fnb_val),
Forecast = as.vector(rev_snaive))
df_snaive %>%
e_charts(x = Date) %>%
e_line(Actual,
legend = T,
symbol = "none") %>%
e_line(Forecast,
legend = T,
symbol = "none") %>%
e_datazoom(
type = "slider",
toolbox = FALSE,
bottom = 10
) %>%
e_title("Actual vs Forecast with Snaive") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(Date, axisPointer = list(show = TRUE))Dan ditutup dengan evaluasi untuk model kita:
forecast::accuracy(rev_snaive,
fnb_val)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -1.797753 9.454801 7.460674 -24.93139 44.61182 0.4120454 0Dapat dilihat dari hasil diatas bahwa nilai MAE yang kita hasilkan sebesar 7.460674 yang dimana lebih besar dari model SNaive pertama kita.
13. FINAL ASSUMPTIONS
Berdasarkan dari hasil-hasil Evaluasi diatas, dapat dilihat bahwa proses transformasi data secara general menurunkan nilai MAE di semua model yang kita pakai. Namun hanya 2 model yang memiliki nilai MAE dibawah 6, yaitu model STLM ARIMA dan TBATS. Maka dari itu kita akan mencoba menggunakan 2 model ini untuk Uji Asumsi kali ini.
13.1 Auto-Correlation
H0: residual tidak berkorelasi.
H1: residual berkorelasi.
Nilai yang diinginkan p-value > 0.05 atau gagal tolak H0 yang berarti residual tidak berkorelasi.
Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time
series bisa melakukan uji Ljung-box dengan menggunakan
fungsi Box.test(residual model, type = "Ljung-Box).
13.1.1 STLM ARIMA
arima_acf <- ggAcf(residuals(recf_arima),
lag.max = 504) +
geom_segment(lineend = "butt",
color = "#c1232b") +
geom_hline(yintercept = 0,
color = "#c1232b") +
geom_hline(yintercept = c(0.0635, -0.0635),
color = "#27727b",
linetype = "dashed") +
theme_minimal() +
theme(
plot.background = element_rect(fill = "White",
color = "White"),
panel.background = element_rect(fill = "White"),
plot.title = element_text(colour = "#27727b",
family = "Cardo",
size = 15),
axis.title.x = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 13),
axis.title.y = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 13),
legend.background = element_rect(fill = "White",
color = "White"),
legend.key = element_rect(fill = "White",
color = "White"),
legend.title = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 11),
legend.text = element_text(color = "Black",
face = "bold",
family = "Cardo")
)
ggplotly(arima_acf)# menggunakan Ljung-Box test
Box.test(recf_arima$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: recf_arima$residuals
## X-squared = 0.00056028, df = 1, p-value = 0.9811Berdasarkan hasil diatas, nilai p-value yang didapatkan adalah sebesar 0.9811 yang dimana > 0.05. Oleh karena itu, dapat disimpulkan bahwa data forecast dengan model STLM ARIMA kita telah lolos uji asumsi auto-korelasi.
13.1.2 TBATS
tbats_acf <- ggAcf(residuals(recf_tbats),
lag.max = 504) +
geom_segment(lineend = "butt",
color = "#c1232b") +
geom_hline(yintercept = 0,
color = "#c1232b") +
geom_hline(yintercept = c(0.0635, -0.0635),
color = "#27727b",
linetype = "dashed") +
theme_minimal() +
theme(
plot.background = element_rect(fill = "White",
color = "White"),
panel.background = element_rect(fill = "White"),
plot.title = element_text(colour = "#27727b",
family = "Cardo",
size = 15),
axis.title.x = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 13),
axis.title.y = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 13),
legend.background = element_rect(fill = "White",
color = "White"),
legend.key = element_rect(fill = "White",
color = "White"),
legend.title = element_text(colour = "Black",
face = "bold",
family = "Cardo",
size = 11),
legend.text = element_text(color = "Black",
face = "bold",
family = "Cardo")
)
ggplotly(tbats_acf)# menggunakan Ljung-Box test
Box.test(recf_tbats$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: recf_tbats$residuals
## X-squared = 0.19209, df = 1, p-value = 0.6612Berdasarkan hasil diatas, nilai p-value yang didapatkan adalah sebesar 0.6612 yang dimana > 0.05. Oleh karena itu, dapat disimpulkan bahwa data forecast dengan model TBATS kita telah lolos uji asumsi auto-korelasi.
13.2 Normality
H0: residual menyebar normal.
H1: residual tidak menyebar normal.
Nilai yang diinginkan p-value > 0.05 atau gagal tolak H0 yang berarti residual menyebar normal.
Untuk mengecek normality residual pada hasil forecasting time series
kita bisa melakukan uji normality (shapiro test) dengan menggunakan
fungsi shapiro.test(residual model)
13.2.1 STLM ARIMA
resarima <- data.frame(residual = recf_arima$residuals)
resarima %>%
e_charts() %>%
e_histogram(residual,
name = "Error",
legend = F) %>%
e_density(residual,
areaStyle = list(opacity = .4,
color = "#27727b"),
itemStyle = list(color = "#27727b",
borderType = "dashed"),
y_index = 1,
legend = F) %>%
e_title(text = "Normality of Residuals STLM-ARIMA") %>%
e_hide_grid_lines() %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(axisLabel = list(fontFamily = "Cardo")) # menggunakan shapiro test
shapiro.test(recf_arima$residuals)##
## Shapiro-Wilk normality test
##
## data: recf_arima$residuals
## W = 0.99581, p-value = 0.01132Berdasarkan hasil diatas, nilai p-value yang didapatkan adalah sebesar 0.01132 yang dimana < 0.05. Oleh karena itu, dapat disimpulkan bahwa data forecast dengan model STLM ARIMA kita tidak lolos uji asumsi normality.
13.2.2 TBATS
restbats <- data.frame(residual = recf_tbats$residuals)
restbats %>%
e_charts() %>%
e_histogram(residual,
name = "Error",
legend = F) %>%
e_density(residual,
areaStyle = list(opacity = .4,
color = "#27727b"),
itemStyle = list(color = "#27727b",
borderType = "dashed"),
y_index = 1,
legend = F) %>%
e_title(text = "Normality of Residuals TBATS") %>%
e_hide_grid_lines() %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)") %>%
e_x_axis(axisLabel = list(fontFamily = "Cardo")) # menggunakan shapiro test
shapiro.test(recf_tbats$residuals)##
## Shapiro-Wilk normality test
##
## data: recf_tbats$residuals
## W = 0.99863, p-value = 0.6866Berdasarkan hasil diatas, nilai p-value yang didapatkan adalah sebesar 0.6866 yang dimana > 0.05. Oleh karena itu, dapat disimpulkan bahwa data forecast dengan model TBATS kita telah lolos uji asumsi normality.
14. ANOTHER PERFORMANCE WITH TEST DATA
Dikarenakan satu-satunya model yang lolos Uji Asumsi adalah model TBATS. Maka kita akan menggunakan model tersebut untuk prediction performance akhir kita. Untuk melakukan prediction performance, kita bisa untuk memasukan hasil prediksi kita ke dalam data test yang telah disediakan.
fnb_test <- read_csv("data-test.csv")
fnb_test$visitor <- rev_tbats
# Lihat apakah bentuk data sudah sesuai atau belum.
paged_table(fnb_test)Setelah data prediksi sudah dimasukkan, save data kembali ke bentuk csv agar bisa disubmit ke Algoritma Capstone Leaderboard.
write.csv(fnb_test,
file = "rangga_submission2.csv")Setelah data csv disubmit, kita bisa melihat apakah nilai MAE dari data prediksi kita sudah memenuhi syarat atau belum.
Berdasarkan dari gambar diatas, data prediksi kita mendapatkan nilai MAE sebesar 15.85 pada data test yang disediakan. Ini berarti dengan model TBATS akhir kita, kita tidak memenuhi salah satu syarat penilaian:
Dapat disimpulkan bahwa model tuning yang dibuat dengan menggunakan Transformasi Data memiliki indikasi Overfitting, yang dimana model tersebut tidak memiliki performa yang baik pada data test kita, namun memiliki performa yang baik pada data validation yang kita miliki.
15. FINAL CONCLUSION
Dari semua proses yang kita jalankan dalam case ini, kita memiliki 2 Final Model yang memiliki kekurangan dan kelebihan masing-masing.
Model STLM-ARIMA: Memiliki performa yang baik, namun tidak dapat lolos Uji Asumsi Normalitas.
Model TBATS w/ Transformed Data: Lolos semua Uji Asumsi, namun performanya cenderung overfit.
Merunut pada diskusi yang dilakukan oleh Pak Rob Hyndman, saya akan memilih model STLM ARIMA untuk submission prediksi akhir saya.
Setelah hasilnya sudah dipilih kita akan melakukan Exploratory Data Analysis kembali untuk hasil forecasting kita.
#Pilih hasil dari model STLM-ARIMA
rangga_forecast <- read_csv("rangga_submission1.csv")
rangga_forecast <- rangga_forecast %>%
select(datetime, visitor)Lalu kita bisa melakukan Analysis Seasonal pada hasil forecast kita.
fnb_decomts1 <- ts(data = rangga_forecast$visitor,
frequency = 13) %>%
decompose(type = "additive")
fnb_dts1 <- rangga_forecast %>%
mutate(Hour = as.factor(hour(datetime)),
Seasonal = fnb_decomts1$seasonal) %>%
distinct(Hour,
Seasonal) %>%
arrange(Hour) %>%
group_by(Hour)
fnb_dts1 %>%
e_charts(Hour) %>%
e_bar(Seasonal,
stack = T,
legend = F) %>%
e_title("Analisis Seasonal berdasarkan Hourly Forecast") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)",
formatter = htmlwidgets::JS("
function(params)
{
return `<strong style='color:Black;'>Time:</strong> ${params.value[0]}
<br/><strong style='color:Black;'>Seasonal:</strong> ${params.value[1]}`
} "))Berdasarkan plot diatas, dapat disimpulkan bahwa:
Peak Season diramalkan akan terjadi di pukul 20:00.
Low Season diramalkan akan terjadi di pukul 10:00.
Trend Seasonal keseluruhannya terus meningkat tiap jamnya hingga mulainya pukul 21:00 yang dimana trend nya mulai menurun.
Dikarenakan objek timeseries kita tidak memiliki Multiple Seasonality, maka kita akan melakukan analisis lanjutan mengenai Perkiraan-perkiraan yang akan terjadi di tiap jam dan harinya menggunakan plot dibawah ini. (Pertanyaannya mengacu pada proses Problem Identification diatas.)
rangga_forecast <- rangga_forecast %>%
mutate(Hour = hour(datetime),
Day = wday(rangga_forecast$datetime, label = TRUE, abbr = FALSE))
rangga_daynhour <- rangga_forecast %>%
select(Day, Hour, visitor) %>%
group_by(Day, Hour) %>%
summarise(visitor = sum(visitor)) %>%
ungroup() %>%
mutate(Day= as.factor(Day),
Hour= as.factor(Hour)) %>%
group_by(Day)
rangga_daynhour %>%
e_charts(Hour) %>%
e_bar(visitor,
stack = T) %>%
e_title("Visitors Forecast",
subtext = "periode Senin, 19 Februari hingga Minggu, 25 Februari 2018") %>%
e_tooltip(backgroundColor = "rgba(255,255,255,0.7)",
formatter = htmlwidgets::JS("
function(params)
{
return `<strong style='color:Black;'>Day:</strong> ${params.seriesName}
<br/><strong style='color:Black;'>Time:</strong> ${params.value[0]}
<br/><strong style='color:Black;'>Visitors:</strong> ${params.value[1]}`
} ")) %>%
e_legend(
orient = "horizontal",
bottom = "2%"
)Berdasarkan plot diatas, dapat disimpulkan bahwa:
Forecast waktu dengan jumlah Visitor terbanyak akan jatuh pada hari Sabtu, 24 Februari pada pukul 20:00 dengan perkiraan 66 Visitors hingga pukul 21:00 dengan perkiraan 62 Visitors di hari yang sama.
Hari dengan jumlah visitor terbanyak juga akan jatuh pada hari Sabtu, 24 Februari dengan perkiraan 460 Visitors dan diikuti dengan hari Minggu, 25 Februari dengan perkiraan 439 Visitors.
Forecast waktu dengan jumlah Visitor paling sedikit akan jatuh pada hari Jumat, 23 Februari pada pukul 10:00 hingga pukul 12:00 dengan perkiraan 1 Visitor pada tiap jamnya.
Hari dengan jumlah visitor paling sedikit juga akan jatuh pada hari Jumat, 23 Februari dengan perkiraan 339 Visitors dan diikuti dengan hari Kamis, 22 Februari dengan perkiraan 352 Visitors.