Latar Belakang & Tujuan

Sepeda adalah alat transportasi yang populer dan serba guna, digunakan oleh banyak orang untuk berbagai keperluan mulai dari alat transportasi menuju tempat kerja atau sekolah dan juga bisa digunakan untuk berolahraga.

Penyewaan sepeda mulai menjadi populer sekitar tahun 1960-an dan 1970-an. Pada saat itu, terutama di beberapa kota di Eropa dan Amerika Utara, muncul kebutuhan akan penyediaan sepeda sewa bagi wisatawan yang ingin menjelajahi kota-kota tersebut dengan cara yang lebih fleksibel dan independen. Konsep penyewaan sepeda ini pun kemudian berkembang dan menyebar ke berbagai negara dan kota di seluruh dunia. Dengan meningkatnya minat akan kegiatan bersepeda dan kesadaran akan manfaatnya, bisnis penyewaan sepeda semakin marak dan berkembang hingga saat ini.

Dimasa lalu penyewaan sepeda dilakukan dengan cara tradisional dimana penyewa harus mengambil dan mengembalikan sepeda ke tempat asal sepeda disewa. Dewasa ini sistem penyewaan sepeda sudah tidak harus mengembalikan ke tempat asal dan bisa dikembalikan ke pool/terminal dari penyedia jasa sewa. Hal ini dimungkinkan karena setiap sepeda yang disewakan sudah terintegrasi dengan pool/terminal dari penyedia jasa sewa, dan setiap aktivitas baik penyewaan maupun pengembalian sudah terekam oleh sistem dan tersimpan menjadi data set.

Hasil rekaman data tersebut biasanya meliputi informasi terkait akun penyewa dan waktu sepeda mulai disewa dan waktu sepeda dikembalikan. Biasanya akan dilakukan rekapan jumlah sepeda yang disewa baik dalam keterangan waktu harian, mingguan, bulanan, hingga tahunan.

Berdasarkan informasi terkait data penyewaan sepeda tersebut, akan dicoba untuk melakukan analisis dan pemodelan menggunakan metode Time Series & Forecasting.

Informasi Dataset

Data Bike sharing ini merupakan data yang digunakan dalam penelitian yang berjudul Event labeling combining ensemble detectors and background knowledge yang telah dilakukan oleh Fanaee-T, Hadi, dan Gama pada tahun 2014 dan telah dipublikasikan pada Jurnal Progress in Artificial Intelligence Published by Springer-Verlag.

Bike sharing dataset berisi informasi mengenai rekap jumlah penggunaan sepeda harian di Washington DC, Amerika Serikat dari tanggal 1 Januari 2011 s.d 31 Desember 2012 yang diperoleh dari sistem Capital bikehare dan diberikan informasi tambahan mengenai cuaca dan musiman yang diambil dari freemeteo.com.

Dataset ini diperoleh melalui kaggle Bike sharing dataset yang diposting oleh lakshmi25npathi. Dimana dari postingan tersebut tersedia dua versi dataset Bike sharing dataset dengan format csv sebagai berikut:

1. day.csv yang berisikan rekap jumlah penyewaan sepeda harian.
2. hour.csv yang berisikan rekap jumlah penyewaan sepeda perjam.

Pada analisis kali ini akan digunakan dataset yang merekap jumlah penyewaan sepeda harianday.csv.

Dataset ini memiliki informasi kolom sebagai berikut:

• instant : indeks data
• dteday : tanggal
• season : musim
  • 1: springer
  • 2: summer
  • 3: fall
  • 4: winter
• yr : tahun (0: 2011 dan 1: 2012)
• mnth : bulan ( 1 s.d 12)
• hr : jam (0 s.d 23)
• holiday : weather day is holiday or not based on Holiday schedule
  • 0: not holiday
  • 1: holiday
• weekday : daftar hari (0 to 6)
• workingday : if day is neither weekend nor holiday is 1, otherwise is 0.
  • 0: Weekend or holiday
  • 1: workingday
• weathersit :
  • 1: Clear, Few clouds, Partly cloudy, Partly cloudy
  • 2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist
  • 3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds
  • 4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog
• temp : Normalized temperature in Celsius. The values are divided to 41 (max)
• atemp: Normalized feeling temperature in Celsius. The values are divided to 50 (max)
• hum: Normalized humidity. The values are divided to 100 (max)
• windspeed: Normalized wind speed. The values are divided to 67 (max)
• casual: count of casual users
• registered: count of registered users
• cnt: count of total rental bike including both casual and registered

