Alasan Memilih Topik Nilai Tukar Rupiah terhadap Dolar AS

Fluktuasi jangka panjang dalam nilai tukar dapat dikaitkan dengan pergeseran fundamental pasar atau faktor ekonomi, termasuk preferensi untuk komoditas domestik atau impor, suku bunga, tingkat produktivitas relatif, tingkat harga, dan hambatan perdagangan. Namun, ada kalanya variasi nilai tukar terlalu signifikan dan impulsif untuk diperhitungkan sepenuhnya oleh variabel tersebut karena terkadang nilainya dapat bervariasi lebih dari 2-3 poin persentase dalam satu hari. Teknik pasar aset yang digunakan untuk menentukan nilai tukar jangka pendek sangat cocok dalam situasi ini, yaitu dalam situasi krisis COVID-19 yang dinilai tidak biasa.

Pada waktu ini, umumnya investor akan mempertimbangkan dua faktor penting, yaitu suku bunga riil relatif dan ekspektasi tentang nilai tukar masa depan untuk mengatasi dilema dalam memilih antara aset domestik dan asing. Namun, investor asing mempertimbangkan faktor-faktor lain selain antisipasi perubahan suku bunga. Sebelum menginvestasikan uang dalam mata uang lain, investor akan didorong oleh ekspektasi tentang arah masa depan nilai tukar.

Perubahan ekspektasi relatif dari pertumbuhan ekonomi masa depan berfungsi sebagai saluran untuk efek COVID-19 pada pasar mata uang. Artinya, kondisi perekonomian yang mengalami kenaikan tajam dalam kasus COVID-19 akan memperkirakan nilai mata uang melemah karena pertumbuhan ekonomi akan dikoreksi secara menurun, serta kekhawatiran bahwa sektor ekonomi tertentu mungkin perlu ditutup. Tentu hal ini akan menurunkan permintaan mata uang lokal dan dana investasi yang diperlukan. Kondisi perekonomian dunia yang berubah drastis karena COVID-19 ini sangat menarik untuk dianalisis lebih dalam, terutama pada nilai tukar rupiah terhadap mata uang internasional.

Berdasarkan latar belakang yang telah disebutkan sebelumnya, kelompok kami memilih topik terkait nilai tukar rupiah dengan nilai tukar dolar amerika yang berhubungan dengan adanya pandemi COVID-19 ini untuk dianalisis lebih lanjut.

Pentingnya Meramal Nilai Tukar Rupiah terhadap dolar AS

Nilai tukar mata uang suatu negara adalah salah satu indikator penting dalam suatu perekonomian. Nilai tukar juga mempunyai implikasi yang luas, baik dalam konteks ekonomi domestik maupun internasional, mengingat hampir semua negara di dunia melakukan transaksi valuta asing (Yanthi dan Artini 2013). Menurut Mulyani (2020) pakar ekonomi dari Universitas Sebelas Maret (UNS) Surakarta menyoroti lemahnya nilai tukar rupiah terhadap dolar Amerika Serikat, sebagai akibat pandemi COVID-19.Nilai rupiah saat ini sudah menyentuh angka Rp16.000, dan menjadi yang terlemah dalam kurun waktu 5 tahun terakhir (Alfira et al. 2021).

Modal asing yang masuk pada suara negara mempengaruhi tingkat suku bunga, tingkat inflasi dan juga nilai tukar rupiah terhadap dolar Amerika. Akibat dari melemahnya Rupiah terhadap USD adalah naiknya harga barang-barang impor, sehingga akan berdampak langsung bagi masyarakat luas (Jana 2015). Mengingat pentingnya menjaga kestabilan nilai tukar Rupiah terhadap USD demi kepentingan ekonomi dalam negeri, dibutuhkan forecasting (peramalan) nilai tukar Rupiah terhadap USD agar para pengambil kebijakan bisa bergerak lebih cepat dalam mengambil keputusan.

Landasan Teori dari Penelitian Sebelumnya

1. Pandemi yang melanda hampir seluruh dunia, menyebabkan ketidakpastian pada sektor ekonomi. Hal ini menyebabkan Rupiah melemah 6,7 persen terhadap dolar hingga bulan September 2020. Kondisi terlemah Rupiah yaitu pada tanggal 23 Maret 2020, dimana nilai tukar Rupiah berada pada nilai tukar Rp.16.575. Sejak awal tahun 2020 hingga pertengahan September 2020 capital outflow mengalami peningkatan menjadi Rp169,22 triliun sebagai akibat dari kepanikan pasar keuangan global pada awal pandemi. Fluktuasi nilai tukar Rupiah juga didorong sejumlah faktor lain diantaranya, penyebaran Covid-19 yang mulai memasuki fase kedua dibeberapa negara, kembali memanasnya perang dagang antara Cina dan Amerika, dan kondisi Brexit yang tidak menentu. Hal tersebut mendorong volatilitas di pasar keuangan domestik dan globalselama September 2020 (Suryaputri & Kurniawati 2021).

2. Mulyani (2020), pakar ekonomi dari Universitas Sebelas Maret (UNS) Surakarta, menyoroti lemahnya nilai tukar rupiah terhadap dolar Amerika Serikat, sebagai akibat pandemi Covid-19. Nilai rupiah saat ini sudah menyentuh angka Rp16.000, dan menjadi yang terlemah dalam kurun waktu 5 tahun terakhir (Alfira et al. 2021).

3. Hasil pengujian membuktikan bahwa betul nilai tukar rupiah mengalami perubahan secara signifkan setelah pemberlakuan PSBB. Namun, perubahan signifikan tersebut cenderung bergerak ke arah negatif atau dengan kata lain perubahan signifikan tersebut merupakan penurunan nilai tukar rupiah terhadap dolar AS. Penurunan nilai tukar rupiah terhadap dolar AS pasca adanya pengumuman PSBB ini sejalan dengan banyaknya kontroversi di masyarakat terkait pembelakuan PSBB dan kontroversi tersebut menumbuhkan sentimen negatif masyarakat terhadap kinerja pemerintah (Made 2020).

Time Series

Analisis time series yaitu jenis analisis yang melibatkan penggunaan data deret waktu atau data time series untuk membuat model yang akan digunakan sebagai dasar peramalan. Model time series stasioner merupakan jenis model time series untuk data stasioner yakni data yang mempunyai tingkat rata-rata dan varian konstan atau tidak berubah signifikan dari waktu ke waktu. Model stokastik dalam mendeskripsikan data time series terdiri dari apa yang disebut model stasioner dan model tidak stasioner (nonstasioner).

Model Stasioner dan Tak Stasioner

Model stasioner merupakan jenis model time series untuk data stasioner yakni data yang mempunyai tingkat rata-rata dan varian konstan atau tidak berubah signifikan dari waktu ke waktu. Adapun model time series untuk data stasioner meliputi model autoregressive (AR), moving average (MA), dan autoregressive moving average (ARMA). Namun, peramalan yang banyak dilakukan dalam bidang industri, bisnis, dan ekonomi, di mana terdapat banyak data deret waktu (time series), sering kali lebih baik direpresentasikan dengan model nonstasioner yakni model yang tidak memiliki rata-rata konstan alami dari waktu ke waktu (Sitanini et al. 2020). Oleh karena itu, tidak mengherankan bahwa banyak dari metode peramalan yang menggunakan model time series nonstasioner seperti rata-rata bergerak tertimbang eksponensial (exponentially weighted moving average).

Beberapa model time series lainnya untuk data yang tidak stasioner meliputi model ARIMA (autoregressive integrated moving average), SARIMA (seasonal autoregressive integrated moving average), ARCH (autoregressive conditional heteroskedasticity), GARCH (generalized autoregressive conditional heteroskedasticity) dan sebagainya.

Dalam analisis time series, terdapat dua konsep penting yang dikenal dengan fungsi autokorelasi (autocorrelation function, ACF) dan fungsi autokorelasi parsial (partial autocorrelation function, PACF). Kedua jenis korelasi tersebut biasanya digunakan dalam spesifikasi model.

Uji Stasioneritas Data Time Series

Stasioneritas merupakan konsep penting dalam analisis time series. Seperti telah dibahas sebelumnya, data time series dikatakan stasioner apabila nilai rata-rata dan variansnya tidak mengalami perubahan yang secara sistematik sepanjang waktu atau dengan kata lain, rata-rata dan variansnya konstan.

Apabila data tidak stasioner, maka sebelum mencari model terbaik, data yang ada perlu distasionerkan terlebih dahulu. Apabila data yang digunakan dalam model ada yang tidak stasioner, maka ada kemungkinan terjadinya spurious regression. Spurious regression adalah regresi yang memiliki \(R^2\) yang tinggi, tetapi tidak mempunyai hubungan yang berarti.

Secara umum, uji untuk mengetahui kestasioneran data time series dapat dikategorikan menjadi tiga, yakni melalui grafik, korelogram, dan uji unit root (ADF Test & Phillips-Perron Test) (Al’afi et al. 2020).

1. Uji stasioneritas data time series menggunakan grafik dapat dilakukan dengan membuat plot antara data observasi dengan variabel waktu (t). Jika dari plot tersebut, terlihat rata-rata dan variansnya konstan, maka data time series dikatakan stasioner. Sebaliknya, jika grafik tidak menunjukkan rata-rata dan varians konstan, maka data time series dikatakan tidak stasioner.

2. Pada dasarnya, korelogram merupakan metode pengujian stasioneritas data time series berdasarkan fungsi autokorelasi (ACF) yang diperoleh dengan memplotkan antara \(ρ_k\) dan \(k\) (lag). Untuk data yang stasioner, korelogram menurun dengan cepat seiring dengan meningkatnya \(k\). Sedangkan untuk data yang tidak stasioner, korelogram cenderung tidak menuju nol (turun lambat).

3. Uji Stasioneritas Augmented Dickey-Fuller (ADF) adalah pengujian yang dilakukan terhadap data deret waktu (time series) untuk mengetahui apakah data deret waktu tersebut stasioner atau tidak. Sebagian analisis deret waktu mensyaratkan agar data stasioner terlebih dahulu sebelum dilakukan analisis lebih lanjut, misalnya analisis data menggunakan Arima. Oleh karena itu, untuk memenuhi syarat tersebut, maka uji stasioneritas perlu dilakukan. Kestasioneran yang dimaksud adalah stasioner data terhadap rata-rata dan stasioner terhadap varian.

   Hipotesis uji stasioner data deret waktu menggunakan ADF adalah:

