1 Time Series and Forecasting

Time series data merupakan data yang berhubungan dengan waktu yang memiliki interval waktu yang sama

Forecasting merupakan suatu metode untuk memprediksi/meramalkan data di masa depan

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.

1.1 Karakteristik data time series

1. Data harus terurut berdasarkan periode waktunya

dat <- data.frame(year = c(2011, 2012, 2014, 2016, 2015, 2013),
                  value = c(1,2,3,4,5,6))
dat

#>   year value
#> 1 2011     1
#> 2 2012     2
#> 3 2014     3
#> 4 2016     4
#> 5 2015     5
#> 6 2013     6

dat %>% 
  arrange(year)

#>   year value
#> 1 2011     1
#> 2 2012     2
#> 3 2013     6
#> 4 2014     3
#> 5 2015     5
#> 6 2016     4

2. Tidak boleh ada waktu atau periode yang terlewat/bolong

dat <- data.frame(date = as.Date(c("2020-02-01","2020-02-03","2020-02-04")),
                                 value = c(23,34,20))
dat

#>         date value
#> 1 2020-02-01    23
#> 2 2020-02-03    34
#> 3 2020-02-04    20

library(padr)
dat <- dat %>% pad()

3. Data tidak boleh ada yang missing

mengisi dengan hasil forecast
mengisi dengan nilai rata-rata yang dekat dengan nilai yang missing

library(zoo)
dat %>% mutate(value = na.fill(value, fill = "extend"))

#>         date value
#> 1 2020-02-01  23.0
#> 2 2020-02-02  28.5
#> 3 2020-02-03  34.0
#> 4 2020-02-04  20.0

Selain menggunakan function na.fill dari package zoo, kita juga dapat menggunakan berbagai fungsi imputasi dari package imputeTS untuk mengisi nilai NA pada objek time series Package imputeTS

1.2 Tahapan membuat object time series

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

ts(data, start, frequency)

data = data yang akan kita prediksi
start = waktu awal mula data yang akan diprediksi
frequency = pola berulang dari data

#>      Data Minute  Hour      Day       Week        Month            Year
#> 1   Daily      -     -        -          7          7*4          7*4*12
#> 2  Hourly      -     -       24       24*7       24*7*4       24*7*4*12
#> 3 Minutes      -    60    60*24    60*24*7    60*24*7*4    60*24*7*4*12
#> 4 Seconds     60 60*60 60*60*24 60*60*24*7 60*60*24*7*4 60*60*24*7*4*12

Dive deeper :

data jam -> pola harian -> freq = 24
data harian -> pola mingguan -> freq = 7
data harian -> pola weekday -> freq = 5
data jam -> pola mingguan -> freq = 24*7
data bulanan -> pola tahunan -> freq = 12
data bulanan -> pola kuartalan -> freq = 12*4
data tahunan -> pola tahunan -> freq = 1

Langkah awal dalam membuat object time series:

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.

read data

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

#>   year CO2.emissions..kt. CO2.emissions..metric.tons.per.capita.
#> 1 1970           35822.92                              0.3119519
#> 2 1971           38987.54                              0.3306215
#> 3 1972           43340.27                              0.3580080
#> 4 1973           49134.13                              0.3954703
#> 5 1974           51260.99                              0.4021567
#> 6 1975           53963.57                              0.4128050
#>   Methane.emissions..kt.of.CO2.equivalent.
#> 1                                   126665
#> 2                                   119293
#> 3                                   195338
#> 4                                   136928
#> 5                                   122570
#> 6                                   121065
#>   Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.
#> 1                                                         51908.88
#> 2                                                         49004.80
#> 3                                                         79261.42
#> 4                                                         56928.71
#> 5                                                         50997.79
#> 6                                                         50922.15
#>   Other.greenhouse.gas.emissions..HFC..PFC.and.SF6..thousand.metric.tons.of.CO2.equivalent.
#> 1                                                                                 130405.80
#> 2                                                                                  94658.48
#> 3                                                                                 422459.50
#> 4                                                                                 174340.90
#> 5                                                                                 110213.86
#> 6                                                                                 109630.88
#>   Total.greenhouse.gas.emissions..kt.of.CO2.equivalent.
#> 1                                              339677.1
#> 2                                              294004.9
#> 3                                              733432.9
#> 4                                              409241.7
#> 5                                              329453.2
#> 6                                              329119.3

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).

Dari data co2 ini, kita akan menggunakan 2 kolom sebagai informasi 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.

mengetahui range atau periode waktu data time series

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

# range waktu
range(co2$year)

#> [1] 1970 2012

mengetahui frequency atau banyak data dalam satu pola musiman

untuk mengetahui frequency data yang dimiliki dapat diinspect berdasarkan:

data yang disusun per periode apa?
pola yang ingin dilihat apa (harian/mingguan/bulanan/kuartalan/tahunan)?

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)
co2_ts

#> Time Series:
#> Start = 1970 
#> End = 2012 
#> Frequency = 1 
#>  [1] 0.3119519 0.3306215 0.3580080 0.3954703 0.4021567 0.4128050 0.4612390
#>  [8] 0.6002978 0.6677802 0.6601457 0.6426495 0.6634071 0.6822242 0.6640976
#> [15] 0.6944021 0.7347680 0.7229173 0.7184145 0.7552095 0.7348064 0.8243417
#> [22] 0.9735380 1.0788747 1.1452296 1.1416286 1.1420774 1.2669930 1.3738789
#> [29] 1.0218517 1.1599925 1.2452416 1.3748183 1.4102338 1.4364045 1.5098982
#> [36] 1.5084806 1.5015768 1.6118554 1.7638951 1.8651654 1.7679079 2.3335854
#> [43] 2.4089201

cek class dari data co2_ts

# class ts object
class(co2_ts)

#> [1] "ts"

memvisualisasikan data time series untuk melihat pola datanya

# inspect pola data
co2_ts %>% autoplot()

Note :

Pergerakan emisi gas co2 di indonesia secara umum naik terus (trend naik). Terjadi penurunan pada tahun 1998, pada saat itu indonesia sedang terjadi krismon (krisis moneter) dimana bisa saja banyak pabrik yang tutup sehingga menyebabkan adanya menurunan emisi gas co2. Dan dari tahun 2010 kenaikan emisi co2 semakin meningkat tajam, bisa saja disebabkan adanya peningkatan jumlah kendaraan bermotor di indonesia yang semakin banyak.

Gunakan function window() untuk subset periode tertentu

window(x = co2_ts, start = 1995, end = 2000) %>% 
  autoplot()

Read data nybirth.csv

birth <- read.csv("data_input/nybirth.csv")
head(birth, 210)

