1 Time Series and Forecasting

Glossary:

  1. Timeseries: Data yang berhubungan dengan waktu.
  2. Forecasting: Prediksi/Peramalan dimasa depan.
  3. Objek ts: untuk analisis time series di R, kita perlu convert dataframe menjadi objek ts.
  • Data: Pilih satu variabel yang ingin diamati
  • start: periode awal data kita.
  • frequency: untuk penentuan pola musiman yang ingin diperhatikan.
  1. Seasonal Adjustment: Penyesuaian musiman, data yang sudah dilakukan proses penyesuaian musiman dapat menunjukkan pola (behaviour) yang riel dimana pengaruh musim sudah dikeluarkan dari data aslinya. (Bps, 2019).

2 Working with Time Series

Apa itu time series? Apa perbedaannya dengan regresi atau yang lainnya?

Time series adalah suatu object dalam statistik dimana object tersebut berhubungan dengan suatu deret waktu tertentu. Objek ini banyak ditemui dalam kehidupan sehari-hari, contoh: harga daging sapi harian, curah hujan bulanan, kuantitas penumpang bulanan, pendapatan tahunan, dll.

Discussion: Apa perbedaan Time Series dengan Regresi?

Perbedaan mendasar antara time series dengan regresi adalah jika pada regresi untuk memprediksi suatu nilai Y dipengaruhi oleh beberapa faktor yaitu x1,x2,..,xn. Sedangkan jika time series, untuk memprediksi suatu nilai Y dipengaruhi oleh nilai Y itu sendiri pada masa lampau (\(Y_{t-1}\)).

Regression

\[y = \beta_0+\beta_1*x_1+\beta_2*x_2+...+\beta_n*x_n\]

Time Series

\[y_t = \beta_0+\beta_1*y_{t-1}+\beta_2*y_{t-2}+...+\beta_n*y_{t-n}\]

Analisis time series berhubungan dengan suatu data yang memiliki nilai numerik pada interval waktu tertentu. Proses untuk memprediksi nilai pada anilisis time series disebut sebagai peramalan atau forecasting. Ide utama dalam melakukan forecasting itu adalah korelasi dari data numerik.

2.1 Karakteristik Data Time Series

Time series data: data yang berhubungan dengan waktu dan memiliki interval waktu yang tetap/sama.

💡 Syarat utama data time series:

  1. Data harus urut sesuai periode waktu dari data terlama sampai ke data terbaru
  2. Interval waktunya harus tetap/sama
  3. Tidak boleh ada data yang terlewat untuk setiap interval
  4. Tidak boleh ada yang missing

❓ Knowledge Check

Apakah data berikut sudah memenuhi syarat data time series yang baik?

  1. Demand product
df_demand <- data.frame(
  date = ymd(c("2021-5-3", "2021-5-4", "2021-5-6", "2021-5-7")),
  demand = c(29, 79, 41, 88)
  )
df_demand

Kesimpulan: data ini belum memenuhi syarat time series. ada periode waktu yang bolong

💡 Pengecekan yang dapat dilakukan:

# inspect periode terlewat
all(seq.Date(from = as.Date("2021-05-03"), to = as.Date("2021-05-07"), by = "day") == df_demand$date)
#> [1] FALSE
seq.Date(from = as.Date("2021-05-03"),
         to = as.Date("2021-05-07"),
         # length.out = 5,
         by = "day")
#> [1] "2021-05-03" "2021-05-04" "2021-05-05" "2021-05-06" "2021-05-07"
  1. Price product
df_price <- data.frame(
  date = ymd(c("2021-5-16", "2021-5-19", "2021-5-18", "2021-5-17", "2021-5-20", "2021-5-21")),
  price = c(1000, 1001, 1002, 1003, 1004, 1005)
  )

df_price

Kesimpulan: belum memenuhi syarat data time series. datanya tidak terurut

💡 Perbaikan yang dapat dilakukan sesuai syarat time series:

  • Mengurutkan data berdasarkan waktu: arrange()
# mengurutkan df_price
df_price %>% 
  arrange(date)
  • Melakukan padding untuk memastikan interval data sama: pad() dari package padr

Secara default, pad() akan menambal periode waktu (tanggal) berdasarkan kolom yang tipe datanya date. Mengisi nilai yang terlewat atau missing (NA), cara yang umum dilakukan dengan package zoo:

  • na.fill(): mengisi NA dengan sebuah nilai, Gunakan fill="extend" untuk mengisi dengan nilai rata-rata dengan nilai yang missing dari rata-rata nilai sebelum dan setelahnya
  • na.aggregate(): nilai aggregasi (mean, median)
  • na.locf(): nilai terakhir sebelum missing

Note: metode untuk mengisi missing value disesuaikan dengan perspektif dari businessnya.

# kembali ke kasus pertama, data df_demand
# case: toko tutup di tanggal 5
anyNA(df_demand)
#> [1] FALSE
# kita isi NA dengan nilai 0
df_demand %>% 
  # pad(start_val = min(df_demand$date), end_val = max(df_demand$date)) %>% 
  pad() %>% 
  mutate(demand = na.fill(demand, fill="extend"))
  # mutate(demand = na.locf(demand))

2.2 Time Series Object

Kita akan menggunakan data emisi CO2 di Indonesia dimana datanya sudah tersimpan dalam folder data_input dengan nama environment_1970f.csv. Dari data co2 ini, kita akan menggunakan 2 kolom yang kita butuhkan yaitu kolom year untuk menunjukkan waktu dan CO2.emissions..metric.tons.per.capita. sebagai nilai Y yang kita amati untuk membuat object ts.

2.2.1 Read data

# read data co
co2 <- read.csv("data_input/environment_1970f.csv")
head(co2)

Data co2 terdiri dari 43 observasi yang mewakili kontribusi gas emisi per tahun terhadap atmosfer Indonesia (43 tahun, 1970-2012). Data ini terdiri dari 7 variabel, yaitu:

  • year: tahun.
  • emisi co2 (kt): emisi yang berasal dari pembakaran bahan bakar fosil dan pembuatan semen, termasuk yang dihasilkan selama konsumsi.
  • emisi co2 (metrik ton per kapita): idem.
  • emisi metana (kt setara co2): emisi yang berasal dari aktivitas manusia (pertanian) dan dari produksi metana industri.
  • emisi nitro oksida (ribu metrik ton setara co2): emisi dari pembakaran biomassa pertanian, kegiatan industri, dan pengelolaan ternak.
  • emisi gas rumah kaca dan lainnya, HFC, PFC dan SF6 (ribu metrik ton setara co2): emisi hasil samping dari hidrofluorokarbon, perfluorokarbon, dan sulfur hexafluoride.
  • total emisi gas rumah kaca (setara dengan co2): total CO2 tidak termasuk pembakaran biomassa siklus pendek (pembakaran limbah pertanian dan savannah), tetapi termasuk pembakaran biomassa lainnya (kebakaran hutan, pembusukan pasca-pembakaran, kebakaran gambut, dan pembusukan lahan gambut yang dikeringkan), semua sumber CH4 antropogenik, sumber N2O dan gas-F (HFC, PFC, dan SF6).

💡 Syarat utama data time series:

  1. Data harus urut sesuai periode waktu dari data terlama sampai ke data terbaru (arrange)
  2. Interval waktunya harus tetap/sama
  3. Tidak boleh ada data yang terlewat untuk setiap interval
  4. Tidak boleh ada yang missing

2.2.2 cek range atau periode waktu

untuk mengetahui periode awal dan akhir data time series kita bisa menggunakan fungsi range().

range(co2$year)
#> [1] 1970 2012
# untuk mengetahui apakah periode sudah terurut
all(1970:2012 == co2$year)
#> [1] TRUE

setelah data terurut, cek missing value data

anyNA(co2$CO2.emissions..metric.tons.per.capita.)
#> [1] FALSE

untuk mengurutkan data berdasarkan periode waktu dan melakukan treatment lain agar sesuai syarat data time series. (pada kasus ini akan error karena data deret waktu dalam tahun yang tidak bertipe date)

# co2 %>% 
#   arrange(year) %>% # mengurutkan waktu
#   pad() %>% # mengeluarkan interval waktu yang bolong
#   mutate(CO2.emissions..metric.tons.per.capita. = 
#            na.fill(CO2.emissions..metric.tons.per.capita., fill = "extend")) # mengisi data kosong

2.2.3 membuat Object TS

Untuk membuat sebuah object time series pada R kita bisa menggunakan function ts() dengan parameter yang diperlukan yaitu :

ts(data, start, frequency)

  • data = kolom data yang akan kita prediksi

  • start = waktu awal mula data yang akan diprediksi (opsional)

  • frequency = pola berulang dari data

  • membuat object ts dengan ts() function

# simpan dalam object co2_ts
co2_ts <- ts(data = co2$CO2.emissions..metric.tons.per.capita.,
             start = 1970,
             frequency = 1)

# cek kelas dari co2_ts
class(co2_ts)
#> [1] "ts"

2.2.4 Cara menentukan frequency

Frequency merupakan banyaknya data dalam satu pola musiman, untuk mengetahui frequency kita perlu mengetahui bagaimana data tersusun dan pola apa yang ingin kita amati.

Latihan :

  1. data jam -> pola harian -> freq = 24
  2. data harian -> pola mingguan -> freq = 7
  3. data harian -> pola weekday -> freq = 5 (data weekend harus dibuang)
  4. data jam -> pola mingguan -> freq = 24*7
  5. data bulanan -> pola tahunan -> freq = 12
  6. data bulanan -> pola kuartalan -> freq = 3
  7. data tahunan -> pola tahunan -> freq = 1 (tidak ada pola musiman)
  8. data tahunan -> pola dekade -> freq = 10
  9. data mingguan -> pola tahunan -> freq = 52
  • memvisualisasikan data time series untuk melihat pola datanya
# inspect pola data
co2_ts %>% 
  autoplot()

Insight : di antara waktu 1996 - 1999 ada penurunan tajam, dan di tahun 2009 turun lalu menukik naik di tahun 2010.

2.2.5 subset object ts

  • dari plot diatas, misalkan kita ingin melakukan subset periode pada object ts dapat menggunakan fungsi window()

window(object_ts, start, end)

Misalkan dari data co2_ts kita ingin mensubset periode dari 1995-2000

# window
window(co2_ts, start = 1996, end = 1999) %>% autoplot()

2.3 Dive Deeper

Kita akan coba menggunakan data births terdiri dari 168 observasi tingkat kelahiran per bulan di New York. Data ini terdiri dari 2 variabel, yaitu:

  • date: tanggal.
  • births: tingkat kelahiran.

data birth data tingkat kelahiran di new york city dari tahun Jan 1946 - Des 1959

  1. Read data nybirth.csv
# read nybirth.csv
birth <- read.csv("data_input/nybirth.csv")
head(birth)
  1. Data wrangling & exploration
  • ubah tipe data pada kolom date
# pak bani
birth$date <- ymd(birth$date)

# cara dplyr
birth %>% 
  mutate(date = ymd(date))
  • cek kelengkapan interval waktu
