Data yang digunakan adalah data historis harga emas mingguan dengan periode januari 2020 - mei 2024 dengan total 226 baris.

library(readxl)
df = read_excel("C:\\Users\\Diva\\Documents\\SEMESTER 6\\PSB\\tugas akhir\\Data Harga Emas2.xlsx")
Y = df[, 2]
t = df[, 1]# store data in Y
df=df[211:437,]
df
## # A tibble: 227 × 2
##    Date                Harga
##    <dttm>              <dbl>
##  1 2020-01-10 00:00:00 1554.
##  2 2020-01-17 00:00:00 1558.
##  3 2020-01-24 00:00:00 1564.
##  4 2020-01-31 00:00:00 1584.
##  5 2020-02-07 00:00:00 1573.
##  6 2020-02-14 00:00:00 1581.
##  7 2020-02-21 00:00:00 1643.
##  8 2020-02-28 00:00:00 1610.
##  9 2020-03-06 00:00:00 1684.
## 10 2020-03-13 00:00:00 1563.
## # ℹ 217 more rows
# Membentuk Objek Time Series Ekspor
ts<-ts(df$`Harga`, frequency=12, start=2017)

## Membentuk Time Series Plot
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
ggplot(data = df, aes(x = Date, y = `Harga`)) +
  geom_line(color="navy", size=1.2) +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Statistika Deskriptif

library(psych)
## Warning: package 'psych' was built under R version 4.3.2
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
knitr::kable(describe(df$Harga))
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 227 1857.604 148.4879 1841.55 1851.17 132.1738 1494.4 2401.5 907.1 0.7737547 1.846555 9.855488

Identifikasi Distribusi Data Menggunakan Uji Kolmogorov-Smirnov

Hipotesis:

H0 = Data mengikuti suatu distribusi

H1 = Data tidak mengikuti suatu distribusi

data <- df$`Harga`
# Uji distribusi normal
result_normal <- ks.test(data, "pnorm", mean = mean(data), sd = sd(data))
## Warning in ks.test.default(data, "pnorm", mean = mean(data), sd = sd(data)):
## ties should not be present for the Kolmogorov-Smirnov test
# Uji distribusi eksponensial
result_exponential <- ks.test(data, "pexp", rate = 1/mean(data))
## Warning in ks.test.default(data, "pexp", rate = 1/mean(data)): ties should not
## be present for the Kolmogorov-Smirnov test
# Uji distribusi gamma
result_gamma <- ks.test(data, "pgamma", shape = 2, rate = 2)
## Warning in ks.test.default(data, "pgamma", shape = 2, rate = 2): ties should
## not be present for the Kolmogorov-Smirnov test
# Uji distribusi Weibull
result_weibull <- ks.test(data, "pweibull", shape = 2, scale = 2)
## Warning in ks.test.default(data, "pweibull", shape = 2, scale = 2): ties should
## not be present for the Kolmogorov-Smirnov test
distribusi <- c("Normal", "Eksponensial","Gamma", "Weibull")
statistik <- c(result_normal$statistic, result_exponential$statistic, result_gamma$statistic, result_weibull$statistic)
nilai_p <- c(result_normal$p.value, result_exponential$p.value, result_gamma$p.value, result_weibull$p.value)
tabel <- data.frame(distribusi,round(statistik,4), round(nilai_p,4))
rownames(tabel) <- NULL
colnames(tabel) <- c("Distribusi", "Nilai Statistik Hitung", "Nilai p")
Keterangan <- c("Tak Tolak H0", "Tolak H0", "Tolak H0", "Tolak H0")
tabel$Keterangan <- Keterangan
knitr::kable(tabel)
Distribusi Nilai Statistik Hitung Nilai p Keterangan
Normal 0.0632 0.3255 Tak Tolak H0
Eksponensial 0.5623 0.0000 Tolak H0
Gamma 1.0000 0.0000 Tolak H0
Weibull 1.0000 0.0000 Tolak H0

Berdasarkan hasil di atas, dapat dikatakan bahwa data harga emas bulanan di Indonesia tahun 2015-2023 mengikuti distribusi Normal pada taraf nyata 5%.

Pendugaan Parameter Distribusi

Penerapan Gibs Sampling

###data
y <- df$`Harga`
mean.y <- mean(y)
var.y <- var(y)
n <- length(y)

mu0 <- 1857 #rata2
t20 <- 21904 #sd^2
s20 <- 1/2 #ini bingung dpt darimana
nu0 <- 1/6 #ini juga

###starting value
S <- 5000 #iterasi
PHI <- matrix(nrow=S, ncol=2)
PHI[1,] <- phi <- c(mean.y,1/var.y)

### Gibbs sampling
set.seed(1)
for(s in 2:S){
  #generate a new theta value from its full conditional
  mun <- (mu0/t20 + n*mean.y*phi[2])/(1/t20 + n*phi[2])
  t2n <- 1/(1/t20 + n*phi[2])
  phi[1] <- rnorm(1,mun,sqrt(t2n))
  
  #generate a new 1/sigma^2 value from its full conditional
  nun <- nu0 + n
  s2n <- (nu0*s20 + (n-1)*var.y + n*(mean.y-phi[1])^2)/nun
  phi[2] <- rgamma(1,nun/2,nun*s2n/2)
  
  PHI[s,] <- phi
}

# Trace Plot
plot(PHI[,1], type="l", ylab = "Theta", xlab="Iterasi ke-")

plot(1/PHI[,2], type="l", ylab = "Sigma^2", xlab="Iterasi ke-")

#Penduga Parameter

# Penduga Bagi Theta
theta.duga <- mean(PHI[,1])

# Penduga Bagi Sigma^2
varians <- 1/PHI[,2]
sigma2.duga <- mean(varians)

cat(paste("theta duga:", theta.duga, "\n" , "sigma^2 duga:", sigma2.duga))
## theta duga: 1857.60829822797 
##  sigma^2 duga: 22259.2686941118

#Pendekatan 95% Credible Interval

#95% credible interval

## Theta Duga
nilai.kuantil <- quantile(PHI[,1], 
                          probs=c(0.025,0.975))
cat(paste("Pendekatan 95% credible interval bagi theta duga: ["
          , round(nilai.kuantil[1], 4), " ; ",
          round(nilai.kuantil[2], 4), "]\n", sep=""))
## Pendekatan 95% credible interval bagi theta duga: [1838.134 ; 1876.7111]
## Sigma2 Duga
nilai.kuantil <- quantile(varians, 
                          probs=c(0.025,0.975))
cat(paste("Pendekatan 95% credible interval bagi sigma^2 duga: ["
          , round(nilai.kuantil[1],4), " ; ",
          round(nilai.kuantil[2],4), "]\n", sep=""))
## Pendekatan 95% credible interval bagi sigma^2 duga: [18451.5063 ; 26831.695]