#>           date births
#> 1   1946-01-01 26.663
#> 2   1946-02-01 23.598
#> 3   1946-03-01 26.931
#> 4   1946-04-01 24.740
#> 5   1946-05-01 25.806
#> 6   1946-06-01 24.364
#> 7   1946-07-01 24.477
#> 8   1946-08-01 23.901
#> 9   1946-09-01 23.175
#> 10  1946-10-01 23.227
#> 11  1946-11-01 21.672
#> 12  1946-12-01 21.870
#> 13  1947-01-01 21.439
#> 14  1947-02-01 21.089
#> 15  1947-03-01 23.709
#> 16  1947-04-01 21.669
#> 17  1947-05-01 21.752
#> 18  1947-06-01 20.761
#> 19  1947-07-01 23.479
#> 20  1947-08-01 23.824
#> 21  1947-09-01 23.105
#> 22  1947-10-01 23.110
#> 23  1947-11-01 21.759
#> 24  1947-12-01 22.073
#> 25  1948-01-01 21.937
#> 26  1948-02-01 20.035
#> 27  1948-03-01 23.590
#> 28  1948-04-01 21.672
#> 29  1948-05-01 22.222
#> 30  1948-06-01 22.123
#> 31  1948-07-01 23.950
#> 32  1948-08-01 23.504
#> 33  1948-09-01 22.238
#> 34  1948-10-01 23.142
#> 35  1948-11-01 21.059
#> 36  1948-12-01 21.573
#> 37  1949-01-01 21.548
#> 38  1949-02-01 20.000
#> 39  1949-03-01 22.424
#> 40  1949-04-01 20.615
#> 41  1949-05-01 21.761
#> 42  1949-06-01 22.874
#> 43  1949-07-01 24.104
#> 44  1949-08-01 23.748
#> 45  1949-09-01 23.262
#> 46  1949-10-01 22.907
#> 47  1949-11-01 21.519
#> 48  1949-12-01 22.025
#> 49  1950-01-01 22.604
#> 50  1950-02-01 20.894
#> 51  1950-03-01 24.677
#> 52  1950-04-01 23.673
#> 53  1950-05-01 25.320
#> 54  1950-06-01 23.583
#> 55  1950-07-01 24.671
#> 56  1950-08-01 24.454
#> 57  1950-09-01 24.122
#> 58  1950-10-01 24.252
#> 59  1950-11-01 22.084
#> 60  1950-12-01 22.991
#> 61  1951-01-01 23.287
#> 62  1951-02-01 23.049
#> 63  1951-03-01 25.076
#> 64  1951-04-01 24.037
#> 65  1951-05-01 24.430
#> 66  1951-06-01 24.667
#> 67  1951-07-01 26.451
#> 68  1951-08-01 25.618
#> 69  1951-09-01 25.014
#> 70  1951-10-01 25.110
#> 71  1951-11-01 22.964
#> 72  1951-12-01 23.981
#> 73  1952-01-01 23.798
#> 74  1952-02-01 22.270
#> 75  1952-03-01 24.775
#> 76  1952-04-01 22.646
#> 77  1952-05-01 23.988
#> 78  1952-06-01 24.737
#> 79  1952-07-01 26.276
#> 80  1952-08-01 25.816
#> 81  1952-09-01 25.210
#> 82  1952-10-01 25.199
#> 83  1952-11-01 23.162
#> 84  1952-12-01 24.707
#> 85  1953-01-01 24.364
#> 86  1953-02-01 22.644
#> 87  1953-03-01 25.565
#> 88  1953-04-01 24.062
#> 89  1953-05-01 25.431
#> 90  1953-06-01 24.635
#> 91  1953-07-01 27.009
#> 92  1953-08-01 26.606
#> 93  1953-09-01 26.268
#> 94  1953-10-01 26.462
#> 95  1953-11-01 25.246
#> 96  1953-12-01 25.180
#> 97  1954-01-01 24.657
#> 98  1954-02-01 23.304
#> 99  1954-03-01 26.982
#> 100 1954-04-01 26.199
#> 101 1954-05-01 27.210
#> 102 1954-06-01 26.122
#> 103 1954-07-01 26.706
#> 104 1954-08-01 26.878
#> 105 1954-09-01 26.152
#> 106 1954-10-01 26.379
#> 107 1954-11-01 24.712
#> 108 1954-12-01 25.688
#> 109 1955-01-01 24.990
#> 110 1955-02-01 24.239
#> 111 1955-03-01 26.721
#> 112 1955-04-01 23.475
#> 113 1955-05-01 24.767
#> 114 1955-06-01 26.219
#> 115 1955-07-01 28.361
#> 116 1955-08-01 28.599
#> 117 1955-09-01 27.914
#> 118 1955-10-01 27.784
#> 119 1955-11-01 25.693
#> 120 1955-12-01 26.881
#> 121 1956-01-01 26.217
#> 122 1956-02-01 24.218
#> 123 1956-03-01 27.914
#> 124 1956-04-01 26.975
#> 125 1956-05-01 28.527
#> 126 1956-06-01 27.139
#> 127 1956-07-01 28.982
#> 128 1956-08-01 28.169
#> 129 1956-09-01 28.056
#> 130 1956-10-01 29.136
#> 131 1956-11-01 26.291
#> 132 1956-12-01 26.987
#> 133 1957-01-01 26.589
#> 134 1957-02-01 24.848
#> 135 1957-03-01 27.543
#> 136 1957-04-01 26.896
#> 137 1957-05-01 28.878
#> 138 1957-06-01 27.390
#> 139 1957-07-01 28.065
#> 140 1957-08-01 28.141
#> 141 1957-09-01 29.048
#> 142 1957-10-01 28.484
#> 143 1957-11-01 26.634
#> 144 1957-12-01 27.735
#> 145 1958-01-01 27.132
#> 146 1958-02-01 24.924
#> 147 1958-03-01 28.963
#> 148 1958-04-01 26.589
#> 149 1958-05-01 27.931
#> 150 1958-06-01 28.009
#> 151 1958-07-01 29.229
#> 152 1958-08-01 28.759
#> 153 1958-09-01 28.405
#> 154 1958-10-01 27.945
#> 155 1958-11-01 25.912
#> 156 1958-12-01 26.619
#> 157 1959-01-01 26.076
#> 158 1959-02-01 25.286
#> 159 1959-03-01 27.660
#> 160 1959-04-01 25.951
#> 161 1959-05-01 26.398
#> 162 1959-06-01 25.565
#> 163 1959-07-01 28.865
#> 164 1959-08-01 30.000
#> 165 1959-09-01 29.261
#> 166 1959-10-01 29.012
#> 167 1959-11-01 26.992
#> 168 1959-12-01 27.897

Data di atas merupakan data persentase kelahiran di New York per bulan. Terdiri dari:

date: tanggal saat dilakukan pencatatan persentase kelahiran
births: persentase kelahiran.

Cek class data birth

class(birth)

#> [1] "data.frame"

Mengetahui range waktu data birth untuk membentuk menjadi object TS

mengubah kolom data date menjadi tipe data date
membuat kolom bulan dari kolom date

# manipulate data
birth_clean <- birth %>% mutate(date = ymd(date),
                                month = month(date, label = T))

?locales
Sys.setlocale("LC_TIME", "English")

#> [1] "English_United States.1252"

# cek range
range(birth_clean$date)

#> [1] "1946-01-01" "1959-12-01"

Eksplorasi data
cek kelengkapan interval waktu

birth_clean %>% pad() %>% is.na() %>% colSums()

#>   date births  month 
#>      0      0      0

visualisasi data

ggplot(data = birth_clean, aes(x = date, y = births))+
  geom_point(col = "blue")+
  geom_point(data = birth_clean %>% 
                    filter(month %in% c("Jan")), col = "red")+
  geom_line()

Insight:

Membuat object time series

membuat object ts dari data birth dan disimpan dalam object birth_ts. Pola yang ingin dilihat adalah pola tahunan.

# answer here
birth_ts <- ts(data = birth_clean$births, start = c(1946,1), frequency = 12)
birth_ts

#>         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
#> 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
#> 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
#> 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
#> 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
#> 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
#> 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
#> 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
#> 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
#> 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
#> 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
#> 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
#> 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
#> 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
#> 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
#>         Nov    Dec
#> 1946 21.672 21.870
#> 1947 21.759 22.073
#> 1948 21.059 21.573
#> 1949 21.519 22.025
#> 1950 22.084 22.991
#> 1951 22.964 23.981
#> 1952 23.162 24.707
#> 1953 25.246 25.180
#> 1954 24.712 25.688
#> 1955 25.693 26.881
#> 1956 26.291 26.987
#> 1957 26.634 27.735
#> 1958 25.912 26.619
#> 1959 26.992 27.897

Lihatlah visualisasi object birth_ts dari Jan 1950-Des 1955

library(TSstudio)
# answer here
birth_ts %>% ts_plot(slider = T)

2 Decomposition

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

Komponen/unsur dalam time series :