# range waktu
range(birth$date)
#> [1] "1946-01-01" "1959-12-01"
# cek apakah interval data sudah sesuai dan tidak ada yang bolong
interval <- seq.Date(from = as.Date("1946-01-01"),
                     to = as.Date("1959-12-01"),
                     by = "month")
all(interval == birth$date)
#> [1] TRUE
  • cek missing value
anyNA(birth)
#> [1] FALSE
  1. Buat object ts buatlah sebuah object ts dari data birth lalu simpanlah dalam object birth_ts dengan data bulanan yang ada. Pola yang ingin dilihat adalah pola tahunan.
# answer here
birth_ts <- ts(data = birth$births,
               start = c(1946, 1),
               frequency = 12)

waktu paling awal : 1950-05-3 (data harian, ingin lihat pola tahunan) frequency = 365 start = c(1950, 123)

  1. Visualisasi object birth_ts menggunakan autoplot()
birth_ts %>% autoplot()

2.3.1 Visualisasi data birth

birth_clean <- birth %>% 
  mutate(date = ymd(date),
         month = month(date, label = T))
ggplot(data = birth_clean, aes(x = date, y = births)) +
  geom_point() +
  geom_point(data = birth_clean %>% 
               filter(month %in% c("Jan", "Feb", "Jul")),
             aes(col = month)) +
  geom_line() +
  scale_color_manual(values = c("red", "blue", "yellow"))

jan: turun ke feb feb: paling dasar jul: paling tinggi

2.4 Decomposition

Decomposition adalah suatu tahapan dalam time series analisis yang digunakan untuk menguraikan beberapa komponen dalam time series data.

💡 Komponen dalam time series :

  • 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

Sebelum melakukan modeling forecasting kita perlu mengamati objek timeseries dari hasil decompose. Ide utama dari decompose adalah untuk menguraikan ketiga komponen dari objek ts (trend, seasonal, residual).

Jika pada hasil decompose, trend masih membetuk sebuah pola maka dapat dicurigai masih ada seasonality yang belum ditangkap. Seharusnya trend cenderung naik atau cendurung turun.

Untuk dapat menguraikan data time series kita menjadi 3 komponen tersebut, kita dapat menggunakan fungsi decompose().

birth_decom <- decompose(birth_ts)

Memvisualisasikan hasil decompose menggunakan autoplot() dari package forecast

# autoplot
birth_decom %>% autoplot()

Pada hasil decompose kita mendapatkan informasi visualisasi:

  1. Data : pola data asli
  2. Seasonal (S) : pola musiman atau pola berulang dari data
  3. Trend (T) : pola data secara global (naik atau turun)
  4. Remainder (E) : pola data yang tidak dapat ditangkap oleh seasonal dan trend

💡 Notes: Jika pada hasil decompose, trend masih membentuk sebuah pola maka dapat dicurigai masih ada seasonality yang belum ditangkap. Seharusnya trend cenderung naik atau cenderung turun secara smooth. Penyebabnya:

  • frequency yang kita tetapkan belum tepat, atau
  • terdapat multiseasonality pada data (pola berulang lebih dari satu, visitor bioskop)

2.5 Knowledge Check

  1. Manakah pernyataan di bawah ini yang salah terkait data time series?
  • Data time series harus disusun secara berurut sesuai waktu dengan interval waktu yang tetap.
  • Data time series dapat diurai menjadi pola trend, seasonal, dan error.
  • Data time series disebut additive karena variasi seasonalnya relatif konstan dari waktu ke waktu.
  • Data time series disebut multiplicative karena menunjukkan trend yang linear.
  1. Misalkan kita memiliki object retail merupakan dataframe penjualan bulanan dari tahun 2010 sampai 2015. Kita ingin menganalisis pergerakan nilai pada kolom sales dari dataframe retail. Cara membuat object ts dengan pola seasonality tahunan adalah …
  • ts(retail, frequency = 1)
  • ts(retail, frequency = 12)
  • ts(retail$sales, frequency = 1)
  • ts(retail$sales, frequency = 12)
  1. Decompose adalah bagian dari eksplorasi data pada analisis time series. Berikut adalah manfaat dari melakukan decompose, kecuali …
  • Mengetahui pola naik atau turun-nya data dari trend.
  • Dapat digunakan untuk analisis pola musiman dari pola seasonality.
  • Mengetahui nilai di waktu yang akan datang berdasarkan pola data historis.
  • Mengetahui apakah frequency yang diatur dalam object time series sudah sesuai atau belum.

Gambar tersebut menunjukkan hasil decomposition untuk data time series terkait tingkat kejahatan harian di Chicago dari tahun 2014 sampai dengan 2019.

  1. Berdasarkan grafik tersebut, bagaimana pola trend yang terdapat pada data?
  • tidak terdapat pola tertentu
  • pola meningkat
  • pola menurun
  • pola naik turun secara berulang

2.6 Additive vs Multiplicative Time Series

Ada 2 jenis model pada data time series, yaitu :

  1. Model Additive : Model time series yang memiliki varians konstan mengikuti trend dan seasonalnya

\[Y_t = T_t + S_t + E_t\]

Data = Trend + Seasonal + Error

  1. Model Multiplicative : Model time series yang memiliki varians semakin tinggi/rendah mengikuti trend dan seasonal yang ada

\[Y_t = T_t * S_t * E_t\]

Data = Trend * Seasonal * Error

2.6.1 Additive Time Series

Jika kita perhatikan lagi pada data birth_ts memiliki pola additive, karena varians dari polanya tetap atau konstan. Untuk melakukan decompose pada pola additive tambahkan parameter type = "additive" dimana secara default fungsi decompose() memiliki type = “additive”.

# additive time series
birth_ts %>% autoplot()

# decompose(birth_ts, type = "additive")

Melakukan inspect komponen time series pada data birth_dc

1. Trend

Trend diperoleh dari hasil perhintungan center moving average (CMA). Tujuan utamanya untuk smoothing data sehingga diperoleh trend yang cenderung naik/ cenderung turun.

Trend = Data - Seasonality - Error

birth_decom$trend %>% autoplot()

Pendekatan manual menggunakan Center Moving Average

birth_trend <- ma(x = birth_ts, order = 12, centre = T)
birth_trend %>% autoplot()

# birth_trend %>% ts_plot()

2. Seasonal

Karena polanya additive maka:

Seasonal + Error = Data - Trend

birth_decom$seasonal %>% autoplot()

Pendekatan manual:

# matematic didalamnya
sea.err_birth <- birth_ts - birth_trend
sea.err_birth %>% autoplot()

# mean of each month
mean.month_birth <- sea.err_birth %>% 
  matrix(ncol =  12, byrow = T) %>% 
  colMeans(na.rm = T)

# mean global
mean.glob_birth <- mean(mean.month_birth)

# Seasonality Calculation
birth_seasonal <- ts(rep(mean.month_birth - mean.glob_birth, 12), start = start(birth_ts), frequency = 12)
 
birth_seasonal %>% autoplot()

3. Error

Untuk memperoleh informasi error pada model additive, dapat menggunakan rumus:

Error = Data - Trend - Seasonal

birth_decom$random %>% autoplot()

# matematis didalamnya

birth_error <- birth_ts - (birth_trend + birth_seasonal)
birth_error %>% autoplot()

2.6.2 Analisis Seasonal

Untuk dapat melakukan analisis seasonal, kita bisa lakukan dengan cara melakukan inputasi komponen seasonal yang ada.

birth %>% 
  mutate(monthly = month(date, label = T, abbr = F),
         seasonality = birth_decom$seasonal) %>% 
  distinct(monthly, seasonality) %>% 
  ggplot(aes(x = monthly, y = seasonality)) +
  geom_col()+
  theme_minimal()

birth_ts %>% seasonplot()

2.6.2.1 Summary day 1

  1. Data time series -> data yang berhubungan dengan deret waktu
  2. Forecasting -> prediksi data time series menggunakan pola data di masa lalu
  3. Syarat data time series:
    • harus urut dari paling lama ke paling baru
    • interval waktu harus sama
    • tidak ada periode waktu yang terlewat
    • tidak ada missing value
  4. Membuat object time series, gunakan fungsi ts(), parameter:
    • data: kolom data yang berisi nilai yang ingin dianalisis atau diforecast
    • start: waktu awal data time series
    • frequency: pola berulang/pola musiman dari data yang ingin kita amati
  5. Proses EDA data time series:
    • fungsi autoplot(), secara umum akan terlihat informasi trend, apakah ada pola musiman atau tidak, tipe data time seriesnya additive atau multiplicative
    • decompose(), akan mengeluarkan nilai berikut ketika di-autoplot:
      • trend: pola general naik atau turun
      • seasonal: pola berulang musiman
      • remainder: data sisanya akan mengetahui frequency di awal sudah benar atau belum (melihat pola data dari trend)

– END OF DAY 1 –

2.6.3 Multiplicative Time Series

Jika kita memiliki pola multiplicative pada data time series dan ingin membuat decomposenya cukup menambahkan type = "multiplicative" pada fungsi decompose().

Ketika kita menemukan pola data kita mengandung multiplikative :

cara 1: membuat data multiplikative tersebut menjadi additive dengan fungsi log. Setelah memperoleh hasil forecast kita dapat mengembalikan nilainya dengan exp.

AirPassengers %>% autoplot()

# transformasi log agar data menjadi tipe additive
log(AirPassengers) %>% autoplot()

(Opsional) Sifat logaritma: perkalian menjadi penjumlahan

Rumus logaritma \[^alog(b) = c <=> b^c = a\]

multiplicative \[y = T * S * E\] \[log(y) = log(T * S * E)\] additive \[log(y) = log(T) + log(S) + log(E)\]

decompose(log(AirPassengers), type = "additive") %>% autoplot()

Cara 2: Tetap menggunakan model multiplicative, kemudian nanti hasil dibandingkan dengan memilih model dengan error yang paling kecil.

Parameter type dalam fungsi decompose(), secara default type = "additive"

multi:

  1. transformasi/scaling -> additive
  2. decompose, dengan catatan selalu tambahkan parameter yang menjelaskan bahwa datanya adalah multiplicative
decompose(AirPassengers, type = "multiplicative") %>% autoplot()

Decompose dapat kita gunakan untuk melakukan analisis seasonal dengan tujuan untuk mendapatkan insight dari efek seasonal yang ada.

Tujuan dan manfaat melakukan decompose: 1. Mengidentifikasi tipe time series yang kita miliki, multiplicative atau additive 2. Kita mengetahui komponen trend, seasonal, dan error pada data kita

2.6.4 Seasonality Adjustment [OPSIONAL]

Objek time series yang telah dihilangkan efek seasonalnya.

Tujuan dari seasonality adjusment adalah untuk mengetahui pola data time series yang dimiliki tanpa memperhitungkan efek seasonal yang ada pada data.

Proses ini biasanya dilakukan ketika ingin melihat trend serta penyimpangan siklus trend pada data time series tanpa memperhitungkan kejadian musiman (seasonal).

  1. Read data AirPassengers
