library(urca)
library(timeSeries)
## Loading required package: timeDate
library(readxl)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
list.files()
##  [1] "~$data_rapi.xlsx"      "~$OLAH DATA FIX.xlsx"  "data_rapi.xlsx"       
##  [4] "Kode-R-Arima-Fix.html" "Kode-R-Arima-Fix.Rmd"  "kode-R-arima.html"    
##  [7] "Kode R Arima Fix.Rmd"  "kode R arima.nb.html"  "kode R arima.Rmd"     
## [10] "OLAH DATA FIX.xlsx"    "rsconnect"
dat <- read_excel("data_rapi.xlsx")

dat <- as.data.frame(dat)

dat
susun_data <- reshape2::melt(dat, id = c("Tahun"))


wilayah <- colnames(dat)

wilayah <- wilayah[-c(1)]
                   
susun_data[,2] = factor(susun_data[,2], levels = wilayah)

susun_data
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.2
ggplot(data = susun_data, aes(x = Tahun, y = value, group = variable, colour = variable)) +
  geom_line()+
  geom_point()

Gambar di atas disajikan jumlah mahasiswa, dari tahun 2017 sampai tahun 2022 untuk setiap wilayah. Terlihat bahwa pada Tahun 2017, jumlah mahasiswa paling banyak berada pada luar pulau papua, yakni sebanyak 850. Terlihat bahwa pada Tahun 2018, jumlah mahasiswa paling banyak berada pada Kabupaten Manokwari, yakni sebanyak 792. Terlihat bahwa pada Tahun 2019, jumlah mahasiswa paling banyak berada pada Kabupaten Manokwari, yakni sebanyak 829. Terlihat bahwa pada Tahun 2020, jumlah mahasiswa paling banyak berada pada Kabupaten Manokwari, yakni sebanyak 640. Terlihat bahwa pada Tahun 2021, jumlah mahasiswa paling banyak berada pada Kabupaten Manokwari, yakni sebanyak 639. Terlihat bahwa pada Tahun 2022, jumlah mahasiswa paling banyak berada pada Kabupaten Manokwari, yakni sebanyak 748. Sementara pada Tahun 2018 dan Tahun 2020, tidak terdapat mahasiswa di Kabupaten Pegunungan Arfak.

#Uji Stasioner pada Level

Mahyus (2014:73) menyatakan untuk dapat menggunakan model ARMA ataupun ARIMA maka sebuah data haruslah stasioner. Apabila data tersebut stasioner pada ordo 0 (nol) atau pada “Level” maka data tersebut dapat menggunakan model ARMA. Sementara itu, apabila data stasioner pada ordo 1 (“first difference”) ataupun ordo 2 (“second difference”) maka data tersebut dapat menggunakan model ARIMA.

Selanjutnya dilakukan uji stasioner dengan pendekatan Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

library(urca)
library(timeSeries)

nilai_pvalue = vector(mode = "numeric")
nama_wilayah <- vector(mode = "character")
hasil_ujistasioner <- vector(mode = "character")
k = 0


for(i in 2 : (length(dat))  )
{

x = dat[,i]
nama <- colnames(dat)[c(i)]

hasil <- kpss.test(x, null="Trend")



k = k + 1

nilai_pvalue[k] <- hasil$p.value
nama_wilayah[k] <- nama

if(  nilai_pvalue[k] >= 0.05  )
{
  hasil_ujistasioner[k] = paste0("Tidak Stasioner pada Level; p >= 0.05")
}
if(  nilai_pvalue[k] < 0.05  )
{
  hasil_ujistasioner[k] = paste0("Stasioner pada Level; p < 0.05")
}



}
## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value greater than printed p-value
#hasil@cval

hasil_ujistasioner <- data.frame(nama_wilayah, nilai_pvalue, hasil_ujistasioner)




colnames(hasil_ujistasioner) = c("Wilayah", "P-Value", "Hasil Uji Stasioner")


hasil_ujistasioner

Berdasarkan hasil pengujian stasioner pada Level, terdapat beberapa wilayah, dengan hasil pengujian tidak stasioner pada level.

#Uji Stasioner pada First Difference

Oleh karena terdapat beberapa wilayah dengan hasil pengujian data tidak stasioner pada level, maka akan dilakukan pengujian stasioner pada tahap first difference, menggunakan uji stasioner dengan pendekatan Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

library(urca)
library(timeSeries)

nilai_pvalue = vector(mode = "numeric")
nama_wilayah <- vector(mode = "character")
hasil_ujistasioner <- vector(mode = "character")
k = 0


