library(tseries)
## Warning: package 'tseries' was built under R version 4.3.2
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(bayesforecast)
## Warning: package 'bayesforecast' was built under R version 4.3.3
## Registered S3 methods overwritten by 'bayesforecast':
## method from
## autoplot.ts forecast
## forecast.ts forecast
## fortify.ts forecast
## print.garch tseries
##
## Attaching package: 'bayesforecast'
## The following object is masked from 'package:tseries':
##
## garch
## The following objects are masked from 'package:base':
##
## beta, gamma
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
library(knitr)
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.2
## Warning: package 'lubridate' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some() masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library("forecast")
## Warning: package 'forecast' was built under R version 4.3.3
##
## Attaching package: 'forecast'
##
## The following objects are masked from 'package:bayesforecast':
##
## fourier, naive
library("TTR")
## Warning: package 'TTR' was built under R version 4.3.2
library("TSA")
## Warning: package 'TSA' was built under R version 4.3.3
## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
##
## The following object is masked from 'package:readr':
##
## spec
##
## The following objects are masked from 'package:stats':
##
## acf, arima
##
## The following object is masked from 'package:utils':
##
## tar
library("graphics")
library("dplyr")
library(MASS)
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
library(bayesplot)
## Warning: package 'bayesplot' was built under R version 4.3.3
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
data1 <- read.csv("D:/Sem 6/Bayes/Project UAS/data bawang merah (2).csv")
data1 <- data1[350:494,]
data1$Tanggal <- as.Date(data1$Tanggal, format = "%d-%b-%y")
head(data1)
## Tanggal Bawang.Merah
## 350 2024-12-15 59833
## 351 2024-12-16 60042
## 352 2024-12-17 61292
## 353 2024-12-18 61917
## 354 2024-12-19 61917
## 355 2024-12-20 61917
str(data1)
## 'data.frame': 145 obs. of 2 variables:
## $ Tanggal : Date, format: "2024-12-15" "2024-12-16" ...
## $ Bawang.Merah: int 59833 60042 61292 61917 61917 61917 62542 62542 62917 NA ...
summary(data1)
## Tanggal Bawang.Merah
## Min. :2024-12-15 Min. :48750
## 1st Qu.:2025-01-20 1st Qu.:54444
## Median :2025-02-25 Median :55870
## Mean :2025-02-25 Mean :56138
## 3rd Qu.:2025-04-02 3rd Qu.:57381
## Max. :2025-05-08 Max. :62917
## NA's :12
set.seed(021)
# Menggunakan interpolasi linier untuk mengisi nilai NA
data1$Bawang.Merah <- na.approx(data1$Bawang.Merah)
# Mengecek apakah masih ada NA
sum(is.na(data1$Bawang.Merah))
## [1] 0
data1.ts <-ts(data1$Bawang.Merah)
plot(data1$Tanggal, data1$Bawang.Merah, type="l", xlab="Periode", ylab="harga bawang merah papua", main="Plot Data harga bawang merah papua", xaxt="n")
axis(1, at=seq(min(data1$Tanggal), max(data1$Tanggal), by="month"), labels=format(seq(min(data1$Tanggal), max(data1$Tanggal), by="month"), "%b %Y"), las=1, cex.axis=0.8)
train<-data1.ts[1:116]
train.ts<-ts(train)
plot(data1$Tanggal[1:116], train.ts, type="l", xlab="Periode", ylab="harga bawang merah papua", main="Plot harga bawang merah papua Train", xaxt="n")
axis(1, at=seq(min(data1$Tanggal[1:116]), max(data1$Tanggal[1:116]), by="month"), labels=format(seq(min(data1$Tanggal[1:116]), max(data1$Tanggal[1:116]), by="month"), "%b %Y"), las=1, cex.axis=0.8)
test<-data1.ts[117:145]
test.ts<-ts(test)
plot(data1$Tanggal[117:145], test.ts, type="l", xlab="Periode", ylab="harga bawang merah papua", main="Plot harga bawang merah papua Test", xaxt="n")
axis(1, at=seq(min(data1$Tanggal[117:145]), max(data1$Tanggal[117:145]), by="month"), labels=format(seq(min(data1$Tanggal[117:145]), max(data1$Tanggal[117:145]), by="month"), "%b %Y"), las=1, cex.axis=0.8)
acf(train.ts)
tseries::adf.test(train.ts)
##
## Augmented Dickey-Fuller Test
##
## data: train.ts
## Dickey-Fuller = -2.6863, Lag order = 4, p-value = 0.2917
## alternative hypothesis: stationary
index <- seq(1:116)
bc <- boxcox(train.ts~index, lambda = seq(-7,7,by=0.1))
train.diff<-diff(train.ts,differences = 1)
plot(data1$Tanggal[2:116], train.diff, type="l", xlab="Periode", ylab="Data Difference 1", main="Plot Difference harga bawang merah papua", xaxt="n")
axis(1, at=seq(min(data1$Tanggal[2:116]), max(data1$Tanggal[2:116]), by="month"), labels=format(seq(min(data1$Tanggal[2:116]), max(data1$Tanggal[2:116]), by="month"), "%b %Y"), las=1, cex.axis=0.8)
mod1_auto1 <- auto.sarima(train.ts,
seasonal=F,
stepwise = FALSE,
stationary = FALSE,
trace = TRUE
)
##
## ARIMA(0,1,0) : 1973.24
## ARIMA(0,1,0) with drift : 1977.845
## ARIMA(0,1,1) : 1963.132
## ARIMA(0,1,1) with drift : 1967.51
## ARIMA(0,1,2) : 1967.647
## ARIMA(0,1,2) with drift : 1972.065
## ARIMA(0,1,3) : 1972.387
## ARIMA(0,1,3) with drift : 1976.771
## ARIMA(0,1,4) : 1968.677
## ARIMA(0,1,4) with drift : 1973.324
## ARIMA(0,1,5) : 1969.698
## ARIMA(0,1,5) with drift : 1974.162
## ARIMA(1,1,0) : 1963.32
## ARIMA(1,1,0) with drift : 1967.813
## ARIMA(1,1,1) : 1967.605
## ARIMA(1,1,1) with drift : 1972.026
## ARIMA(1,1,2) : 1968.028
## ARIMA(1,1,2) with drift : 1972.457
## ARIMA(1,1,3) : 1970.909
## ARIMA(1,1,3) with drift : 1975.198
## ARIMA(1,1,4) : 1971.979
## ARIMA(1,1,4) with drift : 1976.601
## ARIMA(2,1,0) : 1967.952
## ARIMA(2,1,0) with drift : 1972.427
## ARIMA(2,1,1) : 1972.23
## ARIMA(2,1,1) with drift : Inf
## ARIMA(2,1,2) : 1972.138
## ARIMA(2,1,2) with drift : 1976.49
## ARIMA(2,1,3) : 1975.223
## ARIMA(2,1,3) with drift : Inf
## ARIMA(3,1,0) : 1967.267
## ARIMA(3,1,0) with drift : 1971.504
## ARIMA(3,1,1) : 1968.404
## ARIMA(3,1,1) with drift : 1972.711
## ARIMA(3,1,2) : 1969.966
## ARIMA(3,1,2) with drift : 1974.331
## ARIMA(4,1,0) : 1964.487
## ARIMA(4,1,0) with drift : 1969.028
## ARIMA(4,1,1) : 1968.945
## ARIMA(4,1,1) with drift : 1973.521
## ARIMA(5,1,0) : 1968.636
## ARIMA(5,1,0) with drift : 1973.238
##
##
##
## Best model: ARIMA(0,1,1)
##
##
## SAMPLING FOR MODEL 'Sarima' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.005145 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 51.45 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.817 seconds (Warm-up)
## Chain 1: 0.667 seconds (Sampling)
## Chain 1: 1.484 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'Sarima' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7.7e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.77 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.767 seconds (Warm-up)
## Chain 2: 0.742 seconds (Sampling)
## Chain 2: 1.509 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'Sarima' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.733 seconds (Warm-up)
## Chain 3: 0.727 seconds (Sampling)
## Chain 3: 1.46 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'Sarima' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 5e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.829 seconds (Warm-up)
## Chain 4: 0.774 seconds (Sampling)
## Chain 4: 1.603 seconds (Total)
## Chain 4:
Parameter Penjelasan train.ts Data input berupa time series (deret waktu) untuk pelatihan model SARIMA. Harus dalam format ts atau tsibble, tergantung paket yang digunakan.
seasonal = FALSE Menyatakan bahwa model tidak mengikutkan komponen musiman. Artinya model SARIMA dikurangi komponen musiman, menjadi ARIMA biasa.
stepwise = FALSE Proses pencarian model tidak menggunakan algoritma stepwise (lebih cepat, tapi kurang menyeluruh). Dengan FALSE, semua kemungkinan kombinasi parameter diperiksa (comprehensive search), tetapi waktu komputasi bisa lebih lama.
stationary = FALSE Mengizinkan pencarian model dengan asumsi data tidak stasioner (nanti fungsi akan mencari orde diferensiasi d untuk menjadikannya stasioner).
trace = TRUE Menampilkan log progres model selama pencarian, berguna untuk melihat model mana saja yang diuji dan kriteria AIC/BIC-nya.
mod1_auto1$model
##
## y ~ Sarima(0,1,1)
## 116 observations and 1 dimension
## Differences: 1 seasonal Differences: 0
## Current observations: 115
##
## Priors:
## Intercept:
## [1] "mu0 [ ] ~ student ( mu = 0 ,sd = 2.5 ,df = 6 )"
##
## Scale Parameter:
## [1] "sigma0 [ ] ~ student ( mu = 0 ,sd = 1 ,df = 7 )"
##
## [,1]
## [1,] "ma [ 1 ] ~ normal ( mu = 0 ,sd = 0.5 )"
library(MCMCpack)
## Warning: package 'MCMCpack' was built under R version 4.3.3
## Loading required package: coda
## Warning: package 'coda' was built under R version 4.3.3
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2025 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
library(mcmc)
## Warning: package 'mcmc' was built under R version 4.3.3
mcmc_plot(mod1_auto1)
autoplot(mod1_auto1)
#Eksplorasi
sisaan.da <- residuals(mod1_auto1)
par(mfrow=c(2,2))
qqnorm(sisaan.da)
qqline(sisaan.da, col = "blue", lwd = 2)
plot(c(1:length(sisaan.da)),sisaan.da)
acf(sisaan.da)
pacf(sisaan.da)
#1) Sisaan Menyebar Normal
tseries::jarque.bera.test(sisaan.da)
##
## Jarque Bera Test
##
## data: sisaan.da
## X-squared = 2.4117, df = 2, p-value = 0.2994
#2) Sisaan saling bebas/tidak ada autokorelasi
Box.test(sisaan.da, type = "Ljung") #tak tolak H0 > sisaan saling bebas
##
## Box-Ljung test
##
## data: sisaan.da
## X-squared = 0.19453, df = 1, p-value = 0.6592
#3) Sisaan homogen
Box.test((sisaan.da)^2, type = "Ljung") #tak tolak H0 > sisaan homogen
##
## Box-Ljung test
##
## data: (sisaan.da)^2
## X-squared = 1.6527, df = 1, p-value = 0.1986
#4) Nilai tengah sisaan sama dengan nol
t.test(sisaan.da, mu = 0, conf.level = 0.95) #tak tolak h0 > nilai tengah sisaan sama dengan 0
##
## One Sample t-test
##
## data: sisaan.da
## t = -0.55787, df = 115, p-value = 0.578
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -278.0404 155.8424
## sample estimates:
## mean of x
## -61.09898
#---FORECAST---#
ramalan.da <- forecast::forecast(mod1_auto1, data=train , h = 30)
ramalan.da
## Point Forecast Lo 0.8 Hi 0.8 Lo 0.9 Hi 0.9
## 117 55114.12 53633.15 56635.36 53226.18 57054.99
## 118 55193.30 53650.04 56690.07 53301.28 57083.80
## 119 55164.67 53670.50 56666.26 53150.81 57063.18
## 120 55147.86 53726.45 56631.45 53244.96 57045.43
## 121 55144.13 53655.80 56578.08 53205.24 56988.48
## 122 55254.54 53795.46 56782.58 53376.59 57184.78
## 123 55149.96 53610.26 56618.07 53155.44 57048.01
## 124 55123.48 53755.89 56629.76 53202.29 56944.82
## 125 55187.74 53719.60 56698.02 53350.77 57177.38
## 126 55143.02 53619.91 56583.94 53248.34 56935.40
## 127 55117.02 53615.10 56671.43 53228.57 57050.50
## 128 55133.12 53591.13 56602.80 53280.52 57123.70
## 129 55163.62 53693.32 56703.39 53339.05 57152.85
## 130 55186.28 53637.17 56734.30 53304.61 57140.33
## 131 55186.02 53652.53 56653.95 53293.41 57050.47
## 132 55106.07 53681.73 56549.31 53256.47 56896.83
## 133 55181.87 53652.57 56733.41 53237.99 57192.09
## 134 55202.93 53653.69 56704.69 53260.41 57228.47
## 135 55190.77 53695.57 56689.87 53289.23 57094.30
## 136 55190.04 53624.94 56650.80 53205.95 57065.99
## 137 55176.52 53723.20 56622.44 53261.32 57112.19
## 138 55178.56 53659.19 56692.39 53225.64 57312.08
## 139 55194.43 53675.60 56651.84 53052.77 57119.45
## 140 55129.46 53617.27 56641.73 53161.28 57065.83
## 141 55167.28 53667.38 56623.83 53230.99 57050.98
## 142 55197.11 53586.35 56689.79 53119.19 57093.14
## 143 55183.25 53757.85 56634.07 53324.92 57062.11
## 144 55253.14 53733.57 56678.98 53311.58 57232.44
## 145 55155.61 53643.39 56668.86 53220.85 57142.75
## 146 55198.75 53697.63 56715.84 53368.12 57115.42
data.ramalan.da <- ramalan.da$mean
data.ramalan.da
## Time Series:
## Start = 117
## End = 146
## Frequency = 1
## [1] 55114.12 55193.30 55164.67 55147.86 55144.13 55254.54 55149.96 55123.48
## [9] 55187.74 55143.02 55117.02 55133.12 55163.62 55186.28 55186.02 55106.07
## [17] 55181.87 55202.93 55190.77 55190.04 55176.52 55178.56 55194.43 55129.46
## [25] 55167.28 55197.11 55183.25 55253.14 55155.61 55198.75
# Definisikan end_date sebagai tanggal terakhir dari data aktual
end_date <- max(data1$Tanggal)
# Buat rentang tanggal ramalan (pastikan panjangnya sama dengan data.ramalan.da)
forecast_dates <- seq.Date(from = end_date + 1, by = "day", length.out = length(data.ramalan.da))
# Gabungkan tanggal aktual dan ramalan
all_dates <- c(as.Date(data1$Tanggal), forecast_dates)
# Buat data frame gabungan
df <- data.frame(
Tanggal = all_dates,
Value = c(data1.ts, rep(NA, length(data.ramalan.da))),
Type = c(rep("data1.ts", length(data1.ts)), rep("data.ramalan.da", length(data.ramalan.da)))
)
# Masukkan nilai ramalan ke kolom Value
df$Value[df$Type == "data.ramalan.da"] <- data.ramalan.da
# Tentukan batas x-axis sampai tanggal terakhir ramalan
x_range <- c(min(df$Tanggal), max(df$Tanggal))
# Plot data aktual
plot(df$Tanggal[df$Type == "data1.ts"], df$Value[df$Type == "data1.ts"],
type = "l", col = "black",
xlab = "Periode", ylab = "Harga",
xaxt = "n",
xlim = x_range) # Set xlim agar mencakup semua tanggal
# Tambahkan garis ramalan
lines(df$Tanggal[df$Type == "data.ramalan.da"], df$Value[df$Type == "data.ramalan.da"], col = "red")
# Atur label sumbu X agar mencakup semua bulan hingga Juni 2025
axis(1,
at = seq(from = min(df$Tanggal), to = max(df$Tanggal), by = "month"),
labels = format(seq(from = min(df$Tanggal), to = max(df$Tanggal), by = "month"), "%d %b %Y"),
las = 1, # label horizontal
cex.axis = 0.8)
# Tambahkan legenda
legend("bottomleft", legend = c("Data Aktual", "Ramalan"),
col = c("black", "red"), lty = 1, cex = 0.8)
fitted_values_train <- fitted(mod1_auto1)
accuracy(data.ramalan.da[1:30], head(test.ts,n=length(test.ts)))
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -308.8249 1606.785 1138.343 -0.650953 2.119265 -0.08112776 0.6502123
actual_train <- train.ts
mape_train <- mean(abs((actual_train - fitted_values_train) / actual_train)) * 100
mape_train
## [1] 1.626039