AirPassengers
#>      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
#> 1949 112 118 132 129 121 135 148 148 136 119 104 118
#> 1950 115 126 141 135 125 149 170 170 158 133 114 140
#> 1951 145 150 178 163 172 178 199 199 184 162 146 166
#> 1952 171 180 193 181 183 218 230 242 209 191 172 194
#> 1953 196 196 236 235 229 243 264 272 237 211 180 201
#> 1954 204 188 235 227 234 264 302 293 259 229 203 229
#> 1955 242 233 267 269 270 315 364 347 312 274 237 278
#> 1956 284 277 317 313 318 374 413 405 355 306 271 306
#> 1957 315 301 356 348 355 422 465 467 404 347 305 336
#> 1958 340 318 362 348 363 435 491 505 404 359 310 337
#> 1959 360 342 406 396 420 472 548 559 463 407 362 405
#> 1960 417 391 419 461 472 535 622 606 508 461 390 432

Data AirPassengers merupakan data banyaknya penumpang pesawat terbang dari bulan Jan 1949 - Des 1960.

  1. Membuat visualisasi
# prepare dataset
air_df <- data.frame(date = seq(from = ymd("1949-01-01"),
                                to = ymd("1960-12-01"),
                                by = "month"),
                     value = as.data.frame(AirPassengers)$x)

# add monthly column
air_df <- air_df %>%
  mutate(monthly = month(date, label = T))

# visualize
air_df %>% 
  ggplot(aes(x = date, y =value))+
    geom_point(data = air_df %>% filter(monthly == "Jul"), col = "firebrick")+
  geom_line()

  1. Melakukan decomposition

Melakukan decomposition pada data AirPassengers dengan cara 2

air_dc <- AirPassengers %>% decompose(type = "multiplicative")
  1. Menghilangkan efek seasonality dengan menggunakan fungsi seasadj dari package forecast
air_df <- air_df %>% 
  mutate(seasadj = seasadj(air_dc))
air_df %>% 
  ggplot(aes(x = date, y = seasadj))+
  geom_point(data = air_df %>% filter(monthly == "Jul"), col = "firebrick")+
  geom_line()

Notes: Seasonal adjustment hanya digunakan sebagai analisis deskriptif, dan tidak diperhitungkan dalam proses forecasting.

3 Forecasting Model

💡 Forecasting Model:

  1. Simple Moving Average (SMA)
  2. Exponential Smoothing
    • Simple Exponential Smoothing (SES)
    • Double Exponential Smoothing (Holts Exponential Smoothing)
    • Triple Exponential Smoothing (Holts Winter)
  3. Autoregresive Integrated Moving Average (ARIMA)

3.1 Single Moving Average (SMA)

Metode yang menggunakan rataan beregerak untuk melakukan forecasting. Karena menggunakan rataan, bobot yang digunakan sama untuk setiap observasi di masa lalu. Metode ini digunakan untuk data yang tidak mengandung trend dan seasonal (datanya bergerak disekitar rata-rata).

Fungsi yang digunakan untuk forecasting SMA adalah SMA(objek time series, n) dari library TTR. Parameter yang diugunakan, yaitu:

  • x :objek time series yang akan di forecast

  • n: jumlah observasi di masa lalu yang digunakan untuk forecasting

Kita akan coba melakukan forecasting pada data curah hujan tahunan sejak 1813-1912 yang tersimpan pada folder data_input dengan nama file precip1.dat.

  • read data menggunakan scan()
rain <- scan("data_input/precip1.dat", skip = 1)
head(rain)
#> [1] 23.56 26.07 21.86 31.24 23.65 23.88
  • membuat object ts
rain_ts <- ts(data = rain,
              start = 1813,
              frequency = 1)
  • Visualisasikan data
rain_ts %>% autoplot()

  • Melakukan forecasting menggunakan SMA() dengan ordo 3
# buat model
rain_sma3 <- SMA(x = rain_ts, n = 3)
# forecast dari hasil model
rain1913 <- mean(tail(rain_ts, 2), 26.01000)
rain1914 <- mean(tail(rain_ts, 1), 26.01000, rain1913)
rain1915 <- mean(rain1913, rain1914, 26.01000)

rain1913;rain1914;rain1915
#> [1] 26.335
#> [1] 27.88
#> [1] 26.335

Discussions: tes dengan ordo 10

rain_sma10 <- SMA(x = rain_ts, n = 10)

Memvisualisasikan data historis dengan hasil forecast

rain_ts %>% 
  autoplot() + 
  autolayer(rain_sma3) +
  autolayer(rain_sma10)

Kesimpulan: semakin tinggi ordonya maka model yang dihasilkan semakin smooth. Ada pola yang terbentuk tetapi terjadi keterlambatan (lagging) yang tidak kita inginkan, jika ada lagging maka bisa dibilang modelnya belum baik.

💡 Note:

  • SMA tepat digunakan ketika data yang kita miliki tidak memiliki trend maupun seasonal (data stasioner).
  • penentuan ordo tidak ada rule of thumb nya
  • semakin besar ordo, semakin smooth datanya dan hasil forecastnya

3.2 Exponential Smoothing

Metode SMA hanya mempertimbangkan n observasi di masa lampau untuk melakukan forecast baik pola trend maupun seasonal cenderung tidak tertangkap sehingga dibutuhkan metode lain yang memperhitungkan keseluruhan data di masa lampau, yaitu metode eksponensial. Fungsi yang digunakan untuk membuat model exponential smoothing, yaitu:

3.2.1 ets()

dengan fungsi ets() dari library forecast yang merupakan error, trend, dan seasonal. Parameter yang digunakan yaitu:

  • y: objek timeseries yang ingin ditraining
  • model: model time series untuk error, trend, dan seasonal.
    • A: aditif
    • M: multiplikatif
    • N: none
    • Z: auto
  • alpha: bobot untuk smoothing eror, nilainya antara 0 - 1.
    • ~ 1 observasi yang paling baru diberikan bobot yang paling tinggi tinggi dibandingkan dengan data periode lama
    • ~ 0 observasi yang paling baru diberikan bobot yang paling tinggi tinggi namun perbedaannya sangat sedikit dengan data periode lama
  • beta: bobot untuk smoothing trend, nilainya antara 0 - 1.
  • gamma: bobot untuk smoothing seasonal, nilainya antara 0 - 1.

SES: ets(data, model = “*NN”)

Holts exponential: ets(data, model = “**N”)

Holts Winter: ets(data, model = “***“)

3.2.2 HoltWinters()

dengan fungsi HoltWinters(), parameter yang digunakan yaitu: * objek timeseries yang digunakan. * secara default parameter alpha, beta, dan gamma adalah NULL, dimana apabila kita tidak mendefinisikan nilainya, maka model HoltWinters() akan mencari nilai parameternya hingga mendapatkan nilai paling optimum. Sehingga jika objek time series tidak mengandung trend dan seasonal parameter beta dan gamma harus diubah menjadi FALSE.

  • alpha: parameter untuk bobot bobot error
  • beta: parameter untuk bobot trend
  • gamma: parameter untuk bobot seasonal
  • seasonal: untuk menentukan perhitungan berdasarkan tipe data time series, secara default akan digunakan additive, nilai yang dapat dimasukkan yaitu additive dan multiplicative

Simple Eksponential Smoothing: HoltWinters(data, beta = F, gamma = F)

Holt’s Exponential: HoltWinters(data, gamma = F)

Holt’s Winter Explonential: HoltWinters(data)

3.2.3 Simple Exponential Smoothing (SES)

Simple Exponential Smoothing merupakan metode forecasting yang digunakan untuk data time series yang tidak mengandung trend dan seasonal1.

Kita akan coba menggunakan data rain_ts yang kita miliki untuk melakukan forecasting menggunakan SES.

# visualisasikan rain
rain_ts %>% autoplot()

Data rain_ts tidak memiliki pola trend maupun seasonal. Karena datanya memiliki pola tahunan, oleh karena itu data rain tidak memiliki pola seasonal yang dapat dilihat.

SES: ets(data, model = “*NN”)

Melakukan forecasting data rain_ts dengan metode SES (tidak mempertimbangkan pola trend)

# menggunakan fungsi `ets()`
rain_ses <- ets(rain_ts, model = "ZNN")
rain_ses
#> ETS(A,N,N) 
#> 
#> Call:
#>  ets(y = rain_ts, model = "ZNN") 
#> 
#>   Smoothing parameters:
#>     alpha = 1e-04 
#> 
#>   Initial states:
#>     l = 24.8243 
#> 
#>   sigma:  4.2362
#> 
#>      AIC     AICc      BIC 
#> 753.2297 753.4797 761.0452

rata2: 24.8239

# untuk melihat hasil predict menggunakan fitted values
rain_ses$fitted
#> Time Series:
#> Start = 1813 
#> End = 1912 
#> Frequency = 1 
#>   [1] 24.82431 24.82418 24.82431 24.82401 24.82465 24.82454 24.82444 24.82460
#>   [9] 24.82439 24.82507 24.82498 24.82490 24.82566 24.82551 24.82528 24.82510
#>  [17] 24.82541 24.82546 24.82548 24.82577 24.82527 24.82527 24.82480 24.82475
#>  [25] 24.82501 24.82447 24.82415 24.82442 24.82388 24.82451 24.82434 24.82444
#>  [33] 24.82422 24.82401 24.82417 24.82345 24.82395 24.82376 24.82320 24.82278
#>  [41] 24.82384 24.82394 24.82332 24.82315 24.82289 24.82262 24.82202 24.82236
#>  [49] 24.82310 24.82284 24.82312 24.82279 24.82201 24.82247 24.82315 24.82329
#>  [57] 24.82315 24.82321 24.82286 24.82288 24.82378 24.82357 24.82297 24.82333
#>  [65] 24.82346 24.82380 24.82472 24.82562 24.82617 24.82648 24.82671 24.82667
#>  [73] 24.82622 24.82640 24.82662 24.82606 24.82635 24.82625 24.82589 24.82622
#>  [81] 24.82600 24.82550 24.82581 24.82547 24.82534 24.82515 24.82443 24.82420
#>  [89] 24.82405 24.82378 24.82339 24.82471 24.82430 24.82411 24.82405 24.82387
#>  [97] 24.82376 24.82395 24.82400 24.82400
# bandingkan dengan nilai rata-rata
mean(rain_ts)
#> [1] 24.8239
# cek error menggunakan rmse
sqrt(rain_ses$mse)
#> [1] 4.193615

Menggunakan alpha = 0.03

rain_ses2 <- ets(rain_ts, model = "ZNN", alpha = 0.03)
sqrt(rain_ses2$mse)
#> [1] 4.244161
rain_ses3 <- ets(rain_ts, model = "ZNN", alpha = 0.6)
sqrt(rain_ses3$mse)
#> [1] 5.123114

Notes: Untuk membuat model SES menggunakan HoltWinters() kita bisa setup beta = F, dan gamma = F karena dalam model SES tidak memperhitungkan efek trend dan seasonal. (default alpha, beta, dan gamma = NULL)

