Pendahuluan
Latar Belakang
Curah hujan adalah ketinggian air hujan dalam tempat yang datar, tidak menguap, tidak meresap dan tidak mengalir. Curah hujan 1 mm, artinya dalam luasan satu meter persegi pada tempat yang datar tertampung air setinggi satu milimeter atau tertampung air sebanyak satu liter dalam jangka waktu tertentu. Beberapa faktor yang dapat memengaruhi curah hujan antara lain : tekanan udara, kelembaban, udara, dan suhu udara. Tekanan udara merupakan tenaga yang bekerja untuk menggerakkan massa udara dalam setiap satuan luas tertentu (Siswanti, 2011). Kelembaban udara adalah banyaknya uap air yang terkandung dalam udara atau atmosfer (Swarinoto & Sugiyono, 2011). Dari faktor-faktor tersebut akan dianalisis hubungan antara faktor kelembaban dan tekanan udara yang memengaruhi curah hujan dengan menggunakan model regresi dengan peubah lag.
Tujuan
Tujuan dari analisis ini adalah untuk 1. Menganalisis hubungan antara kelembaban dan tekanan udara dengan melakukan permodelan menggunakan metode Koyck, distributed lag model, dan autoregressive. 2. Menentukan model terbaik dari analisis menggunakan metode Koyck, distributed lag model, dan autoregressive.
Data Preparation
Package
Pada praktikum kali ini dibutuhkan beberapa Package yang terdiri dari sebagai berikut:
dLagMdynlmMLmetricslmtestcar
lapply(c("dLagM","dynlm","MLmetrics","lmtest","car",
"car"),
library,character.only=T)Data
Dataset yang digunakan dalam praktikum kali ini adalah data Jena Climate. Jena Climate adalah data timeseries yang direkam setiap sepuluh menit oleh Stasiun Cuaca Institut Biogeochemistry Max Planck di Jena, Jerman. Data ini terdiri dari 420551 observasi yang direkam sejak 1 Januari 2009 hingga 31 Desember 2016. Dari total 15 peubah yang tersedia, diambil 2 peubah yakni Tekanan Udara sebagai peubah penjelas X dan Kelembaban sebagai peubah penjelas Y. Nilai observasi yang diambil adalah observasi ke-1 hingga ke-287 atau data dari tanggal 1 Januari 2009 hingga 2 Januari 2009.
Input Data
dataclimate <- read.csv("D://Kuliah//Semester 5//STA1341 Metode Peramalan Deret Waktu//Tugas Individu//P4//jena_climate_2009_2016.csv")
dataclimate <- dataclimate[(1:287),-c(2:7,9,11:15)]
colnames(dataclimate)[1] <- "Time Period"
colnames(dataclimate)[2] <- "VAct"
colnames(dataclimate)[3] <- "SHum"
head(dataclimate)## Time Period VAct SHum
## 1 01.01.2009 00:10:00 3.11 1.94
## 2 01.01.2009 00:20:00 3.02 1.89
## 3 01.01.2009 00:30:00 3.01 1.88
## 4 01.01.2009 00:40:00 3.07 1.92
## 5 01.01.2009 00:50:00 3.08 1.92
## 6 01.01.2009 01:00:00 3.14 1.96Splitting Data
Splitting data menjadi data train dan data test dengan proporsi sebanyak 80 persen dan 20 persen. Data train terdiri dari data ke-1 hingga ke-230. Sedangkan, data test terdiri dari data ke-231 hingga ke-287.
train <- dataclimate[(1:230),]
test <- dataclimate[(231:287),]Time Series
Mengubah format data menjadi format time series
data.ts <- ts(dataclimate)
train.ts <- ts(train)
test.ts <- ts(test)Eksplorasi Data
Plot Time Series Kelembaban
plot(ts(dataclimate$SHum), main = "Time Series Plot of Humidity", xlab = "Period", ylab="Humidity")
points(dataclimate$SHum)
Terlihat dari plot time series peubah Kelembaban di atas, data cenderung membentuk pola gabungan. Dimana yang paling terlihat secara visual adalah pola tren positif dan negatif.
Plot Time Series Tekanan Udara
plot(ts(dataclimate$VAct), main = "Time Series Plot of Pressure", xlab = "Period", ylab="Pressure")
points(dataclimate$VAct)
Plot peubah Tekanan Udara memiliki pola yang sama dengan plot peubah Kelembaban, yakni gabungan.
Korelasi
cor(dataclimate$VAct,dataclimate$SHum)## [1] 0.9999012Scatter Plot Kelembaban dan Tekanan Udara
plot(dataclimate$VAct,dataclimate$SHum,main="Scatter Plot Kelembaban dan Tekanan Udara",
xlab="Tekanan Udara",ylab ="Kelembaban")
Scatter plot di atas menunjukkan hubungan positif antara peubah respon Kelembaban dan peubah penjelas Tekanan Udara. Kedua peubah tersebut memiliki nilai korelasi yang sangat tinggi yakni 0.9999012.
Analisis Data
Model Koyck
model.koyck <- dLagM :: koyckDlm(x=train$VAct, y=train$SHum)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0132656 -0.0037232 0.0001132 0.0032745 0.0215018
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0003666 0.0026559 -0.138 0.890347
## Y.1 0.2018448 0.0566232 3.565 0.000444 ***
## X.t 0.4979291 0.0353136 14.100 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005753 on 226 degrees of freedom
## Multiple R-Squared: 0.9997, Adjusted R-squared: 0.9997
## Wald test: 3.868e+05 on 2 and 226 DF, p-value: < 2.2e-16
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: -0.0004592727 0.4979291 0.2018448AIC dan BIC
AIC(model.koyck)## [1] -1707.515BIC(model.koyck)## [1] -1693.78Forecasting
fore.koyck <- forecast(model = model.koyck, x=test$VAct, h=57)MAPE dan Akurasi
# Data Testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$SHum)
# Akurasi Training
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
c("MAPE Testing"=mape.koyck,"MAPE Training"=mape_train)## $`MAPE Testing`
## [1] 0.003218369
##
## $`MAPE Training.MAPE`
## [1] 0.002019012Distributed Lag
Lag Optimum
colnames(train)[1] <- "t"
colnames(train)[2] <- "Xt"
colnames(train)[3] <- "Yt"
finiteDLMauto(formula = Yt ~ Xt,
data = data.frame(train),
model.type = "dlm", error.type = "AIC",trace=TRUE)## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 1 1 0.19775 -1841.407 -1827.672 4.34679 0.07565 0.99984 7.351456e-06
## 2 2 0.19645 -1832.443 -1815.296 4.29768 0.07822 0.99984 9.140555e-06
## 3 3 0.19754 -1822.782 -1802.232 4.33971 0.09635 0.99984 1.767389e-05
## 4 4 0.19465 -1815.619 -1791.675 4.12010 0.10318 0.99984 2.417806e-05
## 5 5 0.19565 -1807.794 -1780.465 4.16253 0.10627 0.99984 6.320855e-05
## 6 6 0.19765 -1801.944 -1771.240 4.29266 0.12158 0.99984 1.506751e-04
## 7 7 0.19486 -1797.970 -1763.899 4.38390 0.13902 0.99984 2.779240e-04
## 10 10 0.19168 -1791.518 -1747.401 4.57174 0.12911 0.99986 2.722030e-03
## 8 8 0.19650 -1788.286 -1750.857 4.36786 0.13961 0.99984 5.039835e-04
## 9 9 0.19797 -1785.750 -1744.972 4.62034 0.13213 0.99985 4.616845e-03Didapatkan nilai lag optimum yakni sebesar 10
model.dlm2 = dLagM::dlm(x = train$Xt,y = train$Yt , q = 10)
summary(model.dlm2)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0101152 -0.0027118 0.0002825 0.0025190 0.0100050
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0043281 0.0019141 -2.261 0.024782 *
## x.t 0.6059616 0.0082022 73.878 < 2e-16 ***
## x.1 0.0176441 0.0125427 1.407 0.161002
## x.2 -0.0051066 0.0127405 -0.401 0.688969
## x.3 -0.0108525 0.0128626 -0.844 0.399794
## x.4 0.0029436 0.0126784 0.232 0.816631
## x.5 -0.0006329 0.0126442 -0.050 0.960128
## x.6 0.0047867 0.0126674 0.378 0.705908
## x.7 0.0080955 0.0126853 0.638 0.524058
## x.8 -0.0099409 0.0125867 -0.790 0.430548
## x.9 -0.0155837 0.0121965 -1.278 0.202775
## x.10 0.0274656 0.0076400 3.595 0.000405 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003999 on 208 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 1.392e+05 on 11 and 208 DF, p-value: < 2.2e-16
##
## AIC and BIC values for the model:
## AIC BIC
## 1 -1791.518 -1747.401AIC dan BIC
AIC(model.dlm2)## [1] -1791.518BIC(model.dlm2)## [1] -1747.401Forecasting
fore.dlm2 <- forecast(model = model.dlm2, x=test$VAct, h=57)MAPE dan Akurasi
#Data Testing
mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$SHum)
#Data Training
mape_train <- GoF(model.dlm2)["MAPE"]
c("MAPE Testing" = mape.dlm2, "MAPE Training" = mape_train)## $`MAPE Testing`
## [1] 0.002870965
##
## $`MAPE Training.MAPE`
## [1] 0.001405104Model Autoregressive / Dynamic Regression
Lag Optimum
datareg <- data.frame(dataclimate$`Time Period`,dataclimate$VAct,dataclimate$SHum)
colnames(datareg)[1] <- "t"
colnames(datareg)[2] <- "Xt"
colnames(datareg)[3] <- "Yt"
pq <- ardlBoundOrders(data = data.frame(datareg) , formula = Yt ~ Xt )
c(p=pq$p$Xt, q=pq$q) #p:lag x, q:lag y## p q
## 2 2Didapatkan nilai parameter p dan q optimum yakni sama-sama sebesar 2.
MODEL ARDL
#MODEL AUTOREGRESSIVE
model.ardl = ardlDlm(x = train$Xt, y = train$Yt, p = 2 , q = 2)
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 3, End = 230
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0092376 -0.0027777 -0.0002518 0.0030509 0.0104214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0009726 0.0018069 -0.538 0.59095
## X.t 0.6115006 0.0072488 84.359 < 2e-16 ***
## X.1 -0.1145395 0.0403499 -2.839 0.00495 **
## X.2 -0.1941582 0.0402423 -4.825 2.61e-06 ***
## Y.1 0.2055252 0.0635535 3.234 0.00141 **
## Y.2 0.3092346 0.0631683 4.895 1.89e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003912 on 222 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 3.34e+05 on 5 and 222 DF, p-value: < 2.2e-16AIC dan BIC
AIC(model.ardl)## [1] -1873.018BIC(model.ardl)## [1] -1849.013Forecasting
fore.ardl <- forecast(model = model.ardl, x=test$VAct, h=57)MAPE dan Akurasi
# Data Testing
mape.ardl <- MAPE(fore.ardl$forecasts, test$SHum)
# Data Training
mape_train <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.002346139
##
## $MAPE_training.MAPE
## [1] 0.001421713Perbandingan
Plot Forecasting
par(mfrow=c(1,1))
plot(test$VAct, test$SHum, type="b", col="black",ylab="Kelembaban",xlab="Tekanan Udara")
points(test$VAct, fore.koyck$forecasts,col="red")
points(test$VAct, fore.dlm2$forecasts,col="blue")
points(test$VAct, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "koyck","DLM 2", "Autoregressive"), lty=1, col=c("black","red","blue","green"), cex=0.8)
Akurasi
akurasi <- matrix(c(mape.koyck,mape.dlm2, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi## MAPE
## Koyck 0.003218369
## DLM 0.002870965
## Autoregressive 0.002346139Dari tabel di atas terlihat bahwa model dengan metode Autoregressive merupakan model terbaik daripada model dengan metode Koyck dan DLM karena memiliki nilai akurasi MAPE paling rendah di antara yang lain.
Diagnostik Model
Uji Non Autokorelasi
bgtest(model.ardl$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model.ardl$model
## LM test = 11.314, df = 1, p-value = 0.0007693Karena nilai p-value < alpha atau 0.0007693 < 0.05 maka Tolak H0. Dapat disimpulkan bahwa tidak terdapat cukup bukti untuk menyatakan bahwa sisaan saling bebas atau tidak terdapat autokorelasi.
Uji Homogenitas
bptest(model.ardl$model)##
## studentized Breusch-Pagan test
##
## data: model.ardl$model
## BP = 11.632, df = 5, p-value = 0.04019Karena nilai p-value < alpha atau 0.04019 < 0.05 maka Tolak H0. Dapat disimpulkan bahwa tidak cukup bukti untuk menyatakan bahwa ragam sisaan homogen.
Uji Normalitas Sisaan
norm <- shapiro.test(residuals(model.dlm2))## 1 2 3 4 5
## 1.646165e-03 1.858613e-03 -1.179648e-03 6.113611e-03 2.056491e-03
## 6 7 8 9 10
## -1.031231e-03 8.596576e-03 -3.831278e-04 5.502377e-03 6.174741e-03
## 11 12 13 14 15
## 1.795575e-03 2.101403e-03 2.365584e-03 4.531548e-03 -1.207321e-03
## 16 17 18 19 20
## 5.661446e-03 8.344908e-04 -6.149244e-04 2.463759e-03 -7.049354e-04
## 21 22 23 24 25
## -2.127193e-04 1.121050e-03 -7.298008e-03 5.663113e-03 2.872605e-03
## 26 27 28 29 30
## -5.035680e-04 3.662776e-03 5.113973e-03 1.085558e-03 -2.632285e-03
## 31 32 33 34 35
## 4.434800e-03 -1.386874e-03 2.095582e-03 1.703759e-03 4.071044e-03
## 36 37 38 39 40
## -6.215080e-04 6.825773e-03 1.000502e-02 2.988953e-03 8.100165e-04
## 41 42 43 44 45
## 3.474139e-03 -7.239060e-04 -2.667238e-04 1.858746e-04 -2.212436e-03
## 46 47 48 49 50
## -2.989083e-03 5.461446e-03 5.992042e-03 2.343790e-03 -4.758250e-03
## 51 52 53 54 55
## -2.985780e-03 -4.815591e-03 -8.011307e-03 1.448969e-03 2.684591e-03
## 56 57 58 59 60
## -5.584216e-03 -4.551204e-03 3.489531e-03 -6.242286e-03 -6.388921e-03
## 61 62 63 64 65
## -4.985852e-03 -8.796317e-03 -9.296628e-04 -8.221242e-03 -1.237878e-03
## 66 67 68 69 70
## -3.244911e-03 5.394655e-03 3.588883e-03 -1.644377e-03 4.000870e-03
## 71 72 73 74 75
## -1.979348e-03 -4.299522e-03 -3.245873e-03 -5.262428e-03 -4.964873e-03
## 76 77 78 79 80
## 3.346756e-03 -4.018370e-03 -7.923054e-04 2.819917e-03 -2.854435e-04
## 81 82 83 84 85
## 2.794828e-03 -6.010026e-03 -8.516928e-05 -1.890086e-03 4.756422e-03
## 86 87 88 89 90
## -4.570922e-03 7.991428e-04 -6.392501e-03 2.220822e-03 3.969146e-04
## 91 92 93 94 95
## -2.695850e-03 -7.797488e-04 -9.778673e-03 -7.575431e-03 -1.291375e-03
## 96 97 98 99 100
## -4.115073e-03 1.520615e-03 -5.011907e-03 3.533387e-03 -1.011519e-02
## 101 102 103 104 105
## -6.106641e-03 -5.586062e-03 -5.677898e-03 -5.135206e-03 -2.349479e-03
## 106 107 108 109 110
## -7.553973e-03 -4.118490e-03 -9.987080e-04 -7.068357e-03 -2.029245e-03
## 111 112 113 114 115
## 2.108216e-03 -6.142337e-03 1.793519e-03 -7.827497e-03 -3.625248e-03
## 116 117 118 119 120
## -1.102048e-03 1.154522e-03 -9.140047e-04 -2.913865e-03 1.347275e-03
## 121 122 123 124 125
## 8.535280e-04 -3.061310e-03 -2.817422e-04 -3.279247e-03 -6.130042e-03
## 126 127 128 129 130
## -5.453607e-03 3.218593e-04 -1.657755e-03 -5.500231e-03 6.640902e-03
## 131 132 133 134 135
## 1.441100e-03 1.962269e-03 -1.619761e-03 1.429431e-03 1.095649e-03
## 136 137 138 139 140
## -1.061543e-03 -1.198345e-03 -1.886963e-03 3.586024e-03 -5.484567e-03
## 141 142 143 144 145
## 5.454756e-03 -4.148859e-03 1.889057e-03 1.581409e-03 5.165986e-03
## 146 147 148 149 150
## 1.802109e-03 6.845779e-03 1.873544e-03 6.994610e-03 -2.804445e-03
## 151 152 153 154 155
## -1.346548e-03 -4.585835e-03 1.896549e-03 -1.121315e-03 -3.896408e-03
## 156 157 158 159 160
## 1.934670e-03 1.435107e-03 -2.759516e-03 4.343080e-03 3.241686e-03
## 161 162 163 164 165
## 1.144286e-03 -2.057084e-04 2.994431e-03 3.618455e-03 -2.831882e-03
## 166 167 168 169 170
## -1.708235e-03 3.879733e-03 9.651017e-04 8.342163e-03 3.590979e-03
## 171 172 173 174 175
## 1.873423e-03 3.041456e-04 -1.196616e-03 6.638485e-03 2.087030e-03
## 176 177 178 179 180
## 7.638871e-03 3.878264e-03 8.535915e-04 1.990207e-04 1.181353e-03
## 181 182 183 184 185
## 6.443135e-03 5.631827e-04 1.647020e-03 -8.106723e-04 -4.522514e-03
## 186 187 188 189 190
## 3.133963e-03 -1.904445e-03 -2.098939e-03 1.093184e-03 2.608810e-04
## 191 192 193 194 195
## 3.450989e-03 5.208476e-03 2.411451e-03 4.073688e-03 4.779547e-03
## 196 197 198 199 200
## 1.944695e-03 -5.053098e-03 -7.650467e-04 1.481326e-03 2.458269e-03
## 201 202 203 204 205
## -4.691695e-04 -4.505928e-03 2.768008e-03 3.386425e-03 -1.691680e-04
## 206 207 208 209 210
## 4.914194e-03 1.745889e-03 -8.010411e-04 3.763867e-03 4.615172e-04
## 211 212 213 214 215
## 6.491145e-03 -1.290952e-03 9.322591e-04 -7.184832e-04 -3.237175e-03
## 216 217 218 219 220
## 4.797366e-04 -3.268980e-03 9.227664e-04 -1.499638e-03 -1.596579e-03norm$p.value## [1] 0.3848119Karena nilai p-value > alpha atau 0.3848119 > 0.05 maka Tak Tolak H0. Dapat disimpulkan bahwa terdapat cukup bukti untuk menyatakan bahwa sisaan menyebar normal.
Kesimpulan
Model terbaik yang didapatkan untuk data Jena Climate adalah model yang menggunakan metode Autoregressive. Model ini dipilih karena memiliki nilai MAPE yang paling kecil dibandingkan model dengan metode yang lain. Hasil uji diagnostik menunjukkan bahwa model dengan metode Autoregressive ini memiliki gejala autokorelasi, ragam sisaan yang tidak homogen, dan sisaan yang menyebar normal. Hal ini berarti diperlukan memerlukan penanganan lebih lanjut untuk mengatasi gejala-gejala tersebut.
Daftar Pustaka
https://www.kaggle.com/datasets/mnassrib/jena-climate
Gujarati, D. N. 2004. Basic Econometrics.Fourth Edition. New York: The McGraw-Hill.
Gujarati, D. N. 2006. Essentials of Econometrics. Third Edition. New York: The McGraw-Hill.
Prakoso, Dipa. 2018. Analisis pengaruh tekanan udara kelembaban udara dan suhu udara terhadap tingkat curah hujan di kota Semarang [skripsi]. Semarang (ID) : Universitas Negeri Semarang.