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

Data Peubah Respon

CEK <- read_xlsx("D:/FILE KULIAH SEMESTER 5/MINGGU 4/MPDW (P)/datasuhu.xlsx")
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.4

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("D:/FILE KULIAH SEMESTER 5/MINGGU 4/MPDW (P)/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.7

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

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.

Time Series Plot Data Peubah Respon

g1<-ggplot(CEK, aes(x=Tahun, y=`Rata-rata suhu`)) +
  geom_line(lty=1, lwd=0.8, col="black") + 
  labs(x="Tahun",y = "Rata-rata Suhu (°C)", 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)

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 = "Rata-rata Suhu (°C)", 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)

Gambar di atas menunjukkan time series plot dari data suhu rata-rata Indonesia pada tahun 1950-2000 yang merupakan hasil splitting data awal atau disebut sebagai data training dengan proporsi sebesar 70% dari data awal. Data ini digunakan untuk mencari model yang sesuai melalui serangkaian algoritma tertentu atau digunakan juga untuk tahapan modelling.

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 = "Rata-rata Suhu (°C)", 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)

Gambar di atas menunjukkan time series plot dari data suhu rata-rata Indonesia pada tahun 2001-2021 yang merupakan hasil splitting data awal atau disebut sebagai data testing dengan proporsi sebesar 30% dari data awal. Data ini digunakan untuk mengetahui akurasi performa model yang diperoleh dari data training atau digunakan juga untuk tahapan forecasting.

Scatter Plot Peubah Penjelas dan Peubah Respon

plot(Data$Populasi, Data$`Rata-rata suhu`, col="black", lwd=2, type="p",xlab="Populasi",ylab="Rata-rata Suhu (°C)",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.8704546

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

1. Single Moving Average (SMA)

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

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); df.sma

##       aktual.sma pemulusan.sma ramal.df.sma
##  [1,]      25.91         25.91           NA
##  [2,]      26.04         26.04           NA
##  [3,]      25.90         25.90           NA
##  [4,]      25.92         25.92           NA
##  [5,]      25.93         25.93           NA
##  [6,]      25.82         25.82           NA
##  [7,]      25.84         25.84           NA
##  [8,]      25.71         25.71           NA
##  [9,]      25.98         25.98           NA
## [10,]      26.08         26.08     25.89444
## [11,]      25.88         25.88     25.91333
## [12,]      25.99         25.99     25.89556
## [13,]      26.05         26.05     25.90556
## [14,]      26.04         26.04     25.92000
## [15,]      26.12         26.12     25.93222
## [16,]      26.23         26.23     25.96556
## [17,]      26.03         26.03     26.00889
## [18,]      26.04         26.04     26.04444
## [19,]      26.20         26.20     26.05111
## [20,]      26.18         26.18     26.06444
## [21,]      25.99         25.99     26.09778
## [22,]         NA            NA     26.09778

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

2. Double Moving Average (DMA)

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"

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); df.dma

##       aktual.smd pemulusan.smd ramal.df.dma
##  [1,]      25.91         25.91           NA
##  [2,]      26.04         26.04           NA
##  [3,]      25.90         25.90           NA
##  [4,]      25.92         25.92           NA
##  [5,]      25.93         25.93           NA
##  [6,]      25.82         25.82           NA
##  [7,]      25.84         25.84           NA
##  [8,]      25.71         25.71           NA
##  [9,]      25.98         25.98           NA
## [10,]      26.08         26.08           NA
## [11,]      25.88         25.88           NA
## [12,]      25.99         25.99           NA
## [13,]      26.05         26.05           NA
## [14,]      26.04         26.04           NA
## [15,]      26.12         26.12     25.93500
## [16,]      26.23         26.23     25.95000
## [17,]      26.03         26.03     25.96357
## [18,]      26.04         26.04     25.97286
## [19,]      26.20         26.20     25.98143
## [20,]      26.18         26.18     26.00071
## [21,]      25.99         25.99     26.02643
## [22,]         NA            NA     26.03714

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

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

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