# menggunakan fungsi `HoltWinters()`
rain_ses_hw <- HoltWinters(rain_ts, beta = F, gamma = F)
rain_ses_hw$fitted
#> Time Series:
#> Start = 1814 
#> End = 1912 
#> Frequency = 1 
#>          xhat    level
#> 1814 23.56000 23.56000
#> 1815 23.62054 23.62054
#> 1816 23.57808 23.57808
#> 1817 23.76290 23.76290
#> 1818 23.76017 23.76017
#> 1819 23.76306 23.76306
#> 1820 23.82691 23.82691
#> 1821 23.79900 23.79900
#> 1822 23.98935 23.98935
#> 1823 23.98623 23.98623
#> 1824 23.98921 23.98921
#> 1825 24.19282 24.19282
#> 1826 24.17032 24.17032
#> 1827 24.13171 24.13171
#> 1828 24.10442 24.10442
#> 1829 24.19549 24.19549
#> 1830 24.22261 24.22261
#> 1831 24.24329 24.24329
#> 1832 24.32812 24.32812
#> 1833 24.21938 24.21938
#> 1834 24.23290 24.23290
#> 1835 24.13369 24.13369
#> 1836 24.13867 24.13867
#> 1837 24.21782 24.21782
#> 1838 24.10257 24.10257
#> 1839 24.04293 24.04293
#> 1840 24.12608 24.12608
#> 1841 24.01280 24.01280
#> 1842 24.18448 24.18448
#> 1843 24.15808 24.15808
#> 1844 24.19889 24.19889
#> 1845 24.16153 24.16153
#> 1846 24.12748 24.12748
#> 1847 24.18133 24.18133
#> 1848 24.02499 24.02499
#> 1849 24.16454 24.16454
#> 1850 24.13476 24.13476
#> 1851 24.01621 24.01621
#> 1852 23.93453 23.93453
#> 1853 24.20964 24.20964
#> 1854 24.25018 24.25018
#> 1855 24.11509 24.11509
#> 1856 24.08964 24.08964
#> 1857 24.04430 24.04430
#> 1858 23.99933 23.99933
#> 1859 23.87319 23.87319
#> 1860 23.97780 23.97780
#> 1861 24.17710 24.17710
#> 1862 24.13110 24.13110
#> 1863 24.21405 24.21405
#> 1864 24.15075 24.15075
#> 1865 23.97658 23.97658
#> 1866 24.10933 24.10933
#> 1867 24.29001 24.29001
#> 1868 24.33729 24.33729
#> 1869 24.31468 24.31468
#> 1870 24.34134 24.34134
#> 1871 24.26847 24.26847
#> 1872 24.28659 24.28659
#> 1873 24.51752 24.51752
#> 1874 24.47295 24.47295
#> 1875 24.33660 24.33660
#> 1876 24.43558 24.43558
#> 1877 24.47717 24.47717
#> 1878 24.56625 24.56625
#> 1879 24.79573 24.79573
#> 1880 25.01341 25.01341
#> 1881 25.14045 25.14045
#> 1882 25.20750 25.20750
#> 1883 25.25411 25.25411
#> 1884 25.23351 25.23351
#> 1885 25.11571 25.11571
#> 1886 25.15248 25.15248
#> 1887 25.19729 25.19729
#> 1888 25.05286 25.05286
#> 1889 25.11768 25.11768
#> 1890 25.08710 25.08710
#> 1891 24.99407 24.99407
#> 1892 25.07019 25.07019
#> 1893 25.01085 25.01085
#> 1894 24.88515 24.88515
#> 1895 24.95884 24.95884
#> 1896 24.87469 24.87469
#> 1897 24.84201 24.84201
#> 1898 24.79420 24.79420
#> 1899 24.62284 24.62284
#> 1900 24.57259 24.57259
#> 1901 24.54141 24.54141
#> 1902 24.48421 24.48421
#> 1903 24.39631 24.39631
#> 1904 24.72686 24.72686
#> 1905 24.62852 24.62852
#> 1906 24.58852 24.58852
#> 1907 24.58059 24.58059
#> 1908 24.54271 24.54271
#> 1909 24.52166 24.52166
#> 1910 24.57541 24.57541
#> 1911 24.59433 24.59433
#> 1912 24.59905 24.59905

visualisasi:

rain_ts %>% 
  autoplot() +
  autolayer(rain_ses$fitted)+
  autolayer(rain_ses2$fitted)+
  autolayer(rain_ses3$fitted) +
  autolayer(rain_ses_hw$fitted)

Kesimpulan: Alpha yang menggunakan nilai otomatis (0.0001) akan menghasilkan model dengan garis lurus yang merupakan nilai rata-rata dari data kita, karena nilai tersebut menghasilkan error yang paling kecil (sesuai dengan teori Ordinary Least Square/OLS pada materi regresi). Semakin besar alphanya maka pola data time series semakin terbentuk/semakin mendekati pola aslinya tetapi masih ada lagging.

Model evaluation

Karena kita melakukan pemodelan tanpa cross validation kita tidak perlu melakukan prediksi/forecast sebab hasil forecast data historis sudah terdapat pada output hasil pemodelan. Kita dapat menghitung error yang dihasilkan oleh model dengan menggunakan fungsi accuracy() dari library forecast.

# use accuracy
accuracy(rain_ses2)
#>                       ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -0.02942263 4.244161 3.371061 -2.893919 13.73626 0.6943556
#>                     ACF1
#> Training set -0.07259652

Untuk melihat nilai error menggunakan MAPE(Mean Absolute Percentage Error), kita bisa tahu seberapa persen nilai error dari prediksi kita.

3.2.4 Holt’s Exponential Smoothing

Holt’s exponential (Double Exponential Smoothing) digunakan untuk data time series yang tidak mengandung efek seasonal2.

Melakukan forecasting data co2_ts dengan metode Holt’s Exponential menggunakan fungsi HoltWinters()

co2_hw <- HoltWinters(co2_ts, gamma = F)
co2_hw
#> Holt-Winters exponential smoothing with trend and without seasonal component.
#> 
#> Call:
#> HoltWinters(x = co2_ts, gamma = F)
#> 
#> Smoothing parameters:
#>  alpha: 0.7468631
#>  beta : 0.1106535
#>  gamma: FALSE
#> 
#> Coefficients:
#>         [,1]
#> a 2.38054022
#> b 0.09548633

Buatlah model untuk prediksi co2_ts, nilai alpha = 0.79 & 0.70 nilai beta = 0.15 & 0.01

# alpha 0.79 beta 0.15
co2_hw2 <- HoltWinters(co2_ts, gamma = F, alpha = 0.79, beta = 0.15)
  
# alpha 0.70 beta 0.01
co2_hw3 <- HoltWinters(co2_ts, gamma = F, alpha = 0.7, beta = 0.01)
co2_hw3$fitted
#> Time Series:
#> Start = 1972 
#> End = 2012 
#> Frequency = 1 
#>           xhat     level      trend
#> 1972 0.3492911 0.3306215 0.01866961
#> 1973 0.3741235 0.3553929 0.01873063
#> 1974 0.4079463 0.3890662 0.01888005
#> 1975 0.4227331 0.4038936 0.01883953
#> 1976 0.4345535 0.4157834 0.01877003
#> 1977 0.4721901 0.4532333 0.01895683
#> 1978 0.5817191 0.5618655 0.01985358
#> 1979 0.6624179 0.6419619 0.02045601
#> 1980 0.6812675 0.6608274 0.02044011
#> 1981 0.6744047 0.6542349 0.02016978
#> 1982 0.6867992 0.6667064 0.02009280
#> 1983 0.7036575 0.6835967 0.02006077
#> 1984 0.6957494 0.6759656 0.01978385
#> 1985 0.7145807 0.6948063 0.01977442
#> 1986 0.7486276 0.7287118 0.01991573
#> 1987 0.7503661 0.7306304 0.01973576
#> 1988 0.7475121 0.7280000 0.01951210
#> 1989 0.7724662 0.7529003 0.01956598
#> 1990 0.7654067 0.7461043 0.01930236
#> 1991 0.8263761 0.8066612 0.01971491
#> 1992 0.9501345 0.9293894 0.02074504
#> 1993 1.0618989 1.0402526 0.02164622
#> 1994 1.1424599 1.1202303 0.02222954
#> 1995 1.1641017 1.1418780 0.02222372
#> 1996 1.1707542 1.1486847 0.02206955
#> 1997 1.2608646 1.2381214 0.02274322
#> 1998 1.3635089 1.3399746 0.02353432
#> 1999 1.1454916 1.1243489 0.02114272
#> 2000 1.1768864 1.1556422 0.02124422
#> 2001 1.2464578 1.2247351 0.02172271
#> 2002 1.3589314 1.3363102 0.02262124
#> 2003 1.4178234 1.3948431 0.02298035
#> 2004 1.4539406 1.4308302 0.02311042
#> 2005 1.5166130 1.4931109 0.02350212
#> 2006 1.5343655 1.5109203 0.02344520
#> 2007 1.5346291 1.5114134 0.02321567
#> 2008 1.6124438 1.5886875 0.02375626
#> 2009 1.7432761 1.7184597 0.02481642
#> 2010 1.8542683 1.8285986 0.02566964
#> 2011 1.8188811 1.7938160 0.02506512
#> 2012 2.2078421 2.1791741 0.02866805

Visualisasi data historis dengan hasil forecast

# gunakan autolayer()
co2_ts %>% 
  autoplot()+
  autolayer(co2_hw$fitted[,1], series = "auto")+
  autolayer(co2_hw2$fitted[,1], series = "a: 0.79 b: 0.15") +
  autolayer(co2_hw3$fitted[,1], series = "a: 0.7 b: 0.01")

Model evaluation

# gunakan accuracy()
accuracy(co2_hw$fitted[,1], co2_ts)
#>                  ME      RMSE        MAE      MPE     MAPE        ACF1
#> Test set 0.02267075 0.1162543 0.07175926 1.255456 6.427522 -0.01218677
#>          Theil's U
#> Test set  0.904001
accuracy(co2_hw2$fitted[,1], co2_ts)
#>                  ME      RMSE        MAE     MPE     MAPE       ACF1 Theil's U
#> Test set 0.01951156 0.1166904 0.07164522 1.01725 6.471199 -0.0710277 0.9050397
accuracy(co2_hw3$fitted[,1], co2_ts)
#>                  ME      RMSE        MAE      MPE     MAPE       ACF1 Theil's U
#> Test set 0.03974211 0.1214789 0.07389839 2.660699 6.311941 0.09662694 0.9264768

Kesimpulan: Kalau meninjau RMSE-nya maka model pertama yang paling baik, tetapi kalau meninjau MAPE-nya model ke-3 yang paling baik. Pemilihan metric error tidak ada rule of thumb, yang penting kita bisa menjelaskan error tersebut. Masih ada lagging antara model dengan data aslinya.

rain_ses$mse
#> [1] 17.58641
co2_hw$SSE
#> [1] 0.5541177

