Analisis Time Series pada Renewable Energy Indonesia

Ichlasul Amal

14 Februari 2021

Objektif

Artikel ini akan membahas penerapan analisis time series dalam sektor lingkungan, yaitu renewable energy di Indonesia. Renewable energy adalah penggunaan dari listrik terbarukan (bagian listrik yang dihasilkan oleh pembangkit listrik terbarukan) dalam jumlah listrik yang dihasilkan oleh semua jenis pembangkit listrik.

Metode yang akan dipakai adalah exponential smooothing dan ARIMA. Metode ini berguna untuk melakukan forecasting pada data time series. Analisis ini dapat digunakan untuk melihat karakteristik atau pola pada penggunaan dan pengelolaan energi di Indonesia. Dataset berupa data tahunan yang diambil dari macrotrends dan dapat diakses melalui tautan renewable energy.

Read data

Melakukan read pada dataset yang akan digunakan dalam analisis.

renewable <- read.csv("renewable.csv", stringsAsFactors = F)

Selanjutnya melihat struktur dari tiap dataset, termasuk di dalamnya tipe tiap kolomnya.

str(renewable)
#> 'data.frame':    26 obs. of  3 variables:
#>  $ date                             : chr  "12/31/1990" "12/31/1991" "12/31/1992" "12/31/1993" ...
#>  $ X..of.Electricity.from.Renewables: num  20.9 20.6 23.8 19.8 15.7 ...
#>  $ Annual.Change                    : num  NA -0.36 3.28 -4.06 -4.12 0.81 -0.93 -5.19 5.43 -1.64 ...

Data pre-processing

Melakukan perubahan pada tipe data yang belum sesuai, yaitu kolom date menjadi tipe Date dan menghapus kolom annual change atau selisih nilai di tiap tahun, karena tidak dibutuhkan dalam analisis time series.

library(lubridate)
library(tidyverse)

renewable <- renewable %>% 
  select(-Annual.Change) %>% 
  mutate(date = mdy(date))

Selanjutnya mengecek apakah terdapat missing value atu tidak dalam dataset renewable.

library(padr)

renewable %>% 
  pad() %>% 
  anyNA()
#> [1] FALSE

Hasil output menunjukkan bahwa tidak terdapat missing value pada dataset, sehingga proses analisis time series dapat dilanjutkan ke tahap berikutnya.

Analisis data eksplorasi

Bagian ini akan menjelajah dataset yang lebih dalam, sehingga dapat diketahui karakteristik awal dari data. Fungsi range digunakan untuk mengetahui data paling awal dan akhirnya.

range(renewable$date)
#> [1] "1990-12-31" "2015-12-31"

Data awal dari dataset renewable ini adalah pada tahun 1990 dan berakhir di tahun 2015, dengan pengambilan data pada setiap akhir tahun.

Dataset renewable dibuat menjadi bertipe data time series. Isian start diisi pada data tahun pertama, frequency diisi dengan satu karena data tahunan.

renewable_ts <- ts(renewable$X..of.Electricity.from.Renewables, start = 1990, frequency = 1)

Selanjutnya dilihat plot dari dataset. Plot renewable energy dapat ditampilkan sebagi berikut.

library(forecast)

renewable_ts %>% 
  autoplot()

Berdasarkan hasil plot, data historical renewable energy cenderung mengalami penurunan. Nilai paling rendah terjadi pada tahun 1997. Data secara visual belum stasioner, namun untuk mengetahui secara statistik, stasioneritas data dapat dilihat melalui uji Augmented Dickey-Fuller.

library(tseries)

renewable_ts %>% 
  adf.test()
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -3.1974, Lag order = 2, p-value = 0.1162
#> alternative hypothesis: stationary

Hipotesis uji Augmented Dickey-Fuller:

-\(H_0\): data tidak stasioner

-\(H_1\): data stasioner

Uji Augmented Dickey-Fuller pada data renewable energy menunjukkan hasil yaitu data belum stasioner, karena p-value = 0.1162 masih lebih tinggi dibandingkan alpha = 0.05, sehingga \(H_0\) gagal ditolak.

Data renewable energy dilakukan differencing untuk mengatasi data belum stasioner. Proses differencing dilakukan sebanyak dua kali. Ini karena proses differencing satu kali masih belum dapat membuat data stasioner.

renewable_ts %>% 
  diff() %>% 
  diff() %>% 
  adf.test()
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -5.051, Lag order = 2, p-value = 0.01
#> alternative hypothesis: stationary

