Pengolahan Data Time Series

Pada kesempatan kali dilakukan pengolahan data menggunakan model time series khususnya model ARIMA (Autoregressive Integrated Moving Average). ARIMA cocok jika observasi dari deret waktu (time series) secara statistik berhubungan satu sama lain (dependent).

# Library yang digunakan pada pengolahan data
## Untuk membaca data dengan format Excel
library(readxl)
## Untuk mendapatkan informasi dari data
library(psych)
## Untuk mengubah data menjadi Time Series, Acf, Pacf, Model Arima dan ADF Test
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# Untuk melihat signifikansi dari Koefesien
library(lmtest) 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Untuk memprediksi data nantinya
library(forecast)

Informasi Data

Data penderita DBD di Kota Jakarta Pusat secara bulanan dari Januari 2008 - November 2017 yang diperoleh dari Dinas Kesehatan. Jumlah data adalah 515 data. Data yang digunakan sebagai percobaan dalam pemodelan sebanyak 490, sedangkan sisanya digunakan untuk perbandingan data asli dengan hasil prediksi.

# Input Data
data_asli <- read_excel("D:/SEMESTER 7/Kapseltat/Tugas Kelompok Time Series/Data Mingguan Penderita DBD JakPus.xlsx")
# Mengubah data training yaitu data yang dimodelkan kedalam bentuk time series
data_asli <- ts(data_asli)
data <- ts(data_asli[1:490,])
# Plot Data
plot(data, main = "Plot Data Penderita DBD JakBar 2008-2017", ylab = "Penderita DBD", xlab = "Waktu" , type = 'o')

Berdasarkan plot deret waktu tersebut data penderita DBD di Kota Jakarta Pusat secara bulanan dari Januari 2008 - November 2017 sangatlah beragam. Terlihat bahwa sekilas data tidak stasioner, yaitu terdapat variansi yang cukup besar.

#Statistika Deskriptif Data
describe(data)
##    vars   n  mean    sd median trimmed   mad min max range skew kurtosis   se
## X1    1 490 43.09 34.29     33   38.19 26.69   1 219   218 1.45     2.52 1.55

Rata-rata penderita DBD di Kota Jakarta Pusat secara bulanan dari Januari 2008 - November 2017 sebesar 43 orang dengan banyak penderita terkecil hanya seorang dan terbesar sebanyak 219 orang pada Maret 2016. Penderita DBD di Kota Jakarta Pusat dalam rentang waktu tersebut juga sangat beragam ditandai dengan besarnya standar deviasi. Kemencengan (skewness) positif yang artinya data menumpuk di nilai yang lebih kecil dari rataan. Kemudian keruncingan (kurtosis) < 3 yaitu distribusi lebih landai dari distribusi normal.

Untuk mengecek kestasioneran data lebih lanjut, data akan diuji dengan 3 cara.

# CEK STASIONERITAS
##Plot Data
plot(data, main = "Plot Data Penderita DBD JakBar 2008-2017")
abline(h = mean(ts(data)))

## Plot ACF
acf(data) 

## Uji ADF (Augmented Dickey-Fuller)
adf.test(data)
## Warning in adf.test(data): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data
## Dickey-Fuller = -5.0258, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
#jika p-value < alpha maka stasioner
#jadi asumsi stasioner dulu

Perhatikan dari plot data terlihat bahwa sekilas data tidak stasioner, yaitu terdapat variansi yang cukup besar. Data tidak stasioner juga dibuktikan pada plot ACF (Autocorrelation Function) data yang membentuk pola dan menurun secara perlahan. Namun, pada berdasarkan uji ADF (Augmented Dickey-Fuller), dengan 𝛂 = 0.05 diperoleh p-value = 0,01 < 𝛂 sehingga disimpulkan bahwa data stasioner.

Maka selanjutnya dimodelkan data asli (Model 1) dan data yang sudah didiferensi (Model 2).

Menentukan Modet Deret Waktu (Model 1)