3.2.4.1 Summary day 2

  1. Tipe data time series:
  • Additive, hubungan antar komponennya penjumlahan (range seasonalnya tidak berubah terhadap waktu)
  • Multiplicative, hubungan antar komponennya perkalian (range seasonalnya berubah terhadap waktu)
  1. Fungsi decompose() untuk data tipe multiplicative bisa menggunakan parameter type = "multiplicative" yang defaultnya adalah additive.
  2. Ada 2 cara untuk treatment data dengan tipe multiplicative,
  • diubah dengan transformasi log -> diubah ke bentuk additive
  • tidak usah diubah, sesuaikan parameter yang digunakan dalam men-decompose atau melakukan forecast Kedua cara di atas merupakan proses tuning yang bisa dilakukan ketika bertemu data multiplicative.

  1. Seasonal adjustment, proses analisis deskriptif untuk data time series dengan cara menghilangkan pola seasonal pada data kita.

  2. Tujuan dari pembuatan model time series memprediksi kondisi di masa yang akan datang (untuk antisipasi atau planning) berdasarkan pola yang ditangkap di masa lalu (forecast).

  3. Model yang dapat digunakan untuk forecasting:

  • Simple Moving Average (SMA): model paling dasar, hanya bisa digunakan untuk kasus data tanpa seasonal dan trend
  • Exponential Smoothing: melakukan pemobotan dalam menghitung data di masa lalu
    • Simple Exponential Smoothing - dipakai ketika data tidak punya trend dan seasonal
    • Double Exponential Smoothing - dipakai ketika data tidak punya seasonal
    • Triple Exponential Smoothing - dipakai untuk data dengan seluruh komponen (trend, seasonal dan error)
  • ARIMA (Autoregressive Integrated Moving Average)
  1. Simple Moving Average (SMA)
  • Menggunakan perataan bergerak untuk membentuk model, fungsi yang digunakan yaitu SMA() dari package TTR, parameter yang digunakan:
    • x :objek time series yang akan di forecast
    • n: jumlah observasi di masa lalu yang digunakan untuk forecasting (orde)
  1. Exponential Smoothing
  • Fungsi pertama yang dapat digunakan yaitu ets(), parameter yang dapat diisi:
    • y: objek timeseries yang ditraining
    • model: penentuan komponen yang terdapat pada data time series, urutan penulisannya yaitu error, trend, seasonal dengan input yang digunakan *** tanda bintang diisi oleh nilai:
      • A: aditif
      • M: multiplikatif
      • N: none
      • Z: auto pada fungsi ets tidak terdapat default nilai parameter additive atau multiplicative, harus ditentukan oleh user atau dipilih secara random oleh model jika menggunakan nilai Z pada komponen di parameter model
    • alpha: bobot untuk smoothing error (range nilai 0-1)
    • beta: bobot untuk smoothing trend (range nilai 0-1)
    • gamma: bobot untuk smoothing seasonal (range nilai 0-1)
      • nilai 1 berarti observasi yang paling baru diberikan bobot yang paling tinggi dibandingkan dengan data periode lama namun penurunannya sangat signifikan
      • nilai 0 berarti observasi yang paling baru diberikan bobot yang paling tinggi namun perbedaannya dengan data periode lama tidak terlalu signifikan
  • Fungsi kedua yang dapat digunakan yaitu HoltWinters(), parameter yang dapat diisi:
    • x: objek timeseries yang ditraining
    • alpha: bobot untuk smoothing error (range nilai 0-1)
    • beta: bobot untuk smoothing trend (range nilai 0-1)
    • gamma: bobot untuk smoothing seasonal (range nilai 0-1)
      • secara default parameter alpha, beta, dan gamma adalah NULL, jika tidak didefinisikan nilainya, maka model akan mencari nilai parameternya yang paling optimum. Apabila objek time series tidak mengandung trend atau seasonal, parameter beta atau gamma harus diubah menjadi FALSE.
    • seasonal: untuk menentukan perhitungan berdasarkan tipe data time series, secara default akan digunakan additive, nilai yang dapat dimasukkan yaitu additive dan multiplicative

– END OF DAY 2 –

3.2.5 Holt’s Winters Exponential

Holt’s Winters Exponential (Triple Exponential Smoothing) merupakan metode forecasting yang tepat digunakan untuk data yang memiliki efek trend dan seasonal3.

3.2.5.1 Workflow Analisis Time Series

  1. Read data fancy.dat dan simpan dengan nama object souvenir
souvenir <- scan("data_input/fancy.dat")
souvenir
#>  [1]   1664.81   2397.53   2840.71   3547.29   3752.96   3714.74   4349.61
#>  [8]   3566.34   5021.82   6423.48   7600.60  19756.21   2499.81   5198.24
#> [15]   7225.14   4806.03   5900.88   4951.34   6179.12   4752.15   5496.43
#> [22]   5835.10  12600.08  28541.72   4717.02   5702.63   9957.58   5304.78
#> [29]   6492.43   6630.80   7349.62   8176.62   8573.17   9690.50  15151.84
#> [36]  34061.01   5921.10   5814.58  12421.25   6369.77   7609.12   7224.75
#> [43]   8121.22   7979.25   8093.06   8476.70  17914.66  30114.41   4826.64
#> [50]   6470.23   9638.77   8821.17   8722.37  10209.48  11276.55  12552.22
#> [57]  11637.39  13606.89  21822.11  45060.69   7615.03   9849.69  14558.40
#> [64]  11587.33   9332.56  13082.09  16732.78  19888.61  23933.38  25391.35
#> [71]  36024.80  80721.71  10243.24  11266.88  21826.84  17357.33  15997.79
#> [78]  18601.53  26155.15  28586.52  30505.41  30821.33  46634.38 104660.67

Data souvenir terdiri dari 84 observasi penjualan souvenir per bulan dari tahun 1987.

  1. Buat time series object dan simpan dengan nama object souvenir_ts
souvenir_ts <- ts(data = souvenir,
   start = 1987,
   frequency = 12)
  1. Visualisasikan data souvenir_ts

Apakah tipe data time series yang terbentuk dari data souvenir_ts?

souvenir_ts %>% autoplot()

  1. Decomposition

Tujuannya untuk menguraikan tipe-tipe data pada setiap komponen time series yang ada agar tepat dalam penentuan pemilihan model untuk forecast nya.

# menguraikan komponen data time series
decompose(souvenir_ts, type = "multiplicative") %>% autoplot()

  1. Splitting Data
  • Data train: 1987 - 1992 (6 tahun) simpan dengan nama object souvenir_train
  • Data test: 1993 (1 tahun) simpan dengan nama object souvenir_test
# test menggunakan `tail()`
souvenir_test <- tail(souvenir_ts, 12)
# train menggunakan `head()`
souvenir_train <- head(souvenir_ts, -12)
  1. Melakukan model fitting pada data train souvenir_train
# model multiplicative
souvenir_hw <- HoltWinters(x = souvenir_train, seasonal = "multiplicative")
souvenir_hw$fitted
#>               xhat     level      trend    season
#> Jan 1988  2094.394  5366.664 155.793593 0.3792503
#> Feb 1988  4701.893  5893.382 183.618749 0.7737193
#> Mar 1988  6918.015  6299.594 200.316755 1.0643246
#> Apr 1988  4823.164  6600.038 207.827863 0.7084694
#> May 1988  5932.309  6799.474 207.198341 0.8466657
#> Jun 1988  4721.566  6993.792 206.232122 0.6557709
#> Jul 1988  5979.705  7321.603 215.352469 0.7933847
#> Aug 1988  4965.222  7624.169 221.894860 0.6328296
#> Sep 1988  6716.584  7729.235 213.130889 0.8456653
#> Oct 1988  7922.416  7441.726 175.574975 1.0400555
#> Nov 1988  8472.102  6920.929 123.336088 1.2026950
#> Dec 1988 25826.620  8235.210 212.675721 3.0571698
#> Jan 1989  3640.862  8756.045 235.792565 0.4049075
#> Feb 1989  8206.953  9914.048 304.972873 0.8031056
#> Mar 1989 10125.421  9137.020 223.805781 1.0816803
#> Apr 1989  6740.459  9306.986 219.766898 0.7075295
#> May 1989  7596.137  8822.672 166.949735 0.8449896
#> Jun 1989  5786.580  8536.398 132.950841 0.6674758
#> Jul 1989  7448.388  9108.213 165.872538 0.8031399
#> Aug 1989  5848.269  9231.414 162.671518 0.6225480
#> Sep 1989  8591.861 10691.818 260.021939 0.7845130
#> Oct 1989 10391.545 10943.573 259.401798 0.9275701
#> Nov 1989 15536.913 10940.729 239.729212 1.3896490
#> Dec 1989 35906.260 11084.309 232.516470 3.1728210
#> Jan 1990  5047.371 11115.026 217.378346 0.4453927
#> Feb 1990  8607.881 12013.085 268.440140 0.7008804
#> Mar 1990 11892.596 10898.649 164.702571 1.0749542
#> Apr 1990  7381.389 11233.995 177.503575 0.6468378
#> May 1990  8768.922 10868.835 136.795182 0.7967669
#> Jun 1990  7441.297 10500.547  98.906033 0.7020454
#> Jul 1990  8457.641 10492.426  90.877256 0.7991495
#> Aug 1990  7401.647 10437.232  79.919612 0.7037691
#> Sep 1990  8546.767 10801.931 101.282595 0.7838760
#> Oct 1990  9749.347 10702.379  86.216846 0.9036716
#> Nov 1990 14248.067 10299.936  49.559603 1.3766920
#> Dec 1990 35441.009 11273.631 118.884343 3.1109030
#> Jan 1991  5137.566 10798.396  74.316067 0.4725193
#> Feb 1991  6477.552 10644.390  57.188314 0.6052894
#> Mar 1991 11749.083 10697.381  56.873446 1.0925055
#> Apr 1991  6176.699 10084.011   6.594578 0.6121237
#> May 1991  8846.749 11589.632 119.045224 0.7555720
#> Jun 1991  8169.919 11651.558 114.760410 0.6943479
#> Jul 1991 10214.367 12785.542 191.218243 0.7871277
#> Aug 1991  9894.120 13444.995 226.343250 0.7237126
#> Sep 1991 11726.589 14945.765 321.945357 0.7680647
#> Oct 1991 13332.351 15227.414 318.922447 0.8575880
#> Nov 1991 23944.954 15657.416 327.255170 1.4979948
#> Dec 1991 46196.612 15492.951 290.368443 2.9269262
#> Jan 1992  7353.186 15648.657 280.266637 0.4616248
#> Feb 1992  9935.127 16125.741 295.031066 0.6050341
#> Mar 1992 16903.966 16371.775 291.355452 1.0144532
#> Apr 1992 11219.777 15860.851 231.171887 0.6972260
#> May 1992 12416.152 16274.941 244.893597 0.7515906
#> Jun 1992 11484.322 15096.243 138.101790 0.7538442
#> Jul 1992 13198.630 15969.776 193.270698 0.8165930
#> Aug 1992 14197.653 17664.767 305.923408 0.7900450
#> Sep 1992 16055.564 20470.133 493.421047 0.7658799
#> Oct 1992 21857.055 24532.621 761.157460 0.8641277
#> Nov 1992 39906.066 26712.949 867.617651 1.4468907
#> Dec 1992 79593.943 26649.787 797.794417 2.8998528
  1. Melakukan prediksi pada data test souvenir_test dengan menggunakan fungsi forecast() dari library forecast
  • object = model time series
  • h = berapa banyak data kedepan yang akan di forecast