dari sekian banyak kolom yang tersedia, untuk keperluan analisis kali ini hanya akan digunakan kolom dteday dan cnt.

Konsep Dasar

Analisis dan prediksi kali ini akan menggunakan suatu metode yang disebut Time Series & Forecasting . Dimana Time series & Forecasting merupakan metode analisis dan pengolahan data yang bertujuan untuk memprediksi suatu nilai di masa mendatang berdasarkan data deret waktu.

Untuk melakukan Time series & Forecasting ada beberapa syarat yang harus terpenuhi, antara lain:

1. Data memiliki interval waktu yang lengkap
2. Data terurut berdasarkan waktu terdahulu hingga waktu terkini
3. Tidak terdapat data yang hilang

1. Data Preparation & EDA

Data preparation dilakukan dengan mengubah tipe data menjadi objek ts, hal ini perlu dilakukan karena dalam analisis kali ini akan digunakan packages forecast yang mensyaratkan data harus bertipe objek ts.

Sebuah objek ts dapat diuraikan menjadi 3 komponen utama yang akan dihitung untuk peramalan. antara lain sebagai berikut:

• tren (T): pergerakan rata-rata, secara global, sepanjang interval;
• musiman (S): pola yang ditangkap pada setiap interval musiman;
• kesalahan (E): pola / nilai yang tidak dapat ditangkap oleh tren dan musiman.

Hubungan antara 3 komponen objek ts tersebut dapat berupa aditif atau multiplikatif, yang dapat dilihat pada plot garis dari data asli:

• aditif: T + S + E

Gambar: Ilustrasi Objek ts Aditif

Berdasarkan gambar diatas pola aditif menunjukkan pola gambar puncak dan lembah yang tetap seiring waktu, hal ini mengindikasikan variansi data yang tetap.

• multiplikatif: T * S * E

  

  Gambar: Ilustrasi Objek ts Multiplikatif

  Berdasarkan gambar diatas pola aditif menunjukkan pola gambar puncak dan lembah yang bertumbuh seiring waktu, hal ini mengindikasikan variansi data yang tidak tetap.

2. Pemodelan

Ada beberapa cara peramalan dalam deret waktu:

2.1 Naive Model

Memprediksi nilai masa depan berdasarkan 1 nilai sebelumnya.

2.2 Simple Moving Average (SMA)

Memprediksi nilai titik data dengan menghitung rata-rata n jumlah data sebelum titik data. Metode ini lebih tepat untuk data tanpa tren dan musiman dan digunakan untuk menghaluskan kesalahan data. Kelemahan dari metode ini adalah bahwa ia memberikan bobot yang sama untuk setiap titik data, meskipun pada kenyataannya, titik data yang lebih baru harus memiliki bobot yang lebih tinggi daripada masa lalu.

2.3 Exponential Smoothing

Memprediksi nilai masa depan dengan menghaluskan kesalahan, tren, dan musiman, di mana ia memberikan bobot yang berbeda pada setiap titik data yang digunakan untuk prediksi. Exponential Smoothing dibagi menjadi tiga:

• simple exponential smoothing (SES): smoothing error
• double exponential smoothing (Holt): smoothing error & tren
• triple exponential smoothing (Holt Winters): smoothing error, tren, & musiman

2.4 ARIMA

ARIMA adalah singkatan dari Auto Regresive Integreted Moving Averange. ARIMA sangat bagus untuk memprediksi data deret waktu stasioner. Data deret waktu stasioner memiliki nilai yang bergerak di sekitar rata-ratanya (tidak ada tren atau musiman).

ARIMA memprediksi nilai masa depan dengan menggabungkan metode Auto Regresive dan Moving Averange, sedangkan terintegrasi berarti berapa kali perbedaan diterapkan pada data (cara untuk mengubah data non-stasioner menjadi data deret waktu stasioner). Ada juga pendekatan khusus untuk data stasioner dengan musiman yaitu menggunakan Seasonal ARIMA (SARIMA).

3. Peramalan

Peramalan dapat dilakukan dengan menggunakan fungsi forecast() dari packages forecast.

4. Evaluasi Model

Evaluasi model perlu dilakukan untuk melihat apakah model yang sudah dibuat sudah cukup baik digunakan dalam melakukan prediksi. Untuk melakukan evaluasi model digunakan nilai residual sebagai pembanding, salahsatunya bisa menggunakan function MAPE() dari package MLmetrics

5. Pengujian Asumsi

Asumsi diujikan untuk mengukur apakah residual yang diperoleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Residual data digunakan karena dengan residual dapat diperoleh informasi dari data aktual maupun dari hasil prediksi menggunakan model. terdapat dua asumsi yang diujikan antara lain