Selanjutnya dilakukan estimasi model dengan melihat plot ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function).

#ACF Manual
T= length(data)
maxLag=12
df <- data.frame(1:T,data)
acfResults <- data.frame(k = c(rep(1:maxLag)), acf=NA, stringsAsFactors = F)
for (k in 1: maxLag) {
  N <- nrow(df)-k
  acf <- sum((data[1:N]-mean(ts(data)))%*%(data[(1+k):(N+k)]-mean(data)))/sum((data-mean(data))^2)
  acfResults[k,1] <- k
  acfResults[k,2] <- acf
}
#ACF dengan R
acf(data, plot = FALSE)
## 
## Autocorrelations of series 'data', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.865  0.836  0.781  0.766  0.688  0.615  0.566  0.500  0.464  0.373 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.320  0.257  0.225  0.160  0.115  0.073  0.029  0.004 -0.039 -0.062 -0.094 
##     22     23     24     25     26 
## -0.104 -0.150 -0.166 -0.167 -0.160
#Plot ACF Manual dan dengan R
par(mfrow=c(1,2))
plot(acfResults[,2], type = "h", ylim = c(-.2,1), xlab = "lag", ylab = "ACF", main = 'Plot ACF manual')
abline(h = 0)
abline(h = c(1,-1)*1.96/sqrt(length(data)), lty = 2, col = "blue")

acf(data,main="Plot ACF dengan R")

#PACF Manual
PacfResults <- data.frame(k = c(1:maxLag), pacf = NA, stringsAsFactors = FALSE)
rho = acfResults[,2]
phi = matrix(0, nrow = maxLag, ncol = maxLag)
phi[1,1] = rho[1]
for (i in 2:maxLag){
  phi[i,i] = (rho[i]-(sum(phi[i-1,1:i-1]*rho[(i-1):1])))/(1-sum(phi[i-1,1:i-1]*rho[1:i-1]))
  for (j in 1:maxLag){
    if (j<i){
      phi[i,j] = phi[i-1,j]-(phi[i,i]*phi[i-1,i-j])
    }
    else {0}
  }
  PacfResults[k,1] <- k
  PacfResults[,2] <- diag(phi)
}
#PACF dengan R
pacf(data, plot = FALSE)
## 
## Partial autocorrelations of series 'data', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.865  0.347  0.037  0.136 -0.168 -0.192  0.016 -0.092  0.084 -0.138 -0.064 
##     12     13     14     15     16     17     18     19     20     21     22 
## -0.010  0.041 -0.036  0.012 -0.016 -0.064  0.039 -0.022 -0.012  0.003 -0.008 
##     23     24     25     26 
## -0.120  0.000  0.095  0.085
#Plot PACF Manual dan dengan R
par(mfrow=c(1,2))
plot(PacfResults[,2], type = "h", ylim = c(-.2,1), xlab = "lag", ylab = "PACF", main = 'Plot PACF manual')
abline(h=0)
abline(h=c(1,-1)*1.96/sqrt(length(data)),lty=2,col="blue")

pacf(data,main='Plot PACF dengan R')

Berdasarkan gambar plot ACF dan PACF di atas, terlihat bahwa plot ACF menurun secara eksponensial, sedangkan pada plot PACF terdapat cut-off setelah lag-2, maka didapatkan model ARIMA(2,0,0) yang artinya model hanya mengandung Autoregressive (AR).

Selanjutnya, akan dicari estimasi parameter untuk model pertama ini, yaitu ARIMA(2,0,0).

#Estimasi Parameter
model_1 <- arima(ts(data),order=c(2,0,0))
summary(model_1)
## 
## Call:
## arima(x = ts(data), order = c(2, 0, 0))
## 
## Coefficients:
##          ar1     ar2  intercept
##       0.5614  0.3501    42.2670
## s.e.  0.0424  0.0425     7.9625
## 
## sigma^2 estimated as 257.5:  log likelihood = -2056.09,  aic = 4120.18
## 
## Training set error measures:
##                        ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.002586618 16.04667 10.48246 -22.94941 40.18491 0.9355582
##                      ACF1
## Training set -0.009452142