# forecast menggunakan model multiplicative
souvenir_pred <- forecast(object = souvenir_hw, h = 12)

Visualisasi hasil forecast

# visualisasi hasil forecast dari model multiplicative
souvenir_ts %>% 
  autoplot() +
  autolayer(souvenir_hw$fitted[,1], series = "train") +
  autolayer(souvenir_pred$mean, series = "test")

  1. Model evaluation menggunakan data souvenir_test
# model multiplicative
accuracy(souvenir_hw$fitted[,1], souvenir_train) # data train
#>                ME     RMSE      MAE        MPE     MAPE      ACF1 Theil's U
#> Test set 189.1763 2200.411 1512.406 -0.1299785 13.13863 0.1829319 0.3931926
accuracy(souvenir_pred$mean, souvenir_test) # data test
#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -4060.199 4541.688 4060.199 -21.15116 21.15116 0.5148249 0.6278038

Kesimpulan: model yang dibentuk tidak ada lagging dan selisih MAPE antara training dan testing sebesar 8%-an. Modelnya sudah bisa dikatakan baik.

Tunning Model dimana data di transformasikan log

souvenir_log <- log(souvenir_train)
souvenir_log %>% autoplot()

  • Modeling menggunakan data yang sudah di transformasi
# souvenir_hw_log <- HoltWinters(souvenir_train) # salah memasukkan data
souvenir_hw_log <- HoltWinters(souvenir_log)
  • Forecasting

Note: Hasil forecast dari data yang sudah di transformasikan log harus dikembalikan kembali ke bentuk data aslinya menggunakan exp()

souvenir_pred_log <- forecast(souvenir_hw_log, h = 12)
souvenir_pred_log_exp <- exp(souvenir_pred_log$mean)
  • Model evaluation
accuracy(exp(souvenir_hw_log$fitted[,1]), exp(souvenir_log)) # data train
#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -55.25421 2302.798 1552.995 -1.715425 13.25152 0.1180214  0.404583
accuracy(souvenir_pred_log_exp, souvenir_test) # data test
#>                 ME     RMSE      MAE      MPE    MAPE        ACF1 Theil's U
#> Test set -5769.926 6138.801 5769.926 -25.9165 25.9165 -0.04706218 0.6797808
  • Visualisasi
souvenir_ts %>% 
  autoplot() +
  autolayer(exp(souvenir_hw_log$fitted[,1]), series = "train") +
  autolayer(souvenir_pred_log_exp, series = "test")

3.2.6 Knowledge Check

  1. Model forecasting manakah yang cocok digunakan untuk data yang tidak memiliki komponen trend dan seasonal? Dapat memilih lebih dari satu jawaban.
  • Simple Moving Average
  • Simple Exponential Smoothing
  • Holt Exponential Smoothing
  • Holt-Winter’s Exponential Smoothing
  1. Grafik di bawah menunjukkan data tahunan total penumpang pesawat di Australia dari tahun 1970 sampai 2016. Berdasarkan pola data tersebut, model forecasting manakah yang cocok untuk data ini?
  • Simple Moving Average
  • Simple Exponential Smoothing
  • Holt Exponential Smoothing
  • Holt-Winters Exponential Smoothing
  1. Perhatikan grafik decompose dibawah ini. Berdasarkan grafik tersebut, tipe time series dan model forecasting apakah yang sesuai? Berikan alasannya…

Triple exponential smoothing/Holt-Winters Exponential Smoothing karena data tersebut punya error, trend, dan seasonal

3.3 Dive Deeper Exponential Smoothing

Buatlah satu buah Rmd baru lalu kerjakan kasus berikut ini: Gunakanlah data nybirth.csv. Buatlah model time series yang tepat untuk data tersebut, analisislah pola time series yang ada pada data tersebut, apakah multiplicative atau additive? Lakukanlah forecasting dan lihatlah model evaluationnya. Bandingkan hasil evaluasi serta visualisasinya lalu tulislah kesimpulan Anda. Gunakanlah beberapa nilai bobot sebagai tuning pemodelan Anda. Lampirkan hasil visualisasi dan evaluasinya di Online Class Guide (berikan keterangan parameter bobot yang Anda gunakan)

4 ARIMA (Autoregressive Integrated Moving Average)

ARIMA adalah gabungan antara dua metode, yaitu Auto Regressive (AR) dan Moving Average (MA). I nya menjelaskan Integrated. Tujuan utama dari ARIMA adalah melakukan autocorrelation pada data.

4.1 Stasionarity

Stasionarity time series memiliki arti bahwa pada data time series yang kita miliki tidak memiliki trend maupun seasonal dan memiliki variansi konstan4.

Berdasarkan gambar dibawah ini, manakah yang termasuk data time series yang stasioner?

Apabila data yang kita miliki belum stasioner, maka kita bisa melakukan differencing.

4.2 Differencing

Jadi untuk membuat datanya stasioner, cara yang paling umum digunakan adalah dengan melakukan differencing diff, yaitu mengurangi data saat ini dengan data sebelumnya5. Terkadang, tergantung pada kompleksitas data, jumlah differencing bisa lebih dari 1 kali.

Pemahaman konsep mengapa perlu differencing???

  • data sebelum dilakukan differencing
data.frame(
  xt = c(AirPassengers),
  lag1 = lag(x = as.vector(AirPassengers), n = 1),
  lag2 = lag(x = as.vector(AirPassengers), n = 2),
  lag3 = lag(x = as.vector(AirPassengers), n = 3),
  lag4 = lag(x = as.vector(AirPassengers), n = 4),
  lag5 = lag(x = as.vector(AirPassengers), n = 5) 
) %>% 
  na.omit() %>% 
  GGally::ggcorr(label = T)

  • data setelah dilakukan differencing
data.frame(
  xt = c(diff(AirPassengers)),
  lag1 = lag(x = as.vector(diff(AirPassengers)), n = 1),
  lag2 = lag(x = as.vector(diff(AirPassengers)), n = 2),
  lag3 = lag(x = as.vector(diff(AirPassengers)), n = 3),
  lag4 = lag(x = as.vector(diff(AirPassengers)), n = 4),
  lag5 = lag(x = as.vector(diff(AirPassengers)), n = 5) 
) %>% 
  na.omit() %>% 
  GGally::ggcorr(label = T)

Ide utama melakukan differencing adalah agar ketika melakukan prediksi, tidak ada multicolinearity terhadap data-data sebelumnya.

Model ARIMA: Arima memiliki bentuk model ARIMA(p,d,q) dimana6 :

4.3 AR(p) : Auto Regressive

Auto regressive di ARIMA artinya kita membuat model linear regression berdasarkan lag dari datanya sebagai prediktor7 atau bisa dikatakan melakukan regresi terhadap dirinya sendiri. Nilai p berarti berapa banyak data yang akan dipakai ketika melakukan auto regressive.

Model autoregressive dapat dituliskan sebagai berikut.

\[ Y_t = \alpha + \beta_1Y_{t-1}+\beta_2Y_{t-2}+...+\beta_pY_{t-p}+ \epsilon_1 \]

dimana, \(Y_{t-1}\) adalah lag1, \(\beta_1\) adalah coefficient dari lag1, dan \(\alpha\) adalah intercept.

Untuk mencari order p untuk model AR, kita dapat melihat dari plot PACF (Partial Autocorrelation Function).

autokorelasi pada saat t dan t - k tanpa dipengaruhi oleh korelasi periode ditengahnya (intermediate lags). Misalnya autokorelasi untuk lag(2) maka bisa dikatakan ingin mencari autokorelasi antara hari rabu ke jumat secara langsung.

4.4 I(d) : Integrated

Integrated adalah berapa kali data dilakukan differencing untuk membuat suatu data stationer. Nilai d dapat diketahui dengan mencari tahu berapa kali differencing yang dilakukan pada data.

4.5 MA(q) : Moving Average

Moving Average dalam ARIMA artinya kita melakukan rata-rata berjalan terhadap data time series itu sendiri8. MA digunakan untuk melakukan smoothing error. Nilai q berarti berapa banyak data yang diperlukan untuk smoothing error menggunakan moving average.

Moving Average Model (MA model).

\[ Y_t = \alpha + \epsilon_t + \phi_1\epsilon_{t-1}+\phi_2\epsilon_{t-2}+...+\phi_q\epsilon_{t-q} \]

Untuk mencari order q untuk model MA, kita dapat melihat dari plot ACF (Auto correlation Function)

mencari autokorelasi pada saat t dan t - k, dan dipengaruhi oleh periode ditengahnya. Misal: autokorelasi untuk lag(2) maka bisa dikatakan, akan dihitung nilai autokorelasi di hari jumat ke rabu, namun mempertimbangkan korelasi antara jumat ke kamis dan kamis ke rabu.

Tahap analisis menggunakan metode ARIMA:

  1. Objek time series
  2. Melakukan pengujian stasioner menggunakan adf.test

H0 : data tidak stasioner H1 : data stasioner

kita mengharapkan p-value < alpha, sehingga pola data stasioner

  1. Catat informasi berapa kali dilakukan differencing sampai data menjadi stasioner
  2. Lakukan pembuatan model arima menggunakan auto.arima(objek_ts, seasonal = F).
  3. Tuning model dengan mencari ordo (p) dan (q) secara manual.
  • plot acf vs pacf untuk mencari ordo (p) dan (q) secara manual
  • acf (autocorrelation function), pacf (partial autocorrelation function).
  • acf: mencari autokorelasi pada saat t dan t - k, dan diperngaruhi oleh periode ditengahnya. Misal: autokorelasi untuk lag(2) maka bisa dikatakan, akan dihitung nilai autokorelasi di hari jumat ke rabu, namun mempertimbangkan korelasi antara jumat ke kamus dan kamis ke rabu.
  • pacf: mencari autokorelasi pada saat t dan t - k tanpa dipengaruhi oleh korelasi periode ditengahnya.
  1. Memperhatikan garis keluar dari batas pada lag ke berapa saja. (signifikan pada lag ke berapa aja)
  • membuat model perbandingan dari beberapa kemungkinan ordo nya
  • mencari model terbaik dari nilai AIC yang terkecil
  1. Model evaluation dengan nilai MAPE
  2. Cek asumsi time series
  • normality residual assumption
  • autocorrelation residual assumption

4.6 Non-Seasonal ARIMA(p,d,q)

Kita akan coba melakukan forecast pada data environment_1970f.csv, dimana yang akan kita ubah kedalam bentuk object time series adalah data Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.

  1. Read Data
nitro <- read.csv("data_input/environment_1970f.csv")
head(nitro)
  1. Membuat object time series untuk kolom Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.
nitro_ts <- ts(data = nitro$Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.,
               start = 1970,
               frequency = 1)
  1. Melakukan pengecekan object ts stasioner

Stasioner adalah data yang bergerak disekitar rata-rata globalnya.

  • melihat berdasarkan visual
nitro_ts %>% autoplot()

> justifikasi awal berdasarkan visualisasi, apakah ada indikasi terdapat trend atau seasonal? kurang jelas

  • melakukan pengujian menggunakan adf.test()
adf.test(nitro_ts)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  nitro_ts
#> Dickey-Fuller = -2.7827, Lag order = 3, p-value = 0.2645
#> alternative hypothesis: stationary