Trend : pola data secara general, cenderung untuk naik atau turun. Jika pada trend masih terdapat pola (fluktuatif) artinya masih ada pola yang belum tertangkap dengan baik.
Seasonal : pola musiman yang membentuk pola berulang pada periode waktu yang tetap
Residual/Error : pola yang tidak dapat ditangkap dalam trend dan seasonal

Sebelum melakukan model forecasting kita perlu mengamati objek time series dari hasil decompose. Ide utama dari decomposition 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 object time series kita menjadi 3 komponen tersebut, kita dapat menggunakan fungsi decompose().

# decompose
birth_dc <- birth_ts %>% 
            decompose()

Memvisualisasikan hasil decompose menggunakan autoplot() dari package forecast

# visualize
birth_dc %>% autoplot()

Pada hasil decompose kita mendapatkan informasi visualisasi:

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

Terdapat 2 jenis model pada data time series, yaitu :

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

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

Berikut ini contoh data time series yang bertipe additive :
Berikut ini contoh data time series yang bertipe multiplicative :

AirPassengers %>% autoplot()

2.1 Data dengan 4 frekuensi, ceritanya seasonal

birth_ts1 <- ts(data = birth_clean$births, start = c(1946,1), frequency = 4)
birth_ts1 %>% ts_plot(slider = T)

2.1.1 Decompose untuk data 4 iterasi

birth_dc1 <- birth_ts1 %>% 
            decompose()
birth_dc1 %>% autoplot()

2.2 Summary

Time series adalah suatu data yang berhubungan dengan runtun waktu.
Analisis time series adalah analisis data runtun waktu dimana yang akan kita prediksi adalah data time series itu sendiri dengan data historicalnya.
Untuk membuat object time series dapat menggunakan fungsi ts(data, start, frequency)

data : data yang ingin di forecast
start : awal mula dari data time series
frequency : pola berulang yang ingin diamati

Komponen data time series :

trend : pola data secara global, naik atau turun
seasonal : pola berulang dalam satu musim
random/error : pola tidak tertangkap oleh trend dan seasonal

Tipe pola data time series :

additive : pola time series apabila trend bertambah, maka ukuran seasonalnya tetap

\(X = T+S+E\)

multiplicative : pola time series apabila trend bertambah, maka ukuran seasonal makin besar atau kecil

\(X = T*S*E\)

Manfaat melakukan decompose:

menemukan pola trend, seasonal, dan random pada data time series
bisa melakukan analisis seasonal untuk melihat pola seasonalnya pada data yang nantinya akan berguna untuk pembuatan keputusan yang berdasarkan pola seasonal yang diperoleh
bisa melakukan analisis trend yang dapat digunakan untuk mengetahui pola data yang sesungguhnya pada data tanpa adanya efek seasonal pada data

2.3 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”.

Melakukan inspect komponen time series padaa 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.

birth_dc %>% autoplot()

Trend = Data - Seasonality - Error

birth_dc$trend %>% autoplot()

Pendekatan manual menddungkan Center Moving Average

dalam menentukan ordo berdasarkan frekuensi data yang ada pada object time series

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

birth_dc$trend %>% autoplot()

# contoh
dat <- c(2,3,1,2,1)
dat_ts <- ts(dat,frequency = 2)
dat_ts

#> Time Series:
#> Start = c(1, 1) 
#> End = c(3, 1) 
#> Frequency = 2 
#> [1] 2 3 1 2 1

ma(dat_ts,order = 2,centre = T)

#> Time Series:
#> Start = c(1, 1) 
#> End = c(3, 1) 
#> Frequency = 2 
#> [1]   NA 2.25 1.75 1.50   NA

2. Seasonal

Karena polanya additive maka:

Seasonal + Error = Data - Trend

birth_dc$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, 15), 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_dc$random %>% autoplot()

# matematis didalamnya

birth_error <- birth_ts - (birth_trend + birth_seasonal)
birth_error %>% 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. Mengetahui tipe time series yang kita miliki, multiplicative atau additive 2. Kita mengetahui komponen trend, seasonal, dan error pada data kita

2.4 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.

# menggunakan AirPassenger

AirPassengers %>% autoplot()

log(AirPassengers) %>% autoplot()

decompose(log(AirPassengers)) %>% autoplot()

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

# menggunakan AirPassenger
air_decom <- decompose(AirPassengers)

decompose(AirPassengers, type = "multiplicative") %>% autoplot()

Cara mendapatkan komponen time series pada jenis time series multiplicative

Trend

air_decom$trend %>% autoplot()

# proff manual

air_trend <- ma(AirPassengers, order = 12, centre = T)
air_trend %>% autoplot()

Seasonality

SEASONALITY * ERROR = DATA / TREND

air_decom$seasonal %>% autoplot()

# proff manual 

air_seaserr <- AirPassengers / air_trend

mean_seas <- air_seaserr %>% 
  matrix(ncol = 12, byrow = T) %>% 
  colMeans(na.rm = T)

air_seasonal <- ts(rep(mean_seas, 12), start = c(1949, 1), frequency = 12)
air_seasonal %>% autoplot()

Error ERROR = DATA / (TREND * SEASONALITY)

air_decom$random %>% autoplot()

# proff manual
air_error <- AirPassengers/(air_trend*air_seasonal)
air_error %>% autoplot()

2.4.1 [Optional] Analisis Seasonal

Untuk dapat melakukan analisis seasonal, kita dapat mengamati komponen seasonal secara visual:

# using ts_seasonal
ts_seasonal(ts.obj = birth_ts, type = "all")

3 Forecasting Model

Forecasting Model:

Simple Moving Average (SMA)
Exponential Smoothing
- Simple Exponential Smoothing (SES)
- Double Exponential Smoothing (Holts Exponential Smoothing)
- Triple Exponential Smoothing (Holts Winter)
Autoregresive Integrated Moving Average (ARIMA)

3.1 Simple 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 sering digunakan untuk data yang tidak mengandung trend dan seasonal (datanya bergerak disekitar rata-rata).

SMA tepat digunakan ketika data yang kita miliki tidak memiliki trend maupun seasonal.

Fungsi yang digunakan untuk forecasting SMA adalah SMA(objek time series, n) dari library TTR. Parameter yang digunakan, 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()

# read data rain
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)
rain_ts

#> Time Series:
#> Start = 1813 
#> End = 1912 
#> Frequency = 1 
#>   [1] 23.56 26.07 21.86 31.24 23.65 23.88 26.41 22.67 31.69 23.86 24.11 32.43
#>  [13] 23.26 22.57 23.00 27.88 25.32 25.08 27.76 19.82 24.78 20.12 24.34 27.42
#>  [25] 19.44 21.63 27.49 19.43 31.13 23.09 25.85 22.65 22.75 26.36 17.70 29.81
#>  [37] 22.93 19.22 20.63 35.34 25.89 18.65 23.06 22.21 22.18 18.77 28.21 32.24
#>  [49] 22.27 27.57 21.59 16.93 29.48 31.60 26.25 23.40 25.42 21.32 25.02 33.86
#>  [61] 22.67 18.82 28.44 26.16 28.17 34.08 33.82 30.28 27.92 27.14 24.40 20.35
#>  [73] 26.64 27.01 19.21 27.74 23.85 21.23 28.15 22.61 19.80 27.94 21.47 23.52
#>  [85] 22.86 17.69 22.54 23.28 22.17 20.84 38.10 20.65 22.97 24.26 23.01 23.67
#>  [97] 26.75 25.36 24.79 27.88

Visualisasikan data

rain_ts %>% autoplot()

Data curah hujan yang kita miliki bergerak disekitar rata-rata

Melakukan forecasting menggunakan SMA() dengan ordo 3

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

Memvisualisasikan data historis dengan hasil forecast

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

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.

Simple Exponential Smoothing digunakan untuk data time series yang tidak mengandung trend dan seasonal. Fungsi yang digunakan untuk membuat model exponential smoothing, yaitu:

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