Berdasarkan hasil identifikasi model, diperoleh nilai AIC sebesar 4120.18.

Menentukan Modet Deret Waktu (Model 2)

Data asli didiferensi sekali dan kemudia di cek stasioneritasnya kembali.

# Cara differensi data jika data tidak stasioner
data2 <- diff(ts(data))
plot(data2, main = "Plot Data Differensi 1x")
abline(h = mean(data2))

# CEK STASIONERITAS
##Plot Data
plot(data2, main = "Plot Data Penderita DBD JakBar 2008-2017")
abline(h=mean(ts(data2)))

## Plot ACF
acf(data2) 

## Uji ADF (Augmented Dickey-Fuller)
adf.test(data2)
## Warning in adf.test(data2): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data2
## Dickey-Fuller = -7.8722, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

Dari plot data yang telah didiferensi sekali terlihat data cukup stasioner dibanding data asli sebelumnya. Kemudian, pada plot ACF (Autocorrelation Function) tidak membentuk pola yang signifikan. Begitu pula, pada hasil uji ADF (Augmented Dickey-Fuller), dengan 𝛂 = 0.05 diperoleh p-value = 0,01 < 𝛂 sehingga disimpulkan bahwa data yang telah didiferensi sekali telah stasioner.

Selanjutnya dilakukan estimasi model dengan melihat plot ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function).

acf(data2,main="Plot ACF dengan R")

pacf(data2,main='Plot PACF dengan R')

Dari hasil plot ACF dan PACF di atas, model yang mungkin cocok adalah ARIMA(1,1,1). Selanjutnya, dapat digunakan auto arima di R sebagai pertimbangan dalam menentukan model. Hasilnya adalah ARIMA(5,1,2).

model_auto = auto.arima(data)
summary(model_auto)
## Series: data 
## ARIMA(5,1,2) 
## 
## Coefficients:
##          ar1      ar2      ar3     ar4      ar5      ma1     ma2
##       0.3836  -0.6718  -0.4012  0.0956  -0.1400  -0.8392  0.9452
## s.e.  0.0521   0.0561   0.0575  0.0521   0.0485   0.0258  0.0339
## 
## sigma^2 estimated as 245.8:  log likelihood=-2036.94
## AIC=4089.88   AICc=4090.18   BIC=4123.42
## 
## Training set error measures:
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.09906185 15.55067 10.41434 -13.57487 36.26411 0.9294784
##                   ACF1
## Training set 0.0221461

Berdasarkan ketiga model ARIMA di atas yaitu ARIMA(2,0,0), ARIMA(1,1,1), dan ARIMA(5,1,2) kemudian ditentukan model mana yang cukup baik. Dapat dilihat melalui nilai AIC,MAPE beserta RMSE dari masing-masing model. Diperoleh AIC dan RMSE paling kecil terletak pada ARIMA(5,1,2), sedangkan MAPE paling kecil ada pada ARIMA(1,1,1). Maka, dipilih model ARIMA(5,1,2) sebagai model yang cocok untuk data yang dipunya.

Signifikansi dari Parameter Model

coeftest(model_auto) 
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.383631   0.052145   7.3570 1.880e-13 ***
## ar2 -0.671767   0.056091 -11.9763 < 2.2e-16 ***
## ar3 -0.401218   0.057499  -6.9778 2.999e-12 ***
## ar4  0.095590   0.052087   1.8352  0.066478 .  
## ar5 -0.139983   0.048480  -2.8875  0.003883 ** 
## ma1 -0.839207   0.025794 -32.5347 < 2.2e-16 ***
## ma2  0.945169   0.033936  27.8517 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dapat dilihat bahwa pvalue dari koefisien ar4 > α maka dapat kita simpulkan bahwa hanya koefisien ar4 yang tidak siginifikan pada model yang dimiliki. Sehingga model memiliki semua parameter kecuali ar4. Kemudian, langkah selanjutnya dilakukan uji diagnostik residual data.