Notes H0 : data tidak stasioner H1 : data stasioner

dengan alpha = 0.05 p-value > alpha artinya gagal tolak H0, jadi datanya belum stasioner.

  1. Melakukan differencing data agar stasioner
  • mencoba melakukan pengujian data no_ts dengan satu kali differencing, gunakan fungsi diff(). Beberapa parameter yang dapat digunakan: differences: jumlah differencing dilakukan lag: lompatan differencing
nitro_ts %>% 
  diff() %>% 
  adf.test()
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -6.1019, Lag order = 3, p-value = 0.01
#> alternative hypothesis: stationary

Data kita perlu dilakukan differencing 1x hingga dia stasioner.

artinya nilai d pada ARIMA(p,d,q) = 1

# contoh diff lebih dari 1 kali
nitro_ts %>% 
  diff(differences = 2)

ARIMA(p,1,q)

4.6.0.1 Summary day 3

  1. Triple Exponential Smoothing: digunakan untuk data yang mengandung komponen error, trend, dan seasonal. Masih menggunakan fungsi ets() dan HoltWinters().
  2. Workflow pengolahan data time series:
  • read data
  • buat object time series
  • EDA (lihat bentuknya, tentukan tipe data time seriesnya)
  • Decompose (melihat apakah frequency yang digunakan sudah benar)
  • Splitting data ke train dan test
  • Model fitting dengan data train
  • Forecast data test
  • Evaluasi
  • Kesimpulan
  1. ARIMA (Auto Regressive Integrated Moving Average) -> Bentuknya yaitu ARIMA(p,d,q)
  • melihat apakah datanya stasioner atau tidak, jika tidak stasioner lakukan differencing.
  • pengecekan differencing menggunakan fungsi adf.test() jika p-value > alpha maka belum stasioner
  • pengaplikasian differencing menggunakan fungsi diff()
  • dari jumlah differencing yang dilakukan, kita bisa mendapatkan komponen d dari ARIMA(p,d,q)

– END OF DAY 3 –

  1. Model Fitting
  • membuat model arima secara otomatis menggunakan auto.arima()
nitro_auto <- auto.arima(nitro_ts, trace = T)
#> 
#>  ARIMA(2,1,2) with drift         : Inf
#>  ARIMA(0,1,0) with drift         : 1062.439
#>  ARIMA(1,1,0) with drift         : 1058.421
#>  ARIMA(0,1,1) with drift         : 1044.127
#>  ARIMA(0,1,0)                    : 1060.239
#>  ARIMA(1,1,1) with drift         : Inf
#>  ARIMA(0,1,2) with drift         : Inf
#>  ARIMA(1,1,2) with drift         : Inf
#>  ARIMA(0,1,1)                    : 1042.989
#>  ARIMA(1,1,1)                    : 1045.207
#>  ARIMA(0,1,2)                    : 1045.127
#>  ARIMA(1,1,0)                    : 1056.115
#>  ARIMA(1,1,2)                    : Inf
#> 
#>  Best model: ARIMA(0,1,1)
nitro_auto
#> Series: nitro_ts 
#> ARIMA(0,1,1) 
#> 
#> Coefficients:
#>           ma1
#>       -0.8244
#> s.e.   0.0842
#> 
#> sigma^2 = 3.21e+09:  log likelihood = -519.34
#> AIC=1042.68   AICc=1042.99   BIC=1046.16

ARIMA(p,d,q) -> ARIMA(0,1,1)

  • melakukan tunning model untuk p dan q nya menggunakan acf dan pacf

p = melihat lag yang keluar dari plot PACF (partial autocorrelation function) q = melihat lag yang keluar dari plot ACF (autocorrelation function)

ACF untuk melihat korelasi suatu data dengan mempertimbangkan data tengahnya. misal ingin melihat korelasi hari senin dengan hari rabu, kita harus mempertimbangkan:

  • korelasi hari senin - selasa
  • korelasi hari senin - rabu
  • korelasi hari selasa - rabu

PACF untuk melihat korelasi data tanpa mempertimbangkan data tengahnya. misal ingin melihat korelasi hari senin dengan rabu, cukup melihat nilai korelasi hari senin - rabu saja

melihat plot pacf dan acf untuk menemukan nilai p dan q

# melihat plot acf dan pacf menggunakan `tsdisplay()`
nitro_ts %>% 
  diff() %>% 
  tsdisplay()

nilai p yang diperoleh dari 5 lag pertama yang keluar dari plot pacf

1, 2, 4, 0

nilai q yang diperoleh dari 5 lag pertama yang kelaur dari plot acf

1, 0

Kombinasi model ARIMA yang terbentuk9 : ARIMA(p,d,q)

  • ARIMA(0,1,1) -> threshold dari auto.arima
  • ARIMA(1,1,1)
  • ARIMA(2,1,1)
  • ARIMA(4,1,0)
  • ARIMA(0,1,0)

parameter pada fungsi Arima()

  • y = object time series
  • order = nilai (p,d,q)
nitro_111 <- Arima(y = nitro_ts, order = c(1,1,1))
nitro_211 <- Arima(y = nitro_ts, order = c(2,1,1))
nitro_410 <- Arima(y = nitro_ts, order = c(4,1,0))
  1. Goodness of Fit (Pemilihan model terbaik)

Untuk pemilihan model terbaik kita bisa melihat nilai AIC dari masing-masing model yang paling kecil

nitro_auto$aic;nitro_111$aic;nitro_211$aic;nitro_410$aic
#> [1] 1042.681
#> [1] 1044.575
#> [1] 1044.902
#> [1] 1046.4
  1. Model Evaluation
  • bisa menggunakan fungsi accuracy() dari package forecast untuk melihat MAPE
accuracy(nitro_auto)
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 7917.811 55326.04 33952.52 -2.335425 25.68435 0.7884324
#>                     ACF1
#> Training set 0.003779372

4.7 Seasonal ARIMA (SARIMA(p,d,q)(P,D,Q)[seasonal])

Seasonal arima adalah metode arima dimana object time series yang ada memiliki pola seasonal10. Tahapan dalam melakukan pemodelan menggunakan SARIMA sama seperti saat membuat pemodelan ARIMA.

Order Model pada SARIMA(p,d,q)(P,D,Q)[seasonal]

ARIMA: (p, d, q)

SARIMA: (p, d, q)(P, D, Q)[seasonal]

  • (p, d, q) untuk index utama

  • (P, D, Q) untuk index seasonal

  • seasonal merupakan frequency

Jika kita menggunakan model SARIMA, diperlukan langkah tambahan ketika melakukan differencing untuk seasonal data, yaitu differencing sesuai pola berulang/frequency data time series untuk menghilangkan pola musiman.

Pada bagian ini, kita akan menggunakan data DAUTONSA.csv yaitu data penjualan motor perbulan dari Jan 1967 - Jun 2019.

  1. Read data
sales <- read.csv("data_input/DAUTONSA.csv")
head(sales)
sales %>% 
  mutate(DATE = ymd(DATE)) %>% 
  ggplot(aes(x = DATE, y = DAUTONSA)) +
  geom_line()

  • untuk kebutuhan analisis, datanya akan kita subset dari tahun 2011 keatas
sales <- sales %>% 
  mutate(DATE = ymd(DATE)) %>% 
  filter(DATE >= ymd("2011-01-01"))
sales %>% head()
  1. Membuat object time series, dimana periode bermula dari Januari 2011
# bu dessy
sales_ts <- ts(data = sales$DAUTONSA, 
               start = c(2011,1), 
               frequency = 12)

Melihat visualisasi data secara umum

sales_ts %>% decompose() %>% autoplot()

Informasi dari object sales_ts:

  • dari data sales_ts apakah terdapat seasonality? ada seasonality
  • dari data sales_ts bagaimana kondisi trendnya? ada trend, naik dari 2011 sampai sekitar 2015 lalu turun
  • dari data sales_ts apakah stasioner? tidak stasioner
tail(sales)
  1. Splitting data

Train data : diambil 7 tahun pertama Test data : sisanya

sales_train <- head(sales_ts, 7*12)
sales_test <- tail(sales_ts, -7*12)
  1. Melakukan pengecekan stasionaritas data train
sales_train %>% adf.test()
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -2.9359, Lag order = 4, p-value = 0.192
#> alternative hypothesis: stationary

p-value > 0.05 (alpha), belum stasioner

Differencing based on frequency, cek apakah sudah stasioner atau belum

Notes: Untuk differencing seasonal, kita akan gunakan parameter lag dimana nilainya mengikuti frekuensi data yang kita miliki

Tahapan differencing pada kasus SARIMA dan langkah selanjutnya jika data belum stasioner:

  1. Differencing seasonal 1x, trend 1x
  2. Differencing seasonal 1x, trend 2x
  3. Transformasi log object ts dengan differencing di langkah kedua
# melakukan satu kali differencing data seasonal dan trend
sales_train %>% 
  diff(lag = 12) %>% # differencing untuk komponen seasonal
  diff() %>% # differencing untuk komponen trend
  adf.test()
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -4.7557, Lag order = 4, p-value = 0.01
#> alternative hypothesis: stationary
data.frame(
  xt = c(diff(AirPassengers)),
  lag1 = lag(x = as.vector(diff(AirPassengers)), n = 1),
  lag2 = lag(x = as.vector(diff(AirPassengers)), n = 2),
  lag3 = lag(x = as.vector(diff(AirPassengers)), n = 3),
  lag4 = lag(x = as.vector(diff(AirPassengers)), n = 4),
  lag5 = lag(x = as.vector(diff(AirPassengers)), n = 5) 
)

bentuk umum: SARIMA(p,d,q)(P,D,Q)[seasonal]

d = 1 D = 1

SARIMA(p,1,q)(P,1,Q)[12]

  1. Model Fitting
  • membuat model SARIMA menggunakan auto.arima()

SARIMA(p,d,q)(P,D,Q)[seasonal/freq]

# sales_auto <- auto.arima(sales_train, seasonal = T) # case kalau komponen seasonal ga terdeteksi model
sales_auto <- auto.arima(sales_train)
sales_auto
#> Series: sales_train 
#> ARIMA(0,1,1)(0,1,0)[12] 
#> 
#> Coefficients:
#>           ma1
#>       -0.6468
#> s.e.   0.0900
#> 
#> sigma^2 = 892.7:  log likelihood = -341.71
#> AIC=687.41   AICc=687.59   BIC=691.94
  • tuning model SARIMA secara manual

Menentukan nilai (p,d,q)(P,D,Q)seasonal

Langkah 1: Tentukan D.

Langkah 2: Tentukan d.

Langkah 3: Tentukan P dan Q. Untuk index seasonal, amati berapa banyak lag yang keluar di ACF/PACF ketika lag = freq/kelipatan freq.

# untuk menampilkan banyak lag yang muncul di plot acf dan pacf bisa tambahkan `lag.max`
sales_train %>% 
  diff(lag = 12) %>% 
  diff() %>% 
  tsdisplay(lag.max = 36)

Untuk seasosal menentukan P dan Q kita melihat kelipatan lag berdasarkan frequenci data kita.

  • P : 0
  • D : 1
  • Q : 0