H0 : Data tidak stasioner

H1 : data stasioner

Metode Pemulusan

1. Single Moving Average (SMA)

   Metode Single Moving Average lebih cocok digunakan untuk data yang tidak memiliki trend naik maupun turun dan atau musiman (Astuti et al. 2019). Pada metode SMA, semakin panjang jangka waktu moving average maka pemulusan yang dihasilkan akan semakin halus.

2. Double Moving Average (DMA)

   Metode Double Moving Average merupakan metode peramalan jangka pendek yang dapat mengatasi data tren (Layakana dan Iskandar 2020). Metode ini merupakan pengembangan dari metode Simple Moving Average.

3. Single Exponential Smoothing (SES)

   Metode Exponential Smoothing adalah metode pemulusan dengan melakukan pembobotan menurun secara eksponensial (Mutmainnah 2019). Nilai 𝜆 merupakan parameter pemulusan dengan nilai 0 < λ < 1, yang berarti jika 𝜆 mendekati nilai nol maka pengaruh (pembobotan) periode lalu cukup besar, sementara jika mendekati nilai 1 maka pengaruh periode sekaranglah yang cukup besar.

4. Double Exponential Smoothing (DES)

   Metode Double Exponential Smoothing tidak jauh berbeda dengan metode Single Exponential Smoothing yang telah dijelaskan sebelumnya, hanya saja pada metode ini dilakukan dua tahap metode SES. Metode ini cocok digunakan untuk data dengan pola tren.

Teknik Evaluasi Kesesuaian Model

1. Sum of Square Error (SSE) merupakan penjumlahan dari kuadrat hasil pengurangan antara data aktual dengan hasil peramalan untuk setiap periode waktu.

2. Mean Squared Error (MSE) merupakan penjumlahan kuadrat semua kesalahan peramalan pada setiap periode dibagi dengan jumlah periode peramalan. MSE memperkuat pengaruh dari angka-angka kesalahan peramalan yang besar, tetapi memperkecil pengaruh dari angka-angka kesalahan peramalan yang kecil.

3. Root Mean Squared Error (RMSE) merupakan nilai rata-rata dari jumlah kuadrat kesalahan. Nilai RMSE yang rendah menunjukkan bahwa variasi nilai yang dihasilkan mendekati variasi nilai observasinya.

4. Mean Absolute Percentage Error (MAPE) merupakan rata-rata kesalahan absolut dari peramalan yang dihasilkan dalam bentuk persentasenya terhadap data aktual. MAPE merupakan satu-satunya satuan yang dinyatakan dalam bentuk persen. MAPE juga dapat digunakan untuk membandingkan ketepatan dari teknik yang sama atau berbeda dalam dua deret yang sangat berbeda dan mengukur ketepatan nilai dugaan model yang dinyatakan dalam bentuk rata-rata persentase absolut kesalahan.

5. Mean Absolute Deviation merupakan metode yang digunakan untuk mengevaluasi metode peramalan menggunakan jumlah dari kesalahan-kesalahan yang absolut. Mean Absolute Deviation (MAD) mengukur ketepatan ramalan dengan merata-rata kesalahan dugaan (nilai absolut masingmasing kesalahan).

Model Time Series Stasioner

1. Model autoregressive (AR) adalah model regresi time series yang menghubungkan nilai pengamatan aktual dengan nilai pengamatan sebelumnya. Ini dapat dilakukan ketika sebuah pengamatan tidak lepas dari pengamatan yang terjadi sebelumnya. Konsep dasar model ini yaitu meregresikan pengamatan aktual dengan nilai pengamatan sebelumnya untuk melakukan peramalan nilai ke depan. Plot PACF dapat digunakan untuk mengidentifikasi orde dari model AR(p).

Adapun persamaan model autoregressive dengan order p atau AR(p) dapat dituliskan sebagai berikut:

\[ y_t=ϕ+ϕ_1Y_{t-1}+ϕ_2Y_{t-2}+...+ϕ_pY_{t-p}+\varepsilon_t \] di mana: \(\varepsilon_t∼WN(0,σ^2)\)

Adapun persamaan untuk distributed-lag model dapat dituliskan sebagai berikut:

\[ y_t=ϕ_0+ϕ_1X_t+ϕ_2Y_{t-1}+ϕ_3Y_{t-2}+...+ϕ_pX_{t-p}+\varepsilon_t \]

2. Berbeda dengan model autoregressive (AR) yang menggunakan pengamatan masa lalu (past values) sebagai variabel bebas dalam regresi, pada model moving average (MA) digunakan nilai error dalam model regresi. Fungsi autokorelasi (autocorrelation function, ACF) yang telah dibahas pada artikel sebelumnya, dapat digunakan untuk menentukan atau mengidentifikasi orde dari model MA(q), karena plot ACF diharapkan akan berhenti terpotong setelah lag q (cut off after lag q).

Adapun persamaan model moving average dengan orde q atau MA(q) dapat dituliskan sebagai:

\[ y_t=μ+wt+θ_1w_{t-1}+θ_2w_{t-2}+...+θ_pw_{t-q}+\varepsilon_t \]

3. Model ARMA merupakan campuran antara model autoregressive (AR) yang mengasumsikan bahwa data sekarang dipengaruhi oleh data sebelumnya dan model moving average (MA) yang mengasumsikan bahwa data sekarang dipengaruhi oleh nilai residual data sebelumnya.

Secara umum, model ARMA dengan orde p dan q atau ARMA(p,q) diberikan sebagai berikut:

\[ y_t=γ+ϕ_1Y_{t-1}+ϕ_2Y_{t-2}+...+ϕ_pY_{t-q}+\varepsilon_t-θ_1e_{t-1}-θ_2e_{t-2}-...-θ_qe_{t-q} \]

Model Time Series Tak Stasioner

1. Differencing (pembedaan) digunakan untuk menstasionerkan data yang tidak stasioner. Jika pembedaan pertama (first difference) berhasil membuat data menjadi stasioner, berarti kita peroleh orde d = 1 untuk Arima. Selanjutkan kita menentukan orde p dan q berdasarkan plot ACF dan PACF sampel dari data yang telah distasionerkan tersebut.

   Namun, ada kalanya pembedaan pertama belum berhasil membuat data menjadi stasioner. Jika kondisi demikian terjadi maka kita dapat lakukan pembedaan kedua (second difference), dan jika datanya telah stasioner berarti kita peroleh d = 2, dan selanjutnya menentukan p dan q berdasarkan plot ACF dan PACF sampel dari data yang telah stasioner.

   Perlu diketahui bahwa untuk diferensiasi atau pembedaan yang lebih tinggi tidaklah disarankan. Umumnya, paling banyak kita hanya melakukan pembedaan sebanyak dua kali.

2. Model ARIMA atau dikenal juga dengan model Box-Jenkins menghadirkan kelas model yang sangat kuat dan fleksibel untuk analisis dan peramalan data deret waktu (Schaffer et al. 2021). Selama bertahun-tahun, model tersebut telah dengan sangat berhasil diterapkan pada banyak masalah dalam penelitian dan praktik.

   Model Arima(p,d,q) dapat dituliskan sebagai:

\[ y_t=γ_0+θ_1Y_{t-1}+θ_2Y_{t-2}+...+θ_nY_{t-q}-λ_1e_{t-1}-λ_2e_{t-2}-...-λ_ne_{t-q} \]

3. Autoregressive Conditional Heteroscedasticity (ARCH)

   Model ARCH digunakan untuk mengatasi keheterogenan ragam dengan

   memodelkan fungsi rataan dan fungsi ragam secara simultan (Wijaya & Nugraha,

   2020). Model umumnya sebagai berikut:

   

4. Generalized Autoregressive Conditional Heteroscedasticity (GARCH)

   

Import Data

gs4_deauth()
data <- read_sheet("https://docs.google.com/spreadsheets/d/1gVxOS7VYqzgGWHF4ht2xU_AASX00KgwJkeEjm73_ftA/edit#gid=1554735906")

## v Reading from "Data MPDW".

## v Range 'Copy of data asli'.

data.ts <- ts(data[,2])

Eksplorasi Pola Data

Data pada algoritma ini umumnya perlu dibagi menjadi 2 bagian, yaitu data training dan data testing. Data training nantinya akan digunakan untuk melatih algoritma dalam mencari model yang sesuai, sedangkan data testing akan dipakai untuk menguji dan mengetahui performa model yang didapatkan pada tahapan testing.

data$Date <- as.Date(data$Date)
ggplot(data, aes(x=Date, y=Nilai_tukar)) +
  geom_line() + 
  scale_x_date(date_breaks = "8 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data Nilai Tukar Rupiah Mingguan",
       subtitle = "(Agustus 2017 - Oktober 2022)") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA))

summary(data$Nilai_tukar)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13197   14026   14301   14300   14540   16377

sd(data$Nilai_tukar)

## [1] 497.0225

Berdasarkan plot time series yang ditampilkan di atas, pola data nilai tukar yang cenderung berbentuk campuran (tren naik di awal, cenderung konstan pada periode selanjutnya, kembali berpola tren naik, lalu kembali konstan), tidak dapat disebut berpola stasioner ataupun tren. Rata-rata nilai tukar rupiah terhadap dolar AS mingguan sebesar 14300 dengan nilai standar deviasi (penyebaran data) sebesar 497.02 dari rata-rata. Data nilai tukar rupiah mingguan mencapai nilai tertinggi pada angka 16377 yaitu pada tanggal 30 Maret 2020 dan menyentuh nilai terendah pada angka 13197 pada tanggal 4 September 2017.

