Tugas Pertemuan 6
Fathoni Sabri
2025-10-02
Library
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.1library(tsibble)## Warning: package 'tsibble' was built under R version 4.5.1## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr##
## Attaching package: 'tsibble'## The following objects are masked from 'package:base':
##
## intersect, setdiff, unionlibrary(tseries)## Warning: package 'tseries' was built under R version 4.5.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(MASS)## Warning: package 'MASS' was built under R version 4.5.1library(dplyr)## Warning: package 'dplyr' was built under R version 4.5.1##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(TSA)## Warning: package 'TSA' was built under R version 4.5.1##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary(forecast)## Warning: package 'forecast' was built under R version 4.5.1## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.Arima TSA
## plot.Arima TSAInput Data
data<-read.csv("C:/Users/Fathoni Sabri/Downloads/data_fathoni.csv")
head(data)## Date Close Volume
## 1 10/9/2023 973.5909 22298500
## 2 10/10/2023 1213.2440 13625200
## 3 10/11/2023 1512.8105 38431400
## 4 10/12/2023 1887.2686 219262300
## 5 10/13/2023 2356.5891 39504100
## 6 10/16/2023 2746.0256 227720100str(data)## 'data.frame': 111 obs. of 3 variables:
## $ Date : chr "10/9/2023" "10/10/2023" "10/11/2023" "10/12/2023" ...
## $ Close : num 974 1213 1513 1887 2357 ...
## $ Volume: int 22298500 13625200 38431400 219262300 39504100 227720100 215604400 311694400 150507700 85371300 ...Data Formating
Mengubah Format Tanggal
data$Date <- as.Date(data$Date, format = "%m/%d/%Y")
head(data)## Date Close Volume
## 1 2023-10-09 973.5909 22298500
## 2 2023-10-10 1213.2440 13625200
## 3 2023-10-11 1512.8105 38431400
## 4 2023-10-12 1887.2686 219262300
## 5 2023-10-13 2356.5891 39504100
## 6 2023-10-16 2746.0256 227720100Memilih Kolom
data <- data %>%
select(Date, Close)
head(data)## Date Close
## 1 2023-10-09 973.5909
## 2 2023-10-10 1213.2440
## 3 2023-10-11 1512.8105
## 4 2023-10-12 1887.2686
## 5 2023-10-13 2356.5891
## 6 2023-10-16 2746.0256Splitting Data
Data dibagi menjadi 80% dan 20%
# hitung jumlah data
n <- nrow(data)
n_train <- floor(0.8 * n)
# bagi data
training <- data[1:n_train, ]
testing <- data[(n_train+1):n, ]
# ubah ke time series
train <- ts(training$Close)
test <- ts(testing$Close)length(train)## [1] 88Eksplorasi Data
mean <- mean(data$Close)
mean## [1] 5523.31Plot Data Penuh
dt <- ts(data$Close)
plot <- dt |> as_tsibble() |>
ggplot(aes(x = index, y = value)) + geom_line() + theme_bw() +
geom_hline(yintercept = mean, color = "red", linetype = "dashed") +
xlab("Periode") + ylab("Harga Close")
plot
Plot Data Latih
dt.train <- ts(training$Close)
plot <- dt.train |> as_tsibble() |>
ggplot(aes(x = index, y = value)) + geom_line() + theme_bw() +
geom_hline(yintercept = mean, color = "red", linetype = "dashed") +
xlab("Periode") + ylab("Harga Close")
plot
Plot Data Uji
dt.test <- ts(testing$Close)
plot <- dt.test |> as_tsibble() |>
ggplot(aes(x = index, y = value)) + geom_line() + theme_bw() +
geom_hline(yintercept = mean, color = "red", linetype = "dashed") +
xlab("Periode") + ylab("Harga Close")
plot
Pemeriksaan Kestasioneran Data
Plot ACF
acf(dt.train)
Terlihat bahwa data turun perlahan (tails off), ini mengindikasikan data tidak stasioner dalam rataan.
Uji ADF
tseries::adf.test(dt.train)##
## Augmented Dickey-Fuller Test
##
## data: dt.train
## Dickey-Fuller = -1.8306, Lag order = 4, p-value = 0.6457
## alternative hypothesis: stationaryP-Value gagal menolak H0, artinya data tidak menyebar stasioner dalam rataan. Hal ini sesuai dengan plot ACF sebelumnya.
Kestasioneran Dalam Ragam
index<-seq(1:88)
bc <- boxcox(
dt.train ~ index,
lambda = seq(-4.0, 4.0, by = 0.001)
)
lambda <- bc$x[which.max(bc$y)]
lambda## [1] 1.222ci_lambda <- bc$x[bc$y > max(bc$y) - 0.5 * qchisq(0.95, 1)]
ci_lower <- min(ci_lambda)
ci_upper <- max(ci_lambda)
c(Batas_Bawah = ci_lower, Batas_Atas = ci_upper)## Batas_Bawah Batas_Atas
## 0.768 1.752Hasil analisis Box-Cox menunjukkan bahwa nilai λ=1 berada pada selang kepercayaan 0.768 hingga 1.752. Ini menunjukkan bahwa data menyebar stasioner terhadap varians.
Penanganan Ketidakstasioneran
Differencing
differencing1 <- diff(dt.train, differences = 1)
ts.plot(differencing1,
main = "Differencing Pertama",
col = "blue", lwd = 2)
length(differencing1)## [1] 87Plot ACF Setelah Differencing
acf(differencing1)
Terjadi perubahan, di mana plot ACF pada data yang tadinya tails off berubah menjadi cut off. Namun, untuk memastikan bahwa ketidakstasioneran pada data masih ada atau tidak. perlu dilakukan uji ADF. Selain itu, data cut of pada \(Lag 0\).
Uji ADF Setelah Differencing
adf.test(differencing1)## Warning in adf.test(differencing1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: differencing1
## Dickey-Fuller = -4.2067, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryUji ADF menghasilkan \(P-value\) 0.01 yang lebih kecil dari nilai signifikannya. Sehingga, cukup bukti untuk mengatakan data telah stasioner dalam rataan.
Ketidakstasioneran Dalam Varians Setelah Differencing
length(differencing1)## [1] 87index2 <- seq(1:87)
bc2 <- boxcox(
differencing1 - min(differencing1) + 1 ~ index2,
lambda = seq(0, 5, by = 0.001)
)
lambda <- bc2$x[which.max(bc2$y)]
lambda## [1] 1.236ci_lambda <- bc2$x[bc2$y > max(bc2$y) - 0.5 * qchisq(0.95, 1)]
ci_lower <- min(ci_lambda)
ci_upper <- max(ci_lambda)
c(Batas_Bawah = ci_lower, Batas_Atas = ci_upper)## Batas_Bawah Batas_Atas
## 0.908 1.664Hasil analisis Box-Cox menunjukkan bahwa nilai λ=1 berada pada selang kepercayaan 0.908 hingga 1.664. Ini menunjukkan bahwa setelah dilakukan differencing, data tetap menyebar stasioner terhadap varians.
Identifikasi Model
6.1 Plot ACF untuk MA(q)
acf(differencing1)
Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 0.
6.2 Plot PACF untuk AR(p=)
pacf(differencing1)
Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 0.
eacf(differencing1)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 x o o o o o o o o o o o o o
## 3 o x o o o o o o o o o o o o
## 4 o x o x o o o o o o o o o o
## 5 x x o x o o o o o o o o o o
## 6 o x o x o o o o o o o o o o
## 7 x x x o o o o o o o o o o omaka model tentatifnya adalah ARIMA(0,1,0). Sebagai perbandingan, kita juga akan melihat hasil dari model ARIMA(1,1,1).
Pendugaan Parameter Model Tentatif
ARIMA(0,1,0)
model1 = Arima(differencing1, order=c(0,1,0), method="ML", include.drift = TRUE)
summary(model1)## Series: differencing1
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## -1.9152
## s.e. 50.4047
##
## sigma^2 = 221065: log likelihood = -650.69
## AIC=1305.39 AICc=1305.53 BIC=1310.29
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.002776646 464.7395 330.5924 NaN Inf 0.9883824 -0.5128758lmtest::coeftest(model1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## drift -1.9152 50.4047 -0.038 0.9697Drift tidak signifikan, artinya data Data sudah stasioner dan tidak punya tren deterministik.
ARIMA(1,1,1)
model2 = Arima(differencing1, order=c(1,1,1), method="ML")
summary(model2)## Series: differencing1
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.0148 -0.9280
## s.e. 0.1158 0.0483
##
## sigma^2 = 120766: log likelihood = -625.16
## AIC=1256.32 AICc=1256.61 BIC=1263.68
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -47.52319 341.4691 255.5505 NaN Inf 0.7640271 -0.03083367lmtest::coeftest(model2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.014822 0.115782 0.128 0.8981
## ma1 -0.927991 0.048290 -19.217 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pada Model ARIMA(1,1,1) hanya 1 parameter saja yang signifikan.
Ringkasan AIC
model_tentatif <- data.frame(
Model = c("ARIMA(0,1,0)","ARIMA(1,1,1)"),
AIC = c(AIC(model1),AIC(model2)),
BIC = c(BIC(model1),BIC(model2))
)
model_tentatif## Model AIC BIC
## 1 ARIMA(0,1,0) 1305.386 1310.294
## 2 ARIMA(1,1,1) 1256.319 1263.682Berdasarkan Nilai AIC terkecil, Model ARIMA (1,1,1) dipilih sebagai model.
Analisis Sisaan
Eksplorasi Sisaan
sisaan <- model2$residuals
par(mfrow = c(2,2),
mar = c(4,4,2,1),
oma = c(0,0,2,0))
# Q-Q Plot
qqnorm(sisaan, main = "Normal Q-Q Plot")
qqline(sisaan, col = "blue", lwd = 2)
# Plot Sisaan vs Waktu
plot(sisaan, type = "p",
main = "Plot Sisaan vs Waktu",
xlab = "Waktu", ylab = "Sisaan")
# ACF
acf(sisaan, main = "ACF Sisaan")
# PACF
pacf(sisaan, main = "PACF Sisaan")
Berdasarkan plot diagnostik residual (Q-Q plot, ACF, PACF, dan plot sisaan terhadap waktu), residual model ARIMA(1,1,1) berdistribusi normal dan tidak menunjukkan autokorelasi yang signifikan. Hal ini menunjukkan bahwa model telah sesuai dengan data dan memenuhi asumsi klasik model ARIMA, sehingga layak digunakan untuk peramalan.
Uji Formal
ks.test(sisaan,"pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.56322, p-value < 2.2e-16
## alternative hypothesis: two-sidedBerdasarkan uji KS tersebut, didapat p-value sebesar \(0.0000\) yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.
Box.test(sisaan, type = "Ljung")##
## Box-Ljung test
##
## data: sisaan
## X-squared = 0.085598, df = 1, p-value = 0.7699Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar \(0.7699\) yang lebih besar dari taraf nyata 5% sehingga gagal tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini sesuai dengan hasil eksplorasi.
Box.test((sisaan)^2, type = "Ljung")##
## Box-Ljung test
##
## data: (sisaan)^2
## X-squared = 3.6426, df = 1, p-value = 0.05632Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar \(0.05632\) yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa tak cukup bukti untuk mengatakan ragam sisaan tidak homogen.
t.test(sisaan, mu = 0, conf.level = 0.95)##
## One Sample t-test
##
## data: sisaan
## t = -1.3033, df = 86, p-value = 0.1959
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -120.00968 24.96329
## sample estimates:
## mean of x
## -47.52319Berdasarkan uji t satu sampel terhadap sisaan model ARIMA(1,1,1), diperoleh nilai p = 0.1959 (> 0.05) sehingga gagal menolak hipotesis nol. Hal ini menunjukkan bahwa rata-rata sisaan tidak berbeda signifikan dari nol. Dengan demikian, model tidak mengandung bias sistematis dalam prediksi, dan asumsi nilai tengah sisaan = 0 telah terpenuhi. Hal ini sesuai dengan eksplorasi yang telah dilakukan.
Overfitting
Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(1,1,1) Kandidat model overfitting adalah ARIMA(2,1,1) dan ARIMA(1,1,2).
ARIMA (2,1,1)
model3 = Arima(differencing1, order=c(2,1,1), method="ML")
summary(model3)## Series: differencing1
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 0.0214 0.0527 -0.9355
## s.e. 0.1155 0.1145 0.0463
##
## sigma^2 = 121915: log likelihood = -625.05
## AIC=1258.11 AICc=1258.6 BIC=1267.92
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -48.70135 341.0423 256.0735 NaN Inf 0.7655909 -0.03666833lmtest::coeftest(model3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.021372 0.115497 0.1850 0.8532
## ar2 0.052749 0.114490 0.4607 0.6450
## ma1 -0.935516 0.046301 -20.2052 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pada ARIMA(2,1,1) dengan menaikkan satu ordo AR ternyata tidak lebih baik dari Model ARIMA(1,1,1). Karena nilai AIC dari Model ARIMA(2,1,1) lebih besar serta ada dua parameter yang tidak siginifikan, yaitu ar1 dan ar2.
ARIMA(1,1,2)
model4 = Arima(differencing1, order=c(1,1,2), method="ML")
summary(model4)## Series: differencing1
## ARIMA(1,1,2)
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ar1 ma1 ma2
## -0.3769 -0.5422 -0.3557
## s.e. NaN NaN NaN
##
## sigma^2 = 122243: log likelihood = -625.16
## AIC=1258.33 AICc=1258.82 BIC=1268.15
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -47.09998 341.5008 255.395 NaN Inf 0.7635622 -0.02363055lmtest::coeftest(model4)## Warning in sqrt(diag(se)): NaNs produced##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.37689 NaN NaN NaN
## ma1 -0.54222 NaN NaN NaN
## ma2 -0.35571 NaN NaN NaNMuncul peringatan “NaNs produced”.
model5<-forecast::Arima(differencing1, order=c(1,1,2))lmtest::coeftest(model5)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.81651 0.92596 -0.8818 0.3779
## ma1 -0.11268 0.94671 -0.1190 0.9053
## ma2 -0.75205 0.88028 -0.8543 0.3929ARIMA(1,1,2) dengan menaikkan satu ordo MA ternyata masih tidak lebih baik dari Model ARIMA(1,1,1). Karena nilai AIC dari Model ARIMA(1,1,2) lebih besar serta tidak ada parameter yang signifikan pada alpha 5%.
Karena model hasil overfitting hasilnya tidak lebih baik daripada Model ARIMA(1,1,1), maka model terbaik yang dipilih adalah ARIMA(1,1,1) dilanjutkan dengan peramalan.
Model Terbaik: ARIMA(1,1,1)
model2 = Arima(differencing1, order=c(1,1,1), method="ML")
summary(model2)## Series: differencing1
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.0148 -0.9280
## s.e. 0.1158 0.0483
##
## sigma^2 = 120766: log likelihood = -625.16
## AIC=1256.32 AICc=1256.61 BIC=1263.68
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -47.52319 341.4691 255.5505 NaN Inf 0.7640271 -0.03083367lmtest::coeftest(model2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.014822 0.115782 0.128 0.8981
## ma1 -0.927991 0.048290 -19.217 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Nilai AIC: \(1256.32\)
Nilai AICc: \(1256.61\)
Nilai BIC:
\(1263.68\)
Parameter ma1
signifikan pada taraf nyata 5%
Peramalan
ramalan <- forecast::forecast(model2, h = length(dt.test))
data.ramalan <- ramalan$mean
plot(ramalan)
# Inverse differencing
pt_akhir_train <- tail(dt.test, 1)
forecast <- diffinv(data.ramalan, differences = 1, xi = pt_akhir_train)
hasil <- forecast[-c(0,1)]
hasil <- forecast[1:length(dt.test)]
# Bandingkan aktual vs forecast
perbandingan <- cbind(
Aktual = dt.test,
Forecast = hasil
)
print(head(perbandingan, 10)) ## Time Series:
## Start = 1
## End = 10
## Frequency = 1
## Aktual Forecast
## 1 5470.873 5345.967
## 2 5495.854 5320.577
## 3 5495.854 5293.700
## 4 5495.854 5266.801
## 5 5570.798 5239.902
## 6 5495.854 5213.002
## 7 5495.854 5186.102
## 8 5495.854 5159.203
## 9 5945.515 5132.303
## 10 6045.440 5105.404akurasi_test <- accuracy(hasil, dt.test)
akurasi_test## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 662.6307 761.6958 662.6307 11.31225 11.31225 0.837394 4.054315ts.plot(dt, col = "black", lty = 1,
xlab = "Waktu", ylab = "Harga Saham")
lines(88:(87+length(hasil)), hasil, col = "blue", lty = 2, lwd = 2)
legend("topleft",
legend = c("Data Aktual", "Forecast"),
col = c("black", "blue"),
lty = c(1,2,3), lwd = c(1,1,1), bty = "n")
Berdasarkan hasil peramalan, dapat dilihat bahwa nilai MAPE yang didapatkan sebesar 11,31% yang sebenarnya masih cukup baik untuk model tersebut, tetapi dikarenakan nilai ACF1 dan Theil’s U yang besar menunjukkan model belum menangkap pola autokorelasi residual dengan baik. Peramalan ini hanya cocok untuk jangka pendek.