aopt1<-out1[which.min(out1$MAPE),1]

Ho1<-HoltWinters(training.ts,alpha=aopt1,beta=F,gamma=F)
ses<-data.frame(training.ts,c(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"))

ramalansesopt<- forecast::forecast(Ho1, h=21);ramalansesopt

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 52       25.84924 25.63263 26.06585 25.51797 26.18052
## 53       25.84924 25.54444 26.15404 25.38308 26.31540
## 54       25.84924 25.47656 26.22192 25.27927 26.41921
## 55       25.84924 25.41927 26.27921 25.19165 26.50683
## 56       25.84924 25.36876 26.32972 25.11441 26.58407
## 57       25.84924 25.32308 26.37540 25.04454 26.65394
## 58       25.84924 25.28106 26.41743 24.98028 26.71820
## 59       25.84924 25.24194 26.45655 24.92045 26.77803
## 60       25.84924 25.20519 26.49329 24.86425 26.83424
## 61       25.84924 25.17042 26.52806 24.81108 26.88740
## 62       25.84924 25.13736 26.56112 24.76051 26.93797
## 63       25.84924 25.10576 26.59272 24.71219 26.98630
## 64       25.84924 25.07545 26.62303 24.66583 27.03265
## 65       25.84924 25.04629 26.65219 24.62123 27.07725
## 66       25.84924 25.01814 26.68034 24.57819 27.12029
## 67       25.84924 24.99092 26.70756 24.53656 27.16192
## 68       25.84924 24.96454 26.73394 24.49621 27.20227
## 69       25.84924 24.93892 26.75956 24.45703 27.24145
## 70       25.84924 24.91401 26.78448 24.41892 27.27956
## 71       25.84924 24.88974 26.80875 24.38181 27.31668
## 72       25.84924 24.86606 26.83242 24.34560 27.35288

plot(ramalansesopt,xlab="periode  waktu", ylab="Yt", main="Peramalan Dengan Metode SES")
lines(datasuhu.ts,lty=1,col="red")
points(datasuhu.ts,cex=0.5)

ses1 <- HoltWinters(training.ts, gamma=FALSE, beta=FALSE, alpha=0.99)
ramalansesopt<- forecast::forecast(ses1, h=21)
ramalansesopt <- as.data.frame(ramalansesopt); ramalansesopt

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 52       25.84924 25.63263 26.06585 25.51797 26.18052
## 53       25.84924 25.54444 26.15404 25.38308 26.31540
## 54       25.84924 25.47656 26.22192 25.27927 26.41921
## 55       25.84924 25.41927 26.27921 25.19165 26.50683
## 56       25.84924 25.36876 26.32972 25.11441 26.58407
## 57       25.84924 25.32308 26.37540 25.04454 26.65394
## 58       25.84924 25.28106 26.41743 24.98028 26.71820
## 59       25.84924 25.24194 26.45655 24.92045 26.77803
## 60       25.84924 25.20519 26.49329 24.86425 26.83424
## 61       25.84924 25.17042 26.52806 24.81108 26.88740
## 62       25.84924 25.13736 26.56112 24.76051 26.93797
## 63       25.84924 25.10576 26.59272 24.71219 26.98630
## 64       25.84924 25.07545 26.62303 24.66583 27.03265
## 65       25.84924 25.04629 26.65219 24.62123 27.07725
## 66       25.84924 25.01814 26.68034 24.57819 27.12029
## 67       25.84924 24.99092 26.70756 24.53656 27.16192
## 68       25.84924 24.96454 26.73394 24.49621 27.20227
## 69       25.84924 24.93892 26.75956 24.45703 27.24145
## 70       25.84924 24.91401 26.78448 24.41892 27.27956
## 71       25.84924 24.88974 26.80875 24.38181 27.31668
## 72       25.84924 24.86606 26.83242 24.34560 27.35288

error1<-ramalansesopt[,1] - testing.ts
mape1<-sum(abs(error1/testing.ts)*100)/length(testing.ts);mape1

## [1] 0.6213737

4. Double Exponential Smoothing

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

aopt<-out[which.min(out$MAPE),1]
bopt<-out[which.min(out$MAPE),2]

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

ramalandesopt<- forecast::forecast(Ho, h=21); 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.98920
## 72       26.01198 23.92854 28.09542 22.82564 29.19832

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

des1 <- HoltWinters(training.ts, gamma=FALSE, beta=0.43, alpha=0.36)
ramalandesopt<- forecast::forecast(des1, h=21)
ramalandesopt <- as.data.frame(ramalandesopt); 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.98920
## 72       26.01198 23.92854 28.09542 22.82564 29.19832

error2<-ramalandesopt[,1] - testing.ts
mape2<-sum(abs(error2/testing.ts)*100)/length(testing.ts);mape2

## [1] 0.349017

5. Metode Terbaik

PEMULUSAN <- c("SMA", "DMA", "SES", "DES")
MAPE <- round(c(mape.sma, mape.dma, mape1, mape2),2)
data.frame(PEMULUSAN,MAPE)

##   PEMULUSAN MAPE
## 1       SMA 0.46
## 2       DMA 0.56
## 3       SES 0.62
## 4       DES 0.35

Berdasarkan ukuran performa menggunakan MAPE diperoleh MAPE terkecil yaitu metode DES dengan nilai 0.35. Jika dilihat dari plot peramalan 21 tahun kedepan dan dibandingkan dengan data testing maka plot yang lebih tepat menggambarkan data testing adalah metode DES. Jadi metode terbaik yang digunakan adalah metode Double Exponential Smoothing.

Packages dan Data

lapply(c("dLagM","dynlm","MLmetrics","readxl","car","lmtest","caTools","ggplot2"),library,character.only=T)[[1]]

##  [1] "dLagM"     "dynlm"     "zoo"       "nardl"     "MLmetrics" "dplyr"    
##  [7] "TSA"       "TTR"       "tseries"   "expsmooth" "fma"       "ggplot2"  
## [13] "fpp2"      "smooth"    "greybox"   "Mcomp"     "forecast"  "readxl"   
## [19] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [25] "base"

pop <- read_excel("D:/FILE KULIAH SEMESTER 5/MINGGU 4/MPDW (P)/pop.xlsx")
pop$Suhu<-pop$Suhu
pop$Populasi<-pop$Populasi/10^6
train<-pop[1:51,];train

## # A tibble: 51 x 3
##    Tahun Populasi  Suhu
##    <dbl>    <dbl> <dbl>
##  1  1950     69.6  25.5
##  2  1951     71.0  25.5
##  3  1952     72.6  25.5
##  4  1953     74.2  25.6
##  5  1954     75.9  25.4
##  6  1955     77.7  25.4
##  7  1956     79.7  25.3
##  8  1957     81.7  25.5
##  9  1958     83.8  25.7
## 10  1959     86.0  25.5
## # ... with 41 more rows
## # i Use `print(n = ...)` to see more rows

test<-pop[-1:-51,];test

## # A tibble: 21 x 3
##    Tahun Populasi  Suhu
##    <dbl>    <dbl> <dbl>
##  1  2001     217.  25.9
##  2  2002     220.  26.0
##  3  2003     223.  25.9
##  4  2004     226.  25.9
##  5  2005     229.  25.9
##  6  2006     232.  25.8
##  7  2007     235.  25.8
##  8  2008     238.  25.7
##  9  2009     241.  26.0
## 10  2010     244.  26.1
## # ... with 11 more rows
## # i Use `print(n = ...)` to see more rows

a. Model Regresi

plot.reg <- lm(Suhu ~ Populasi, data = pop)
summary(plot.reg)

## 
## Call:
## lm(formula = Suhu ~ Populasi, data = pop)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36786 -0.08803  0.00311  0.09288  0.33276 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.512e+01  4.244e-02  591.86   <2e-16 ***
## Populasi    3.501e-03  2.366e-04   14.79   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1283 on 70 degrees of freedom
## Multiple R-squared:  0.7577, Adjusted R-squared:  0.7542 
## F-statistic: 218.9 on 1 and 70 DF,  p-value: < 2.2e-16

b. Uji Autokorelasi

dwtest(plot.reg$model)

## 
##  Durbin-Watson test
## 
## data:  plot.reg$model
## DW = 1.5085, p-value = 0.0122
## alternative hypothesis: true autocorrelation is greater than 0

p-value < 0.05 maka tolak H0, artinya cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada model dengan taraf nyata 5%.

1. Model Koyck

Metode Koyck didasarkan asumsi bahwa semakin jauh jarak lag peubah penjelas dari periode sekarang, maka semakin kecil pengaruh peubah lag terhadap peubah respon.

summary(lm(Suhu~Populasi,train))

## 
## Call:
## lm(formula = Suhu ~ Populasi, data = train)
## 
## 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(train$Populasi,train$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.2291461

e<-forecast(model = k, x=test$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 = test$Populasi, h = 21)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

Model koyck yang diperoleh adalah Yt = 1.938e+01 + 2.291e-01(Yt-1) + 2.55e-03(Xt). Interpretasi dari model tersebut adalah setiap kenaikan populasi satu satuan akan menyebabkan mean temperatur meningkat sebesar 2.55e-03 satuan. Nilai 2.291e-01 berarti jumlah nilai mean temperatur periode ke t dipengaruhi oleh nilai mean temperatur pada periode sebelumnya yaitu sebesar 2.291e-01 kali lipat dari sebelumnya.

#mape data testing
mapektest <- MLmetrics::MAPE(e$forecasts, test$Suhu);mapektest

## [1] 0.003405158

#akurasi data training
mapektrain <- GoF(k)[1,"MAPE"];mapektrain

## [1] 0.003830219

2. Distributed Lag Model with lag optim

qopt<-finiteDLMauto(formula=Suhu~Populasi,data=data.frame(train),model.type = "dlm",error.type = "AIC")
dlmod<-dlm(x=train$Populasi,y=train$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.78956

f<-forecast(dlmod,x=test$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 = test$Populasi, h = 21)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

forecast

## function (model, x, h = 1, interval = FALSE, level = 0.95, nSim = 500) 
## UseMethod("forecast")
## <bytecode: 0x0000000025315038>
## <environment: namespace:dLagM>

Hasil output software R menunjukkan bahwa Xt dan X(t-1) memiliki pengaruh yang signifikan terhadap model. Sehingga berdasarkan model distributed lag, rata-rata suhu di Indonesia dipengaruhi oleh jumlah populasi saat itu dan jumlah populasi pada periode sebelumnya.

#akurasi testing
mapedtest<-MLmetrics::MAPE(f$forecasts,test$Suhu);mapedtest

## [1] 0.00527501

#akurasi training
mapedtrain<-GoF(dlmod)[1,"MAPE"];mapedtrain

## [1] 0.003757701

3. Autoregressive Model With Lag Optimum

pqopt<-ardlBoundOrders(data.frame(train),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.8481

ARD<-ardlDlm(Suhu~Populasi,train,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=test$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 = test$Populasi, h = 21)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

#akurasi testing
mapeatest<-MLmetrics::MAPE(g$forecasts,test$Suhu);mapeatest

## [1] 0.008669551

#akurasi training
mapeatrain<-GoF(ARD)[1,"MAPE"];mapeatrain

## [1] 0.001428362

plot(test$Populasi,test$Suhu,col="black",type="b",ylim = c(25,27),lty=3)
points(test$Populasi,f$forecasts,col="coral")
lines(test$Populasi,f$forecasts,col="coral",lty=1)
points(test$Populasi,g$forecasts,col="blue")
lines(test$Populasi,g$forecasts,col="blue",lty=1)
points(test$Populasi,e$forecasts,col="green")
lines(test$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"))

Perbandingan Akurasi Training

#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

Dilihat dari pemodelan mengunakan data training maka metode yang terbaik adalah Autoregresive Distributed Model dengan nilai mape yang terkecil dibandingkan dengan metode yang lain yaitu mapenya sebesar 0.00143

Perbandingan Akurasi Testing

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

Dalam melakukan peramalan setelah dilihat performa setiap metode terhadap data testing maka metode terbaik adalah Model Koyck yang memiliki mape terkecil yaitu 0.0034

Pengecekan Autokorelasi

dwtest(k$model)

## 
##  Durbin-Watson test
## 
## data:  k$model
## DW = 2.0127, p-value = 0.4421
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(ARD$model)

## 
##  Durbin-Watson test
## 
## data:  ARD$model
## DW = 2.2578, p-value = 0.1216
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(dlmod$model)

## 
##  Durbin-Watson test
## 
## data:  dlmod$model
## DW = 1.8534, p-value = 0.2042
## alternative hypothesis: true autocorrelation is greater than 0

Autokol tertangani menggunakan model Distributed Lag dan Koyck Model dan ARDLM

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.

plot(train$Tahun, train$Suhu, main="Time Series Plot Data Training", 
     xlab="Tahun", ylab="Rata-rata Suhu")
lines(train$Tahun, train$Suhu, col="red")

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$Suhu)

Dilihat dari plot ACF juga tail off (slowly) yang artinya data tersebut tidak stasioner. Secara ekploratif bisa diketahui bahwa data yang digunakan tidaklah stasioner, namun untuk lebih meyakinkan kita akan melakukan uji signifikansi kestasioneran data menggunakan uji ADF Test.

adf.test(train$Suhu)#ga stasioner

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train$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$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: stationary

ts.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(2,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. Menggunakan Data Training

model1. <- stats::arima(train$Suhu, order = c(0,1,1), method = "ML") #ARIMA(0,1,1)
model2. <- stats::arima(train$Suhu, order = c(2,1,0), method = "ML") #ARIMA(2,1,0)
model3. <- stats::arima(train$Suhu, order = c(1,1,2), method = "ML") #ARIMA(1,1,2)
model4. <- stats::arima(train$Suhu, order = c(1,1,1), method = "ML") #ARIMA(1,1,1)
model5. <- stats::arima(train$Suhu, order = c(0,1,2), method = "ML") #ARIMA(0,1,1)
aic1.<-model1.$aic
aic2.<-model2.$aic
aic3.<-model3.$aic
aic4.<-model4.$aic
aic5.<-model5.$aic
AIC <- c(aic1.,aic2.,aic3.,aic4.,aic5.)
Model <- c("ARIMA (0,1,1)","ARIMA (2,1,0)","ARIMA (1,1,2)","ARIMA (1,1,1)","ARIMA (0,1,2)")
Akurasi <- data.frame(Model,AIC);Akurasi

##           Model       AIC
## 1 ARIMA (0,1,1) -48.26892
## 2 ARIMA (2,1,0) -41.92151
## 3 ARIMA (1,1,2) -45.04274
## 4 ARIMA (1,1,1) -46.88146
## 5 ARIMA (0,1,2) -46.99767

Akurasi[which.min(Akurasi$AIC),]

##           Model       AIC
## 1 ARIMA (0,1,1) -48.26892

Berdasarkan nilai AIC diperoleh model terbaik sementara yaitu ARIMA(0,1,1) dengan nilai AIC -48.27

5. Pendugaan Parameter Model

#ARIMA(0,1,1)
coeftest(model1.) 

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.706773   0.096699  -7.309 2.691e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(2,1,0)
coeftest(model2.) 

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.47111    0.13468 -3.4980 0.0004688 ***
## ar2 -0.29592    0.14095 -2.0995 0.0357734 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(1,1,2)
coeftest(model3.) 

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.15775    0.73863 -0.2136   0.8309
## ma1 -0.45158    0.71625 -0.6305   0.5284
## ma2 -0.22137    0.50196 -0.4410   0.6592

#ARIMA(1,1,1)
coeftest(model4.) 

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.136262   0.174744  0.7798    0.4355    
## ma1 -0.756608   0.098593 -7.6741 1.666e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(0,1,2)
coeftest(model5.) 

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)    
## ma1 -0.60555    0.14342 -4.2221 2.42e-05 ***
## ma2 -0.11319    0.13064 -0.8664   0.3863    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Selain dilihat dari nilai AIC terkecil, perlu juga dilihat dari signifikasi pendugaan parameter tiap model,ada dua model tentatif dengan semua parameter yang berpengaruh signifikan yaitu ARIMA(2,1,0) dan ARIMA (0,1,1). Dengan demikian, pemilihan dugaan model terbaik perlu dilihat dari model yang signifikan dan memiliki nilai AIC terkecil dan diperoleh model ARIMA (0,1,1)

auto.arima(train$Suhu)

## Series: train$Suhu 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##           ma1   drift
##       -0.8003  0.0089
## s.e.   0.0888  0.0043
## 
## sigma^2 = 0.01991:  log likelihood = 27.49
## AIC=-48.98   AICc=-48.46   BIC=-43.24

Jika di cek dengan fungsi auto arima ternyata sesuai yaitu model yang diperoleh adalah ARIMA(0,1,1)

6. Overfitting Model

Dilihat dari nilai AIC dan signifikansi parameter model, maka model yang dipilih sebelumnya adalah ARIMA(0,1,1). Selanjutnya kita akan melakukan Overfitting model ARIMA(0,1,1) dan melihat kembali model terbaik dari hasil overfitting-nya

#ARIMA(1,1,1)
model1.ov <- stats::arima(train$Suhu, order = c(1,1,1), method = "ML")
summary(model1.ov)

## 
## Call:
## stats::arima(x = train$Suhu, order = c(1, 1, 1), method = "ML")
## 
## Coefficients:
##          ar1      ma1
##       0.1363  -0.7566
## s.e.  0.1747   0.0986
## 
## sigma^2 estimated as 0.02007:  log likelihood = 26.44,  aic = -46.88
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.02604102 0.1403187 0.1098393 0.09884066 0.4287943 0.8436195
##                     ACF1
## Training set -0.03737112

coeftest(model1.ov) 

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.136262   0.174744  0.7798    0.4355    
## ma1 -0.756608   0.098593 -7.6741 1.666e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(0,1,2)
model2.ov <- stats::arima(train$Suhu, order = c(0,1,2), method = "ML")
summary(model2.ov)

## 
## Call:
## stats::arima(x = train$Suhu, order = c(0, 1, 2), method = "ML")
## 
## Coefficients:
##           ma1      ma2
##       -0.6055  -0.1132
## s.e.   0.1434   0.1306
## 
## sigma^2 estimated as 0.02002:  log likelihood = 26.5,  aic = -47
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE      MAPE      MASE
## Training set 0.02619196 0.1401519 0.1096855 0.0994333 0.4282045 0.8424385
##                     ACF1
## Training set -0.05343932

coeftest(model2.ov) 

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)    
## ma1 -0.60555    0.14342 -4.2221 2.42e-05 ***
## ma2 -0.11319    0.13064 -0.8664   0.3863    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan hasil uji signifikansi di atas, kedua model perluasan overfitting ARIMA (0,1,1) , yaitu ARIMA (1,1,1) dan ARIMA (0,1,2) tidak signifikan, sehingga model yang dipilih untuk tahapan analisis selanjutnya adalah ARIMA (0,1,1).

7. Diagnostik Model ARIMA(0,1,1)

sisaan1 <- residuals(model1.)
qqnorm(sisaan1) 
qqline(sisaan1)

Terlihat bahwa titik-titik data menyebar di sekitar garis quantile-quantile. Sehingga dapat disimpulkan bahwa tidak ada pelanggaran syarat kenormalan sisaan pada kasus ini.

tsdiag(model1., gof = 16, omit.initial =F, main="Uji Ljung")

Dari hasil uji Ljung Box di atas, gambar kedua dan ketiga memberikan gambaran terkait asumsi kebebasan sisaan. Pada gambar kedua terlihat bahwa tidak ada lag yang cut-off kecuali pada lag nol. Pada gambar ketiga terlihat bahwa titik-titik semua berada di atas 0.05 (alpha). Hal ini mengindikasikan bahwa sisaan saling bebas, sehingga dapat disimpulkan bahwa tidak ada pelanggaran syarat kebebasan sisaan pada kasus ini.

Box.test(sisaan1, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan1
## X-squared = 0.10691, df = 1, p-value = 0.7437

Berdasarkan LJung-Box test, diperoleh p-value (0.7437) > α (0.05), maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan antara lag saling bebas atau dapat dikatakan tidak ada autokorelasi antara sisaan lag pada taraf nyata 5%.

8. Peramalan ARIMA(0,1,1)

model1. <- stats::arima(training.ts, order = c(0,1,1), method = "ML")
dugaan <- fitted(model1.) 
plot.ts(training.ts, xlab = "tahun", ylab = "Data Yt", main = "Fitted-Value ARIMA(0,1,1)") 
points(training.ts) 
par(col = "red") 
lines(dugaan) 

ramalan <- forecast::forecast(Arima(training.ts, order=c(0,1,1),method="ML",include.drift = TRUE), h=21) 

plot(ramalan,lwd=2,main = "Peramalan Dengan ARIMA(0,1,1)")
box(lwd=2)

plot(ramalan)

plot(ramalan,xlab="periode  waktu", ylab="Yt", main="Peramalan Dengan Metode ARIMA(0,1,1)")
lines(datasuhu.ts,lty=1,col="red")
points(datasuhu.ts,cex=0.5)

## 9. Akurasi ARIMA(0,1,1)

ramalan. <- as.data.frame(ramalan); ramalan

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 52       25.89896 25.71811 26.07980 25.62238 26.17553
## 53       25.90782 25.72341 26.09223 25.62579 26.18985
## 54       25.91669 25.72878 26.10459 25.62930 26.20407
## 55       25.92555 25.73421 26.11689 25.63292 26.21818
## 56       25.93441 25.73970 26.12913 25.63662 26.23221
## 57       25.94328 25.74524 26.14131 25.64041 26.24614
## 58       25.95214 25.75084 26.15344 25.64429 26.26000
## 59       25.96100 25.75650 26.16551 25.64824 26.27377
## 60       25.96987 25.76220 26.17754 25.65227 26.28747
## 61       25.97873 25.76795 26.18951 25.65637 26.30109
## 62       25.98760 25.77375 26.20144 25.66055 26.31465
## 63       25.99646 25.77959 26.21333 25.66478 26.32814
## 64       26.00532 25.78547 26.22518 25.66909 26.34156
## 65       26.01419 25.79139 26.23699 25.67345 26.35493
## 66       26.02305 25.79735 26.24875 25.67787 26.36823
## 67       26.03192 25.80335 26.26049 25.68235 26.38148
## 68       26.04078 25.80938 26.27218 25.68688 26.39468
## 69       26.04964 25.81545 26.28384 25.69147 26.40782
## 70       26.05851 25.82154 26.29547 25.69610 26.42091
## 71       26.06737 25.82768 26.30707 25.70079 26.43396
## 72       26.07624 25.83384 26.31863 25.70552 26.44695

error.arima<-ramalan.[,1] - testing.ts
mape.arima<-sum(abs(error.arima/testing.ts)*100)/length(testing.ts);mape.arima

## [1] 0.2945791

MAPE nya 0.294