Pengujian Asumsi pada time series hanya dilakukan apabila menggunakan model ARIMA dan SARIMA karena kedua model tersebut berdasar pada Regresi.

5.1. Normalitas Residual

Diharapkan residual yang dihasilkan dari model yang telah dibuat berdistribusi normal, untuk melakukan pengujian normalitas pada residual dilakukan dengan uji Shapiro-Wilk menggunakan fungsi shapiro.test. Dengan hipotesis pengujian sebagai berikut:

• H0: residual berdistribusi normal
• H1: residual tidak berdistribusi normal

kriteria pengujian: tolak H0 apabila nilai p-value < alpha (0.05)

5.2. Autocorelation

Diharapkan tidak terjadi autokorelasi residual dari model yang telah dibuat, untuk melakukan pengujian autokorelasi pada residual dilakukan dengan uji Ljung Box menggunakan fungsi Box.test. Dengan hipotesis pengujian sebagai berikut:

• H0: tidak terdapat autokorelasi pada residual
• H1: terdapat autokorelasi residual

kriteria pengujian: tolak H0 apabila nilai p-value < alpha (0.05)

Objek TS

Dikarenakan data yang digunakan merupakan rekapan jumlah penyewa sepeda harian, maka pada tahapan ini akan dicoba untuk membuat objek ts dengan mengatur frekuensi = 7 untuk pola mingguan.

Setelah dilakukan pembuatan objek ts, juga akan dilakukan tahapan visualisasi untuk melihat apakah terdapat pola dari objek ts yang sudah dibuat menggunakan fungsi autoplot().

#membuat objek ts 1
ts(data = bike$cnt,      #mengambil kolom nilai pada objek ts
          start = range(bike$dteday)[[1]], #menginformasikan keterangan waktu mulai pada data
          frequency = 7) %>%  #pola musiman diindikasikan berulang setiap 7 hari

#visualisasi objek ts
  autoplot()

dari visualisasi objek ts diatas terlihat bahwa

• data fluktuatif
• tren dari data terlihat menaik meskpun belum terlihat begitu jelas karena terlihat ada penurunan

berdasarkan pada informasi data yang disediakan oleh peneliti sebelumnya, yakni adanya kemungkinan pengaruh pola musim di wilayah eropa terhadap jumlah penyewaan sepeda, maka akan coba kita bandingkan rata-rata jumlah penyewa sepeda berdasarkan musim.