#Uji Diagnostik
qqnorm(residuals(model_auto))
qqline(residuals(model_auto))

hist(residuals(model_auto))

shapiro.test(residuals(model_auto))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_auto)
## W = 0.90254, p-value < 2.2e-16
Box.test(residuals(model_auto),type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(model_auto)
## X-squared = 0.24179, df = 1, p-value = 0.6229
checkresiduals(residuals(model_auto))
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

Berdasarkan QQ Plot, terlihat pada awal dan akhir data nilainya cukup jauh dari garis. Sedangkan pada histogram sekilas data residual berdistribusi normal dan terdapat nilai ekstrim tinggi data. Untuk menguji kenormalan digunakan Uji Shapiro-Wilk diperoleh p-value < 0.5. Artinya untuk uji kenormalan Shapiro-Wilk 𝐻0 ditolak, sehingga tidak berdistribusi normal. Selanjutnya, pada grafik residual ACF masih terdapat lag sehingga dapat diasumsikan adanya korelasi antar residual. Untuk pembuktian lebih lanjut, gunakan Uji Ljung-Box diperoleh p-value < 0.05, sehingga 𝐻0 ditolak. Artinya, ada korelasi antar residual. Adanya korelasi dan data tidak berdistribusi normal mengakibatkan tidak terpenuhinya asumsi residual dalam pemodelan.

Pemeriksaan Data pada data Test

#grafik observasi
plot(data_asli,col="black", ylab="Penderita DBD JakPus", xlab="waktu",type='l') 
#grafik model
lines(fitted.values(model_auto), type="l",lty=2,col="red")

Prediksi

Selanjutnya diperiksa prediksi berdasarkan model terpilih.

#Prediksi
prediksi = forecast(model_auto, h = 5)
autoplot(prediksi)

summary(prediksi)
## 
## Forecast method: ARIMA(5,1,2)
## 
## Model Information:
## Series: data 
## ARIMA(5,1,2) 
## 
## Coefficients:
##          ar1      ar2      ar3     ar4      ar5      ma1     ma2
##       0.3836  -0.6718  -0.4012  0.0956  -0.1400  -0.8392  0.9452
## s.e.  0.0521   0.0561   0.0575  0.0521   0.0485   0.0258  0.0339
## 
## sigma^2 estimated as 245.8:  log likelihood=-2036.94
## AIC=4089.88   AICc=4090.18   BIC=4123.42
## 
## Error measures:
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.09906185 15.55067 10.41434 -13.57487 36.26411 0.9294784
##                   ACF1
## Training set 0.0221461
## 
## Forecasts:
##     Point Forecast      Lo 80    Hi 80     Lo 95    Hi 95
## 491       10.84053  -9.253166 30.93422 -19.89012 41.57118
## 492       10.51339 -12.365174 33.39196 -24.47636 45.50315
## 493       14.93132 -11.343940 41.20658 -25.25322 55.11586
## 494       14.15296 -14.637697 42.94362 -29.87855 58.18448
## 495       12.49798 -20.242388 45.23834 -37.57409 62.57004

KESIMPULAN

Model dari pengolahan data penderita DBD di Kota Jakarta Pusat secara bulanan dari Januari 2008 - November 2017 diperoleh bahwa model yang terpilih adalah ARIMA(5,1,2) dan memiliki RMSE sebesar ≈15% sehingga model ini dapat digunakan untuk memprediksi data ke depannya.

Namun, perhatikan pada uji diagnostik model ARIMA(5,1,2) sebelumnya menunjukkan adanya korelasi antar residual dan residual model tidak berdistribusi normal sehingga asumsi residual dalam pemodelan tidak dipenuhi. Maka, akan lebih baik data dimodelkan dengan data yang memperhatikan korelasi antar residual datanya.