Hasil eksplorasi data nilai tukar rupiah dari bulan Agustus 2017 sampai Oktober 2022 menunjukkan bahwa pada bulan September 2018 sampai November 2018 nilai tukar rupiah di Indonesia mengalami pelemahan terhadap nilai dolar Amerika Serikat. Hal tersebut disebabkan oleh beberapa faktor seperti terjadi bencana alam yang melanda di daerah Sulawesi, Yogyakarta, dan beberapa daerah lainnya. Tentu hal ini menyebabkan perlambatan pertumbuhan ekonomi Indonesia. Pada awal tahun 2019 sampai akhir tahun 2019, nilai tukar rupiah Indonesia mengalami sedikit penguatan. Hal tersebut dipengaruhi oleh beberapa faktor, salah atunya adalah adanya pemilu presiden dan penetapannya. Tentu kebijakan yang dirumuskan oleh presiden memberikan rasa percaya diri pada para ekonom dan investor untuk menanamkan modalnya di dalam negeri. Namun pada awal tahun 2020, tepatnya bulan Januari sampai bulan Juni, nilai tukar rupiah terhadap dolar Amerika Serikat kembali mengalami pelemahan, bahkan mencapai nilai tertinggi sejak 5 tahun terakhir, yaitu mencapai nilai Rp16.000,- untuk 1 dolar USD. Pelemahan tersebut disebabkan oleh adanya pandemi Covid-19 yang membuat semua negara di dunia kewalahan dan belum siap menangani pandemi ini. Perubahan dalam tatanan ekonomi yang sangat signifikan dan serba mendadak membuat sektor ekonomi menjadi salah satu sektor yang paling terdampak dengan mewabahnya pandemi Covid-19. Hal ini tentu berdampak pada nilai tukar rupiah mengalami pelemahan terhadap nilai dolar USD. Kemudian, ketegangan antara Rusia dengan UNI EROPA dan Amerika Serikat, ketegangan antara Amerika Serikat dengan Tiongkok dan Arab Saudi, serta adanya indikasi perang antara Rusia dengan Ukraina di pertengahan tahun 2022 kembali mendorong pelemahan nilai tukar rupiah terhadap dolar USD yang cukup signifikan. Hal ini disebabkan peristiwa-peristiwa di atas berdampak pada peningkatan krisis energi dan krisis pangan terhadap seluruh dunia serta memperburuk keadaan resesi dan inflasi di dunia, khususnya negara berkembang.

Selanjutnya data yang memiliki pola cenderung mirip akan dibagi 2 menjadi data training dan testing.

Eksplorasi Pembagian Data

data$Date <- as.Date(data$Date)
p=(230/271)
f=as.integer(p*271)
jumlah.train <- as.integer(p*nrow(data))
ggplot(data, aes(x=Date, y=Nilai_tukar)) +
  geom_line() + 
  scale_x_date(date_breaks = "8 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data Nilai Tukar Rupiah Mingguan",
       subtitle = "(Agustus 2017 - Oktober 2022)") +
  geom_vline(aes(xintercept = Date[jumlah.train], 
                 col="Pemisah_data_train"), lty=2, lwd=.7) +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA)) +
  scale_color_manual(name = "", values = c(Pemisah_data_train="blue"))

Berdasarkan hasil eksplorasi dan pertimbangan terkait pola data, didapatkan titik potong yang dinilai tepat dalam membagi data training dan data testing, yaitu data training dari periode 14 Agustus 2017 sampe 3 Januari 2022 dan data testing 10 Januari 2022 smpe 17 Oktober 2022. Titik pemisah data training dan data testing tersebut bertepatan pada tanggal 3 Januari 2022.

data[jumlah.train,1]

Pembagian Data Training dan Data Testing

train=1:f
training.ts <- ts(data$Nilai_tukar[train])
testing.ts <- ts(data$Nilai_tukar[-train], start = f+1)

Berikut adalah plot data training nilai tukar rupiah:

ggplot(data[1:230,], aes(x=Date, y=Nilai_tukar)) +
  geom_line() + 
  scale_x_date(date_breaks = "6 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data Training Nilai Tukar Rupiah Mingguan",
       subtitle = "(Agustus 2017 - Januari 2022)") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA))

Berikut adalah plot data testing nilai tukar rupiah:

ggplot(data[231:271,], aes(x=Date, y=Nilai_tukar)) +
  geom_line() + 
  scale_x_date(date_breaks = "2 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data Testing Nilai Tukar Rupiah Mingguan",
       subtitle = "(Januari 2022 - Oktober 2022)") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA))

Plot Time Series dan Korelogram

acf(training.ts,lag.max = 20)

Terlihat dari plot ACF di atas bahwa adanya tails off (meluruh menjadi nol secara asimptotik) yang mengindikasikan bahwa data tidak stasioner.

Augmented Dickey-Fuller Test

Kestasioneran data diuji menggunakan ADF-Test dengan hipotesis sebagai berikut:

H0 : Data tidak stasioner

H1 : Data stasioner

adf.test(training.ts) 

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.ts
## Dickey-Fuller = -3.0781, Lag order = 6, p-value = 0.123
## alternative hypothesis: stationary

Hasil di atas menunjukkan bahwa \(P-value\) hasil uji di atas lebih besar dari \(alpha=0.05\), sehingga terima H0. Artinya, pada taraf nyata 5% tidak cukup bukti untuk menyatakan bahwa data nilai tukar rupiah terhadap dolar AS stasioner.

1. Double Moving Average (DMA)

Sebelum melakukan pemulusan dan peramalan data Nilai Tukar, terlebih dahulu perlu ditentukan nilai n yang menghasilkan nilai eror minimum.

akurasi = data.frame()
for(n in 1:50){
  data.sma<-SMA(training.ts, n=n)
  dma <- SMA(data.sma, n=n)
  At <- 2*data.sma - dma
  Bt <- 2/(n-1)*(data.sma - dma)
  data.dma<- At+Bt
  data.ramal2<- c(NA, data.dma)
  
  error.dma = training.ts-data.ramal2[1:length(training.ts)]
  SSE.dma = sum(error.dma[(2*n):length(training.ts)]^2, na.rm = T)
  MSE.dma = mean(error.dma[(2*n):length(training.ts)]^2, na.rm = T)
  MAPE.dma = mean(abs((error.dma[(2*n):length(training.ts)]/training.ts[(2*n):length(training.ts)])*100), na.rm = T)
  
  akurasi.dma <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
  row.names(akurasi.dma)<- c("SSE", "MSE", "MAPE")
  colnames(akurasi.dma) <- c(paste("Akurasi m =", n))
  akurasi = rbind(akurasi, c(n=n, akurasi=akurasi.dma[3]))
}
colnames(akurasi) = c("n", "MAPE")
akurasi[order(akurasi[,2]),]

Berdasarkan output di atas, didapat bahwa nilai n=3 menghasilkan nilai MAPE paling minimum. Namun selain berdasarkan nilai eror kita juga perlu melihat plot hasil pemulusan data yang dihasilkan.

par(mfrow=c(3,2))
for(n in seq(3,30,5)) {
  data.sma<-SMA(training.ts, n=n)
  data.ramal<-c(NA,data.sma)
  data.gab<-cbind(aktual=c(training.ts,rep(NA,5)),pemulusan=c(data.sma,rep(NA,5)),ramalan=c(data.ramal,rep(data.ramal[length(data.ramal)],4)))
  ts.plot(data.gab[,1], xlab="Time Period", ylab="Data Nilai Tukar", main= paste0("SMA N=",n ," Data Close"))
          points(data.gab[,1])
          lines(data.gab[,2],col="green",lwd=2)
          lines(data.gab[,3],col="red",lwd=2)
          legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8,col=c("black","green","red"), cex=0.8)
}

Berdasarkan plot pemulusan grafik di atas, terlihat bahwa plot pemulusan ketika n=3 masih kurang baik karena plot yang dihasilkan belum cukup mulus. Sementara plot pemulusan dengan n=8 dan n=13 telah terlihat cukup mulus dan cukup optimum. Nilai eror untuk n=13 juga cukup kecil dibandingkan n lainnya. Oleh karena itu, akan dilakukan tiga pemulusan dengan metode DMA yang dipilih berdasarkan nilai eror dan grafik paling optimum, yaitu ketika menggunakan n=3, n=8, dan n=13.

Berikut adalah hasil peramalan data training menggunakan metode DMA untuk 41 periode ke depan beserta plot pemulusan yang dihasilkan.

#Single Moving  Average dengan n=3
df_ts_sma <- TTR::SMA(training.ts,  n=3)
#Double Moving  Average dengan n=3
df_ts_dma <- TTR::SMA(df_ts_sma, n=3)
At <- 2*df_ts_sma-df_ts_dma
Bt <- 2/(3-1)*(df_ts_sma-df_ts_dma)
pemulusan_dma <- At+Bt
data.ramal1<- c(NA, pemulusan_dma)

t = 1:41
f = c()

for (i in t) {
  f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}

data.gab1 <- cbind(aktual = c(training.ts,rep(NA,41)), pemulusan1 = c(df_ts_sma,rep(NA,41)),pemulusan2 = c(pemulusan_dma, rep(NA,41)),At = c(At, rep(NA,41)), Bt = c(Bt,rep(NA,41)),ramalan = c(data.ramal1, f[-1]))
tail(data.gab1)

##        aktual pemulusan1 pemulusan2 At Bt  ramalan
## [266,]     NA         NA         NA NA NA 13125.59
## [267,]     NA         NA         NA NA NA 13095.44
## [268,]     NA         NA         NA NA NA 13065.30
## [269,]     NA         NA         NA NA NA 13035.15
## [270,]     NA         NA         NA NA NA 13005.01
## [271,]     NA         NA         NA NA NA 12974.86

#Plot time series
ts.plot(data.gab1[,1], xlab="Weekly", ylab="Close", main="Plot Pemulusan dan Peramalan dengan Metode DMA n=3")
points(data.gab1[,1])
lines(data.gab1[,3],col="green",lwd=2)
lines(data.gab1[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)

# Ukuran Keakuratan
error.dma <- data.gab1[, 1] - data.gab1[, 3]

## SSE (Sum Square Error)
SSE.dma <- sum(error.dma^2, na.rm = T)

## MSE (Mean Squared Error)
MSE.dma <- mean(error.dma^2, na.rm = T)

## RMSE (Root Mean Square Error)
RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))

## MAD (Mean Absolute Deviation)
MAD.dma <- mean(abs(error.dma), na.rm = T)

## MAPE (Mean Absolute Percentage Error)
r.error.dma <- (error.dma/data.gab1[, 1])*100 # Relative Error
MAPE.dma <- mean(abs(r.error.dma), na.rm = T)


akurasi1 <- data.frame(
          "Akurasi Data Training DMA" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
          "Double Moving Average n=3" = c(SSE.dma, MSE.dma, MAPE.dma, RMSE.dma, MAD.dma))