objek timeseries yang digunakan.
model: model time series untuk error, trend, dan seasonal.
- A: aditif
- M: multiplikatif
- N: none
- Z: auto
alpha: bobot untuk smoothing, 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 = "***")

HoltWinters(), parameter yang digunakan, yaitu:

objek time series yang digunakan.
secara default parameter alpha, beta, dan gamma adalah NULL, dimana apabila kita tidak mendefinisikan nilainya, maka model HolWinters() 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.

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

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

Holt’s Winter Explonential: HoltWinters(data)

3.2.1 Simple Exponential Smoothing (SES)

Simple Exponential Smoothing merupakan metode forecasting yang tepat digunakan untuk data yang tidak memiliki trend maupun seasonal¹.

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

# visualisasikan rain
rain_ses <- ets(rain_ts, model = "ANN") # ZZZ : Aoto
rain_ses

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

Data rain_ts tidak memiliki pola trend maupun seasonal. Karena datanya polanya 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 = "ZZZ")

# untuk memperoleh hasil prediksi
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

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

# HoltWinters(x = rain_ts, alpha = T, beta = F, gamma = F)
rain_ses2 <- ets(y = rain_ts, model = "ANN", alpha = .6)
rain_ses2

#> ETS(A,N,N) 
#> 
#> Call:
#>  ets(y = rain_ts, model = "ANN", alpha = 0.6) 
#> 
#>   Smoothing parameters:
#>     alpha = 0.6 
#> 
#>   Initial states:
#>     l = 24.3141 
#> 
#>   sigma:  5.1751
#> 
#>      AIC     AICc      BIC 
#> 791.2695 791.3932 796.4799

Kesimpulan: > Ketika alpha mendekati 1: yang terjadi adalah model akan memberikan bobot yang paling besar untuk data terbaru, namun perubahan bobotnya sangat signifikan

Ketika alpha mendekati 0: yang terjadi adalah model akan memberikan bobot yang paling besar untuk data terbaaru, namun perubahan bobotnya secara perlahan

Visualisasi data historis dengan hasil forecast

# gunakan autolayer()
rain_ts %>% 
  autoplot() +
  autolayer(rain_ses$fitted, series = "SES alpha = 0.0001") + # fitted adalah data hasil prediksi dari data yang digunakan
  autolayer(rain_ses2$fitted, series = "SES alpha = 0.06")

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(object = rain_ses2$fitted, x = rain_ts)

#>                  ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set 0.04047917 5.123114 4.056663 -2.692735 16.35586 -0.331274 0.8526713

accuracy(object = rain_ses$fitted, x = rain_ts)

#>                     ME     RMSE      MAE       MPE     MAPE        ACF1
#> Test set -0.0004906603 4.193615 3.339691 -2.758248 13.61432 -0.06601945
#>          Theil's U
#> Test set 0.6795944

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

3.2.2 Holt’s Exponential Smoothing

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

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

co2_ts %>% autoplot()

co2_holt <- HoltWinters(x = co2_ts, gamma = F)
co2_holt2 <- HoltWinters(x = co2_ts, alpha = .9, beta = .1, gamma = F)

Alpha = 0.9 karena juka diperhatikan data terbaru itu nilainya jauh berbeda dengan data terdahulu (1970) itulah kenapa kita ingin memberikan bobot yang sangat besar untuk data terbaru, sedangkan data terdahulu diberikan bobot yang sangat kecil. Data terdahulu kurang relevan dengan kondisi saat ini.

Visualisasi data historis dengan hasil forecast

# gunakan autolayer()
co2_ts %>% autoplot() +
  autolayer(co2_holt$fitted[,1]) +
  autolayer(co2_holt2$fitted[,1])

Model evaluation

# gunakan accuracy()
accuracy(co2_holt$fitted[,1], x = 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_holt2$fitted[,1], x = co2_ts)

#>                 ME      RMSE        MAE      MPE     MAPE       ACF1 Theil's U
#> Test set 0.0199014 0.1174784 0.06970157 1.105452 6.279417 -0.1491025 0.8963148

MAPE : 6.27% artinya kesalahan model kita terhadap

3.2.3 Holt’s Winters Exponential

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

3.2.3.1 Workflow Analisis Time Series

Read data fancy.dat dan simpan dengan nama object souvenir

# read data souvenir
souvenir <- scan("data_input/fancy.dat")
head(souvenir)

#> [1] 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74

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

Buat time series object dan simpan dengan nama object souvenir_ts

# object souvenir_ts
souvenir_ts <- ts(data = souvenir, start = c(1987,1), frequency = 12)

Visualisasikan data souvenir_ts

Apakah tipe data time series yang terbentuk dari data souvenir_ts?