Uji Augmented Dickey-Fuller pada data renewable energy dengan data telah dilakukan differencing seanyak dua kali menunjukkan bahwa data telah stasioner, karena p-value = 0.01 lebih kecil dibandingkan alpha = 0.05, sehingga \(H_0\) ditolak.

renewable_ts %>% 
  diff() %>% 
  diff() %>% 
  autoplot()

Plot yang dihasilkan setelah data dilakukan differencing dua kali telah stabil, sehingga secara visual data telah stasioner.

Cross validation

Dataset renewable energy akan dilakukan cross validation, yaitu pembagian dataset menjadi data train dan data test. Data test diambil pada 10 tahun terakhir dari dataset dan sisa datanya adalah menjadi data train.

train_renewable <- head(renewable_ts, n = length(renewable_ts) - 1*10)
test_renewable <- tail(renewable_ts, n = 1*10)

Pembuatan dan evaluasi model

Bagian ini membahas pembuatan model menggunakan exponential smoothing dan ARIMA. Model-model yang dibuat kemudian dibandingkan bagaimana kinerjanya. Kinerja yang paling baik yaitu yang memiliki error (menggunakan Mean Absolute Percentage Error (MAPE)) paling kecil. MAPE adalah salah satu cara dalam statistika untuk mengukur kinerja model. MAPE dihitung berdasarkan rata-rata dari nilai absolut dari selisih nilai aktual dan nilai prediksi yang lalu dibagi dengan nilai aktual. Model pertama yang akan dibangun adalah Holt’s Exponential Smoothing. Model ini mengasumsikan data memiliki pola trend.

library(MLmetrics)

model_hw_r <- HoltWinters(train_renewable, gamma = F)
forecast_hw_r <- forecast(model_hw_r, h = 1*10)
mape_hw_r <- MAPE(y_pred = forecast_hw_r$mean, y_true = test_renewable)

Model selanjutnya yang dicobakan adalah ARIMA. Sebelum masuk ke pemodelan ARIMA, akan dicari nilai p, d, dan q yang mungkin menggunakan plot Autocorrelation Function (ACF) dan Partial Autocorrelation Function (PACF). Fungsi yang bisa digunakan adalah tsdisplay().

train_renewable %>% 
  diff() %>% 
  diff() %>% 
  tsdisplay()

  • Dies down/decay: mengalami penurunan atau kenaikan secara lambat

  • cuts off lag: mengalami penurunan atau kenaikan secara cepat

Plot ACF cuts off setelah lag 1, sehingga opsi order q bisa diisi dengan 1. Pada plot PACF, cuts off setelah lag 2, sehingga opsi order p bisa diisi dengan 1 dan 2. Maka kombinasi model ARIMA(p,d,q) yang dapat dibentuk sebagai berikut:

  1. ARIMA(1,2,1)
  2. ARIMA(2,2,1)
arima1_r <- Arima(y = train_renewable, order = c(1,2,1), seasonal = F)
arima2_r <- Arima(y = train_renewable, order = c(2,2,1), seasonal = F)

Lalu dilihat nilai AIC masing-masing model ARIMA.

arima1_r$aic
#> [1] 76.8748
arima2_r$aic
#> [1] 78.05855

Model arima1_r memiliki nilai AIC terkecil, sehingga model ARIMA ini yang akan digunakan untuk melakukan Forecasting dengan model ARIMA. Forecasting dilakukan pada 10 data (tahun) ke depan.

forecast_arima_r <- forecast(arima1_r, h=1*10)
mape_arima_r <- MAPE(forecast_arima_r$mean, test_renewable)

Model ARIMA lain yang bisa dilakukan adalah menggunakan fungsi auto.arima(). Fungsi ini akan secara otomatis memilih model ARIMA yang terbaik berdasarkan nilai AIC, AICc atau BIC.

auto_arima_r <- auto.arima(train_renewable)
forecast_auto_arima_r <- forecast(auto_arima_r, h=1*10)
mape_auto_arima_r <- MAPE(y_pred = forecast_auto_arima_r$mean, y_true = test_renewable)

Setelah dilakukan forecast pada semua model, dilakukan pembuatan plot forecast dari setiap model dan plot data test.

test_renewable %>% 
  autoplot(series = "test") +
  autolayer(forecast_hw_r$mean, series = "forecast exponential smoothing") +
  autolayer(forecast_arima_r$mean, series = "forecast arima") +
  autolayer(forecast_auto_arima_r$mean, series = "forecast auto arima")

Berdasarkan plot forecast, model exponential smoothing (Holt’s Exponential Smoothing) lebih mendekati pada nilai aslinya. Sehingga secara visual, model exponential smoothing adalah model yang terbaik untuk melakukan forecast pada dataset renewable energy.