akurasi1

selisihdma1<-data.gab1[231:271,6]-testing.ts
selisihdma1

## Time Series:
## Start = 231 
## End = 271 
## Frequency = 1 
##  [1]  -133.8558  -179.2005  -261.2445  -310.1898  -284.6341  -335.9788
##  [7]  -334.4232  -412.0677  -386.0120  -432.9564  -477.2015  -522.3460
## [13]  -551.6898  -577.9349  -694.8792  -777.6234  -796.4681  -997.7132
## [19] -1034.3577  -913.0023  -853.7464 -1061.5911 -1321.3353 -1329.6800
## [25] -1517.9242 -1549.8695 -1599.4134 -1607.6585 -1441.3031 -1638.0472
## [31] -1359.1917 -1500.9359 -1609.5804 -1681.4251 -1668.5697 -1868.0140
## [37] -1993.6585 -2237.8031 -2249.9476 -2457.1927 -2645.1370

SSEtestingdma1<-sum(selisihdma1^2)
MSEtestingdma1<-SSEtestingdma1/length(testing.ts)
MAPEtestingdma1<-sum(abs(selisihdma1/testing.ts)*100)/length(testing.ts)

akurasitestingdma <- matrix(c(SSEtestingdma1,MSEtestingdma1,MAPEtestingdma1),nrow=3,ncol=1)
row.names(akurasitestingdma)<- c("SSE", "MSE", "MAPE")
colnames(akurasitestingdma) <- "Akurasi Data Testing DMA n=3"
akurasitestingdma

##      Akurasi Data Testing DMA n=3
## SSE                  6.963651e+07
## MSE                  1.698451e+06
## MAPE                 7.473236e+00

#Single Moving  Average dengan n=8
df_ts_sma2 <- TTR::SMA(training.ts,  n=8)
#Double Moving  Average dengan N=8
df_ts_dma2 <- TTR::SMA(df_ts_sma2, n=8)
At <- 2*df_ts_sma2-df_ts_dma2
Bt <- 2/(8-1)*(df_ts_sma2-df_ts_dma2)
pemulusan_dma2 <- At+Bt
data.ramal2<- c(NA, pemulusan_dma2)

t = 1:41
f = c()

for (i in t) {
  f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}

data.gab2 <- cbind(aktual = c(training.ts,rep(NA,41)), pemulusan1 = c(df_ts_sma2,rep(NA,41)),pemulusan2 = c(pemulusan_dma2, rep(NA,41)),At = c(At, rep(NA,41)), Bt = c(Bt,rep(NA,41)),ramalan = c(data.ramal2, f[-1]))
tail(data.gab2)

##        aktual pemulusan1 pemulusan2 At Bt  ramalan
## [266,]     NA         NA         NA NA NA 14725.57
## [267,]     NA         NA         NA NA NA 14735.75
## [268,]     NA         NA         NA NA NA 14745.93
## [269,]     NA         NA         NA NA NA 14756.12
## [270,]     NA         NA         NA NA NA 14766.30
## [271,]     NA         NA         NA NA NA 14776.48