souvenir_ts %>% 
  autoplot()

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
souvenir_ts %>% 
  decompose(type = "multiplicative") %>%
  autoplot()

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()`
test <- tail(souvenir_ts, 12)

# train menggunakan `head()`
train <- head(souvenir_ts, 6*12)
# train <- head(souvenir_ts, length(souvenir_ts)-length(test))

Melakukan model fitting pada data train train

# ets(model = "MMM")
model_souvenir <- HoltWinters(x = train, seasonal = "multiplicative")
model_souvenir

#> Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
#> 
#> Call:
#> HoltWinters(x = train, seasonal = "multiplicative")
#> 
#> Smoothing parameters:
#>  alpha: 0.3469842
#>  beta : 0.07501578
#>  gamma: 0.5711478
#> 
#> Coefficients:
#>             [,1]
#> a   2.758252e+04
#> b   8.079173e+02
#> s1  4.676809e-01
#> s2  6.030877e-01
#> s3  9.592971e-01
#> s4  7.056491e-01
#> s5  6.754072e-01
#> s6  7.911595e-01
#> s7  8.912119e-01
#> s8  8.937350e-01
#> s9  8.856460e-01
#> s10 9.134738e-01
#> s11 1.392572e+00
#> s12 2.915102e+00

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_forecast <- forecast(object = model_souvenir, h = 12)
souvenir_forecast

#>          Point Forecast    Lo 80     Hi 80    Lo 95     Hi 95
#> Jan 1993       13277.67 11237.89  15317.44 10158.10  16397.24
#> Feb 1993       17609.17 15156.03  20062.32 13857.41  21360.94
#> Mar 1993       28784.94 25268.64  32301.23 23407.23  34162.65
#> Apr 1993       21744.01 18610.83  24877.19 16952.23  26535.79
#> May 1993       21357.80 18034.65  24680.95 16275.49  26440.12
#> Jun 1993       25657.32 21616.65  29698.00 19477.65  31837.00
#> Jul 1993       29622.05 24836.99  34407.12 22303.93  36940.18
#> Aug 1993       30427.98 25296.39  35559.56 22579.90  38276.06
#> Sep 1993       30868.11 25438.75  36297.48 22564.62  39171.61
#> Oct 1993       32576.03 26643.50  38508.56 23503.01  41649.05
#> Nov 1993       50786.55 41491.30  60081.80 36570.70  65002.41
#> Dec 1993      108667.82 88629.46 128706.18 78021.79 139313.85

h = 12 artinya kita akan melakukan forecasting untuk 12 bulan ke depan. Jumlah h yang kita tentukan, akan sesuai dengan data test yang kita siapkan.

Point forecast adalah forecast yang dihasilkan. 80 dan 95 ini adalah confidence interval * Lo : menandakan batas bawah * Hi : menandakan batas atas

Summary:

Perkenalan entang model-model yang ada. ada SMA yaitu data yang tidak memiliki train dan seasonal. liat ordo nya, jika 10, berarti liat rata2 10 periode terakhir, semua data dianggpap penting, ga semua dihitung. penentuan ordonya user yang menentukan.

Exponential smoothing ada 3 model, SES, Holt dan Holt Winter, komponen yang akan digunakan untuk melakukan smoothing adalah

alpha -> smoothing error beta -> smoothing trend theta -> smoothing seasonal

Simple Exponential memberikan bobot terbesar untuk data terbaru, seberapa kecil atau besar bisa dilihat dari paramete alpha theta gamma nya. alpha makin besar berarti data kondisi sekarang makin berpengaruh dan data makin ke belakan makin tidak berpengaruh.
Holt Exponential Smoothing digunakan untuk data yang memiliki trend
Holt winters Smoothing digunakan pada data yang memiliki trend dan juga seasonal yang semuanya multiplikatif

Intrepretasi hasil forecast :

Visualisasi hasil forecast

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

Model evaluation menggunakan data souvenir_test

# model multiplicative
accuracy(object = souvenir_forecast$mean, 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

MAPE: Mean Absolute percentage error

Dive Deeper

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.

Read Data

birth <- read.csv("data_input/nybirth.csv")

Buatlah object ts

birth_ts <- ts(data = birth$births, start = c(1946,1), frequency = 12)
birth_ts %>% autoplot()

Visualisasikan object ts

birth_ts %>% decompose(type = "additive") %>% autoplot()

Splitting train test, untuk data train 12 tahun, data test 2 tahun

# test
birth_test <- tail(birth_ts, 2*12)

# train menggunakan `head()`
birth_train <- head(birth_ts, 12*12)

Lakukan decompose untuk data train dan visualisasikan untuk menentukan additive atau multiplicative

decompose(x = birth_train, type = "additive") %>% autoplot()

Modelkan dengan metode yang tepat yang sudah dipelajari

model dengan nilai smoothing nya auto
model dengan nilai smoothing dengan mendefinisikan nilai alpha beta dan gamma secara manual

# a
model_birth_auto <- HoltWinters(x = birth_train, seasonal = "additive")
# b
model_birth_manual <- HoltWinters(x = birth_train, alpha = 0.5, beta = 0.01, gamma = 0.005, seasonal = "additive")
model_birth_hanung <- HoltWinters(x = birth_train, alpha = 0.1, beta = 0.4, gamma = 0.1, seasonal = "additive")

Lakukanlah forecasting untuk 2 tahun kedepan (24 bulan)

birth_forecast_auto <- forecast(object = model_birth_auto, h = 24)
birth_forecast_manual <- forecast(object = model_birth_manual, h = 24)
birth_forecast_hanung <- forecast(object = model_birth_hanung, h = 24)

Visualisasikan data train, test dan hasil forecastnya

birth_train %>% 
  autoplot()+
  autolayer(birth_test, series = "test")+
  autolayer(birth_forecast_manual$mean, series = "forecast_manual")+
  autolayer(birth_forecast_auto$mean, series = "forecast_auto")+
  autolayer(birth_forecast_hanung$mean, series = "forecast Hanung")+
  autolayer(model_birth_manual$fitted[,1], series = "fitted")

Model Evaluation

accuracy(birth_forecast_manual$mean, birth_test)

#>                  ME     RMSE       MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -0.2765909 1.090862 0.9364968 -1.172147 3.456039 0.6976645 0.6511813

accuracy(birth_forecast_auto$mean, birth_test)

#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -1.231122 1.579416 1.249804 -4.612144 4.674416 0.7004142  0.965476

accuracy(birth_forecast_hanung$mean, birth_test)

#>                  ME      RMSE       MAE         MPE     MAPE      ACF1
#> Test set 0.01246057 0.7761363 0.6616637 -0.06608348 2.418264 0.6171999
#>          Theil's U
#> Test set 0.4560637

3.3 Summary

Pola data time series

Additive, lebih mudah diprediksi/diforecast karena perubahan pada data time series tersebut cenderung tetap/statis.
Multiplicative, lebih sulit diprediksi/diforecast karena perubahan pada data time series cenderung dinamis/berubah-ubah biasanya kelipatannya/multiplikatif (perubahan yang terjadi bisa semakin membesar/mengecil). Sehingga, umumnya dilakukan transformasi data dengan pola multiplicative menjadi additive dengan transformasi `log()

Forecasting model:

SMA(Single moving average) : memperhitungkan beberapa data terakhir saja tanpa melihat pola data secara keseluruhan. Cara kerja dari SMA selain dia merata-ratakan datanya, sebenarnya setiap datanya diberikan bobot yang sama sesuai ordonya, misalnya ordo 3, maka setiap data diberikan bobot 1/3.
Simple exponential smoothing (SES) : cocok untuk data yang tidak memiliki pola trend dan seasonal. Parameter smoothing didalamnya yaitu alpha (untuk eror).
Holt Exponential : cocok untuk data yang tidak memiliki pola seasonal. Parameter smoothing didalamnya yaitu alpha (untuk eror) dan beta (untuk trend).

Untuk membuat model exponential smoothing bisa menggunakan dua fungsi yaitu ets() dan HoltWinters().

fungsi ets(object_ts, model = "ZZZ")
fungsi HoltWinters(object_ts, alpha = NULL/FALSE, beta = NULL/FALSE, gamma = NULL/FALSE, seasonal = "additive/multiplicative")

Untuk melakukan evalusi model bisa menggunakan nilai error :

RMSE(Root Mean Square Error) : range nilainya cenderung seperti data aslinya
MAPE(Mean Absolute Percentage Error) : range nilainya 0%-100%

—–End Of day 2———-

Knowledge check:

Bagaimana karakteristik data time series?

data waktu terurut
interval waktu sama
tidak boleh ada data yang missing

Berapa minimal observasi yang perlu dimiliki untuk time series forecasting?

untuk pemodelan minimal data yang kita miliki setidaknya 2 periode

untuk validasi data test setidaknya siapkan 1 periode

Kenapa kita harus melakukan decompose data?

untuk mengidentifikasi pola error, trend, dan seasonal untuk mempermudah memilih model yang tepat
Ketika nilai freq = 1 tidak bisa dilakukan decompose, karena dia tidak memiliki pembanding dengan periode sebelumnya

Apa perbedaan dari data additive dan multiplicative?

pola additive cenderung kenaikan dan penurunannya konstan
pola multiplicative cenderung kenaikan dan penurunan lembar semakin lama semakin besar atau kecil

Bagaimana kita lakukan cross validation untuk data time series?

tidak boleh random/acak, harus terurut
data test yang digunakan adalah data terbaru
jumlah data menyesuaikan periode

Kapan menggunakan simple moving average?
Bagaimana karakteristik simple moving average?
Kapan menggunakan model simple exponential smoothing, holt exponential smoothing, dan holt winter?
Bagaimana peran alpha beta dan gamma pada modelling?

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 konstan. Berdasarkan gambar dibawah ini, manakah yang termasuk data time series yang stasioner?

knitr::include_graphics("assets/TS_Variation.PNG")

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 sebelumnya⁴. Terkadang, tergantung pada kompleksitas data, jumlah differencing bisa lebih dari 1 kali.

x <- c(5, -3, 6, 7)
x

#> [1]  5 -3  6  7

diff(x) # differencing 1 kali

#> [1] -8  9  1

diff(diff(x)) # differencing 2 kali

#> [1] 17 -8

Pemahaman konsep mengapa perlu differencing???

AirPassengers %>% 
  diff() %>%
  autoplot()

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()

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()

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

summary:

model arima > melakukan regresi > antar lag tidak boleh memiliki hubungan yang kuat > data harus stasioner > ketika data kita tidak stasioner > kita haru membuat data menjadi stasioner dengan cara lakukan differencing > differencing akan menghilangkan efek multikolinearitas

Quick Summary

ARIMA memerlukan model dalam keadaan stasioner. Teknik yang bisa dilakukan untuk membuat data menjadi stasioner adalah dengan menggunakan differencing. Kestasioneran data diperlukan karena model regresi memerlukan kondisi korelasi antar variabel tidak kuat. Jika dibiarkan kuat, maka akan terjadi masalah multikolinearitas.

4.3 AR(p) : Auto Regressive