Semua nilai MAPE yang telah diperoleh dari setiap model, dijadikan dalam satu dataframe, kemudian diurutkan dari nilai terkecil ke terbesar.

perf <- as.data.frame(rbind(mape_hw_r, mape_arima_r, mape_auto_arima_r))
colnames(perf) <- c("MAPE")
perf %>% arrange(MAPE)
#>                         MAPE
#> mape_hw_r         0.08864083
#> mape_auto_arima_r 0.12730579
#> mape_arima_r      0.14863042

Nilai MAPE terkecil adalah model exponential smoothing (Holt’s Exponential Smoothing), dengan nilai MAPE sebesar 0.08864083 atau 8.86%. Nilai ini kurang dari 10%, sehingga model memiliki tingkat akurasi peramalan yang tinggi. Berikut ditampilkan hasil forecast, juga interval konfidensi 80% dan interval konfidensi 95%.

forecast_hw_r
#>      Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
#> 2006      13.142469 9.656627 16.62831 7.8113329 18.47361
#> 2007      12.782369 8.870683 16.69406 6.7999620 18.76478
#> 2008      12.422269 8.126751 16.71779 5.8528416 18.99170
#> 2009      12.062169 7.414411 16.70993 4.9540358 19.17030
#> 2010      11.702069 6.726946 16.67719 4.0932754 19.31086
#> 2011      11.341969 6.059732 16.62421 3.2634838 19.42045
#> 2012      10.981869 5.409417 16.55432 2.4595392 19.50420
#> 2013      10.621769 4.773487 16.47005 1.6775929 19.56595
#> 2014      10.261669 4.149992 16.37335 0.9146658 19.60867
#> 2015       9.901569 3.537390 16.26575 0.1683963 19.63474

Hasil forecast menggunakan exponential smoothing (Holt’s Exponential Smoothing) lebih detail dapat ditampilkan dalam plot berikut ini.

plot(forecast_hw_r)

Hasil forecast ditampilkan dalam garis warna biru, sedangkan interval konfidensi 80% ditampilkan dalam area yang berwarna abu-abu tua, dan interval konfidensi 95% ditampilkan dalam area berwarna abu-abu muda.

Uji asumsi

Bagian ini membahas uji asumsi pada model. Uji asumsi pada analisis time series ada 2, yaitu uji Ljung-Box dan normality residual.

Hipotesis uji Ljung-Box sebagai berikut:

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

Uji Ljung-Box digunakan untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time series.

Box.test(forecast_hw_r$residuals, type = "Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  forecast_hw_r$residuals
#> X-squared = 0.018311, df = 1, p-value = 0.8924

Hasil uji Ljung-Box menunjukkan bahwa model tidak memiliki autokorelasi pada residualnya karena p-value lebih dari alpha = 0.05, sehingga \(H_0\) gagal ditolak. Uji asumsi berikutnya adalah normality residual.

Hipotesis uji normality residual sebagai berikut:

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

Pengecekan normality residual pada hasil forecasting time series bisa dilakukan menggunakan shapiro test.

shapiro.test(forecast_hw_r$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  forecast_hw_r$residuals
#> W = 0.9356, p-value = 0.3649

P-value yang dihasilkan adalah 0.3649, maka lebih dari alpha = 0.05, sehingga \(H_0\) gagal ditolak. Ini berati residual pada model berdistribusi normal.

Simpulan

Peramalan pada persentase penggunaan listrik terbarukan di Indonesia dapat dilakukan menggunakan metode exponential smoothing dan ARIMA. Metode exponential smoothing yang digunakan adalah Holt’s Exponential Smoothing, karena pada dataset hanya terdapat efek trend dan tidak terdapat efek seasonal. ARIMA yang digunakan dapat dibagi menjadi dua, pertama dengan cara mencari parameter p, d, dan q secara manual menggunakan plot ACF dan PACF, dan yang kedua menggunakan fungsi auto.arima() yang akan mencari model ARIMA terbaik berdasarkan nilai AIC.

Model terbaik yang didapatkan adalah saat menggunakan metode exponential smoothing, ini dilihat dari nilai MAPE yang paling kecil, yaitu sebesar 8.86%. Model exponential smoothing memiliki akurasi peramalan yang tinggi karena memiliki nilai MAPE yang kurang dari 10%. Uji asumsi pada model exponential smoothing, yaitu independensi residual dan normalitas residual semuanya signifikan. Sehingga tidak terdapat korelasi residual antar lag dan residual berdistribusi normal.