#menghitung rata-rata jumlah penyewa berdasarkan musim dan tahun
a <- bike %>% 
  group_by(season, yr) %>% 
  summarise(mean_cnt = mean(cnt)) %>% 
  mutate(label = glue("Musim = {season}
                      Tahun = {yr}
                      Jumlah Penyewa = {round(mean_cnt,2)}"))

#membuat visualisasi diagram batang
plot1 <- ggplot(data = a, aes(y = mean_cnt, 
                              x = reorder(season,mean_cnt), 
                              text = label))+ 
  geom_col(aes(fill = yr))+
  labs(title = "Rata-Rata Jumlah Penyewa Sepeda berdasarkan Musim dan Tahun",
       x = "Musim",
       y = "Rata-Rata Jumlah Penyewa")+
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none")+
  facet_wrap(facets = ~yr)

#membuat plot interaktif
ggplotly(plot1, tooltip = "text")

dari visualisasi rata-rata jumlah penyewa berdasarkan musim dan tahun diatas terlihat bahwa:

• Jumlah penyewa sepeda terbanyak terjadi di musim gugur (Fall) yang biasa terjadi pada bulan Agustus hingga Oktober.
• Jumlah penyewa paling sedikit terjadi di musim semi (Springer) yang biasa terjadi pada bulan Februari hingga April.

untuk melihat lebih detail pola dari objek ts yang sudah kita buat maka parlu dilakukan dekomposisi yang akan memecah data menjadi tiga bagian utama yakni trend, seasonal dan remainder/error menggunakan fungsi decompose().

ts(data = bike$cnt,      #mengambil kolom nilai pada objek ts
          start = range(bike$dteday)[[1]], #menginformasikan keterangan waktu mulai pada data
          frequency = 7) %>%  #pola musiman diindikasikan berulang setiap 7 hari
decompose() %>%  #melakukan dekomposisi
  autoplot() #visualisasi objek ts 1 

dari visualisasi hasil dekomposisi objek ts diatas terlihat bahwa:

• pada panel pertama merupakan plot objek ts yang menunjukan bahwa objek ts tidak memiliki efek multiplikatif.
• pada panel kedua merupakan plot tren yang menunjukan bahwa komponen estimasi tren menunjukkan sebuah pola. Pola tren ini mungkin bersumber dari musiman ekstra yang tidak ditangkap dari periode alami yang lebih tinggi dalam kasus ini, sehingga dapat dianggap sebagai data multi-musiman. Untuk mengatasi musiman yang kompleks ini, kita perlu mengubah data menjadi objek msts() yang menerima pengaturan frekuensi ganda.
• pada panel ketiga merupakan plot seasonal, yang menunjukan bahwa objek ts memiliki seasonal

Objek msts ke-1

Dari hasil pembuatan objek ts dan dekomposisi sebelumnya kita dapati bahwa ada kemungkinan terjadi multi musiman, yang sebelumnya kita gunakan pola musiman mingguan kini kita akan coba tambahkan informasi pola musiman yang lebih tinggi yakni bulanan untuk membuat objek msts() dengan menambahkan frekuensi = 7*4.

Untuk Objek msts, dekomposisi dilakukan menggunakan fungsi mstl()

#membuat objek msts 1
bike$cnt %>% 
  msts(seasonal.periods = c(7,7*4)) %>% # multiseasonal ts (weekly, monthly)
  mstl() %>% # multiseasonal ts decomposition
  autoplot() 

hasil pembuatan objek msts kedua masih terlihat adanya pola pada tren, sehingga akan dicoba untuk membuat objek msts dengan menambahkan informasi pola musiman yang lebih tinggi.

Objek msts ke-2

Pada tahapan ini akan dicoba untuk membuat objek msts yang menangkap pola musiman bulanan dan kuartal.

#membuat objek ts 3
bike$cnt %>% 
  msts(seasonal.periods = c(7*4,7*4*3)) %>% # multiseasonal ts (Weekly, monthly, quartly)
  mstl() %>% # multiseasonal ts decomposition
  autoplot() 

hasil pembuatan objek msts kedua pola pada tren sudah nampak lebih smooth, dan menunjukkan tren yang naik.

Dari ketiga objek ts yang sudah coba buat, Objek ts terakhir menunjukkan dekomposisi terbaik, terlihat pada bagian tren yang sudah terlihat lebih smooth dan tidak lagi menunjukan adanya indikasi pola. Maka akan digunakan objek msts 2 untuk tahapan selanjutnya.

#menyimpan objek ts final
bike_msts2 <- bike$cnt %>% 
  msts(seasonal.periods = c(7*4,7*4*3))

Model ETS Holt-Winters

1. Pemodelan

#pembuatan model ets Holt-Winters
bike_ets <- stlm(bike_train,  # menggunakan data train
                 s.window = c(7*4,7*4*3),
                 method = "ets",   # metode yang digunakan adalah ETS holts-winter
                 lambda = NULL)    # model stlm akan mengidentifikasi apakah data additive/multiplikatif

#menampilkan ringkasan model yang terbentuk
summary(bike_ets$model)

#> ETS(A,N,N) 
#> 
#> Call:
#>  ets(y = x, model = etsmodel, allow.multiplicative.trend = allow.multiplicative.trend) 
#> 
#>   Smoothing parameters:
#>     alpha = 0.2111 
#> 
#>   Initial states:
#>     l = 1210.1622 
#> 
#>   sigma:  805.9655
#> 
#>      AIC     AICc      BIC 
#> 12851.11 12851.14 12864.52 
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
#> Training set 39.12542 804.7189 598.7306 -4.83282 18.86044 0.3768627 0.1941238

dari hasil pemodelan diatas diketahui bahwa model ETS yang terbentuk adalah ETS(A,N,N) yang menunjukkan

• A = Pola data Aditif
• N pertama = tidak terdapat tren
• N kedua = tidak terdapat seasonal

dan nilai MAPE sebesar 18.86044 menunjukkan bahwa rata-rata kesalahan prediksi model tersebut adalah sebesar 18.86044% dari nilai aktual pada data train

2. Forecasting & Evaluasi model

Untuk melakukan forecasting akang digunakan function forecast(model, h = banyaknya data yg akan di prediksi)

dan untuk melakukan evaluasi model akan digunakan nilai error antara nilai aktual dari data test dengan hasil prediksi model. Untuk melakukan evaluasi model akan digunakan function MAPE(nilai prediksi, nilai aktual) dari packages MLmetrics

# forecasting
bike_ets_f <- forecast(bike_ets,  #model yang digunakan untuk prediksi
                       h = 7*4*3) #banyaknya nilai yang akan diprediksi

# model evaluation: Mean Absolute Percentage Error (MAPE)
MAPE(bike_ets_f$mean, # nilai hasil prediksi
     bike_test)       # nilai asli data test

#> [1] 3.968679

dari hasil evaluasi model diatas diperoleh nilai MAPE sebesar 3.968679 menunjukkan bahwa rata-rata kesalahan prediksi model tersebut adalah sebesar 3.968679% dari nilai aktual pada data test.

Model SARIMA

1. Pemodelan

#pembuatan model ets Holt-Winters
bike_sarima <- stlm(bike_train,  # menggunakan data train
                 s.window = c(7*4,7*4*3),
                 method = "arima",   # metode yang digunakan adalah arima
                 lambda = NULL)    # model stlm akan mengidentifikasi apakah data additive/multiplikatif

#menampilkan ringkasan model yang terbentuk
summary(bike_sarima$model)

#> Series: x 
#> ARIMA(0,1,3) with drift 
#> 
#> Coefficients:
#>           ma1      ma2      ma3   drift
#>       -0.6016  -0.1869  -0.0616  8.8597
#> s.e.   0.0390   0.0479   0.0399  4.6530
#> 
#> sigma^2 = 611067:  log likelihood = -5218.5
#> AIC=10447   AICc=10447.09   BIC=10469.35
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 1.194816 778.6812 569.4619 -5.654257 18.09717 0.3584399
#>                      ACF1
#> Training set 0.0003992559

dari hasil pemodelan diatas diketahui bahwa model ARIMA yang terbentuk adalah ARIMA(0,1,3) yang menunjukkan bahwa

• model tidak memiliki komponen AR (nilai AR = 0),
• memiliki satu tingkat differencing (nilai I = 1),
• dan memiliki tiga komponen MA (nilai MA = 3)

dan nilai MAPE sebesar 18.09717 menunjukkan bahwa rata-rata kesalahan prediksi model tersebut adalah sebesar 18.09717% dari nilai aktual dari data train

2. Forecasting & Evaluasi model

Untuk melakukan forecasting akang digunakan function forecast(model, h = banyaknya data yg akan di prediksi)

dan untuk melakukan evaluasi model akan digunakan nilai error antara nilai aktual dari data test dengan hasil prediksi model. Untuk melakukan evaluasi model akan digunakan function MAPE(nilai prediksi, nilai aktual) dari packages MLmetrics

# forecasting
bike_sarima_f <- forecast(bike_sarima,  #model yang digunakan untuk prediksi
                       h = 7*4*3) #banyaknya nilai yang akan diprediksi

# model evaluation: Mean Absolute Percentage Error (MAPE)
MAPE(bike_sarima_f$mean, # nilai hasil prediksi
     bike_test)       # nilai asli data test

#> [1] 4.44462

dari hasil evaluasi model diatas diperoleh nilai MAPE sebesar 4.44462 menunjukkan bahwa rata-rata kesalahan prediksi model tersebut adalah sebesar 4.44462% dari nilai aktual pada data test.

3. Pengujian Asumsi

3.1. Normality

Pengujian normalitas untuk sisaan/error dilakukan dengan uji Shapiro-Wilk menggunakan fungsi shapiro.test

• H0: sisaan/error berdistribusi normal
• H1: sisaan/error tidak berdistribusi normal

kriteria pengujian: tolak H0 apabila nilai p-value < alpha (0.05)

shapiro.test(bike_sarima_f$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  bike_sarima_f$residuals
#> W = 0.95109, p-value = 8.202e-14

berdasarkan hasil pengujian Shapiro-Wilk kita peroleh nilai p-value < alpha, maka kita menolak H0 yang artinya sisaan/error yang kita peroleh tidak berdistribusi normal

3.2. Autocorelation

Pengujian normalitas untuk sisaan/error dari model dengan uji Ljung Box menggunakan fungsi Box.test

• H0: tidak terdapat autokorelasi pada sisaan/error
• H1: terdapat autokorelasi pada sisaan/error

kriteria pengujian: tolak H0 apabila nilai p-value < alpha (0.05)

Box.test(bike_sarima_f$residuals, type = "Ljung-Box")

#> 
#>  Box-Ljung test
#> 
#> data:  bike_sarima_f$residuals
#> X-squared = 0.00010361, df = 1, p-value = 0.9919

berdasarkan hasil pengujian Ljung-Box kita peroleh nilai p-value > alpha, maka kita menerima H0 yang artinya tidak terdapat autokorelasi pada sisaan/error