Auto regressive di ARIMA artinya kita membuat model linear regression berdasarkan lag dari datanya sebagai prediktor⁵. Linear regression models, baik ketika prediktornya tidak ada korelasi dan independen satu dengan lainnya. 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.

4.4 MA(q) : Moving Average

Moving Average dalam ARIMA artinya kita melakukan rata-rata berjalan terhadap data time series itu sendiri⁶.

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}\]

Arima memiliki bentuk model (p,d,q) dimana⁷ :

Auto Regressive -> AR(p)
AR ini mirip dengan regresi linear, hanya saja dia melakukan regresi terhadap dirinya sendiri. Nilai p berarti berapa banyak data yang akan dipakai ketika melakukan auto regressive.

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.

Integrated -> I(d)
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.
Moving Average -> MA(q)
MA digunakan untuk melakukan smoothing error. Nilai q berarti berapa banyak data yang diperlukan untuk smoothing error menggunakan moving average.

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:

Objek time series
Melakukan pengujian stasioner menggunakan adf.test

H0 : data tidak stasioner H1 : data stasioner

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

Catat informasi berapa kali dilakukan differencing sampai data menjadi stasioner
Lakukan pembuatan model arima menggunakan auto.arima(objek_ts, seasonal = F).
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.

Memperhatikan pola dari plot acf dan pacf. Dimana ada dua tipe pola yaitu: tails off atau cuts of lag.

tails off: garis saat lag ke-1 mengalami penurunan yang lambat.
cuts off lag: garis saat lag ke-1 mengalami penurunan yang cepat.

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

Model evaluation dengan nilai MAPE
Cek asumsi time series

normality residual assumption
autocorrelation residual assumption

4.5 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.

Read Data

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

#>   year CO2.emissions..kt. CO2.emissions..metric.tons.per.capita.
#> 1 1970           35822.92                              0.3119519
#> 2 1971           38987.54                              0.3306215
#> 3 1972           43340.27                              0.3580080
#> 4 1973           49134.13                              0.3954703
#> 5 1974           51260.99                              0.4021567
#> 6 1975           53963.57                              0.4128050
#>   Methane.emissions..kt.of.CO2.equivalent.
#> 1                                   126665
#> 2                                   119293
#> 3                                   195338
#> 4                                   136928
#> 5                                   122570
#> 6                                   121065
#>   Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.
#> 1                                                         51908.88
#> 2                                                         49004.80
#> 3                                                         79261.42
#> 4                                                         56928.71
#> 5                                                         50997.79
#> 6                                                         50922.15
#>   Other.greenhouse.gas.emissions..HFC..PFC.and.SF6..thousand.metric.tons.of.CO2.equivalent.
#> 1                                                                                 130405.80
#> 2                                                                                  94658.48
#> 3                                                                                 422459.50
#> 4                                                                                 174340.90
#> 5                                                                                 110213.86
#> 6                                                                                 109630.88
#>   Total.greenhouse.gas.emissions..kt.of.CO2.equivalent.
#> 1                                              339677.1
#> 2                                              294004.9
#> 3                                              733432.9
#> 4                                              409241.7
#> 5                                              329453.2
#> 6                                              329119.3

Membuat object time series untuk kolom Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.

no_ts <- ts(data = no$Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent., start = 1970,  frequency = 1)

Melakukan pengecekan object ts stasioner

Stasioner adalah data yang bergerak disekitar rata-rata globalnya.

melihat berdasarkan visual

no_ts %>% autoplot()

melakukan pengujian menggunakan adf.test()

adf.test(no_ts)

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  no_ts
#> Dickey-Fuller = -2.7827, Lag order = 3, p-value = 0.2645
#> alternative hypothesis: stationary

Menguji stasioneritas data H0: data tidak stasioner H1: data stasioner

pvalue>alpha, artinya data belum stasioner

Melakukan differencing data agar stasioner

mencoba melakukan pengujian data no_ts dengan satu kali differencing

# melihat secara visual
no_ts %>% diff %>% autoplot()

Untuk mengetahui data stasioner atau tidak kita akan menggunakan adf.test()

# lakukan uji statistik
no_ts %>% diff() %>% adf.test()

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

H0 : Data tidak stasioneer H1 : Data Stasioner

Data stasioner menunjukkan bahwa data lag(prediktor) saling berkorelasi.

Kesimpulan:

Model Fitting