for(i in 2 : (length(dat))  )
{

x = dat[,i]

x <- diff(x)

nama <- colnames(dat)[c(i)]

hasil <- kpss.test(x, null="Trend")



k = k + 1

nilai_pvalue[k] <- hasil$p.value
nama_wilayah[k] <- nama

if(  nilai_pvalue[k] >= 0.05  )
{
  hasil_ujistasioner[k] = paste0("Tidak Stasioner pada First Difference; p >= 0.05")
}
if(  nilai_pvalue[k] < 0.05  )
{
  hasil_ujistasioner[k] = paste0("Stasioner pada First Difference; p < 0.05")
}



}
## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value
#hasil@cval

hasil_ujistasioner <- data.frame(nama_wilayah, nilai_pvalue, hasil_ujistasioner)




colnames(hasil_ujistasioner) = c("Wilayah", "P-Value", "Hasil Uji Stasioner")


hasil_ujistasioner

Berdasarkan pengujian stasioner pada first difference, masih terdapat beberapa wilayah, yang memberikan hasil tidak stasioner.

#Uji Stasioner pada Second Difference

Oleh karena terdapat beberapa wilayah dengan hasil pengujian data tidak stasioner pada first difference, maka akan dilakukan pengujian stasioner pada tahap second difference, menggunakan uji stasioner dengan pendekatan Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

nilai_pvalue = vector(mode = "numeric")
nama_wilayah <- vector(mode = "character")
hasil_ujistasioner <- vector(mode = "character")
k = 0


for(i in 2 : (length(dat))  )
{

x = dat[,i]

x <- diff(x)

x <- diff(x)

nama <- colnames(dat)[c(i)]

hasil <- kpss.test(x, null="Trend")



k = k + 1

nilai_pvalue[k] <- hasil$p.value
nama_wilayah[k] <- nama

if(  nilai_pvalue[k] >= 0.05  )
{
  hasil_ujistasioner[k] = paste0("Tidak Stasioner pada First Difference; p >= 0.05")
}
if(  nilai_pvalue[k] < 0.05  )
{
  hasil_ujistasioner[k] = paste0("Stasioner pada First Difference; p < 0.05")
}



}
## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value

## Warning in kpss.test(x, null = "Trend"): p-value smaller than printed p-value
#hasil@cval

hasil_ujistasioner <- data.frame(nama_wilayah, nilai_pvalue, hasil_ujistasioner)




colnames(hasil_ujistasioner) = c("Wilayah", "P-Value", "Hasil Uji Stasioner")


hasil_ujistasioner

Berdasarkan hasil uji stasioner pada second difference, diketahui seluruh nilai p < 0.05 untuk semua wilayah.

#ARIMA

Berikut disajikan hasil peramalan JUMLAH MAHASISWA untuk tahun 2023, 2024 dan 2025 menggunakan metode ARIMA.

nama_wilayah <- vector(mode = "character")

simpan_hasil <- 0

k = 0




for(i in 2 : length(dat) )
{

x <- dat[,i]

x <- diff(x)
x <- diff(x)


hasil_arima <- arima(x, order = c(0, 1, 0), method = c( "CSS")  )



prediksi <- predict(hasil_arima, 3)

prediksi <- prediksi$se
prediksi <- unlist(prediksi)

prediksi <- round(prediksi, digits = 0)


k = k + 1
nama_wilayah[k] = colnames(dat)[c(i)]

if(i == 2)
{
  simpan_hasil = prediksi
}
if(i > 2)
{
  simpan_hasil = rbind(simpan_hasil, prediksi)
}


}

rownames(simpan_hasil) = nama_wilayah
colnames(simpan_hasil) = c("2023","2024","2025")

print(simpan_hasil)
##                        2023 2024 2025
## KAB. FAK-FAK             25   35   43
## KAB. KAIMANA             10   14   17
## KAB. MANOKWARI          263  372  456
## KAB. MANOKWARI SELATAN   24   33   41
## KAB. PEGUNUNGAN ARFAK    10   14   18
## KAB. TELUK BINTUNI       64   90  110
## KAB. TELUK WONDAMA       65   92  112
## PAPUA                   212  299  367
## PAPUA BARAT DAYA        122  173  212
## PAPUA PEGUNUNGAN         45   64   78
## PAPUA TENGAH             84  118  145
## PAPUA SELATAN             4    6    7
## LUAR PULAU PAPUA        263  372  455