untuk keseluruhan data menentukan p dan q kita melihat 5 lag pertama yang keluar

  • p: 1, 2, 0
  • d: 1
  • q: 1, 0

Kombinasi model yang mungkin terbentuk adalah : SARIMA(p,d,q)(P,D,Q)[seasonal]

  • SARIMA(0,1,1)(0,1,0)[12] -> auto.arima
  • SARIMA(0,1,0)(0,1,0)[12] (pa rizky)
  • SARIMA(2,1,0)(0,1,0)[12] (bu dessy)
  • SARIMA(1,1,1)(0,1,0)[12] (pa eri)
  • SARIMA(1,1,1)(2,1,0)[12] (pa jeremy) -> kurang tepat karena nilai P tidak dalam opsi kita
  • SARIMA(2,1,1)(0,1,0)[12] (bu cherisa)
rizky <- Arima(sales_train, order = c(0,1,0), seasonal = c(0,1,0))
dessy <- Arima(sales_train, order = c(2,1,0), seasonal = c(0,1,0))
eri <- Arima(sales_train, order = c(1,1,1), seasonal = c(0,1,0))
# jeremy <- Arima(sales_train, order = c(1,1,1), seasonal = c(2,1,0))
cherisa <- Arima(sales_train, order = c(2,1,1), seasonal = c(0,1,0))

Goodness of fit menggunakan nilai AIC

sales_auto$aic; rizky$aic; dessy$aic; eri$aic; cherisa$aic
#> [1] 687.4144
#> [1] 709.6523
#> [1] 688.9844
#> [1] 689.4116
#> [1] 689.9545
  1. Melakukan forecasting untuk 18 bulan kedepan
# forecast, h = 18
sales_fc <- forecast(sales_auto, h = 18)

Visualisasi

# menggunakan autoplot dan autolayer
sales_ts %>% 
  autoplot()+
  autolayer(sales_auto$fitted, series = "training") +
  autolayer(sales_fc, series = "forecast")

  1. Model evaluation
  • menggunakan accuracy() dari package forecast dan melihat nilai MAPEnya
# accuracy
# accuracy(sales_auto)
# accuracy(sales_fc$mean, sales_test)
accuracy(sales_fc, sales_test)
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -4.11363 27.27477 19.30697 -1.132620 4.466261 0.4645816
#> Test set     14.06787 27.57037 23.20688  4.370293 7.204645 0.5584247
#>                     ACF1 Theil's U
#> Training set -0.02948427        NA
#> Test set      0.08457512  0.515082

5 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 ini11:

  1. Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
  2. Residual memiliki rata-rata 0.

Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.

5.1 No-autocorrelation residual

Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time series bisa menggunakan uji Ljung-box dengan menggunakan fungsi Box.test(residual model, type = "Ljung-Box) yang akan menghasilkan sebuah nilai p-value. Hipotesis yang digunakan yaitu:

\(H_0\): residual has no-autocorrelation \(H_1\): residual has autocorrelation

yang diinginkan p-value > 0.05 (alpha), no-autocorrelation

# menggunakan Ljung-Box test
Box.test(sales_auto$residuals, type = "Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  sales_auto$residuals
#> X-squared = 0.075662, df = 1, p-value = 0.7833

p-value (0.7833) > 0.05 (alpha) berarti gagal tolak H0 -> residual tidak berautokorelasi

Apabila dari hasil test tersebut masih terdapat autokorelasi, maka kita akan melakukan tuning kembali terhadap model. Pada kasus model exponential smoothing, update nilai bobot dapat dilakukan untuk tuning, sedangkan pada kasus ARIMA dan SARIMA, kita akan mengganti nilai P, Q, p, atau q. Pada kasus ARIMA/SARIMA ini, pemilihan nilai yang diupdate dapat dicek menggunakan visualisasi plot residual acf dengan fungsi acf(residual model) dan pacf(residual model)

acf(sales_auto$residuals, lag.max = 36)

pacf(sales_auto$residuals, lag.max = 36)

Note: untuk p dan q cukup melihat sejauh 4-5 lag apakah ada garis yang keluar atau tidak, sedangkan untuk P dan Q lihat pada 3x kelipatan frequency. Jika ada nilai yang keluar pada area yang kita tinjau tersebut, maka nilai tersebut yang kita update.

5.2 Normality of residual

Untuk mengecek normality residual pada hasil forecasting time series kita bisa melakukan uji normality (shapiro test) dengan menggunakan fungsi shapiro.test(residual model). Hipotesis yang digunakan yaitu:

\(H_0\): residual menyebar normal \(H_1\): residual tidak menyebar normal

yang diinginkan p-value > 0.05 (alpha), residual menyebar normal

shapiro.test(sales_auto$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  sales_auto$residuals
#> W = 0.96645, p-value = 0.02739

p-value < 0.05 sehingga tolak H0, residual tidak menyebar normal -> tuning model

Note : Untuk mengecek asumsi cukup di model yang memang paling baik.

6 [Additional] STLM (Seasonal Trend with Loess Model)

Apabila dalam decompose biasa, dalam mendapatkan komponen trend dengan cara central moving average (CMA) dimana secara konsep setiap data yang ingin dirata-ratakan diberikan bobot yang sama sesuai ordo yang ditetapkan. Karena merata-ratakan data tengahnya hasilnya kita kehilangan data awal dan data akhirnya, sehingga ada beberapa informasi yang hilang. Ada salah satu cara untuk mendapatkan decompose data namun tetap mempertahankan informasi dari seluruh data yang kita miliki yaitu dengan menggunakan STL(Seasonal Trend with Loess). 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. Kekurangan dari STL hanya bisa melakukan decompose pada additive data, apabila terdapat multiplicative data dapat menggunakan transformasi log()12.

Untuk memodelkan hasil STL, kita bisa menerapkan STLM(Seasonal Trend with Loess Model) dimana kita bisa menerapkan metode exponential smoothing (ETS) dan ARIMA. Selain itu, STLM dapat digunakan sebagai alternative cara untuk menangkap seasonal yang belum bisa ditangkap oleh metode ETS dan ARIMA biasa.

6.1 Decompose menggunakan stl()

Dengan menggunakan data rata-rata penjualan souvenir per bulan yang kita miliki, kita kan coba melakukan decompose dengan stl() dan memodelkan dengan menggunakan fungsi stlm() dengan menerapkan metode arima.

Perlu mengatur parameter s.window untuk mendefinisikan pola seasonal yang ingin ditangkap dengan melakukan transformasi log pada data train dan test.

souvenir_train %>% decompose(type = "multiplicative") %>% autoplot()

log(souvenir_train) %>% stl(s.window = 12) %>% autoplot()

6.2 Modeling

Dengan menggunakan informasi data train sales_train dan test sales_test, kita akan coba untuk membuat model dengan menggunakan metode stlm() dengan mengatur parameter method = "ets".

Parameter stlm() yang harus di atur:

  • y : object time series

  • s.window : pola seasonal yang ingin ditangkap

  • method : metode forecast yang akan digunakan, tersedia ets dan arima

  • modeling

model_stlm <- stlm(y = log(souvenir_train), s.window = 12, method = "ets")

6.3 Forecasting

Melakukan forecast untuk 12 bulan mendatang

forecast_stlm <- forecast(model_stlm, h = 12)

Visualisasi hasil forecast

souvenir_ts %>% autoplot(series = "actual") +
  autolayer(souvenir_test, series = "test data") +
  autolayer(exp(forecast_stlm$mean), series = "forecast") +
  scale_color_manual(values = c("black", "blue", "firebrick"))

6.4 Model Evaluation

Melakukan evaluasi model dengan menggunakan fungsi MAPE() dari package MLmetrics

library(MLmetrics)
MAPE(y_pred = exp(forecast_stlm$mean), y_true = souvenir_test)*100
#> [1] 16.54777
accuracy(exp(forecast_stlm$mean), souvenir_test)
#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -2916.257 4342.346 3665.434 -13.97002 16.54777 0.5238158 0.4981687
accuracy(souvenir_pred$mean, souvenir_test)
#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -4060.199 4541.688 4060.199 -21.15116 21.15116 0.5148249 0.6278038

6.4.0.1 Summary day 4

  1. Workflow forecasting dengan metode non seasonal ARIMA:
  • read data
  • cek karakteristik (syarat) data time series
  • buat object time series dan EDA (visualisasikan data)
  • decomposition
  • splitting data
  • lakukan pengujian stasioner dengan adf.test()
  • lakukan differencing jika dibutuhkan, catat berapa kali differencing dilakukan
  • buat model auto arima
  • tuning model dengan mencari nilai p dan q menggunakan plot ACF dan PACF
  • bandingkan AIC antara beberapa model
  • forecast model
  • visualisasi hasil dan evaluasi
  • cek asumsi
  1. Metode ARIMA untuk data time series dengan pola seasonal menggunakan Seasonal ARIMA (SARIMA). Tahapan pada metode SARIMA hampir sama dengan non seasonal ARIMA, tetapi diperlukan langkah tambahan yaitu melakukan differencing komponen seasonal data.
  2. Pengecekan asumsi model time series menggunakan residual dari hasl modeling. Nilai residual yang diharapkan:
  • Residual yang tidak ber-autokorelasi, jika ada korelasi maka artinya ada informasi yang tertinggal. Pengujian menggunakan Box.test(residual_model, type = "Ljung-Box)
  • Nilai rata-rata residual bernilai 0 (Berdistribusi normal). Pengujian menggunakan shapiro.test(residual_model)
  1. Evaluasi model menggunakan error, berikut beberapa contoh metrics yang dapat digunakan beserta nilai yang ditolerir (subjektif):
  • MAPE (kurang lebih 10%)
  • RMSE (tergantung nilai data kita) Setelah dievaluasi, bandingkan metrics yang sama antara metrics model yang ditraining dengan metrics forecast dan cek overfittingnya. Batas overfit subjektif, tetapi umumnya jika beda metrics antara training dan testing <10% dapat dikatakan tidak overfit.

8 Inclass Questions

  1. Cara buang data weekend?
library(lubridate)
library(padr)
library(zoo)
library(dplyr)

the_date <- seq.Date(from = as.Date("2022-11-01"),
                     to = as.Date("2022-11-30"),
                     by = "day")
value <- 1:30

df <- data.frame(date = the_date, value = value)

df_no_weekend <- df %>%
  mutate(day = wday(date, label = T)) %>% # memunculkan nama hari
  filter(!(day %in% c("Sat", "Sun"))) %>% # memilih selain sabtu dan minggu
  select(-day) # menghilangkan kolom yang tidak dibutuhkan

df_no_weekend

  1. Simple Exponential Smoothing↩︎

  2. Trend Methods↩︎

  3. Holt’s Exponential Smoothing↩︎

  4. Stationarity and Differencing↩︎

  5. Explanation diff() in R↩︎

  6. Non Seasonal Model↩︎

  7. Auto Regression Model↩︎

  8. Moving Average Model↩︎

  9. ARIMA in R↩︎

  10. Seasonal ARIMA↩︎

  11. Residuals Diagnostic↩︎

  12. STL Decomposition↩︎