melakukan tunning model untuk p, d, 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()`
no_ts %>% diff() %>% tsdisplay()

Dies down/decay: mengalami penurunan atau kenaikan secara lambat, ordo yang dipakai adalah 0
cuts off lag: mengalami penurunan atau kenaikan secara cepat

PACF untuk ordo p ACF untuk ordo q kalo dies down diisi dengan 0 kalo cut-off liat dimana, sampai Lag yang masih signifikan.

Asumsi pertama : ACF itu cut off, PACF itu dies down PACF -> dies down berarti q = 0 ACF -> cut-off, terakhir di Lag 1 berarti q = 1 Kalau seandainya, ACF punya 2 yang signifikan misalnya Lag 2 juga signifikan, maka ACF nya = 2

ARIMA(p,d,q) = ARIMA(0,1,1)

Kombinasi model ARIMA yang terbentuk ⁸:

ARIMA(p,d,q)

ARIMA(p,d,q)

Asumsi kedua : Keduanya cut-off PACF = q = 1 ACF -> p = 1,2,4 (karena cut-off di lag 1,2,4) ARIMA(p,d,q) 1. ARIMA(1,1,1) 2. ARIMA(2,1,1) 3. ARIMA(4,1,1)

parameter pada fungsi Arima()

y = object time series
order = nilai (p,d,q)

# create model
no_arima1 <- Arima(y = no_ts, order = c(0,1,1), seasonal = F)
no_arima2 <- Arima(y = no_ts, order = c(1,1,1), seasonal = F)
no_arima3 <- Arima(y = no_ts, order = c(2,1,1), seasonal = F)
no_arima4 <- Arima(y = no_ts, order = c(4,1,1), seasonal = F)
no_arima3

#> Series: no_ts 
#> ARIMA(2,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2      ma1
#>       -0.0038  -0.2197  -0.7690
#> s.e.   0.1780   0.1634   0.1224
#> 
#> sigma^2 estimated as 3.227e+09:  log likelihood=-518.45
#> AIC=1044.9   AICc=1045.98   BIC=1051.85

cek aic masing-masing model

no_arima1$aic

#> [1] 1042.681

no_arima2$aic

#> [1] 1044.575

no_arima3$aic

#> [1] 1044.902

no_arima4$aic

#> [1] 1048.008

Kita ingin mendapatkan nilai AIC terendah karena AIC menandakan banyaknya informasi yang hilang.

Membuat model arima secara otomatis menggunakan auto.arima()

auto.arima(no_ts, seasonal = F)

#> Series: no_ts 
#> ARIMA(0,1,1) 
#> 
#> Coefficients:
#>           ma1
#>       -0.8244
#> s.e.   0.0842
#> 
#> sigma^2 estimated as 3.21e+09:  log likelihood=-519.34
#> AIC=1042.68   AICc=1042.99   BIC=1046.16

Goodness of Fit (Pemilihan model terbaik)

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

accuracy(no_arima1)

#>                    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

accuracy(no_arima2)

#>                   ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
#> Training set 8204.36 55267.96 33751.68 -2.00415 25.36724 0.7837685 -0.03222488

Model Evaluation

bisa menggunakan fungsi accuracy() dari package forecast untuk melihat MAPE

no_ts %>% autoplot() + 
  autolayer(no_arima1$fitted, series = "fitted 1") +
  autolayer(no_arima2$fitted, series = "fitted 2")

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

Seasonal arima adalah metode arima dimana object time series yang ada memiliki pola seasonal⁹. 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)[m]

(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

Read data

# read data DAUTONSA.csv
sales <- read.csv("data_input/DAUTONSA.csv") %>% 
        mutate(DATE = ymd (DATE))
tail(sales, 20)

#>           DATE DAUTONSA
#> 611 2017-11-01  340.959
#> 612 2017-12-01  362.059
#> 613 2018-01-01  280.479
#> 614 2018-02-01  328.988
#> 615 2018-03-01  416.105
#> 616 2018-04-01  326.924
#> 617 2018-05-01  394.895
#> 618 2018-06-01  372.417
#> 619 2018-07-01  321.267
#> 620 2018-08-01  338.095
#> 621 2018-09-01  329.564
#> 622 2018-10-01  329.158
#> 623 2018-11-01  306.410
#> 624 2018-12-01  342.587
#> 625 2019-01-01  277.863
#> 626 2019-02-01  285.093
#> 627 2019-03-01  372.433
#> 628 2019-04-01  295.505
#> 629 2019-05-01  346.105
#> 630 2019-06-01  327.864

Data di atas merupakan data rata-rata penjualan motor di suatu dealer dari tahun 1967-2019.

DATE: tanggal saat dilakukan pencatatan penjualan
DAUTONSA: rata-rata penjualan.
untuk kebutuhan analisis, datanya akan kita subset dari tahun 2011 keatas

# filter tahun 2011 keatas
sales_clean <- sales %>% filter(DATE >= "2011-01-01")
sales_clean

#>           DATE DAUTONSA
#> 1   2011-01-01  258.576
#> 2   2011-02-01  338.615
#> 3   2011-03-01  445.315
#> 4   2011-04-01  405.418
#> 5   2011-05-01  362.214
#> 6   2011-06-01  348.950
#> 7   2011-07-01  329.517
#> 8   2011-08-01  340.459
#> 9   2011-09-01  327.466
#> 10  2011-10-01  321.021
#> 11  2011-11-01  304.949
#> 12  2011-12-01  363.464
#> 13  2012-01-01  323.424
#> 14  2012-02-01  430.209
#> 15  2012-03-01  528.054
#> 16  2012-04-01  431.834
#> 17  2012-05-01  491.050
#> 18  2012-06-01  459.854
#> 19  2012-07-01  393.252
#> 20  2012-08-01  452.837
#> 21  2012-09-01  411.972
#> 22  2012-10-01  369.820
#> 23  2012-11-01  382.692
#> 24  2012-12-01  444.846
#> 25  2013-01-01  383.879
#> 26  2013-02-01  446.362
#> 27  2013-03-01  544.516
#> 28  2013-04-01  466.742
#> 29  2013-05-01  519.361
#> 30  2013-06-01  499.660
#> 31  2013-07-01  446.915
#> 32  2013-08-01  516.827
#> 33  2013-09-01  383.073
#> 34  2013-10-01  398.003
#> 35  2013-11-01  406.615
#> 36  2013-12-01  421.205
#> 37  2014-01-01  344.189
#> 38  2014-02-01  409.921
#> 39  2014-03-01  534.342
#> 40  2014-04-01  470.222
#> 41  2014-05-01  564.577
#> 42  2014-06-01  499.801
#> 43  2014-07-01  479.974
#> 44  2014-08-01  556.870
#> 45  2014-09-01  414.495
#> 46  2014-10-01  430.949
#> 47  2014-11-01  428.752
#> 48  2014-12-01  475.786
#> 49  2015-01-01  386.010
#> 50  2015-02-01  419.495
#> 51  2015-03-01  525.027
#> 52  2015-04-01  477.314
#> 53  2015-05-01  563.176
#> 54  2015-06-01  491.397
#> 55  2015-07-01  481.665
#> 56  2015-08-01  496.815
#> 57  2015-09-01  446.569
#> 58  2015-10-01  447.743
#> 59  2015-11-01  387.789
#> 60  2015-12-01  472.123
#> 61  2016-01-01  361.397
#> 62  2016-02-01  427.270
#> 63  2016-03-01  502.115
#> 64  2016-04-01  447.760
#> 65  2016-05-01  468.254
#> 66  2016-06-01  452.263
#> 67  2016-07-01  437.189
#> 68  2016-08-01  433.380
#> 69  2016-09-01  419.560
#> 70  2016-10-01  377.266
#> 71  2016-11-01  374.906
#> 72  2016-12-01  444.215
#> 73  2017-01-01  307.619
#> 74  2017-02-01  368.826
#> 75  2017-03-01  446.580
#> 76  2017-04-01  411.024
#> 77  2017-05-01  433.909
#> 78  2017-06-01  399.619
#> 79  2017-07-01  376.212
#> 80  2017-08-01  394.215
#> 81  2017-09-01  406.634
#> 82  2017-10-01  345.309
#> 83  2017-11-01  340.959
#> 84  2017-12-01  362.059
#> 85  2018-01-01  280.479
#> 86  2018-02-01  328.988
#> 87  2018-03-01  416.105
#> 88  2018-04-01  326.924
#> 89  2018-05-01  394.895
#> 90  2018-06-01  372.417
#> 91  2018-07-01  321.267
#> 92  2018-08-01  338.095
#> 93  2018-09-01  329.564
#> 94  2018-10-01  329.158
#> 95  2018-11-01  306.410
#> 96  2018-12-01  342.587
#> 97  2019-01-01  277.863
#> 98  2019-02-01  285.093
#> 99  2019-03-01  372.433
#> 100 2019-04-01  295.505
#> 101 2019-05-01  346.105
#> 102 2019-06-01  327.864

Membuat object time series, dimana periode bermula dari Januari 2011

sales_ts <- ts(data = sales_clean$DAUTONSA, start = c(2011,1), frequency = 12)

Melihat visualisasi data secara umum

sales_ts %>% autoplot()

Splitting data

Test data : sisanya Train data : diambil 7 tahun pertama

# data train selama 7 tahun
sales_train <- head(sales_ts, 7*12)
# data test 1.5 tahun (12+6)
sales_test <- tail(sales_ts, 18)

Melakukan pengecekan stasionaritas data

# mencoba dengan seasonal period (frequency)
sales_train %>% adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -2.9359, Lag order = 4, p-value = 0.192
#> alternative hypothesis: stationary

H0: data tidak stasioner H1: data stasioner

Data tidak stasioner

Differencing based on frequency, cek apakah sudah stasioner atau belum

sales_train %>% diff() %>% autoplot() # lansung hilangin efek trend

sales_train %>% diff(lag = 12) %>% diff() %>% autoplot() # hilangin efek seasonal dulu

sales_train %>% diff(lag = 12) %>% diff() %>% adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -4.7557, Lag order = 4, p-value = 0.01
#> alternative hypothesis: stationary

Pertama, buang data seasonal. D menyatakan nilai differencing seasonal. disini diff bisa dua kali, ada data yg tahunan dua kali Lalu data trend di buang. d menyatakan differencing trend

Tentukan nilai differrence untuk index utama dan seasonal:

d = 1

D = 1

sales_train %>% diff(lag = 12) %>% diff(lag = 1) %>% 
  tsdisplay()

Untuk seasosal, AR murni
* P : 0 * D : 1 * Q : 0

untuk keseluruhan data, MA murni * p: 0 * d: 1 * q: 1

Kombinasi model yang mungkin terbentuk adalah :

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

SARIMA (0,1,1)(0,1,0)[m]
SARIMA (0,1,1)(2,1,0)

Model Fitting

membuat model SARIMA menggunakan auto.arima()

# create model using auto.arima()
sarima1 <- Arima(y = sales_train, order = c(0,1,1), seasonal = c(0,1,0))
sarima2 <- Arima(y = sales_train, order = c(0,1,1), seasonal = c(2,1,0))
sarima.auto <- auto.arima(y = sales_train, seasonal = T)

# cek nilai AIC
sarima1$aic

#> [1] 687.4144

sarima2$aic

#> [1] 690.8078

tunning model SARIMA menggunakan nilai residual

model forecast yang baik akan menghasilkan residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.

Menentukan nilai (p,d,q)(P,D,Q)[m]

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.

Kombinasi model yang mungkin terbentuk adalah :

ARIMA ( )( )[ ]
SARIMA

# create model

Melakukan forecasting untuk 18 bulan kedepan

sales_forecast <- forecast(sarima1, h = 18)

kalau data di log(), lakukan exp() disini

Visualisasi

sales_train %>% autoplot() + 
  autolayer(sales_test, series = "data test") +
  autolayer(sarima1$fitted, series = "fitted") +
  autolayer(sales_forecast$mean, series = "forecast")

Model evaluation

menggunakan accuracy() dari package forecast dan melihat nilai MAPEnya

# accuracy
accuracy(sales_forecast$mean, sales_test)

#>                ME     RMSE      MAE      MPE     MAPE       ACF1 Theil's U
#> Test set 14.06787 27.57037 23.20688 4.370293 7.204645 0.08457512  0.515082

Kesimpulan:

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 ini¹⁰:

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

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

5.1 No-autocorrelation residual

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

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

Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time series bisa menggunakan beberapa cara, yaitu:

menggunakan visualisasi plot residual acf dengan fungsi acf(residual model)
melakukan uji Ljung-box dengan menggunakan fungsi Box.test(residual model, type = "Ljung-Box)

# menggunakan visualisasi
acf(sarima1$residuals)

# menggunakan Ljung-Box test
Box.test(sarima1$residuals, type = "Ljung-Box")

#> 
#>  Box-Ljung test
#> 
#> data:  sarima1$residuals
#> X-squared = 0.075662, df = 1, p-value = 0.7833

\(H_0\): residual has no-autocorrelation > 0.05 \(H_1\): residual has autocorrelation < 0.05

Kesimpulan: nilai p-value (0.78) > (0.75) artinya residual kita tidak terjadi autokorelasi

jika terjadi autokorelasi 1. diff() 2. log 3. cek data, harusnya ada yang outlier -> hist. boxplot 4. ## Normality residual

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

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

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

shapiro.test(sarima1$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  sarima1$residuals
#> W = 0.96645, p-value = 0.02739

p-value > alpha maka data tidak normally distributed

hist(sarima1$residuals)

Mau mentolerir error ini ga, kemampuan model kita hanya bisa memprediksi kurleb 50 observasi saja, antara -50 sampai 50 saja.

Note : Untuk mengecek asumsi cukup di model akhir yang akan dipilih

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

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()¹¹.

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 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

sales_train %>% 
  stl(s.window = 12) %>% 
  autoplot()

6.2 Modeling

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

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

ets_stlm <- stlm(y = sales_train, s.window = 12, method = "ets")
arima_stlm <- stlm(y = sales_train, s.window = 12, method = "arima")

arima_stlm$model

#> Series: x 
#> ARIMA(2,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2    ma1
#>       -1.3573  -0.5765  0.970
#> s.e.   0.0903   0.0897  0.046
#> 
#> sigma^2 estimated as 562.9:  log likelihood=-380.07
#> AIC=768.15   AICc=768.66   BIC=777.82

cek nilai aic

ets_stlm$model$aic

#> [1] 920.6818

arima_stlm$model$aic

#> [1] 768.1463

6.3 Forecasting

Melakukan forecast untuk 18 bulan mendatang

forecast_stlm_arima <- forecast(arima_stlm , h =18)

Visualisasi hasil forecast

sales_train %>% 
  autoplot()+
  autolayer(sales_test, series = "test")+
  autolayer(forecast_stlm_arima$mean, series = "forecast")

7 Dive Deeper

Gunakan data euretail untuk forecasting, data euretail ini sudah dalam bentuk object time series dengan frequency = 4

data("euretail")
class(euretail)

#> [1] "ts"

frequency(euretail)

#> [1] 4

inspect data

euretail %>% autoplot()

split data train dan data test

#data test (kita akan melakukan prediksi 3 tahun terakhir)
eur_test <- tail(euretail, 3*4)
#data train (sisa dari data test)
eur_train <- head(euretail, 13*4)

eur_train %>% adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -0.40927, Lag order = 3, p-value = 0.9826
#> alternative hypothesis: stationary

eur_train %>% diff(lag = 4) %>% adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -3.8857, Lag order = 3, p-value = 0.02218
#> alternative hypothesis: stationary

eur_train %>% diff(lag = 4) %>% diff() %>% adf.test()

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

Differencing pada data seasonal membuat data menjadi stasioner

eur_train %>% diff(lag = 4) %>% tsdisplay()

p: 1 d: 0 q: 0

P: 0 D: 1 Q: 3

Buatlah model dari kombinasi di atas

eur_sarima1 <- Arima(y = eur_train, order = c(1,0,0), seasonal = c(0,1,3))
eur_sarima2 <- Arima(y = eur_train, order = c(1,0,0), seasonal = c(0,1,0))

Cek aic model

eur_sarima1$aic

#> [1] 68.18751

eur_sarima2$aic

#> [1] 81.93948

forecast data test

eur_forecast <- forecast(object = eur_sarima1, h = 12)

visualisasi hasil forecast dan data actual

eur_train %>% 
  autoplot() +
  autolayer(object = eur_sarima1$fitted, series = "fitted") +
  autolayer(object = eur_test, series = "test") +
  autolayer(object = eur_forecast$mean, series = "forecast")

Evaluasi model

accuracy(eur_forecast$mean, eur_test)

#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -3.985693 4.194059 3.985693 -4.129324 4.129324 0.6565622  9.408637

8 Reference

More resource to learn more :

Time Series : In-class materials

Davel Zhafran

September 12, 2021

1 Time Series and Forecasting

1.1 Karakteristik data time series

1.2 Tahapan membuat object time series

2 Decomposition

2.1 Data dengan 4 frekuensi, ceritanya seasonal

2.1.1 Decompose untuk data 4 iterasi

2.2 Summary

2.3 Additive Time Series

2.4 Multiplicative Time Series

2.4.1 [Optional] Analisis Seasonal

3 Forecasting Model

3.1 Simple Moving Average (SMA)

3.2 Exponential Smoothing

3.2.1 Simple Exponential Smoothing (SES)

3.2.2 Holt’s Exponential Smoothing

3.2.3 Holt’s Winters Exponential

3.2.3.1 Workflow Analisis Time Series

3.3 Summary

4 ARIMA (Autoregressive Integrated Moving Average)

4.1 Stasionarity

4.2 Differencing

4.3 AR(p) : Auto Regressive

4.4 MA(q) : Moving Average

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

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

5 Assumption

5.1 No-autocorrelation residual

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

6.1 Decompose menggunakan `stl()`

6.2 Modeling

6.3 Forecasting

7 Dive Deeper

8 Reference

Time Series : In-class materials

Davel Zhafran

September 12, 2021

1 Time Series and Forecasting

1.1 Karakteristik data time series

1.2 Tahapan membuat object time series

2 Decomposition

2.1 Data dengan 4 frekuensi, ceritanya seasonal

2.1.1 Decompose untuk data 4 iterasi

2.2 Summary

2.3 Additive Time Series

2.4 Multiplicative Time Series

2.4.1 [Optional] Analisis Seasonal

3 Forecasting Model

3.1 Simple Moving Average (SMA)

3.2 Exponential Smoothing

3.2.1 Simple Exponential Smoothing (SES)

3.2.2 Holt’s Exponential Smoothing

3.2.3 Holt’s Winters Exponential

3.2.3.1 Workflow Analisis Time Series

3.3 Summary

4 ARIMA (Autoregressive Integrated Moving Average)

4.1 Stasionarity

4.2 Differencing

4.3 AR(p) : Auto Regressive

4.4 MA(q) : Moving Average

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

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

5 Assumption

5.1 No-autocorrelation residual

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

6.1 Decompose menggunakan stl()

6.2 Modeling

6.3 Forecasting

7 Dive Deeper

8 Reference

6.1 Decompose menggunakan `stl()`