Eksplorasi dan Pemulusan Data Rata-rata suhu di Indonesia pada tahun 1950-2021
Kelompok 2 Paralel 1
M. Nachnoer Novatron F.A. (G1401201014)
Herdian Partawijaya (G1401201022)
Fani Srinoverti (G1401201023)
Hana Faiza Amalina (G1401201051)
Hanaa Budinur Izdihaar (G1401201076)
A. Latar Belakang
Perubahan iklim merupakan suatu kondisi yang ditandai dengan berubahnya pola iklim dunia yang mengakibatkan fenomena cuaca yang tidak menentu. Perubahan iklim terjadi karena adanya perubahan faktor iklim, seperti suhu udara dan curah hujan yang terjadi secara terus menerus dalam jangka waktu yang panjang (Kementerian Lingkungan Hidup, 2004). Perubahan iklim juga dipengaruhi oleh kondisi cuaca yang tidak stabil, seperti curah hujan yang tidak menentu, badai yang sering terjadi, suhu udara yang ekstrim, serta arah angin yang berubah drastis (Ratnaningayu, 2013).
Laporan yang dikeluarkan pada tahun 2001 oleh Intergovernmental Panel on Climate Change menyimpulkan bahwa temperatur udara global telah meningkat sebesar 0.6 derajat Celcius sejak 1861. Pemanasan tersebut disebabkan oleh aktivitas manusia yang menambah gas-gas rumah kaca ke atmosfer. Pemenang hadiah Nobel, Paul Crutzen, berpendapat bahwa ada periode geologis baru dalam sejarah manusia, yaitu antroposen. Pada periode baru ini, aktivitas manusialah yang memberikan dampak besar pada hampir semua aspek bumi.
Masyarakat saat ini dihadapkan pada perubahan iklim dunia. Salah satu dampak yang paling terasa adalah efek rumah kaca, yakni peningkatan kadar gas rumah kaca di atmosfer. Hansen (2011) menyebutkan bahwa terdapat peningkatan suhu di darat dan laut serta peningkatan emisi memiliki pengaruh atas peningkatan suhu ini. Jika suhu terus meningkat saat emisi meningkat, maka kehidupan manusia, hewan, dan tumbuhan akan berubah.
B. Metadata
Data Temperatur Rata-Rata Indonesia (Y)
Data yang digunakan pada penulisan ini adalah data rata-rata suhu di Indonesia pada tahun 1950-2021 dengan satuan °C. Data tersebut diperoleh dari: climateknowledgeportal.wordbank.org.
Data Populasi Indonesia (X)
Kemudian data yang digunakan sebagai peubah penjelas atau independen adalah data populasi penduduk tahunan di Indonesia pada tahun 1950-2021. Data tersebut diperoleh dari: climateknowledgeportal.wordbank.org
C. Call Package
lapply(c("readxl","readr","Mcomp","smooth","fpp2","tseries","lmtest",
"dLagM","TTR","TSA","dplyr","MLmetrics"),library, character.only=T)[[1]]## [1] "readxl" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"D.Data Awal
1. Data Peubah Respon
CEK <- read_csv2("C:/Users/falco/Downloads/datasuhu.csv")## i Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.## Rows: 121 Columns: 3
## -- Column specification --------------------------------------------------------
## Delimiter: ";"
## dbl (3): Tahun, Rata-rata suhu, 5-yr smooth
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.CEK <- CEK[,-3]
CEK <- CEK[-1:-49,]
glimpse(CEK)## Rows: 72
## Columns: 2
## $ Tahun <dbl> 1950, 1951, 1952, 1953, 1954, 1955, 1956, 1957, 1958,~
## $ `Rata-rata suhu` <dbl> 25.50, 25.52, 25.49, 25.55, 25.42, 25.38, 25.28, 25.5~head(CEK)## # A tibble: 6 x 2
## Tahun `Rata-rata suhu`
## <dbl> <dbl>
## 1 1950 25.5
## 2 1951 25.5
## 3 1952 25.5
## 4 1953 25.6
## 5 1954 25.4
## 6 1955 25.42. Data Peubah Penjelas
Data yang digunakan sebagai peubah penjelas adalah data populasi Indonesia dengan periode 1950 - 2021. Penggunaan data populasi Indonesia digunakan untuk menjawab tujuan yaitu penentuan persamaan model regresi deret waktu antara suhu rata-rata dengan variabel penjelas di Indonesia menggunakan metode metode terbaik. Berdasarkan studi literatur dan eksplorasi data, data Rata-rata suhu dan populasi di Indonesia mempunyai hubungan linear yang bagus dengan korelasi yang tinggi.
pop <- read_xlsx("C:/Users/falco/Downloads/pop.xlsx")
pop <- pop[,-3]
pop$Populasi <- pop$Populasi/10^6
glimpse(pop)## Rows: 72
## Columns: 2
## $ Tahun <dbl> 1950, 1951, 1952, 1953, 1954, 1955, 1956, 1957, 1958, 1959, 1~
## $ Populasi <dbl> 69.56762, 71.01935, 72.57083, 74.20799, 75.92533, 77.74150, 7~head(pop)## # A tibble: 6 x 2
## Tahun Populasi
## <dbl> <dbl>
## 1 1950 69.6
## 2 1951 71.0
## 3 1952 72.6
## 4 1953 74.2
## 5 1954 75.9
## 6 1955 77.73. Penggabungan Data Peubah Respon dan Peubah Penjelas
Data <- merge(CEK,pop)
Data <- as.data.frame(Data)
head(Data)## Tahun Rata-rata suhu Populasi
## 1 1950 25.50 69.56762
## 2 1951 25.52 71.01935
## 3 1952 25.49 72.57083
## 4 1953 25.55 74.20799
## 5 1954 25.42 75.92533
## 6 1955 25.38 77.741504. Data Training dan Testing dari Peubah Respon
Splitting data dilakukan untuk membagi data awal menjadi data training dan data testing dengan proporsi data training 70% dan data testing 30%. Data training merupakan data yang digunakan untuk mencari model yang sesuai dengan data melalui serangkaian algoritma tertentu atau digunakan juga untuk tahapan modelling.Data testing merupakan data yang digunakan untuk mengetahui akurasi performa model yang diperoleh dari data training atau digunakan juga untuk tahapan forecasting.
train<-CEK[1:51,]
test<-CEK[-1:-51,]
#data time series
datasuhu.ts<-ts(CEK$`Rata-rata suhu`)
datapop.ts<-ts(pop$Populasi)
training.ts<-ts(train$`Rata-rata suhu`)
testing.ts<-ts(test$`Rata-rata suhu`)- Awalnya saat menggunakan data awal Rata-rata Suhu periode 1901-2021 splitting data dilakukan dengan proporsi 80% training dan 20% testing. Jika splitting di 70% dan 30%, data training isinya cenderung stasioner semua sedangkan data testing polanya trend. Jadi untuk data awal (1901-2021) agar peramalannya akurat maka menggunakan splitting 80% training dan 20% testing.
- Tapi setelah menggunakan data yg sekarang, yaitu data Rata-rata suhu periode 1950-2021 (menyesuaikan dengan periode data peubah penjelas “populasi”) proporsi 80% train dan 20% test itu justru menurunkan rsquare dari model regresi time series dan tidak sesuai dengan pola datanya.
- Oleh sebab itu kembali menggunakan 70% train dan 30% test, dan setelah dicek untuk data yg sekarang, trend data training tidak monoton hanya pola stasioner saja, sehingga akan memperbesar peluang mendapatkan hasil peramalan yang bagus karena diduga akan mendekati pola data testing (20 tahun terakhir) sekaligus mendapatkan model regresi dengan performa yang tinggi.
E. Eksplorasi Data
1. Time Series Plot Data Peubah Respon dan Peubah Penjelas
g1<-ggplot(CEK, aes(x=Tahun, y=`Rata-rata suhu`)) +
geom_line(lty=1, lwd=0.8, col="black") +
labs(x="Tahun",y = "Temperatur", title="Time Series Plot Data Peubah Respon",
subtitle = "Periode 1950 - 2021") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5)) +
geom_point(size=1.2)
g2<-ggplot(pop, aes(x=Tahun, y=Populasi)) +
geom_line(lty=1, lwd=0.8, col="black") +
labs(x="Tahun",y = "Populasi (Juta jiwa)", title="Time Series Plot Data Peubah Penjelas",
subtitle = "Periode 1950 - 2021") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5)) +
geom_point(size=1.2)
gridExtra::grid.arrange(g1,g2,ncol=2)
- Plot di atas menunjukkan data deret waktu dari temperatur rata-rata tahunan di Indonesia dalam rentang waktu 1950 sampai 2021. Dapat dilihat bahwa data ini memiliki pola gabungan yaitu pola data horizontal dan pola data trend positif.
- Pada rentang periode 30 tahun pertama, terlihat pola data yang cenderung horizontal karena data berfluktuasi di sekitar rata-rata. Setelahnya, terjadi peningkatan temperatur rata-rata yang cukup signifikan sampai tahun 2021. Hal tersebut terjadi karena kegiatan industri yang meningkat dengan begitu pesat dan kegiatan manusia yang banyak menghasilkan gas rumah kaca menjadi penyebab utama pemanasan global (Ramlan 2002).
- Plot data populasi di atas menunjukkan data deret waktu dari populasi tahunan di Indonesia dalam rentang waktu 1950 sampai 2021. Dapat dilihat bahwa data ini memiliki pola trend yang positif dari awal hingga akhir periode.
2. Time Series Plot Data Training dari Peubah Respon
ggplot(train, aes(x=Tahun, y=`Rata-rata suhu`)) +
geom_line(lty=1, lwd=0.8, col="black") +
labs(x="Tahun",y = "Temperatur", title="Time Series Plot Data Training",
subtitle = "Periode 1950 - 2000") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5)) +
geom_point(size=1.2)
3. Time Series Plot Data Testing dari Peubah Respon
ggplot(test, aes(x=Tahun, y=`Rata-rata suhu`)) +
geom_line(lty=1, lwd=0.8, col="black") +
labs(x="Tahun",y = "Temperatur", title="Time Series Plot Data Testing",
subtitle = "Periode 2001 - 2021") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5)) +
geom_point(size=1.2)
Scatter Plot Peubah Penjelas dan Peubah Respon
plot(Data$Populasi, Data$`Rata-rata suhu`, col="navyblue", lwd=2, type="p",xlab="Populasi",ylab="Temperatur",main="Scatter Plot Populasi dan Temperatur")
abline(lm(Data$`Rata-rata suhu` ~ Data$Populasi), col = "red", lwd = 2)
Dari scatter plot di atas dapat dilihat bahwa peubah Populasi dan peubah Rata-rata Suhu mempunyai hubungan linear positif, artinya semakin meningkat jumlah populasi di Indonesia maka rata-rata suhunya juga akan meningkat.
Nilai Korelasi
cor(Data$Populasi, Data$`Rata-rata suhu`)## [1] 0.8704546Scatter plot di atas menunjukkan bahwa peubah penjelas (populasi) dan peubah respon (temperatur) memiliki hubungan linear positif yang kuat dengan nilai koefisien korelasi sebesar 0,87.
F. Pemulusan
1. Single Moving Average (SMA)
a. Mencari parameter n optimum
#fungsi untuk mencari parameter optimum SMA
nopsma<-function(training.ts, min, max, metode=c("SSE","MSE","RMSE","MAD","MAPE"),axis.y.label=" "){
par(mfrow=c(2,ceiling((max-min+1)/2)))
d<-data.frame()
for(n in min:max){
sma<-SMA(training.ts,n)
ramal<-c(NA,sma)
df<-cbind(aktual=c(training.ts,NA),mulus=c(sma,NA),ramal)
error<-df[,1]-df[,3]
sse<-sum(error^2,na.rm=T)
mse <- mean(error^2,na.rm=T)
rmse<-sqrt(mse)
mad<-mean(abs(error),na.rm=T)
rerror<-error/df[,1]*100
mape<-mean(abs(rerror), na.rm = T)
ak<-data.frame(n=n,SSE=sse,MSE=mse,RMSE=rmse,MAD=mad,MAPE=mape)
d<-rbind(d,ak)
ts.plot(training.ts, lty=3, col="black", main=paste("n =",n), ylab=axis.y.label);points(training.ts);lines(sma,
col="red");lines(ramal, col="green", lty=1)
legend("topleft",c("Aktual","Pemulusan","Ramal"),
lty =c(3,1,1) ,col=c("black","red","green"))
}
opt<-d$n[which.min(d[,metode])]
return(list(Akurasi=d,n.optimum=paste("n optimum adalah",opt,"dengan",metode,"sebesar",round(d[which(d$n==opt),metode],3))))
}
#parameter optimum
SMA<-nopsma(training.ts,1,10,"MAPE","Rataan Suhu");SMA
## $Akurasi
## n SSE MSE RMSE MAD MAPE
## 1 1 1.4127000 0.02825400 0.1680893 0.1302000 0.5088765
## 2 2 1.2395750 0.02529745 0.1590517 0.1278571 0.4994362
## 3 3 1.1593778 0.02415370 0.1554146 0.1236111 0.4827585
## 4 4 1.0926438 0.02324774 0.1524721 0.1222872 0.4775863
## 5 5 1.0665920 0.02318678 0.1522721 0.1223478 0.4777157
## 6 6 1.0368972 0.02304216 0.1517964 0.1228519 0.4793815
## 7 7 0.9837714 0.02235844 0.1495274 0.1194156 0.4656674
## 8 8 0.9549937 0.02220916 0.1490274 0.1182558 0.4609680
## 9 9 0.9032099 0.02150500 0.1466458 0.1160847 0.4523759
## 10 10 0.9073360 0.02213015 0.1487620 0.1203415 0.4689074
##
## $n.optimum
## [1] "n optimum adalah 9 dengan MAPE sebesar 0.452"Dengan melakukan fungsi pengulangan untuk metode SMA diperoleh n optimum yaitu 9 dengan ukuran performa model MAPE sebesar 0.425
b. Evaluasi Performa SMA Menggunakan Data Testing
df.ts.sma <- TTR::SMA(testing.ts, n=9)
ramal.df.sma <- c(NA,df.ts.sma)
df.sma <- cbind(aktual.sma=c(testing.ts, NA), pemulusan.sma=c(testing.ts, NA), ramal.df.sma)
error.sma <- df.sma[,1]-df.sma[,3]
rerror.sma <- (error.sma/df.sma[,1])*100
mape.sma <- mean(abs(rerror.sma), na.rm=T);mape.sma## [1] 0.45571662. Double Moving Average (DMA)
Parameter optimum
nopdma<-function(training.ts, min, max, metode=c("SSE","MSE","RMSE","MAPE","MAD"),axis.y.label=" ",just.graph=F,just.output=F){
d<-data.frame()
par(mfrow=c(2,ceiling((max-min+1)/2)))
for(n in min:max){
sma<-SMA(training.ts,n)
ddma <- SMA(sma, n)
At <- 2*sma - ddma
Bt <- 2/(n-1)*(sma - ddma)
dma<- At+Bt
ramal<- c(NA, dma)
data.gab2 <- cbind(aktual = c(training.ts,NA), pemulusan1 = c(sma,NA),
pemulusan2 = c(dma, NA),At = c(At, NA), Bt = c(Bt,NA),
ramalan = ramal )
data.gab2
error = data.gab2[,1]-data.gab2[,6]
sse<-sum(error^2,na.rm=T)
mse <- mean(error^2,na.rm=T)
rmse<-sqrt(mse)
mad<-mean(abs(error),na.rm=T)
rerror<-error/data.gab2[,1]*100
mape<-mean(abs(rerror), na.rm = T)
ak<-data.frame(n,SSE=sse,MSE=mse,RMSE=rmse,MAD=mad,MAPE=mape)
d<-rbind(d,ak)
if(just.output==F){
ts.plot(training.ts, lty=3, col="black", main=paste("n =",n),ylab=axis.y.label)
points(training.ts);lines(sma, col="red");lines(ramal, col="green", lty=1)
legend("topleft",c("Aktual","Pemulusan","Ramal"),
lty =c(3,1,1) ,col=c("black","red","green"))
}
}
opt<-d$n[which.min(d[,metode])];opt
if(just.graph==F){
return(list(DMA=dma,Akurasi=d,n.optimum=paste("n optimum adalah",opt,"dengan",metode,"sebesar",round(d[which(d$n==opt),metode],3))))
}
}
# Parameter n Optimum
DMA<-nopdma(training.ts,2,15,"MAPE","Rataan Suhu",just.output = T);DMA## $DMA
## Time Series:
## Start = 1
## End = 51
## Frequency = 1
## [1] NA NA NA NA NA NA NA NA
## [9] NA NA NA NA NA NA NA NA
## [17] NA NA NA NA NA NA NA NA
## [25] NA NA NA NA 25.47098 25.49435 25.52844 25.54611
## [33] 25.57747 25.62608 25.57457 25.56704 25.60406 25.63715 25.63197 25.65390
## [41] 25.70190 25.76370 25.80372 25.82415 25.83310 25.84797 25.85147 25.88037
## [49] 25.91657 25.96028 25.97956
##
## $Akurasi
## n SSE MSE RMSE MAD MAPE
## 1 2 2.4931938 0.05194154 0.2279069 0.1668229 0.6520165
## 2 3 1.9074025 0.04146527 0.2036302 0.1670531 0.6528499
## 3 4 1.5607548 0.03547170 0.1883393 0.1475616 0.5765849
## 4 5 1.2975208 0.03089335 0.1757650 0.1356810 0.5300116
## 5 6 1.1945537 0.02986384 0.1728116 0.1400681 0.5468422
## 6 7 1.0639007 0.02799739 0.1673242 0.1387290 0.5415854
## 7 8 0.9383179 0.02606439 0.1614447 0.1365048 0.5327141
## 8 9 0.8438414 0.02481887 0.1575400 0.1271550 0.4958286
## 9 10 0.7517040 0.02349075 0.1532669 0.1232559 0.4803702
## 10 11 0.6708112 0.02236037 0.1495339 0.1196926 0.4664474
## 11 12 0.6716029 0.02398582 0.1548736 0.1185811 0.4617090
## 12 13 0.5619866 0.02161487 0.1470200 0.1137013 0.4423279
## 13 14 0.5090805 0.02121169 0.1456423 0.1107897 0.4300857
## 14 15 0.4859327 0.02208785 0.1486198 0.1177945 0.4572236
##
## $n.optimum
## [1] "n optimum adalah 14 dengan MAPE sebesar 0.43"Dengan melakukan fungsi pengulangan untuk metode DMA diperoleh n optimum yaitu 14 dengan ukuran performa model MAPE sebesar 0.43
Evaluasi Performa DMA Menggunakan Data Testing
df.ts.dma <- TTR::SMA(testing.ts, n=14)
ramal.df.dma <- c(NA,df.ts.dma)
df.dma <- cbind(aktual.smd=c(testing.ts,NA), pemulusan.smd=c(testing.ts,NA), ramal.df.dma)
error.dma <- df.dma[,1]-df.dma[,3]
rerror.dma <- (error.dma/df.dma[,1])*100
mape.dma <- mean(abs(rerror.dma), na.rm=T);mape.dma## [1] 0.5640028Plot DMA
nopdma(training.ts,2,5,"MAPE","Rataan Suhu",just.graph = T)
nopdma(training.ts,6,10,"MAPE","Rataan Suhu",just.graph=T)
nopdma(training.ts,11,15,"MAPE","Rataan Suhu",just.graph=T)
3. Single Exponential Smoothing
Parameter optimum SES
out1<-data.frame()
for(i in seq(0.01,0.99,by=0.01)){
Ho1<-HoltWinters(training.ts,alpha=i,beta=F,gamma=F)
error1<-c(training.ts)-c(Ho1$fitted[,1])
sse1<-Ho1$SSE
mse1<-sse1/length(training.ts)
mape1<-mean(abs(error1/training.ts*100))
rmse1<-sqrt(mse1)
akurasi1<-cbind("Alpha"=i,"SSE"=sse1,
"MSE"=mse1,"RMSE"=rmse1,"MAPE"=mape1)
out1<-rbind(out1,akurasi1)
}
out1[which.min(out1$MAPE),]## Alpha SSE MSE RMSE MAPE
## 99 0.99 1.402324 0.02749655 0.1658208 0.0314612aopt1<-out1[which.min(out1$MAPE),1]Dengan melakukan fungsi pengulangan untuk metode SES diperoleh parameter optimum yaitu alpha 0.99 dengan ukuran performa model MAPE sebesar 0.03146
Plot SES Dengan Parameter Optimum
Ho1<-HoltWinters(training.ts,alpha=aopt1,gamma=F)
ses<-data.frame(training.ts,c(NA,NA,Ho1$fitted[,1]));colnames(ses)<-c("y","yhat")
plot(1950:2000,training.ts,xlab="Tahun",ylab="Suhu Udara (°C)");lines(1950:2000,training.ts,lty=3,lwd=3,col="green");points(training.ts)
lines(1950:2000,y=ses$yhat,col="#FFB200",lwd=3);
legend("topleft",c("Aktual","Fitted SES"),lty=c(3,1),col=c("green","#FFB200"))
Hasil ramal SES 20 tahun kedepan
ramalansesopt<- forecast::forecast(Ho1, h=20);ramalansesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 52 25.86940 25.65055 26.08825 25.53470 26.20410
## 53 25.88936 25.58139 26.19732 25.41837 26.36035
## 54 25.90931 25.53276 26.28587 25.33342 26.48521
## 55 25.92927 25.49481 26.36373 25.26482 26.59372
## 56 25.94923 25.46372 26.43474 25.20670 26.69175
## 57 25.96918 25.43749 26.50088 25.15603 26.78234
## 58 25.98914 25.41497 26.56331 25.11102 26.86726
## 59 26.00910 25.39537 26.62282 25.07049 26.94771
## 60 26.02906 25.37817 26.67994 25.03362 27.02450
## 61 26.04901 25.36298 26.73505 24.99981 27.09822
## 62 26.06897 25.34949 26.78845 24.96862 27.16932
## 63 26.08893 25.33749 26.84037 24.93970 27.23815
## 64 26.10888 25.32679 26.89098 24.91277 27.30500
## 65 26.12884 25.31724 26.94044 24.88760 27.37008
## 66 26.14880 25.30872 26.98887 24.86401 27.43358
## 67 26.16876 25.30114 27.03637 24.84185 27.49566
## 68 26.18871 25.29440 27.08303 24.82097 27.55645
## 69 26.20867 25.28843 27.12891 24.80128 27.61606
## 70 26.22863 25.28316 27.17409 24.78266 27.67459
## 71 26.24858 25.27855 27.21862 24.76505 27.73212plot(ramalansesopt,xlab="periode waktu", ylab="Yt", main="Peramalan Dengan Metode SES")
lines(datasuhu.ts,lty=1,col="red")
points(datasuhu.ts,cex=0.5)
Terilhat bahwa plot hasil nilai peramalan bergerak naik secara positif.
Evaluasi Performa SES Menggunakan Data Testing
ses1 <- HoltWinters(testing.ts, gamma=FALSE, beta=FALSE, alpha=0.99); ses1## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = testing.ts, alpha = 0.99, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.99
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 25.9919error1<-c(testing.ts)-c(ses1$fitted[,1])## Warning in c(testing.ts) - c(ses1$fitted[, 1]): longer object length is not a
## multiple of shorter object lengthmape1<-mean(abs(error1/testing.ts*100));mape1## [1] 0.018103814. Double Exponential Smoothing
parameter Optimum DES
out<-data.frame()
for(i in seq(0.01,0.99,by=0.01)){
for(j in seq(0.01,0.99,by=0.01)){
Ho<-HoltWinters(training.ts,alpha=i,beta=j,gamma=F)
error<-c(training.ts)-c(Ho$fitted[,1])
sse<-Ho$SSE
mse<-sse/length(training.ts)
mape<-mean(abs(error/training.ts*100))
rmse<-sqrt(mse)
akurasi<-cbind("Alpha"=i,"Beta"=j,"SSE"=sse,
"MSE"=mse,"RMSE"=rmse,"MAPE"=mape)
out<-rbind(out,akurasi)
}
}
out[which.min(out$MAPE),]## Alpha Beta SSE MSE RMSE MAPE
## 3508 0.36 0.43 1.324848 0.0259774 0.1611751 0.3334745aopt<-out[which.min(out$MAPE),1]
bopt<-out[which.min(out$MAPE),2]Dengan melakukan fungsi pengulangan untuk metode DES diperoleh parameter optimum yaitu alpha 0.36 dan beta 0.43 dengan ukuran performa model MAPE sebesar 0.3335
Plot DES dengan parameter optimum
Ho<-HoltWinters(training.ts,alpha=aopt,beta=bopt,gamma=F)
des<-data.frame(training.ts,c(NA,NA,Ho$fitted[,1]));colnames(des)<-c("y","yhat")
plot(1950:2000,training.ts,xlab="Tahun",ylab="Suhu Udara (°C)");lines(1950:2000,training.ts,col="green",lty=3,lwd=3);points(training.ts)
lines(1950:2000,des$yhat,col="#FFB200",lwd=3);
legend("topleft",c("Aktual","Fitted Des"),lty=c(3,1),col=c("green","#FFB200"))
Hasil ramal DES 20 taun kedepan
ramalandesopt<- forecast::forecast(Ho, h=20); ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 52 25.93552 25.72263 26.14841 25.60993 26.26111
## 53 25.93934 25.69990 26.17879 25.57314 26.30555
## 54 25.94317 25.66450 26.22184 25.51698 26.36935
## 55 25.94699 25.61766 26.27632 25.44332 26.45066
## 56 25.95081 25.56105 26.34058 25.35472 26.54691
## 57 25.95464 25.49616 26.41311 25.25345 26.65582
## 58 25.95846 25.42415 26.49277 25.14131 26.77561
## 59 25.96228 25.34589 26.57867 25.01959 26.90497
## 60 25.96610 25.26202 26.67019 24.88930 27.04291
## 61 25.96993 25.17303 26.76683 24.75118 27.18868
## 62 25.97375 25.07930 26.86820 24.60580 27.34170
## 63 25.97757 24.98113 26.97402 24.45364 27.50151
## 64 25.98140 24.87877 27.08403 24.29507 27.66772
## 65 25.98522 24.77242 27.19802 24.13040 27.84004
## 66 25.98904 24.66226 27.31583 23.95990 28.01818
## 67 25.99287 24.54844 27.43729 23.78380 28.20193
## 68 25.99669 24.43109 27.56229 23.60231 28.39107
## 69 26.00051 24.31032 27.69070 23.41559 28.58543
## 70 26.00434 24.18625 27.82242 23.22381 28.78486
## 71 26.00816 24.05896 27.95736 23.02712 28.98920plot(ramalandesopt)
plot(ramalandesopt,xlab="periode waktu", ylab="Yt", main="Peramalan Dengan Metode DES")
lines(datasuhu.ts,lty=1,col="red")
points(datasuhu.ts,cex=0.5)
Evaluasi Performa DES Menggunakan Data Testing
Ho<-HoltWinters(testing.ts,alpha=aopt,beta=bopt,gamma=F)
error<-c(testing.ts)-c(Ho$fitted[,1])## Warning in c(testing.ts) - c(Ho$fitted[, 1]): longer object length is not a
## multiple of shorter object lengthmape<-mean(abs(error/testing.ts*100));mape## [1] 0.38008875. Metode Terbaik
PEMULUSAN <- c("SMA", "DMA", "SES", "DES")
MAPE <- round(c(mape.sma, mape.dma, mape1, mape),2)
data.frame(PEMULUSAN,MAPE)## PEMULUSAN MAPE
## 1 SMA 0.46
## 2 DMA 0.56
## 3 SES 0.02
## 4 DES 0.38Berdasarkan ukuran performa menggunakan MAPE diperoleh MAPE terkecil yaitu metode SES dengan nilai 0.02. Jika dilihat dari plot peramalan 20 tahun kedepan dan dibandingkan dengan data testing maka plot yang lebih tepat menggambarkan data testing adalah metode SES. Jadi metode terbaik yang digunakan adalah metode Single Exponential Smoothing.
G. Pemodelan Regresi Time Series
1. Model Koyck
trainpop<- Data[1:51,-1]
testpop<- Data[-c(1:51),-1]
colnames(trainpop)[1]<-colnames(testpop)[1]<-"Suhu"
#Model regresi sederhana
summary(lm(Suhu~Populasi,trainpop))##
## Call:
## lm(formula = Suhu ~ Populasi, data = trainpop)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35429 -0.09536 0.01146 0.09056 0.36120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.516e+01 6.119e-02 411.18 < 2e-16 ***
## Populasi 3.174e-03 4.301e-04 7.38 1.71e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1374 on 49 degrees of freedom
## Multiple R-squared: 0.5264, Adjusted R-squared: 0.5167
## F-statistic: 54.47 on 1 and 49 DF, p-value: 1.714e-09# model koyck
k<-koyckDlm(trainpop$Populasi,trainpop$Suhu);summary(k)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.387798 -0.069832 0.003516 0.079177 0.328643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.938e+01 3.525e+00 5.497 1.54e-06 ***
## Y.1 2.291e-01 1.401e-01 1.635 0.108652
## X.t 2.550e-03 6.153e-04 4.144 0.000141 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1352 on 47 degrees of freedom
## Multiple R-Squared: 0.5581, Adjusted R-squared: 0.5393
## Wald test: 29.79 on 2 and 47 DF, p-value: 4.401e-09
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: 25.14053 0.002549488 0.2291461e<-forecast(model = k, x=testpop$Populasi, h=21);e## $forecasts
## [1] 25.85663 25.86580 25.87546 25.88496 25.89445 25.90425 25.91430 25.92445
## [9] 25.93454 25.94459 25.95475 25.96504 25.97519 25.98504 25.99460 26.00382
## [17] 26.01268 26.02126 26.02964 26.03736 26.04396
##
## $call
## forecast.koyckDlm(model = k, x = testpop$Populasi, h = 21)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"#mape data testing
mapektest <- MLmetrics::MAPE(e$forecasts, testpop$Suhu);mapektest## [1] 0.003405158#akurasi data training
mapektrain <- GoF(k)[1,"MAPE"];mapektrain## [1] 0.0038302192. Distributed Lag Model with lag optim
qopt<-finiteDLMauto(formula=Suhu~Populasi,data=trainpop,model.type = "dlm",error.type = "AIC")
dlmod<-dlm(x=trainpop$Populasi,y=trainpop$Suhu,q=qopt$`q - k`);summary(dlmod)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.297902 -0.081588 -0.006217 0.072111 0.298897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.36253 0.09228 274.831 < 2e-16 ***
## x.t -0.14181 0.04777 -2.969 0.00469 **
## x.1 0.14659 0.04826 3.037 0.00389 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1272 on 47 degrees of freedom
## Multiple R-squared: 0.6091, Adjusted R-squared: 0.5925
## F-statistic: 36.62 on 2 and 47 DF, p-value: 2.586e-10
##
## AIC and BIC values for the model:
## AIC BIC
## 1 -59.43765 -51.78956f<-forecast(dlmod,x=testpop$Populasi,h=21);f## $forecasts
## [1] 25.95275 25.97255 25.99221 26.02147 26.03396 26.02981 26.03436 26.04649
## [9] 26.06593 26.08186 26.08944 26.09855 26.12335 26.15201 26.17910 26.20750
## [17] 26.23619 26.26027 26.27989 26.32607 26.39079
##
## $call
## forecast.dlm(model = dlmod, x = testpop$Populasi, h = 21)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"#akurasi testing
mapedtest<-MLmetrics::MAPE(f$forecasts,testpop$Suhu);mapedtest## [1] 0.00527501#akurasi training
mapedtrain<-GoF(dlmod)[1,"MAPE"];mapedtrain## [1] 0.0037577013. Autoregressive Model With Lag Optimum
pqopt<-ardlBoundOrders(trainpop,Suhu~Populasi);pqopt## $p
## Populasi
## 1 15
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7
## p = 1 -53.52151 -50.22716 -46.89920 -44.00399 -40.75964 -40.10707 -35.95724
## p = 2 -49.45381 -48.54938 -45.25847 -42.65339 -39.42943 -39.29839 -35.16105
## p = 3 -49.96497 -49.96497 -48.62345 -46.12269 -43.74500 -43.12085 -39.14161
## p = 4 -44.02036 -44.20249 -44.20249 -44.15993 -41.83745 -41.58383 -37.62724
## p = 5 -40.39291 -41.45981 -39.96107 -39.96107 -39.88554 -39.65785 -35.87085
## p = 6 -40.34474 -43.31440 -43.09734 -43.37309 -43.37309 -43.67724 -40.23342
## p = 7 -35.23841 -36.53590 -36.70804 -38.70290 -42.35565 -42.35565 -42.71604
## p = 8 -42.99408 -45.56223 -44.21684 -42.25448 -43.41034 -43.98601 -43.98601
## p = 9 -37.23821 -39.35156 -38.59670 -37.08767 -39.40860 -41.13501 -39.61728
## p = 10 -38.78071 -40.43788 -38.55031 -37.81094 -39.72019 -38.42981 -36.44399
## p = 11 -30.74271 -32.65509 -31.11916 -29.83299 -32.56621 -33.47502 -31.48360
## p = 12 -31.00565 -31.52797 -29.57507 -28.36835 -33.74612 -34.96087 -32.96703
## p = 13 -22.11990 -23.32810 -21.79189 -21.20033 -24.21286 -24.27101 -22.54580
## p = 14 -18.57945 -21.47701 -19.52345 -18.37144 -23.40155 -24.12385 -22.24909
## p = 15 -18.67010 -18.97921 -16.99125 -15.56849 -20.88249 -23.54451 -21.54452
## q = 8 q = 9 q = 10 q = 11 q = 12 q = 13 q = 14
## p = 1 -34.16125 -30.75218 -28.42318 -23.93919 -19.67190 -15.42558 -20.22429
## p = 2 -32.94613 -29.00109 -26.44889 -21.96983 -17.71156 -13.47936 -18.22648
## p = 3 -37.22657 -33.17977 -30.52573 -27.67783 -23.18446 -18.88321 -27.11032
## p = 4 -35.63473 -31.44691 -28.66131 -27.06733 -22.87376 -18.43297 -27.14983
## p = 5 -33.83642 -29.61359 -26.96791 -25.20460 -21.43762 -17.39266 -25.45420
## p = 6 -39.50361 -34.99882 -32.50505 -33.35687 -29.42482 -35.43739 -34.36394
## p = 7 -43.09057 -38.13485 -35.39494 -38.21917 -36.23872 -42.83377 -41.48466
## p = 8 -44.56932 -39.60586 -36.30339 -39.25261 -36.12895 -41.23228 -39.96213
## p = 9 -39.61728 -37.67121 -34.52120 -37.87496 -34.51356 -42.08182 -41.55638
## p = 10 -35.96671 -35.96671 -34.24255 -36.42528 -32.88486 -41.92583 -40.29351
## p = 11 -32.17133 -31.17799 -31.17799 -35.99314 -32.26115 -42.15757 -43.75174
## p = 12 -37.20065 -35.22232 -33.65187 -33.65187 -32.61349 -44.67428 -44.85044
## p = 13 -22.98668 -22.62984 -20.67540 -27.51501 -27.51501 -43.73991 -43.88063
## p = 14 -28.46812 -29.32093 -27.34389 -43.44771 -42.73301 -42.73301 -46.38568
## p = 15 -22.12406 -23.69245 -22.86138 -31.08742 -29.18499 -58.31335 -58.31335
## q = 15
## p = 1 -18.78773
## p = 2 -17.86550
## p = 3 -24.20120
## p = 4 -23.94476
## p = 5 -21.96407
## p = 6 -33.28841
## p = 7 -48.21143
## p = 8 -51.34203
## p = 9 -64.78299
## p = 10 -63.92515
## p = 11 -62.38612
## p = 12 -60.64745
## p = 13 -64.31767
## p = 14 -106.94753
## p = 15 -216.84814
##
## $min.Stat
## [1] -216.8481ARD<-ardlDlm(Suhu~Populasi,trainpop,p=pqopt$p$Populasi,q=pqopt$q);summary(ARD)##
## Time series regression with "ts" data:
## Start = 16, End = 51
##
## Call:
## dynlm(formula = as.formula(model.text), data = data)
##
## Residuals:
## 16 17 18 19 20 21 22
## -8.031e-03 -3.199e-02 -1.975e-02 -1.095e-02 2.971e-02 -2.439e-03 1.828e-02
## 23 24 25 26 27 28 29
## -5.350e-04 1.999e-02 1.600e-02 3.918e-02 1.445e-02 -2.946e-04 -1.079e-02
## 30 31 32 33 34 35 36
## -1.898e-05 -6.163e-02 7.964e-02 -1.158e-02 3.173e-02 -8.312e-02 -1.308e-01
## 37 38 39 40 41 42 43
## 6.465e-02 7.098e-02 -1.320e-02 2.161e-03 -3.631e-03 -4.909e-02 8.575e-02
## 44 45 46 47 48 49 50
## 2.060e-02 -1.307e-02 1.665e-02 -1.276e-01 -2.150e-02 8.980e-02 -6.070e-02
## 51
## 6.112e-02
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 301.0499 140.8507 2.137 0.0994 .
## Populasi.t -0.0439 0.5568 -0.079 0.9409
## Populasi.1 -0.8323 1.4458 -0.576 0.5957
## Populasi.2 1.7186 1.9995 0.860 0.4385
## Populasi.3 -0.7827 2.4597 -0.318 0.7662
## Populasi.4 -1.2129 2.6471 -0.458 0.6706
## Populasi.5 2.1619 2.6478 0.816 0.4601
## Populasi.6 -1.4026 2.7418 -0.512 0.6359
## Populasi.7 -1.4001 2.9439 -0.476 0.6592
## Populasi.8 2.6812 2.8274 0.948 0.3967
## Populasi.9 -0.7938 2.6593 -0.299 0.7802
## Populasi.10 -1.4216 2.6683 -0.533 0.6224
## Populasi.11 2.6175 2.6164 1.000 0.3737
## Populasi.12 -1.3771 2.7784 -0.496 0.6461
## Populasi.13 -0.3859 2.8145 -0.137 0.8976
## Populasi.14 1.2817 2.2922 0.559 0.6059
## Populasi.15 -0.7612 1.0425 -0.730 0.5057
## Suhu.1 -0.5009 0.4535 -1.105 0.3313
## Suhu.2 -0.5256 0.3953 -1.330 0.2543
## Suhu.3 -1.4581 0.6534 -2.232 0.0895 .
## Suhu.4 -0.7833 0.6704 -1.168 0.3075
## Suhu.5 -0.7045 0.4456 -1.581 0.1890
## Suhu.6 -1.4266 0.5969 -2.390 0.0752 .
## Suhu.7 -0.8251 0.6007 -1.373 0.2416
## Suhu.8 -0.5314 0.5156 -1.031 0.3609
## Suhu.9 -1.1054 0.5602 -1.973 0.1197
## Suhu.10 -0.3851 0.5595 -0.688 0.5291
## Suhu.11 -0.3392 0.4296 -0.790 0.4739
## Suhu.12 -0.8023 0.4090 -1.961 0.1214
## Suhu.13 -0.3131 0.4458 -0.702 0.5211
## Suhu.14 -0.6863 0.4087 -1.679 0.1684
## Suhu.15 -0.4451 0.5398 -0.825 0.4559
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1524 on 4 degrees of freedom
## Multiple R-squared: 0.9381, Adjusted R-squared: 0.4584
## F-statistic: 1.955 on 31 and 4 DF, p-value: 0.2722#forecasting
g<-forecast(ARD,x=testpop$Populasi,h=21);g## $forecasts
## [1] 25.36576 25.87952 26.06048 26.14031 26.11562 25.94481 26.00797 26.19236
## [9] 25.92240 25.99651 25.87186 25.94542 26.14045 26.16460 26.36106 25.83063
## [17] 25.81979 25.87449 26.24704 26.63630 26.74700
##
## $call
## forecast.ardlDlm(model = ARD, x = testpop$Populasi, h = 21)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"#akurasi testing
mapeatest<-MLmetrics::MAPE(g$forecasts,testpop$Suhu);mapeatest## [1] 0.008669551#akurasi training
mapeatrain<-GoF(ARD)[1,"MAPE"];mapeatrain## [1] 0.0014283624. Plot Perbandingan
plot(testpop$Populasi,testpop$Suhu,col="black",type="b",ylim = c(25,27),lty=3)
points(testpop$Populasi,f$forecasts,col="coral")
lines(testpop$Populasi,f$forecasts,col="coral",lty=1)
points(testpop$Populasi,g$forecasts,col="blue")
lines(testpop$Populasi,g$forecasts,col="blue",lty=1)
points(testpop$Populasi,e$forecasts,col="green")
lines(testpop$Populasi,e$forecasts,col="green",lty=1)
legend("topleft",c("Aktual","Distributed Lag","Autoregressive Lag","Koyck"),
lty =c(3,1,1) ,col=c("black","coral","blue","green"))
5. Perbandingan Akurasi
#perbandingan akurasi training
data.frame(Metode=c("Autoregressive","Distributed","Koyck"),MAPE=c(mapeatrain,mapedtrain,mapektrain))## Metode MAPE
## 1 Autoregressive 0.001428362
## 2 Distributed 0.003757701
## 3 Koyck 0.003830219#perbandingan akurasi testing
data.frame(Metode=c("Autoregressive","Distributed","Koyck"),MAPE=c(mapeatest,mapedtest,mapektest))## Metode MAPE
## 1 Autoregressive 0.008669551
## 2 Distributed 0.005275010
## 3 Koyck 0.003405158Berdasarkan ukuran performa menggunakan MAPE maka model terbaik adalah model Koyck
6. Pengecekan Autokorelasi
bgtest(k$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: k$model
## LM test = 0.14064, df = 1, p-value = 0.7076bgtest(ARD$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: ARD$model
## LM test = 14.87, df = 1, p-value = 0.0001152bgtest(dlmod$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: dlmod$model
## LM test = 0.23425, df = 1, p-value = 0.6284Autokol tertangani menggunakan model Distributed Lag dan Koyck Model
7. Kesimpulan
SIMPULAN: Metode pemulusan terbaik yang dipakai adalah SES dengan nilai mape 0.02 PEMODELAN TERBAIK MENGGUNAKAN MODEL KOYCK
H. Pemodelan ARIMA
1. Cek kestasioneran
Sebelum melakukan pemodelan, data yang digunakan harus sudah stasioner dari nilai tengah atau ragam. kestasioneran data dapat dilihat secara eksploratif, yaitu menggunakan plot time series data atau plot ACF. Untuk memastikan kestasioneran data pada nilai tengah, dilakukan uji Dicky-Fuller. Apabila p-value yang dihasilkan lebih besar dari alpha (0.05) maka tak tolak H0 yang berarti bahwa data tersebut tidak stasioner.
ts.plot(train$`Rata-rata suhu`)
Dari plot time series terhadap data awal terlihat bahwa data memiliki pola trend positif, dimana data tersebut tidak stasioner baik dalam ragam ataupun rataan.
acf(train$`Rata-rata suhu`)
Dilihat dari plot ACF juga tail off (slowly) yang artinya data tersebut tidak stasioner.
adf.test(train$`Rata-rata suhu`)#ga stasioner##
## Augmented Dickey-Fuller Test
##
## data: train$`Rata-rata suhu`
## Dickey-Fuller = -2.9646, Lag order = 3, p-value = 0.1866
## alternative hypothesis: stationary- Kemudian berdasarkan uji diagnostik menggunakan adf test diperoleh p-value > 0.05 maka tak tolak H0 artinya data tersebut tidak stasioner.
- Karena data tidak stasioner maka kita akan melakukan Differencing untuk men-stasionerkan data.
2. Differencing
d1<-diff(train$`Rata-rata suhu`,differences = 1)
adf.test(d1)##
## Augmented Dickey-Fuller Test
##
## data: d1
## Dickey-Fuller = -5.2743, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryts.plot(d1)
- Setelah melihat pola data yang cenderung trend dan tidak stasioner, maka dilakukan penanganannya dengan melakukan differencing agar pola data menjadi stasioner. differencing dilakukan sekali menggunakan fungsi diff.
- Setelah melakukan differencing sekali, dapat dilihat bahwa pola data sudah stasioner. Berdasarkan uji diagnostik menggunakan adf test diperoleh p-value < 0.05 maka tolak H0 artinya data tersebut sudah stasioner.
3. Kandidat Model
Pemodelan Arima dapat dilakukan dengan mengetahui nilai p dan q menggunakan plot ACF, PACF, atau ARMA. Setelah dilakukan differencing sekali, data kemudian digambarkan plot ACF, PACF, dan tabel ARMA. Berikut adalat plot ACF, PACF, dan tabel ARMA pada data yang telah di differencing sekali.
par(mfrow=c(1,2))
acf(d1)#ARIMA(0,1,1)
pacf(d1)#ARIMA(2,1,0)
eacf(d1)#ARIMA(1,1,2), ARIMA(1,1,1), ARIMA(0,1,2)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 x o o o o o o o o o o o o o
## 3 x x x o o o o o o o o o o o
## 4 x o o x o o o o o o o o o o
## 5 x x o x o o o o o o o o o o
## 6 x o o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o o- Model tentatif yang diperoleh berdasarkan plot ACF adalah ARIMA(0,1,1) dan berdasarkan plot PACF adalah ARIMA(3,1,0).
- Diperoleh tiga model tentatif berdasarkan pola segitiga-nol yang terbentuk pada tabel EACF yaitu ARIMA(1,1,2), ARIMA(1,1,1), dan ARIMA(0,1,2).
4. Evaluasi Performa ARIMA menggunakan Data Testing
Setelah didapatkan model model terpilih dengan menggunakan plot ACF, PACF, dan tabel ARMA, setiap model dicari nilai AIC nya
model1<-arima(test$`Rata-rata suhu`,c(0,1,1));aic1<-model1$aic
model2<-arima(test$`Rata-rata suhu`,c(2,1,0));aic2<-model2$aic
model3<-arima(test$`Rata-rata suhu`,c(1,1,2));aic3<-model3$aic
model4<-arima(test$`Rata-rata suhu`,c(1,1,1));aic4<-model4$aic
model5<-arima(test$`Rata-rata suhu`,c(0,1,2));aic5<-model5$aic
d<-data.frame(model=paste("model",1:5),ar=c(0,2,1,1,0),i=rep(1,5),
ma=c(1,0,2,1,2),aic=round(c(aic1,aic2,aic3,aic4,aic5),2));d## model ar i ma aic
## 1 model 1 0 1 1 -28.23
## 2 model 2 2 1 0 -29.54
## 3 model 3 1 1 2 -27.26
## 4 model 4 1 1 1 -26.33
## 5 model 5 0 1 2 -27.02d[which.min(d$aic),]## model ar i ma aic
## 2 model 2 2 1 0 -29.54Model-model yang akan dibandingkan nilai AIC-nya adalah model ARIMA(0,1,1), ARIMA(2,1,0), ARIMA(1,1,2), ARIMA(1,1,1), dan ARIMA(0,1,2). Model yang dipilih adalah model dengan nilai AIC terendah. Berdasarkan model-model tersebut, model dengan nilai AIC terendah adalah model ARIMA(2,1,0) dengan nilai AIC sebesar -29.54.
5. Kesimpulan
Kandidat model terbaik berdasarkan AIC adalah ARIMA(2,1,0). Namun, ketika menggunakan auto.arima di R, model yang didapatkan ARIMA(0,1,0). Maka, yang masuk akal, model terbaiknya adalah model hasil perhitungan mandiri yaitu ARIMA(2,1,0)
auto.arima(test$`Rata-rata suhu`)## Series: test$`Rata-rata suhu`
## ARIMA(0,1,0)
##
## sigma^2 = 0.01632: log likelihood = 12.79
## AIC=-23.59 AICc=-23.36 BIC=-22.59I. Perbandingan 2 Model Terbaik
plot(CEK$Tahun,CEK$`Rata-rata suhu`,xlab="Tahun",ylab="Rata-rata Suhu (°C)",main="Perbandingan Hasil Prediksi Data Testing")
f<-c(predict(model2,21)$pred);f## [1] 26.08537 26.14763 26.06748 26.06879 26.11233 26.09228 26.07721 26.09494
## [9] 26.09536 26.08542 26.08960 26.09322 26.08931 26.08905 26.09132 26.09045
## [17] 26.08959 26.09045 26.09054 26.09003 26.09021lines(CEK$Tahun,CEK$`Rata-rata suhu`,col="blue",lwd=2)
lines(CEK$Tahun,c(rep(NA,nrow(train)),f),col="red",lwd=2)
lines(CEK$Tahun,c(rep(NA,nrow(train)),e$forecasts),col="green",lwd=2)
legend("topleft",c("Aktual","ARIMA (2, 1, 0)","Koyck"),col = c("blue","red","green"),lwd=2)
dfm<-data.frame(MAPE=c(mapektest,mean(abs(model2$residuals/training.ts*100))));rownames(dfm)<-c("Koyck","Arima (2,1,0)");dfm## MAPE
## Koyck 0.003405158
## Arima (2,1,0) 0.323423091J. Kesimpulan
Berdasarkan grafik dengan garis prediksi paling mendekati data aktual adalah Model Koyck. Artinya, pemodelan terbaik untuk memprediksi data rataan suhu di Indonesia yaitu ketika tidak hanya mempertimbangkan pengaruh rataan suhu di masa kini dan periode sebelumnya, tetapi juga jumlah populasi di masa kini dan periode sebelumnya, mengingat jumlah populasi memiliki korelasi hubungan positif kuat terhadap rataan suhu di Indonesia selama 72 tahun terakhir