#Plot time series
ts.plot(data.gab2[,1], xlab="Weekly", ylab="Close", main="Plot Pemulusan dan Peramalan dengan Metode DMA N=8")
points(data.gab2[,1])
lines(data.gab2[,3],col="green",lwd=2)
lines(data.gab2[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)

# Ukuran Keakuratan
error.dma2 <- data.gab2[, 1] - data.gab2[, 3]

## SSE (Sum Square Error)
SSE.dma2 <- sum(error.dma2^2, na.rm = T)

## MSE (Mean Squared Error)
MSE.dma2 <- mean(error.dma2^2, na.rm = T)

## RMSE (Root Mean Square Error)
RMSE.dma2 <- sqrt(mean(error.dma2^2, na.rm = T))

## MAD (Mean Absolute Deviation)
MAD.dma2 <- mean(abs(error.dma2), na.rm = T)

## MAPE (Mean Absolute Percentage Error)
r.error.dma2 <- (error.dma2/data.gab2[, 1])*100 # Relative Error
MAPE.dma2 <- mean(abs(r.error.dma2), na.rm = T)


akurasi2 <- data.frame(
          "Akurasi Data Training DMA" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
          "Double Moving Average n=13" = c(SSE.dma2, MSE.dma2, MAPE.dma2, RMSE.dma2, MAD.dma2))
akurasi2

selisihdma2<-data.gab2[231:271,6]-testing.ts
selisihdma2

## Time Series:
## Start = 231 
## End = 271 
## Frequency = 1 
##  [1]   54.7282737   49.7093400    7.9911874   -0.6283333   65.2531241
##  [6]   54.2341894   96.1156468   58.7969081  125.1783645  118.5598208
## [11]  114.6404962  109.8217575  120.8038009  134.8844762   58.2659326
## [16]   15.8475849   37.3286503 -123.5906744 -119.9094130   41.7718483
## [21]  141.3535007  -26.1654330 -245.5837816 -213.6027153 -361.5210630
## [26] -353.1405836 -362.3587363 -330.2780609 -123.5967996 -280.0151472
## [31]   39.1661141  -62.2522335 -130.5709722 -162.0899068 -108.9086455
## [36] -268.0271881 -353.3459268 -557.1646654 -528.9834041 -695.9027287
## [41] -843.5212724

SSEtestingdma2<-sum(selisihdma2^2)
MSEtestingdma2<-SSEtestingdma2/length(testing.ts)
MAPEtestingdma2<-sum(abs(selisihdma2/testing.ts)*100)/length(testing.ts)

akurasitestingdma2 <- matrix(c(SSEtestingdma2,MSEtestingdma2,MAPEtestingdma2),nrow=3,ncol=1)
row.names(akurasitestingdma2)<- c("SSE", "MSE", "MAPE")
colnames(akurasitestingdma2) <- "Akurasi Data Testing DMA N=8"
akurasitestingdma2

##      Akurasi Data Testing DMA N=8
## SSE                  2.909009e+06
## MSE                  7.095143e+04
## MAPE                 1.252382e+00

#Single Moving  Average dengan n=13
df_ts_sma3 <- TTR::SMA(training.ts,  n=13)
#Double Moving  Average dengan n=13
df_ts_dma3 <- TTR::SMA(df_ts_sma3, n=13)
At <- 2*df_ts_sma3-df_ts_dma3
Bt <- 2/(13-1)*(df_ts_sma3-df_ts_dma3)
pemulusan_dma3 <- At+Bt
data.ramal3<- c(NA, pemulusan_dma3)

t = 1:41
f = c()

for (i in t) {
  f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}

data.gab3 <- cbind(aktual = c(training.ts,rep(NA,41)), pemulusan3 = c(df_ts_sma3,rep(NA,41)),pemulusan3 = c(pemulusan_dma3, rep(NA,41)),At = c(At, rep(NA,41)), Bt = c(Bt,rep(NA,41)),ramalan3 = c(data.ramal3, f[-1]))
tail(data.gab3)

##        aktual pemulusan3 pemulusan3 At Bt ramalan3
## [266,]     NA         NA         NA NA NA 14312.37
## [267,]     NA         NA         NA NA NA 14313.19
## [268,]     NA         NA         NA NA NA 14314.02
## [269,]     NA         NA         NA NA NA 14314.84
## [270,]     NA         NA         NA NA NA 14315.67
## [271,]     NA         NA         NA NA NA 14316.49

#Plot time series
ts.plot(data.gab3[,1], xlab="Weekly", ylab="Close", main="Plot Pemulusan dan Peramalan dengan Metode DMA n=13")
points(data.gab3[,1])
lines(data.gab3[,3],col="green",lwd=2)
lines(data.gab3[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)

# Ukuran Keakuratan
error.dma3 <- data.gab3[, 1] - data.gab3[, 3]

## SSE (Sum Square Error)
SSE.dma3 <- sum(error.dma3^2, na.rm = T)

## MSE (Mean Squared Error)
MSE.dma3 <- mean(error.dma3^2, na.rm = T)

## RMSE (Root Mean Square Error)
RMSE.dma3 <- sqrt(mean(error.dma3^2, na.rm = T))

## MAD (Mean Absolute Deviation)
MAD.dma3 <- mean(abs(error.dma3), na.rm = T)

## MAPE (Mean Absolute Percentage Error)
r.error.dma3 <- (error.dma3/data.gab3[, 1])*100 # Relative Error
MAPE.dma3 <- mean(abs(r.error.dma3), na.rm = T)


akurasi3 <- data.frame(
          "Akurasi Data Training DMA" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
          "Double Moving Average n=3" = c(SSE.dma3, MSE.dma3, MAPE.dma3, RMSE.dma3, MAD.dma3))
akurasi3

selisihdma3<-data.gab3[231:271,6]-testing.ts
selisihdma3

## Time Series:
## Start = 231 
## End = 271 
## Frequency = 1 
##  [1]   -31.00911   -45.38424   -96.45859  -114.43431   -57.90905   -78.28419
##  [7]   -45.75893   -92.43386   -35.40860   -51.38335   -64.65887   -78.83381
## [13]   -77.20796   -72.48348  -158.45822  -210.23277  -198.10790  -368.38343
## [19]  -374.05836  -221.73330  -131.50784  -308.38298  -537.15752  -514.53265
## [25]  -671.80720  -672.78292  -691.35727  -668.63279  -471.30773  -637.08227
## [31]  -327.25721  -438.03176  -515.70669  -556.58182  -512.75676  -681.23150
## [37]  -775.90644  -989.08138  -970.25631 -1146.53183 -1303.50658

SSEtestingdma3<-sum(selisihdma3^2)
MSEtestingdma3<-SSEtestingdma3/length(testing.ts)
MAPEtestingdma3<-sum(abs(selisihdma3/testing.ts)*100)/length(testing.ts)

akurasitestingdma3 <- matrix(c(SSEtestingdma3,MSEtestingdma3,MAPEtestingdma3),nrow=3,ncol=1)
row.names(akurasitestingdma3)<- c("SSE", "MSE", "MAPE")
colnames(akurasitestingdma3) <- "Akurasi Data Testing DMA n=13"
akurasitestingdma3

##      Akurasi Data Testing DMA n=13
## SSE                   1.076415e+07
## MSE                   2.625403e+05
## MAPE                  2.605386e+00

akurasitestingdma <-matrix(c(SSEtestingdma1,MSEtestingdma1,MAPEtestingdma1,SSEtestingdma2,MSEtestingdma2,MAPEtestingdma2,SSEtestingdma3,MSEtestingdma3,MAPEtestingdma3),nrow=3,ncol=3)
row.names(akurasitestingdma)<- c("SSE", "MSE", "MAPE")
colnames(akurasitestingdma) <- c("dma ske1","dma ske2","dma ske3")
akurasitestingdma

##          dma ske1     dma ske2     dma ske3
## SSE  6.963651e+07 2.909009e+06 1.076415e+07
## MSE  1.698451e+06 7.095143e+04 2.625403e+05
## MAPE 7.473236e+00 1.252382e+00 2.605386e+00

par(mfrow=c(1,1))
plot(1:length(training.ts), training.ts, type="b", col="black", xlim=c(1,length(training.ts)+41),
     ylab="Data Nilai Tukar",xlab="Mingguan", main="Perbandingan Metode Ramalan vs Aktual\nDilihat setelah data training")
lines((length(training.ts)+1):(length(training.ts)+41), testing.ts,col="black")
lines((length(training.ts)+1):(length(training.ts)+41), data.gab1[231:271,6],col="blue")
lines((length(training.ts)+1):(length(training.ts)+41), data.gab2[231:271,6],col="green")
lines((length(training.ts)+1):(length(training.ts)+41), data.gab3[231:271,6],col="red")
legend("bottomleft",c("data aktual", "Ramalan DMA 1", "Ramalan DMA 2", "Ramalan DMA 3"), lty=1, 
       col=c("black", "blue", "green", "red"))

2. Double Exponential Smoothing (DES)

Sebelum melakukan pemulusan dan peramalan data Nilai Tukar Rupiah, terlebih dahulu perlu ditentukan nilai alpha dan beta yang menghasilkan peramalan optimum dengan nilai eror kecil.

#Pengulangan Melihat Nilai Keakuratan
a = seq(0.1,0.9,0.1)
b = seq(0.1,0.9,0.1)
output = data.frame()
for (i in a) {
  for (j in b) {
    des <- HoltWinters(training.ts, alpha = i, beta=j, gamma=F)
    sse <- des$SSE
    mse <- sse/length(training.ts)
    rmse <-sqrt(mse)
    akurasi <- cbind("SSE"=sse, "MSE"=mse, "RMSE"=rmse)
    output <- rbind(output, akurasi)
  }
}

output <- cbind("Alpha" = rep(a,each=9), "Beta"=b, output)
dplyr::arrange(.data=as.data.frame(output), MSE)

Berdasarkan perhitungan nilai eror menggunakan RMSE, SSE, MSE diperoleh nilai parameter paling optimum adalah ketika alpha=0.9 dan beta=0.1. Hasil ini juga menunjukkan bahwa semakin tinggi nilai alpha maka nilai eror semakin kecil dan semakin tinggi nilai beta maka nilai eror semakin besar.

Selanjutnya akan dilakukan analisis pemulusan secara grafik (subjektif).

a = seq(0.1, 0.9, 0.3)
b = seq(0.1, 0.9, 0.3)
output = data.frame()
for (i in a) {
  for (j in b) {
  df_des <- HoltWinters(training.ts, alpha = i, beta=j, gamma=F)
  datades <- data.frame(training.ts, c(NA, NA, df_des$fitted[,1]))
  colnames(datades) = c("y","yhat")
  
  ts.plot(training.ts,xlab="periode  waktu", ylab="Yt",  col="blue", lty=3)
  points(training.ts)
  lines (datades[,2], col="red",lwd=2)
  title(paste0("Double Exponential Smoothing (Alpha=",i ," dan Beta=", j, ")"), cex.main=1, font.main=4 ,col.main="black")
  legend("topleft", c("Data aktual","Fitted DES"), lty=1:3,col=c ("blue","red"))
  }
}

Berdasarkan grafik plot sebelumnya dengan mempertimbangkan nilai eror yang dihasilkan, pemulusan dengan alpha=0.7 dan beta=0.7 dapat dipilih sebagai parameter pemulusan DES yang baik.

Oleh karena itu, kelompok kami memilih 2 model dengan nilai alpha=0.9;beta=0.1 dan alpha=0.7;beta=0.7 yang menghasilkan plot grafik cukup mulus dan nilai eror yang kecil.

Berikut adalah hasil peramalan data training menggunakan metode DES untuk 41 periode ke depan beserta plot pemulusan yang dihasilkan.

#Gabungan Data Aktual, Pemulusan, dan Ramalan
df_des <- HoltWinters(training.ts, alpha = 0.9, beta=0.1, gamma=F)
ramal_des <- forecast::forecast(df_des,h=41)
(df_ramal_des <- data.frame(ramal_des))

data.des <- cbind(aktual=c(training.ts, rep(NA,41)),
            pemulusan=c(NA, df_des$fitted[,2], as.numeric(df_des$coefficients[1]+df_des$coefficients[2]), rep(NA,41)),
            ramalan = c(NA, NA, df_des$fitted[,1], df_ramal_des$Point.Forecast))
data.des <- ts(data.des)
tail(data.des)

## Time Series:
## Start = 266 
## End = 271 
## Frequency = 1 
##     aktual pemulusan  ramalan
## 266     NA        NA 14334.99
## 267     NA        NA 14335.77
## 268     NA        NA 14336.54
## 269     NA        NA 14337.32
## 270     NA        NA 14338.10
## 271     NA        NA 14338.87

plot(ramal_des, xlab="Mingguan", ylab="Nilai Tukar Rupiah", main="Plot Pemulusan Metode DES Ramalan 1")
lines((length(training.ts)+1):(length(training.ts)+41), testing.ts,col="black")

sse <- df_des$SSE
mse <- sse/length(training.ts)
rmse <-sqrt(mse)
akurasi.des.train <- rbind("SSE"=sse, "MSE"=mse, "RMSE"=rmse)
colnames(akurasi.des.train ) <- "Akurasi Data Training 1 DES"
akurasi.des.train 

##      Akurasi Data Training 1 DES
## SSE                 7461495.9561
## MSE                   32441.2868
## RMSE                    180.1146

selisihdes1<-ramal_des$mean-testing.ts
selisihdes1

## Time Series:
## Start = 231 
## End = 271 
## Frequency = 1 
##  [1]    -6.692601   -21.116122   -72.238862   -90.262970   -33.786100
##  [6]   -54.209622   -21.732752   -68.456078   -11.479210   -27.502341
## [11]   -40.826253   -55.049579   -53.472123   -48.796035  -134.819166
## [16]  -186.642101  -174.565623  -344.889536  -350.612862  -198.336188
## [21]  -108.159123  -285.082644  -513.905580  -491.329101  -648.652036
## [26]  -649.676144  -668.298885  -645.622797  -448.346123  -614.169058
## [31]  -304.392384  -415.215319  -492.938645  -533.862167  -490.085493
## [36]  -658.608624  -753.331950  -966.555276  -947.778602 -1124.102514
## [41] -1281.125645

SSEtestingdes1<-sum(selisihdes1^2)
MSEtestingdes1<-SSEtestingdes1/length(testing.ts)
MAPEtestingdes1<-sum(abs(selisihdes1/testing.ts)*100)/length(testing.ts)

#Gabungan Data Aktual, Pemulusan, dan Ramalan
df_des2 <- HoltWinters(training.ts, alpha = 0.7, beta=0.7, gamma=F)
ramal_des2 <- forecast::forecast(df_des2,h=41)
(df_ramal_des2 <- data.frame(ramal_des2))

data.des2 <- cbind(aktual=c(training.ts, rep(NA,41)),
            pemulusan=c(NA, df_des2$fitted[,2], as.numeric(df_des2$coefficients[1]+df_des2$coefficients[2]), rep(NA,41)),
            ramalan = c(NA, NA, df_des2$fitted[,1], df_ramal_des2$Point.Forecast))
data.des2 <- ts(data.des2)
tail(data.des2)

## Time Series:
## Start = 266 
## End = 271 
## Frequency = 1 
##     aktual pemulusan  ramalan
## 266     NA        NA 15185.02
## 267     NA        NA 15210.77
## 268     NA        NA 15236.52
## 269     NA        NA 15262.27
## 270     NA        NA 15288.02
## 271     NA        NA 15313.77

plot(ramal_des2, xlab="Mingguan", ylab="Nilai Tukar Rupiah", main="Plot Pemulusan Metode DES Ramalan 2")
lines((length(training.ts)+1):(length(training.ts)+41), testing.ts,col="black")

sse2 <- df_des2$SSE
mse2 <- sse/length(training.ts)
rmse2 <-sqrt(mse2)
akurasi.des.train2 <- rbind("SSE"=sse2, "MSE"=mse2, "RMSE"=rmse2)
colnames(akurasi.des.train2) <- "Akurasi Data Training 2 DES"
akurasi.des.train2

##      Akurasi Data Training 2 DES
## SSE                 8788440.9269
## MSE                   32441.2868
## RMSE                    180.1146

selisihdes2<-ramal_des2$mean-testing.ts
selisihdes2

## Time Series:
## Start = 231 
## End = 271 
## Frequency = 1 
##  [1]  -30.785258  -20.233964  -46.381888  -39.431181   42.020504   46.571797
##  [7]  104.023483   82.274972  164.226656  173.178340  184.829244  195.580733
## [13]  222.133004  251.783908  190.735592  163.887472  200.938765   55.589669
## [19]   74.841158  252.092647  367.244527  215.295822   11.447701   58.998995
## [25]  -73.349124  -49.398417  -43.046342    4.604561  226.856051   86.007931
## [31]  420.759420  334.911300  282.162790  266.214083  334.965572  191.417257
## [37]  121.668747  -66.579764  -22.828275 -174.177371 -306.225687

SSEtestingdes2<-sum(selisihdes2^2)
MSEtestingdes2<-SSEtestingdes2/length(testing.ts)
MAPEtestingdes2<-sum(abs(selisihdes2/testing.ts)*100)/length(testing.ts)

des.opt<- HoltWinters(training.ts, gamma = FALSE)
des.opt

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0.003852476
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 14318.500000
## b    -4.977974

#Gabungan Data Aktual, Pemulusan, dan Ramalan
df_des3 <- HoltWinters(training.ts, alpha = 1, beta=0.00385, gamma=F)
ramal_des3 <- forecast::forecast(df_des3,h=41)
(df_ramal_des3 <- data.frame(ramal_des3))

data.des3 <- cbind(aktual=c(training.ts, rep(NA,41)),
            pemulusan=c(NA, df_des3$fitted[,2], as.numeric(df_des3$coefficients[1]+df_des3$coefficients[2]), rep(NA,41)),
            ramalan = c(NA, NA, df_des3$fitted[,1], df_ramal_des3$Point.Forecast))
data.des3 <- ts(data.des3)
tail(data.des3)

## Time Series:
## Start = 266 
## End = 271 
## Frequency = 1 
##     aktual pemulusan  ramalan
## 266     NA        NA 14139.15
## 267     NA        NA 14134.17
## 268     NA        NA 14129.19
## 269     NA        NA 14124.21
## 270     NA        NA 14119.22
## 271     NA        NA 14114.24

plot(ramal_des3, xlab="Mingguan", ylab="Nilai Tukar Rupiah", main="Plot Pemulusan Metode DES Ramalan 3")
lines((length(training.ts)+1):(length(training.ts)+41), testing.ts,col="black")

ssedes.train.opt<-des.opt$SSE
msedes.train.opt<-ssedes.train.opt/length(training.ts)
sisaandes.opt<-ramal_des3$residuals
head(sisaandes.opt)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1]         NA         NA    8.00000 -118.03080   71.42362   70.14864

mapedes.train.opt<- sum(abs(sisaandes.opt[3:length(training.ts)]/training.ts[3:length(training.ts)])
*100)/length(training.ts)
akurasides.opt <- matrix(c(ssedes.train.opt,msedes.train.opt,mapedes.train.opt))
row.names(akurasides.opt)<- c("SSE", "MSE", "MAPE")
colnames(akurasides.opt) <- c("Akurasi alpha=0.1 dan beta=0.2")
akurasides.opt

##      Akurasi alpha=0.1 dan beta=0.2
## SSE                    6.931466e+06
## MSE                    3.013681e+04
## MAPE                   7.809192e-01

selisihdes.opt<-ramal_des3$mean-testing.ts
selisihdes.opt

## Time Series:
## Start = 231 
## End = 271 
## Frequency = 1 
##  [1]    -0.9818981   -21.1639911   -78.0453032  -101.8279832   -51.1096853
##  [6]   -77.2917793   -50.5734814  -103.0553794   -51.8370825   -73.6187855
## [11]   -92.7012696  -112.6831676  -116.8642837  -117.9467677  -209.7284708
## [16]  -267.3099778  -260.9920719  -437.0745559  -448.5564540  -302.0383520
## [21]  -217.6198591  -400.3019521  -634.8834602  -618.0655532  -781.1470603
## [26]  -787.9297403  -812.3110524  -795.3935364  -603.8754345  -775.4569415
## [31]  -471.4388396  -588.0203466  -671.5022447  -718.1843387  -680.1662368
## [36]  -854.4479389  -954.9298369 -1173.9117350 -1160.8936330 -1342.9761171
## [41] -1505.7578201

SSEtestingdes.opt<-sum(selisihdes.opt^2)
MSEtestingdes.opt<-SSEtestingdes.opt/length(testing.ts)
MAPEtestingdes.opt<-sum(abs(selisihdes.opt/testing.ts)*100)/length(testing.ts)

akurasitestingdes <-matrix(c(SSEtestingdes1,MSEtestingdes1,MAPEtestingdes1,SSEtestingdes2,MSEtestingdes2,MAPEtestingdes2,SSEtestingdes.opt,MSEtestingdes.opt,MAPEtestingdes.opt),nrow=3,ncol=3)
row.names(akurasitestingdes)<- c("SSE", "MSE", "MAPE")
colnames(akurasitestingdes) <- c("des ske1","des ske2","des ske3")
akurasitestingdes

##          des ske1     des ske2     des ske3
## SSE  1.005351e+07 1.430348e+06 1.568083e+07
## MSE  2.452076e+05 3.488655e+04 3.824592e+05
## MAPE 2.446280e+00 1.029864e+00 3.182277e+00

par(mfrow=c(1,1))
plot(1:length(training.ts), training.ts, type="b", col="black", xlim=c(1,length(training.ts)+41),
     ylab="Data Nilai Tukar",xlab="Mingguan", main="Perbandingan Metode Ramalan vs Aktual\nDilihat setelah data training")
lines((length(training.ts)+1):(length(training.ts)+41), testing.ts,col="black")
lines((length(training.ts)+1):(length(training.ts)+41), ramal_des$mean,col="blue")
lines((length(training.ts)+1):(length(training.ts)+41), ramal_des2$mean,col="green")
lines((length(training.ts)+1):(length(training.ts)+41), ramal_des3$mean,col="red")
legend("bottomleft",c("data aktual","Ramalan DES 1","Ramalan DES 2","Ramalan DES 3"), lty=1, 
       col=c("black", "blue","green", "red"))

Kebaikan Model

Berikut adalah plot perbandingan data testing aktual, data peramalan DMA paling optimum, dan data peramalan DES paling optimum:

par(mfrow=c(1,1))
plot(1:length(training.ts), training.ts, type="b", col="black", xlim=c(1,length(training.ts)+41),
     ylab="Data Nilai Tukar",xlab="Mingguan", main="Perbandingan Metode Ramalan vs Aktual\nDilihat setelah data training")
lines((length(training.ts)+1):(length(training.ts)+41), testing.ts,col="black")
lines((length(training.ts)+1):(length(training.ts)+41), data.gab2[231:271,6],col="blue")
lines((length(training.ts)+1):(length(training.ts)+41), ramal_des2$mean,col="red")
legend("bottomleft",c("data aktual", "Ramalan DMA", "Ramalan DES"), lty=1, 
       col=c("black", "blue", "red"))

Berdasarkan perbandingan grafik di atas, dapat terlihat bahwa kedua garis ramalan tidak jauh berbeda dengan garis aktualnya. Hal ini mengindikasikan bahwa pemodelan data training telah cukup baik untuk digunakan dalam peramalan data testing. Selanjutnya akan dicari nilai eror kedua metode untuk menentukan metode mana yang terbaik1

MSEfull <-matrix(c(SSEtestingdma2,MSEtestingdma2,MAPEtestingdma2,SSEtestingdes2,MSEtestingdes2,MAPEtestingdes2),nrow=3,ncol=2)
row.names(MSEfull)<- c("SSE","MSE","MAPE")
colnames(MSEfull) <- c("DMA","DES")
MSEfull

##               DMA          DES
## SSE  2.909009e+06 1.430348e+06
## MSE  7.095143e+04 3.488655e+04
## MAPE 1.252382e+00 1.029864e+00

Berdasarkan perhitungan nilai eror di atas, telah didapat model pemulusan dan peramalan terbaik adalah metode DES ramalan 2 dengan alpha=0.7 dan beta=0.7.

Differencing

Model yang tidak stasioner dalam rataan memerlukan differencing (pembedaan) untuk menstasionerkan datanya terlebih dahulu. Setelah dilakukan differencing, kita akan memperoleh nilai d untuk model Arima.

Differencing dilakukan untuk menstasionerkan data nilai tukar rupiah yang telah terindentifikasi tak stasioner pada tahap sebelumnya. Data time series \(y_t\) sebagai homogeneous nonstationary tidak stasioner, tetapi dengan pembedaan pertama (first difference) yaitu \(wt=y_t−y_{t−1}=(1−B)y_t\) atau pembedaan order yang lebih tinggi, \(w_t=(1−B)^dy_t\), menghasilkan data time series yang stasioner.

Maka selanjutnya akan dilakukan differencing pertama pada data nilai tukar rupiah.

train.diff<-diff(training.ts, differences = 1) 

plot.ts(train.diff, lty=1, xlab="Mingguan", ylab="Data Difference 1 Nilai Tukar", main="Plot Difference Nilai Tukar")
points(train.diff)

acf(train.diff, lag.max=20, main="data difference nilai tukar Indonesia")

adf.test(train.diff) #stasioner

## Warning in adf.test(train.diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff
## Dickey-Fuller = -5.9883, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

#Plot Box-Cox
b <- boxcox(data$Nilai_tukar[train]~data$Date[train], lambda = seq(2,-10))

#Nilai Rounded Lambda
lambda <- b$x[which.max(b$y)] 
lambda

## [1] -6.484848

A. Penentuan Ordo Model Arima

Penentuan ordo model ARIMA dapat dilakukan dengan menampilkan plot ACF, PACF, dan EACF sebagai berikut:

#Identifikasi Model ARIMA dengan ACF, PACF, dan EACF
acf(train.diff, lag.max=20, main="ACF data nilai tukar Indonesia")

pacf(train.diff, lag.max=20, main="PACF data nilai tukar Indonesia")

eacf(train.diff)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o x o o x o x  o  o  o 
## 1 x o o o o x o o x o x  o  o  o 
## 2 x o o o o x o o o o x  o  o  o 
## 3 x x x o o x x o x o o  o  o  o 
## 4 x x o x o x x o x o x  o  o  o 
## 5 x x o x o x o o o x o  o  o  o 
## 6 x o x x x o o o o o o  o  o  o 
## 7 x o o x x o x o o o o  o  o  o

1. Plot ACF dan PACF di atas tidak menunjukkan adanya pola cuts off ataupun tails off, sehingga model yang dapat digunakan hanyalah Arima(0,1,0).
2. Berdasarkan plot EACF di atas, kita dapat menggunakan orde 0,1;0,2;1,1;1,2;2,2 untuk parameter p,q pada model.
3. Apabila lag 0 diabaikan, terlihat dari plot ACF adanya cuts off after lag ke 6 dan terdapat pula cuts off after lag ke 4 pada plot PACF , sehingga model lain yang dapat digunakan adalah Arima(4,1,0) dan Arima(0,1,6).

B. Pemilihan Model Terbaik

Spesifikasi model yang dapat digunakan adalah Arima(0,1,0), Arima(0,1,1), Arima(0,1,2), Arima(1,1,1), Arima(1,1,2), Arima(2,1,2), Arima(4,1,0), dan Arima(0,1,6). Selanjutnya akan dipilih model terbaik berdasarkan nilai AIC.

#Melihat summary dan signifikansi seluruh pemodelan
model1 <- Arima(training.ts, order=c(0,1,0), method="ML")
summary(model1)

## Series: training.ts 
## ARIMA(0,1,0) 
## 
## sigma^2 = 29970:  log likelihood = -1505.2
## AIC=3012.39   AICc=3012.41   BIC=3015.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE      MAPE      MASE      ACF1
## Training set 4.247196 172.7426 111.5669 0.0238215 0.7729705 0.9961706 0.0532247

model2 <- Arima(training.ts, order=c(0,1,1), method="ML")
summary(model2)

## Series: training.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.0480
## s.e.  0.0627
## 
## sigma^2 = 30024:  log likelihood = -1504.9
## AIC=3013.8   AICc=3013.86   BIC=3020.67
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE      MAPE     MASE
## Training set 4.075828 172.5199 112.0505 0.02291826 0.7766009 1.000489
##                     ACF1
## Training set 0.002473135

lmtest::coeftest(model2)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ma1 0.047954   0.062688   0.765   0.4443

model3 <- Arima(training.ts, order=c(0,1,2), method="ML")
summary(model3)

## Series: training.ts 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1     ma2
##       0.0376  0.0753
## s.e.  0.0663  0.0716
## 
## sigma^2 = 30019:  log likelihood = -1504.38
## AIC=3014.77   AICc=3014.87   BIC=3025.07
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE      MAPE     MASE
## Training set 3.882657 172.1263 112.3354 0.02197603 0.7787627 1.003033
##                     ACF1
## Training set 0.006203937

lmtest::coeftest(model3)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ma1 0.037636   0.066259  0.5680   0.5700
## ma2 0.075312   0.071606  1.0517   0.2929

model4 <- Arima(training.ts, order=c(1,1,1), method="ML")
summary(model4)

## Series: training.ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.4908  -0.4240
## s.e.  0.3191   0.3275
## 
## sigma^2 = 30056:  log likelihood = -1504.52
## AIC=3015.04   AICc=3015.15   BIC=3025.34
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE      MAPE     MASE
## Training set 3.801792 172.2313 112.1421 0.02151112 0.7774392 1.001307
##                     ACF1
## Training set -0.01854928

lmtest::coeftest(model4)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1  0.49081    0.31909  1.5382   0.1240
## ma1 -0.42400    0.32755 -1.2945   0.1955

model5 <- Arima(training.ts, order=c(1,1,2), method="ML")
summary(model5)

## Series: training.ts 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.3274  -0.2851  0.0777
## s.e.  0.3365   0.3340  0.0763
## 
## sigma^2 = 30071:  log likelihood = -1504.08
## AIC=3016.16   AICc=3016.34   BIC=3029.89
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE      MAPE     MASE
## Training set 3.678911 171.8953 112.3112 0.02094799 0.7786863 1.002817
##                     ACF1
## Training set 0.002814549

lmtest::coeftest(model5)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.327387   0.336523  0.9728   0.3306
## ma1 -0.285100   0.333972 -0.8537   0.3933
## ma2  0.077734   0.076323  1.0185   0.3084

model6 <- Arima(training.ts, order=c(2,1,2), method="ML")
summary(model6)

## Series: training.ts 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       0.6488  -0.7414  -0.6503  0.8662
## s.e.  0.1842   0.0968   0.1386  0.0711
## 
## sigma^2 = 29336:  log likelihood = -1500.88
## AIC=3011.75   AICc=3012.02   BIC=3028.92
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 4.059663 169.4051 111.9665 0.02338373 0.776337 0.9997386
##                    ACF1
## Training set 0.03617494

lmtest::coeftest(model6)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.648805   0.184170  3.5229 0.0004269 ***
## ar2 -0.741404   0.096777 -7.6609 1.846e-14 ***
## ma1 -0.650268   0.138594 -4.6919 2.707e-06 ***
## ma2  0.866249   0.071123 12.1795 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model7 <- Arima(training.ts, order=c(4,1,0), method="ML")
summary(model7)

## Series: training.ts 
## ARIMA(4,1,0) 
## 
## Coefficients:
##          ar1     ar2     ar3      ar4
##       0.0559  0.0687  0.0901  -0.1355
## s.e.  0.0653  0.0651  0.0652   0.0652
## 
## sigma^2 = 29506:  log likelihood = -1501.44
## AIC=3012.89   AICc=3013.16   BIC=3030.06
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE     MASE
## Training set 3.970938 169.8968 112.0075 0.02256611 0.776964 1.000105
##                     ACF1
## Training set 0.009643107

lmtest::coeftest(model7)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ar1  0.055904   0.065342  0.8556  0.39224  
## ar2  0.068718   0.065135  1.0550  0.29142  
## ar3  0.090075   0.065184  1.3819  0.16701  
## ar4 -0.135497   0.065240 -2.0769  0.03781 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model8 <- Arima(training.ts, order=c(0,1,6), method="ML")
summary(model8)

## Series: training.ts 
## ARIMA(0,1,6) 
## 
## Coefficients:
##          ma1     ma2     ma3      ma4      ma5      ma6
##       0.1200  0.0351  0.0660  -0.2025  -0.0617  -0.3265
## s.e.  0.0634  0.0682  0.0678   0.0667   0.0805   0.0728
## 
## sigma^2 = 26828:  log likelihood = -1489.94
## AIC=2993.89   AICc=2994.39   BIC=3017.92
## 
## Training set error measures:
##                    ME    RMSE      MAE        MPE      MAPE     MASE
## Training set 6.766619 161.281 114.5097 0.04002731 0.7949972 1.022447
##                     ACF1
## Training set -0.01607387

lmtest::coeftest(model8)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1  0.119984   0.063413  1.8921  0.058475 .  
## ma2  0.035112   0.068196  0.5149  0.606641    
## ma3  0.065984   0.067840  0.9726  0.330733    
## ma4 -0.202518   0.066718 -3.0354  0.002402 ** 
## ma5 -0.061681   0.080518 -0.7661  0.443645    
## ma6 -0.326489   0.072773 -4.4864 7.245e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan output di atas, model ARIMA(2,1,2) yang menghasilkan seluruh parameter signifikan, sementara model ARIMA(4,1,0) dan model ARIMA(0,1,6) menghasilkan beberapa parameter signifikan, serta pada model lainnya tidak ada parameter yang signifikan.

#Melihat summary dan signifikansi seluruh pemodelan
aic_model <- data.frame(
          "Nama Model" = c("Model 1", "Model 2","Model 3", "Model 4","Model 5", "Model 6","Model 7", "Model 8"), 
          "AIC" = c(model1$aic,model2$aic,model3$aic,model4$aic,model5$aic,model6$aic,model7$aic,model8$aic))
aic_model

dplyr::arrange(.data=aic_model, AIC)

Berdasarkan hasil di atas, model 8 Arima(0,1,6) memiliki nilai AIC terkecil dari seluruh model yang diujikan, kemudian dilanjut oleh model 6.

model9 <- auto.arima(training.ts)
model9

## Series: training.ts 
## ARIMA(0,1,0) 
## 
## sigma^2 = 29970:  log likelihood = -1505.2
## AIC=3012.39   AICc=3012.41   BIC=3015.83

Dengan menggunakan fungsi auto.ARIMA untuk mendapatkan parameter optimum pada R, didapat bahwa model ARIMA yang optimum adalah Arima(0,1,0) yang juga merupakan model 1.

Model tentatif dipilih berdasarkan nilai AIC minimum dengan pertimbangan bahwa seluruh parameter signifikan. Oleh karena itu, model 6 Arima(2,1,2) yang menghasilkan seluruh parameter signifikan serta nilai AIC yang cukup minimum akan dianalisis lebih lanjut sebagai model tentatif serta akan digunakan untuk melakukan peramalan terhadap data testing.

C. Overfitting Model 6 Arima(2,1,2)

model6a <- Arima(training.ts, order=c(1,1,2), method="ML")
summary(model6a)

## Series: training.ts 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.3274  -0.2851  0.0777
## s.e.  0.3365   0.3340  0.0763
## 
## sigma^2 = 30071:  log likelihood = -1504.08
## AIC=3016.16   AICc=3016.34   BIC=3029.89
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE      MAPE     MASE
## Training set 3.678911 171.8953 112.3112 0.02094799 0.7786863 1.002817
##                     ACF1
## Training set 0.002814549

lmtest::coeftest(model6a)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.327387   0.336523  0.9728   0.3306
## ma1 -0.285100   0.333972 -0.8537   0.3933
## ma2  0.077734   0.076323  1.0185   0.3084

model6b <- Arima(training.ts, order=c(3,1,2), method="ML")
summary(model6b)

## Series: training.ts 
## ARIMA(3,1,2) 
## 
## Coefficients:
##           ar1      ar2     ar3     ma1    ma2
##       -1.1747  -0.2649  0.2017  1.2970  0.429
## s.e.   0.1825   0.1948  0.0735  0.1781  0.185
## 
## sigma^2 = 27842:  log likelihood = -1494.48
## AIC=3000.96   AICc=3001.34   BIC=3021.56
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE      MAPE    MASE
## Training set 3.425523 164.6682 112.6778 0.01983489 0.7829424 1.00609
##                     ACF1
## Training set 0.003247792

lmtest::coeftest(model6b)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -1.174743   0.182512 -6.4365 1.222e-10 ***
## ar2 -0.264883   0.194803 -1.3597  0.173910    
## ar3  0.201723   0.073493  2.7448  0.006055 ** 
## ma1  1.296952   0.178147  7.2802 3.333e-13 ***
## ma2  0.429041   0.185031  2.3188  0.020408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model6c <- Arima(training.ts, order=c(2,1,1), method="ML")
summary(model6c)

## Series: training.ts 
## ARIMA(2,1,1) 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.2850  0.0597  -0.2366
## s.e.  0.3793  0.0711   0.3755
## 
## sigma^2 = 30110:  log likelihood = -1504.22
## AIC=3016.45   AICc=3016.63   BIC=3030.18
## 
## Training set error measures:
##                    ME    RMSE      MAE        MPE      MAPE     MASE
## Training set 3.710449 172.006 112.2243 0.02108375 0.7780404 1.002041
##                      ACF1
## Training set -0.002297497

lmtest::coeftest(model6c)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.284957   0.379348  0.7512   0.4525
## ar2  0.059725   0.071064  0.8404   0.4007
## ma1 -0.236583   0.375504 -0.6300   0.5287

model6d <- Arima(training.ts, order=c(2,1,3), method="ML")
summary(model6d)

## Series: training.ts 
## ARIMA(2,1,3) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2     ma3
##       -1.4611  -0.7437  1.6043  0.9782  0.2369
## s.e.   0.1178   0.1510  0.1347  0.2236  0.0793
## 
## sigma^2 = 27709:  log likelihood = -1493.95
## AIC=2999.9   AICc=3000.28   BIC=3020.5
## 
## Training set error measures:
##                   ME    RMSE      MAE        MPE      MAPE     MASE        ACF1
## Training set 3.53787 164.275 113.3376 0.02050304 0.7876653 1.011981 -0.01183227

lmtest::coeftest(model6d)

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -1.461062   0.117793 -12.4036 < 2.2e-16 ***
## ar2 -0.743698   0.150983  -4.9257 8.405e-07 ***
## ma1  1.604327   0.134697  11.9106 < 2.2e-16 ***
## ma2  0.978231   0.223646   4.3740 1.220e-05 ***
## ma3  0.236869   0.079273   2.9880  0.002808 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan output di atas, hanya model 6b Arima(3,1,2) dan 6d Arima(2,1,3) yang menghasilkan parameter signifikan, sementara pada model lainnya tidak ada parameter yang signifikan.

#Melihat summary dan signifikansi pemodelan overfitting
aic_model_ov <- data.frame(
          "Nama Model" = c("Model 6", "Model 6a","Model 6b", "Model 6c","Model 6d"), 
          "AIC" = c(model6$aic,model6a$aic,model6b$aic,model6c$aic,model6d$aic))
aic_model_ov

dplyr::arrange(.data=aic_model_ov, AIC)

#Model 6d
lmtest::coeftest(model6d)

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -1.461062   0.117793 -12.4036 < 2.2e-16 ***
## ar2 -0.743698   0.150983  -4.9257 8.405e-07 ***
## ma1  1.604327   0.134697  11.9106 < 2.2e-16 ***
## ma2  0.978231   0.223646   4.3740 1.220e-05 ***
## ma3  0.236869   0.079273   2.9880  0.002808 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model6d$aic

## [1] 2999.9

D. Analisis Sisaan Model Tentatif 6d Arima(2,1,3)

sisaan <- model6d$residuals
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(sisaan, type="o", 
     ylab = "Sisaan", xlab = "Order", main = "Sisaan M6d vs Order")
abline(h = 0, col='red')
acf(sisaan)
pacf(sisaan)

ks.test(sisaan,"pnorm") 

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.4826, p-value < 2.2e-16
## alternative hypothesis: two-sided

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

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 0.032622, df = 1, p-value = 0.8567

t.test(sisaan, mu = 0, conf.level = 0.95) 

## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.32598, df = 229, p-value = 0.7447
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -17.84680  24.92254
## sample estimates:
## mean of x 
##   3.53787

E. Identifikasi Efek ARCH

Identifikasi adanya proses ARCH pada sisaan model ARIMA(2,1,3) dapat dilakukan dengan uji Lagrange Multiplier (LM). Hipotesis yang digunakan adalah sebagai berikut:

H0 : Tidak terjadi heteroskedastisitas

H1 : Terjadi Heteroskedastisitas

Hasil uji LM test yang diperoleh sebagai berikut:

model6d <- arima(training.ts, order=c(2,1,3), method="ML")
aTSA::arch.test(model6d)

## ARCH heteroscedasticity test for residuals 
## alternative: heteroscedastic 
## 
## Portmanteau-Q test: 
##      order   PQ  p.value
## [1,]     4 58.9 5.00e-12
## [2,]     8 91.2 2.22e-16
## [3,]    12 96.0 3.44e-15
## [4,]    16 96.3 1.67e-13
## [5,]    20 96.5 5.32e-12
## [6,]    24 97.6 7.68e-11
## Lagrange-Multiplier test: 
##      order     LM  p.value
## [1,]     4 102.50 0.00e+00
## [2,]     8  44.37 1.81e-07
## [3,]    12  23.82 1.35e-02
## [4,]    16  15.71 4.02e-01
## [5,]    20  11.65 9.00e-01
## [6,]    24   9.09 9.96e-01

for (i in 1:15) {
  ArchTest <- ArchTest(model6d$residuals, lags=i, demean=TRUE)
  cat("P Value LM Test lag ke", i,"adalah" , ArchTest$p.value, "\n") }

## P Value LM Test lag ke 1 adalah 0.04837992 
## P Value LM Test lag ke 2 adalah 0.1116185 
## P Value LM Test lag ke 3 adalah 8.253398e-12 
## P Value LM Test lag ke 4 adalah 3.772793e-11 
## P Value LM Test lag ke 5 adalah 1.551907e-10 
## P Value LM Test lag ke 6 adalah 6.066259e-11 
## P Value LM Test lag ke 7 adalah 2.224219e-10 
## P Value LM Test lag ke 8 adalah 7.867362e-10 
## P Value LM Test lag ke 9 adalah 1.321873e-10 
## P Value LM Test lag ke 10 adalah 4.192927e-10 
## P Value LM Test lag ke 11 adalah 4.041537e-10 
## P Value LM Test lag ke 12 adalah 7.107633e-10 
## P Value LM Test lag ke 13 adalah 1.953893e-09 
## P Value LM Test lag ke 14 adalah 2.819025e-09 
## P Value LM Test lag ke 15 adalah 3.41956e-09

Berdasarkan hasil uji diperoleh nilai \(P-Value\) dari lag ke-3 sampai lag ke-15 yang didapat lebih kecil dari α = 5%, maka keputusan yang diambil adalah Tolak H0. Dengan kata lain, dapat disimpulkan bahwa terdapat unsur heteroskedastisitas pada sisaan model ARIMA(2,1,3). Kondisi yang menolak H0 ini terjadi sampai dengan lag ke-15 sehingga dapat dikatakan bahwa indikasi pemodelan sisaan ini lebih cocok menggunakan model GARCH dibandingkan model ARCH. Hal tersebut bersesuaian dengan pernyataan Desvina & Meijer (2018) yang menyebutkan bahwa ketika pemeriksaan sisaan yang dilakukan dari lag 1 sampai 12 signifikan, maka model ARCH lebih cocok digunakan. Akan tetapi, jika pemeriksaan sisaan dilakukan hingga lebih dari lag 12, maka model GARCH lebih cocok digunakan. Oleh karena itu, selanjutnya akan dilakukan pendugaan model dengan menggunakan GARCH.

F. Pendugaan Model GARCH

Pendugaan model GARCH untuk mendapat model terbaik ditentukan dengan cmempertimbangkan nilai \(p-value\) ARCH-LM, nilai AIC terkecil, dan signifikansi parameternya.

Hasil trial and error menghasilkan model ARIMA(2,1,3)-GARCH(3,5) dengan \(p-value\) uji ARCH-LM sebesar 0,976, artinya efek ARCH sudah teratasi. Selanjutnya akan dilakukan overfitting pada model GARCH tersebut yang hasilnya ditunjukkan pada tabel di bawah:

Berdasarkan tabel di atas, didapat bahwa model ARIMA(2,1,3)-GARCH(3,6) hasil overfitting telah menghasilkan nilai AIC paling minimum sehingga dipilih menjadi model terbaik. Selanjutnya, sebelum melakukan peramalan, dilakukan uji diagnostik kembali pada model ARIMA(2,1,3)-GARCH(3,6) dengan uji Shapiro.Wilk dan uji LM ARCH sebagai berikut:

Tabel di atas menunjukkan hasil uji diagnostik model ARIMA(2,1,3) GARCH(3,6) untuk normalitas dan heterokesdatisitas. Berdasarkan nilai \(P-value\) yang didapatkan dari kedua uji tersebut, model sudah memenuhi asumsi kenormalan dan kehomogenan ragam karena memiliki p-value yang lebih besar daripada \(α = 1%\). Hasil tersebut menunjukkan bahwa model GARCH sudah mengatasi heteroskesdastisitas dan ketidaknormalan sisaan hasil uji formal pada model ARIMA tanpa GARCH.

G. Peramalan

Peramalan dilakukan dengan menggunakan model terbaik yang didapatkan dari pendugaan model GARCH. Hasil ramalan nilai tukar rupiah periode mingguan menggunakan model ARIMA (2,1,3)-GARCH (3,6) selama 41 periode ke depan ditunjukkan pada plot di bawah:

# Peramalan
garchSpec <- ugarchspec(
  variance.model=list(model="sGARCH",
                      garchOrder=c(3,6)),
  mean.model=list(armaOrder=c(2,3)),
  distribution.model="std")
garchFitt <- ugarchfit(spec=garchSpec, data=train.diff)
forc<- ugarchforecast(fitORspec = garchFitt, data = data$Nilai_tukar, n.ahead = 41, n.roll = 0)
plot(forc, which= 1)

Hasil forecasting di atas masih perlu ditransformasi balik karena adanya prosess differencing. Hasil transformasi balik serta perbandingan data ramalan dengan data aktual dapat dilihat di pada tabel berikut.

pt_1 <- data$Nilai_tukar[231] #nilai akhir data latih
hasil.forc.Diff <- forc@forecast$seriesFor[,1]
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1
perbandingan <- data.frame("Aktual"= testing.ts,
                           "Ramalan" = hasil[-1])
head(perbandingan,10)

Date<-c(data$Date[231:271]) 
dataframe <- data.frame(Date, perbandingan) 
plot(dataframe$Ramalan, type="b", col="black")

T <- nrow(dataframe) 
MAPE <- 1/T*sum(abs((dataframe$Aktual-dataframe$Ramalan)/dataframe$Aktual)*100)
MAPE

## [1] 1.987857

Menurut Marikar (2019), nilai MAPE yang kurang dari 10% dikategorikan sebagai hasil peramalan yang sangat baik. Hasil akurasi peramalan data nilai tukar rupiah menggunakan nilai MAPE menunjukkan angka 1.99%. Hal ini berarti peramalan data nilai tukar rupiah terhadap dolar AS dengan menggunakan model ARIMA(2,1,3)-GARCH(3,1) sudah tepat digunakan karena cukup mendekati data aktual.