1 PENDAHULUAN

Laporan ini menyajikan analisis komprehensif menggunakan tiga metode utama dalam Time Series Analysis. Dokumen ini disusun untuk memenuhi kebutuhan UAS Analisis Data Deret Waktu 2 dengan alur sebagai berikut:

  1. Fungsi Transfer (Transfer Function): Digunakan untuk memodelkan hubungan kausalitas satu arah (Input \(\rightarrow\) Output) dengan memperhitungkan jeda waktu (lag).
  2. Vector Autoregression (VAR): Digunakan untuk memodelkan sistem variabel yang saling mempengaruhi (endogen) secara simultan.
  3. Singular Spectrum Analysis (SSA): Metode non-parametrik untuk membedah struktur sinyal (tren dan osilasi) dari deret waktu tunggal.

BAGIAN 1: ANALISIS FUNGSI TRANSFER

Konsep Dasar: Model Fungsi Transfer (Transfer Function Noise Model) menggabungkan model dinamis (hubungan antara Input \(X_t\) dan Output \(Y_t\)) dengan model gangguan (noise). Tujuannya adalah melihat bagaimana variasi pada input mempengaruhi output di masa depan setelah melewati mekanisme delay (penundaan).

1.1 1. Input dan Persiapan Data

Langkah pertama adalah memuat data dan memisahkan variabel menjadi Input (Impor) dan Output (Harga). Data dibagi menjadi training (pembentukan model) dan testing (validasi).

data <- read.csv(file.choose(), header=T, sep=";", dec=",")
head(data)
harga = data$Rerata_Harga_Beras
impor_ton = data$Kuantitas_Impor/1000

beras <- data.frame(harga=harga, impor=impor_ton)
head(beras)

1.2 2. Statistik Deskriptif

Melihat karakteristik dasar data seperti rata-rata dan penyebaran data (standar deviasi).

##STATISTIK DESKRIPTIF ----------

# SUMMARY
summary(beras)
##      harga           impor       
##  Min.   :10860   Min.   :  2000  
##  1st Qu.:11519   1st Qu.: 22081  
##  Median :12304   Median : 38788  
##  Mean   :12452   Mean   :103310  
##  3rd Qu.:13272   3rd Qu.:113208  
##  Max.   :17086   Max.   :567211
# STANDAR DEVIASI
sd(beras$harga)
## [1] 1166.619
sd(beras$impor)
## [1] 131011.8

Interpretasi: Output summary memberikan gambaran rentang data (Min, Max, Mean). Standar deviasi (sd) menunjukkan seberapa besar fluktuasi harga dan volume impor dari rata-ratanya. Nilai SD yang tinggi menandakan volatilitas data yang tinggi.

1.3 3. Partisi Data

Data dibagi menjadi 81% untuk pelatihan model (training) dan 19% untuk pengujian (testing).

##PARTISI DATA ----------
datatraining <- head(beras, round(0.81*nrow(data)))
datatesting <- tail(beras, round(0.19*nrow(data)))
nrow(datatraining)
## [1] 103
nrow(datatesting)
## [1] 24

1.4 4. Plot Data Time Series

Visualisasi data untuk melihat pola tren atau musiman secara kasat mata.

##PLOT DATA ----------

# HARGA
Yt <- ts(datatraining$harga, start = c(2014,1), freq=12)
plot(Yt, ylab="Rupiah (Rp)", xlab="Tahun", main="Harga Beras (Training)")

# IMPOR
Xt <- ts(datatraining$impor, start = c(2014,1), freq=12)
plot(Xt, ylab="Ton", xlab="Tahun", main="Kuantitas Impor Beras (Training)")

Interpretasi: Plot deret waktu membantu kita mendeteksi ketidakstasioneran. Jika grafik menunjukkan tren menanjak/menurun atau varians yang melebar, maka data perlu ditangani (transformasi/differencing) sebelum pemodelan.

1.5 5. Identifikasi Kestasioneran (Awal)

Mengecek kestasioneran menggunakan plot ACF/PACF dan uji formal Augmented Dickey-Fuller (ADF).

##PLOT ACF DAN PACF ----------
par(mfrow=c(1,2))

# HARGA
acf(Yt, lag.max=36)
pacf(Yt, lag.max=36)

# IMPOR
acf(Xt, lag.max=36)
pacf(Xt, lag.max=36)

Penanganan Stasioneritas: Dilakukan transformasi Box-Cox (untuk menstabilkan varians) dan Differencing (untuk menstabilkan rata-rata/menghilangkan tren).

##UJI STASIONER ----------
library(forecast)
library(tseries)

# HARGA
adf.test(Yt)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Yt
## Dickey-Fuller = -2.1998, Lag order = 4, p-value = 0.4938
## alternative hypothesis: stationary
diff_Yt <- diff(Yt, 1); adf.test(diff_Yt)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_Yt
## Dickey-Fuller = -4.1857, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
# IMPOR
lambdaX <- BoxCox.lambda(Xt); lambdaX
## [1] -0.4086458
tf_Xt <- ((Xt^lambdaX)-1)/lambdaX
lambdaX_2 <- BoxCox.lambda(tf_Xt); lambdaX_2
## [1] 1.999924
adf.test(tf_Xt)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  tf_Xt
## Dickey-Fuller = -4.0883, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
diff_tf_Xt <- diff(tf_Xt, 1); adf.test(diff_tf_Xt)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_tf_Xt
## Dickey-Fuller = -7.3699, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
tf_Xt %>% tsdisplay(main="tf_Xt")

diff_tf_Xt %>% tsdisplay(main="diff_tf_Xt")

Interpretasi: Pada uji ADF, jika P-value < 0.05, maka \(H_0\) (data tidak stasioner) ditolak, yang berarti data sudah stasioner. Pastikan data input dan output stasioner sebelum lanjut ke tahap pre-whitening.

1.6 6. Plot ACF dan PACF Setelah Penanganan

Melihat pola autokorelasi pada data yang sudah di-differencing.

##PLOT ACF DAN PACF SETELAH PENANGANAN STASIONERITAS ----------
par(mfrow=c(1,2))

# HARGA
acf(diff_Yt, lag.max=36)
pacf(diff_Yt, lag.max=36)

# IMPOR
acf(diff_tf_Xt, lag.max=36)
pacf(diff_tf_Xt, lag.max=36)

1.7 7. Pemodelan ARIMA untuk Input (Xt) - Pre-whitening Tahap 1

Konsep Pre-whitening: Kita harus “memutihkan” input (menghilangkan pola autokorelasi) agar korelasi yang terukur nanti murni hubungan antara Input dan Output, bukan dipengaruhi oleh masa lalu input itu sendiri.

##IDENTIFIKASI MODEL ARIMA UNTUK Xt ----------
library(lmtest)

modelI1 = Arima(tf_Xt, order=c(1,1,0), method="ML")
coeftest(modelI1) #signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ar1 -0.295761   0.099824 -2.9628 0.003048 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI2 = Arima(tf_Xt, order=c(2,1,0), method="ML")
coeftest(modelI2) #signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.356352   0.097481 -3.6556 0.0002566 ***
## ar2 -0.321410   0.102598 -3.1327 0.0017320 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI3 = Arima(tf_Xt, order=c(0,1,1), method="ML")
coeftest(modelI3) #signifikan
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.79271    0.16448 -4.8195 1.439e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI4 = Arima(tf_Xt, order=c(1,1,1), method="ML")
coeftest(modelI4) #signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.348494   0.093584   3.7239 0.0001962 ***
## ma1 -0.999999   0.028085 -35.6059 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI5 = Arima(tf_Xt, order=c(2,1,1), method="ML")
coeftest(modelI5)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.380490   0.099600   3.8202 0.0001333 ***
## ar2 -0.096427   0.106481  -0.9056 0.3651563    
## ma1 -1.000000   0.028980 -34.5064 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI6 = Arima(tf_Xt, order=c(0,1,2), method="ML")
coeftest(modelI6) #signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.611411   0.101919 -5.9990 1.985e-09 ***
## ma2 -0.388588   0.098339 -3.9515 7.766e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI7 = Arima(tf_Xt, order=c(1,1,2), method="ML")
coeftest(modelI7)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ar1  0.025741   0.346905  0.0742  0.94085  
## ma1 -0.635713   0.343360 -1.8514  0.06411 .
## ma2 -0.364233   0.342262 -1.0642  0.28724  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelI8 = Arima(tf_Xt, order=c(2,1,2), method="ML")
coeftest(modelI8)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.45212    0.21893 -2.0651   0.03891 *  
## ar2  0.20677    0.14677  1.4088   0.15891    
## ma1 -0.15544    0.18724 -0.8302   0.40644    
## ma2 -0.84455    0.18662 -4.5254 6.027e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# MODEL ARIMA YANG MEMENUHI ASUMSI UNTUK Xt
model_Xt <- data.frame(modelI1$aic, modelI2$aic, modelI3$aic, modelI4$aic, 
                       modelI6$aic); model_Xt
min(model_Xt)
## [1] -527.9554

Interpretasi: Model dengan nilai AIC terkecil dan koefisien yang signifikan dipilih sebagai model pre-whitening terbaik untuk seri Input.

1.8 8. Uji Asumsi Residual Input (Xt)

Memastikan sisaan model ARIMA input bersifat acak (white noise).

##UJI ASUMSI RESIDUAL Xt ----------
rX = residuals(modelI6)

# UJI NORMALITAS Xt 
nX = length(rX)
meanX = mean(rX)
sdX = sd(rX)
resX = rnorm(nX,meanX,sdX)
ks.test(rX, resX)
## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  rX and resX
## D = 0.13592, p-value = 0.2973
## alternative hypothesis: two-sided
# UJI AUTOKORELASI Xt
Box.test(rX, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  rX
## X-squared = 0.027134, df = 1, p-value = 0.8692
# UJI HETEROSKEDASTISITAS Xt
Box.test(rX^2, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  rX^2
## X-squared = 0.9984, df = 1, p-value = 0.3177

Interpretasi: * KS Test (Normalitas): P-value > 0.05 berarti residual berdistribusi normal. * Ljung-Box (Autokorelasi): P-value > 0.05 berarti tidak ada autokorelasi pada residual (sudah white noise).

1.9 9. Pre-Whitening (Cek Ljung-Box Kandidat Lain)

##PRE-WHITENING ----------
Box.test(residuals(modelI1), type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  residuals(modelI1)
## X-squared = 0.37154, df = 1, p-value = 0.5422
Box.test(residuals(modelI2), type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  residuals(modelI2)
## X-squared = 0.21803, df = 1, p-value = 0.6405
Box.test(residuals(modelI3), type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  residuals(modelI3)
## X-squared = 6.2105, df = 1, p-value = 0.0127
Box.test(residuals(modelI4), type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  residuals(modelI4)
## X-squared = 0.065996, df = 1, p-value = 0.7973
Box.test(residuals(modelI6), type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  residuals(modelI6)
## X-squared = 0.027134, df = 1, p-value = 0.8692

1.10 10. Korelasi Silang (Cross-Correlation) - Pre-whitening Tahap 2

Menghitung Cross Correlation Function (CCF) antara residual input dan output yang telah difilter. Ini digunakan untuk mengidentifikasi lag (\(b, r, s\)).

##KORELASI SILANG DERET OUTPUT DAN INPUT
library(tfarima)

# ANTARA VARIABEL HARGA DAN IMPOR
umXt <- um(tf_Xt, ar=0, i=1, ma=2)
umYt_Xt <- fit(umXt, diff_Yt)

par(mfrow=c(1,1))
ccf(diff_tf_Xt, diff_Yt, um.x = umXt, um.y = umYt_Xt, lag.max=36)

pccf(diff_tf_Xt, diff_Yt, um.x = umXt, um.y = umYt_Xt, lag.max=36)

#r=0, s=1, b=22

Interpretasi: Lihat lag negatif di mana batang grafik melewati garis batas biru (signifikan). * Lag signifikan pertama menunjukkan \(b\) (delay). * Pola setelahnya menunjukkan orde \(r\) dan \(s\). * Dalam kasus ini, diidentifikasi \(b=22\), \(r=0\), \(s=1\).

1.11 11. Identifikasi dan Estimasi Fungsi Transfer

Membangun model lengkap yang terdiri dari komponen fungsi transfer (input-output) dan komponen noise (ARIMA).

##IDENTIFIKASI FUNGSI TRANSFER AWAL ----------
library(MTS)
library(forecast)

pre_tf <- tfm1(Yt, tf_Xt, orderN=c(0,1,2), orderX=c(0,1,22))
## ARMA coefficients & s.e.: 
##             ma1   ma2
## coef.arma 0.385 0.259
## se.arma   0.110 0.116
## Transfer function coefficients & s.e.: 
##      intercept    X     
## v        -17.0 1258 -238
## se.v      43.1 1182 1058
pre_tf$residuals %>% ggtsdisplay()

pre_tf$residuals %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 7.5638, df = 10, p-value = 0.6714
## 
## Model df: 0.   Total lags used: 10
# MODEL SISAAN
arima_Nt <- auto.arima(pre_tf$residuals, d=1, seasonal=FALSE)
arima_Nt
## Series: pre_tf$residuals 
## ARIMA(4,1,0) 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4
##       -0.7880  -0.5945  -0.4105  -0.2336
## s.e.   0.1095   0.1334   0.1330   0.1081
## 
## sigma^2 = 69828:  log likelihood = -544.07
## AIC=1098.15   AICc=1098.98   BIC=1109.93
lmtest::coeftest(arima_Nt)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.78803    0.10948 -7.1976 6.126e-13 ***
## ar2 -0.59446    0.13342 -4.4556 8.364e-06 ***
## ar3 -0.41054    0.13296 -3.0877  0.002017 ** 
## ar4 -0.23363    0.10808 -2.1616  0.030647 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# GABUNGAN FUNGSI TRANSFER AWAL DENGAN MODEL SISAAN
tf_impor <- tfest(Yt, tf_Xt, um.x = umXt, um.y = umYt_Xt)
noise_um <- um(Yt, ar=c(4), i=c(1), ma=c(0))
tfm_harga <- tfarima::tfm(Yt, inputs=list(tf_impor), noise=noise_um)

1.12 12. Persamaan Model

Menampilkan polinomial lag dari model akhir.

##PERSAMAAN MODEL FUNGSI TRANSFER ----------
printLagpol(tfm_harga$inputs[[1]]$theta)
## 450 + 1300B - 3000B^2
printLagpol(tfm_harga$inputs[[1]]$phi)
## 1 - 0.12B
printLagpol(tfm_harga$noise$phi)
## 1 - 0.31B - 0.1B^2 + 0.037B^3 + 0.081B^4

1.13 13. Uji Diagnostik Model Akhir

Memastikan sisaan model fungsi transfer akhir (\(a_t\)) bersifat white noise.

##UJI DIAGNOSTIK Y ----------

# WHITE NOISE
rY = residuals(tfm_harga)
Box.test(rY, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  rY
## X-squared = 9.5897e-05, df = 1, p-value = 0.9922
# KORELASI SILANG
pccf(rY, tf_Xt, lag.max=36)

Interpretasi: Jika P-value Box-Ljung > 0.05, maka model sudah layak digunakan karena sisaannya acak dan tidak mengandung informasi yang tertinggal.

1.14 14. Validasi Model (Data Testing)

Melakukan peramalan pada data testing dan menghitung akurasi dengan MAPE (Mean Absolute Percentage Error).

##HASIL FORECAST ----------
Ytest <- ts(datatesting$harga)
Xtest <- ts(datatesting$impor)

# FUNGSI MULTI TRANSFER
umx_test <- um(Xtest, ar=0, i=1, ma=2)
umy_test <- fit(umx_test, Ytest)
tf_test <- tfest(Ytest, Xtest, um.x = umx_test, um.y = umy_test)
noise_um_test <- um(Ytest, ar=c(4), i=c(1), ma=c(0))
tfmy_test <- tfarima::tfm(Ytest, inputs=list(tf_test), noise=noise_um_test)
ramal_tfm <- Ytest - residuals(tfmy_test)


##PEMILIHAN MODEL TERBAIK (MAPE) ----------

# PLOT NILAI TESTING DAN PERAMALAN
plot(Ytest, type="l", col="black", xlab=" ", ylab=" "); par(new=TRUE)
plot(ramal_tfm, type="l", col="red", lty=2, main="Harga Beras (Testing)",
     ylab="Rupiah (Rp)", xaxt="n", yaxt="n")
legend("bottomright", 
       c("Data Aktual (Testing)", "Data Peramalan Fungsi Transfer"),
       col=c("black", "red"), lty=c(1,2), cex=0.7)

mape <- function(Y_aktual, Y_forecast){
  rumus = mean(abs(Y_aktual-Y_forecast)/Y_aktual)*100
  print(rumus)
}
mape(Ytest, ramal_tfm)
## [1] 0.9435526

Interpretasi: Nilai MAPE menunjukkan rata-rata kesalahan prediksi dalam persen. Semakin kecil nilai MAPE, semakin akurat model tersebut.

1.15 15. Peramalan Periode Berikutnya

Menggunakan seluruh data untuk memprediksi nilai masa depan (12 periode ke depan).

## PERIODE BERIKUTNYA ----------
outputY <- ts(beras$harga, start = c(2014,1), freq=12)
inputX <- ts(beras$impor, start = c(2014,1), freq=12)

umx_data <- um(inputX, ar=0, i=1, ma=2)
umy_data <- fit(umx_data, outputY)
tf_data <- tfest(outputY, inputX, um.x = umx_data, um.y = umy_data)
noise_um_data <- um(outputY, ar=c(4), i=c(1), ma=c(0))
tfmy_data <- tfarima::tfm(outputY, inputs=list(tf_data), noise=noise_um_data)

HASIL <- predict.tfm(tfmy_data, n.ahead=12)
HASIL
##          Forecast      RMSE   95% LB   95% UB
## Aug 2024 17621.98  237.9332 17155.64 18088.32
## Sep 2024 18152.69  419.1119 17331.25 18974.14
## Oct 2024 18337.37  598.0999 17165.11 19509.62
## Nov 2024 18485.05  753.2710 17008.66 19961.43
## Dec 2024 18583.93  896.9316 16825.97 20341.88
## Jan 2025 18657.18 1029.2330 16639.92 20674.44
## Feb 2025 18702.43 1152.6165 16443.34 20961.52
## Mar 2025 18733.83 1267.4365 16249.70 21217.96
## Apr 2025 18754.17 1374.9586 16059.31 21449.04
## May 2025 18768.27 1476.0114 15875.34 21661.20
## Jun 2025 18777.53 1571.4771 15697.49 21857.57
## Jul 2025 18783.82 1662.0221 15526.32 22041.33

BAGIAN 2: ANALISIS VECTOR AUTOREGRESSION (VAR)

Konsep Dasar: VAR adalah model sistem persamaan di mana setiap variabel diperlakukan sebagai variabel endogen. Model ini sangat berguna untuk melihat hubungan timbal balik antar variabel. Di sini kita menggunakan studi kasus hubungan antara Income dan Savings.

1.16 1. Persiapan Pustaka dan Data

Memuat pustaka vars dan menyiapkan data multivariat.

library(tidyverse)  # read_csv, dplyr, tidyr, ggplot2
library(zoo)        # as.yearqtr, as.Date untuk Quarter
library(tseries)    # adf.test
library(vars)       # VAR, VARselect, IRF, FEVD, uji diagnostik
library(readr)
a <- read_csv(file.choose())
head(a)
quarter_seq <- seq(from = as.yearqtr("1970 Q1"), to = as.yearqtr("2016 Q3"), by = 0.25)

# Gabungkan ke data
a$Quarter <- quarter_seq

# Lihat hasilnya
head(a)

Mengambil subset data agar fokus pada periode tertentu.

# Ambil hanya kolom Quarter, Income, dan Savings
data_var <- a %>%
  dplyr::select(Quarter, Income, Savings)

# Lihat hasilnya
head(data_var)
# Subsetting Data
data_var <- a[50:110, c("Quarter", "Income", "Savings")]

# Cek hasil subset
head(data_var)
tail(data_var)

1.17 2. Eksplorasi Data (Visualisasi)

Visualisasi untuk melihat korelasi visual atau pola pergerakan bersama (co-movement) antara Income dan Savings.

# Plot sederhana: 1 gambar, 2 garis (Income dan Savings)

plot_df <- data_var %>%
mutate(Date = as.Date(Quarter)) %>%
pivot_longer(cols = c(Income, Savings),
names_to = "Variable", values_to = "Value")

ggplot(plot_df, aes(x = Date, y = Value, linetype = Variable)) +
geom_line() +
labs(title = "Income dan Savings per Kuartal",
x = "Tahun",
y = "Pertumbuhan (%)",
linetype = NULL) +
scale_x_date(date_breaks = "2 years", date_labels = "%Y") +
theme_minimal(base_size = 11)

1.18 3. Uji Stasioneritas Data

VAR mensyaratkan semua variabel dalam sistem harus stasioner.

adf_result <- adf.test(data_var$Income)
print(adf_result)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_var$Income
## Dickey-Fuller = -4.2512, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
adf_result <- adf.test(data_var$Savings)
print(adf_result)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_var$Savings
## Dickey-Fuller = -4.4248, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Interpretasi: Jika P-value < 0.05, maka data stasioner. Jika tidak, perlu dilakukan differencing sebelum masuk ke pemodelan VAR.

1.19 4. Penentuan Panjang Lag Optimal

Menentukan berapa periode masa lalu yang berpengaruh signifikan. Kriteria yang umum digunakan adalah AIC (Akaike), HQ (Hannan-Quinn), dan SC (Schwarz).

a <- data.frame(Income = data_var$Income, Savings = data_var$Savings)
a
lag <- VARselect(a, lag.max = 15, type="const")
lag$criteria
##                1         2         3         4         5         6         7
## AIC(n)  2.437840  2.509547  2.485245  2.611974  2.643660  2.672360  2.788128
## HQ(n)   2.527190  2.658464  2.693729  2.880026  2.971278  3.059546  3.234881
## SC(n)   2.676358  2.907078  3.041788  3.327530  3.518227  3.705940  3.980721
## FPE(n) 11.452527 12.320586 12.061391 13.765925 14.332711 14.939153 17.073951
##                8         9        10        11        12        13        14
## AIC(n)  2.940213  2.949924  2.957955  2.762453  2.761289  2.663648  2.694729
## HQ(n)   3.446533  3.515810  3.583408  3.447473  3.505876  3.467802  3.558450
## SC(n)   4.291818  4.460541  4.627584  4.591094  4.748942  4.810314  5.000407
## FPE(n) 20.361185 21.218026 22.275818 19.291750 20.567055 20.244083 23.139452
##               15
## AIC(n)  2.810358
## HQ(n)   3.733646
## SC(n)   5.275048
## FPE(n) 29.553906
lag$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      1      1      1      1

Interpretasi: Pilih lag yang paling banyak direkomendasikan oleh kriteria informasi (misal SC menyarankan lag 1).

1.20 5. Estimasi Model VAR

Mengestimasi parameter model VAR dengan lag terpilih (misal \(p=1\)).

model1 <- vars::VAR(a, p = 1, type="const")
summary(model1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Income, Savings 
## Deterministic variables: const 
## Sample size: 60 
## Log Likelihood: -240.765 
## Roots of the characteristic polynomial:
## 0.4129 0.2434
## Call:
## vars::VAR(y = a, p = 1, type = "const")
## 
## 
## Estimation results for equation Income: 
## ======================================= 
## Income = Income.l1 + Savings.l1 + const 
## 
##            Estimate Std. Error t value Pr(>|t|)   
## Income.l1   0.56743    0.18013   3.150  0.00260 **
## Savings.l1 -0.03911    0.01313  -2.980  0.00424 **
## const       0.35742    0.15778   2.265  0.02731 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.5926 on 57 degrees of freedom
## Multiple R-Squared:  0.16,   Adjusted R-squared: 0.1305 
## F-statistic: 5.428 on 2 and 57 DF,  p-value: 0.006952 
## 
## 
## Estimation results for equation Savings: 
## ======================================== 
## Savings = Income.l1 + Savings.l1 + const 
## 
##            Estimate Std. Error t value Pr(>|t|)  
## Income.l1    3.2042     2.6048   1.230   0.2237  
## Savings.l1  -0.3979     0.1898  -2.096   0.0405 *
## const       -1.9647     2.2815  -0.861   0.3928  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 8.57 on 57 degrees of freedom
## Multiple R-Squared: 0.07498, Adjusted R-squared: 0.04253 
## F-statistic:  2.31 on 2 and 57 DF,  p-value: 0.1085 
## 
## 
## 
## Covariance matrix of residuals:
##         Income Savings
## Income  0.3512   3.765
## Savings 3.7649  73.436
## 
## Correlation matrix of residuals:
##         Income Savings
## Income  1.0000  0.7414
## Savings 0.7414  1.0000

1.21 6. Uji Stabilitas Model

Sistem VAR harus stabil agar analisis selanjutnya valid.

roots(model1)
## [1] 0.4128588 0.2433829

Interpretasi: Semua akar polinomial karakteristik harus memiliki modulus < 1 (berada di dalam lingkaran satuan) agar model stabil.

1.22 7. Uji Kausalitas Granger

Menguji apakah satu variabel berguna untuk memprediksi variabel lain.

causality(model1, cause = "Income")
## $Granger
## 
##  Granger causality H0: Income do not Granger-cause Savings
## 
## data:  VAR object model1
## F-Test = 1.5132, df1 = 1, df2 = 114, p-value = 0.2212
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Income and Savings
## 
## data:  VAR object model1
## Chi-squared = 21.281, df = 1, p-value = 3.967e-06
causality(model1, cause = "Savings")
## $Granger
## 
##  Granger causality H0: Savings do not Granger-cause Income
## 
## data:  VAR object model1
## F-Test = 8.8775, df1 = 1, df2 = 114, p-value = 0.003529
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Savings and Income
## 
## data:  VAR object model1
## Chi-squared = 21.281, df = 1, p-value = 3.967e-06

Interpretasi: Jika P-value < 0.05, maka tolak \(H_0\). Artinya, terdapat hubungan kausalitas (prediktif) dari variabel penyebab ke variabel akibat.

1.23 8. Uji Diagnostik Residual

Memastikan residual model memenuhi asumsi: 1. Serial Test: Tidak ada autokorelasi. 2. ARCH Test: Ragam residual konstan (homoskedastis). 3. Normality Test: Residual berdistribusi normal.

serial.test(model1, lags.pt = 12, type = "PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object model1
## Chi-squared = 43.612, df = 44, p-value = 0.4882
arch.test(model1, lags.multi = 12, multivariate.only = TRUE)
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object model1
## Chi-squared = 127.48, df = 108, p-value = 0.09713
normality.test(model1, multivariate.only = TRUE)
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object model1
## Chi-squared = 3.0729, df = 4, p-value = 0.5457
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object model1
## Chi-squared = 1.0472, df = 2, p-value = 0.5924
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object model1
## Chi-squared = 2.0256, df = 2, p-value = 0.3632

1.24 9. Impulse Response Function (IRF)

Konsep: IRF menunjukkan respons variabel sistem terhadap guncangan (shock) satu standar deviasi pada salah satu variabel dari waktu ke waktu.

# IRF untuk shock pada Income dan Savings, respon keduanya
plot(irf(model1, n.ahead = 10))

Interpretasi: * Grafik menunjukkan bagaimana variabel merespon guncangan. * Jika garis nol berada di luar interval kepercayaan (garis merah), maka respon tersebut signifikan secara statistik.

1.25 10. Forecast Error Variance Decomposition (FEVD)

Konsep: FEVD menjelaskan seberapa besar persentase varians error suatu variabel yang dijelaskan oleh dirinya sendiri vs variabel lain dalam sistem.

fevd_model1 <- fevd(model1, n.ahead = 10)
fevd_model1    # hasil numerik
## $Income
##          Income   Savings
##  [1,] 1.0000000 0.0000000
##  [2,] 0.8764410 0.1235590
##  [3,] 0.8750261 0.1249739
##  [4,] 0.8733945 0.1266055
##  [5,] 0.8732717 0.1267283
##  [6,] 0.8732361 0.1267639
##  [7,] 0.8732315 0.1267685
##  [8,] 0.8732306 0.1267694
##  [9,] 0.8732305 0.1267695
## [10,] 0.8732305 0.1267695
## 
## $Savings
##          Income   Savings
##  [1,] 0.5496192 0.4503808
##  [2,] 0.5154674 0.4845326
##  [3,] 0.5169585 0.4830415
##  [4,] 0.5167083 0.4832917
##  [5,] 0.5167274 0.4832726
##  [6,] 0.5167251 0.4832749
##  [7,] 0.5167253 0.4832747
##  [8,] 0.5167253 0.4832747
##  [9,] 0.5167253 0.4832747
## [10,] 0.5167253 0.4832747
# Plot FEVD

plot(fevd_model1)

Interpretasi: Jika grafik batang suatu variabel (misal Savings) didominasi oleh warna variabel lain (Income), berarti variabel tersebut sangat dipengaruhi oleh variabel lain dalam sistem (bersifat endogen).

1.26 11. Peramalan (Forecasting)

Memprediksi nilai masa depan (8 kuartal ke depan).

# Misal meramalkan 8 kuartal ke depan

h <- 8
fcst <- predict(model1, n.ahead = h, ci = 0.95)
fcst
## $Income
##           fcst      lower    upper       CI
## [1,] 0.5987849 -0.5627170 1.760287 1.161502
## [2,] 0.7634897 -0.4907230 2.017703 1.254213
## [3,] 0.7660564 -0.4988190 2.030932 1.264875
## [4,] 0.7830414 -0.4838127 2.049895 1.266854
## [5,] 0.7861778 -0.4809796 2.053335 1.267157
## [6,] 0.7884161 -0.4787954 2.055628 1.267212
## [7,] 0.7891106 -0.4781099 2.056331 1.267220
## [8,] 0.7894532 -0.4777688 2.056675 1.267222
## 
## $Savings
##            fcst     lower    upper       CI
## [1,] -1.6952881 -18.49122 15.10064 16.79593
## [2,]  0.6285564 -16.79974 18.05686 17.42830
## [3,]  0.2315312 -17.23187 17.69493 17.46340
## [4,]  0.3977513 -17.07003 17.86553 17.46778
## [5,]  0.3860275 -17.08215 17.85420 17.46817
## [6,]  0.4007428 -17.06750 17.86899 17.46825
## [7,]  0.4020586 -17.06620 17.87032 17.46826
## [8,]  0.4037603 -17.06450 17.87202 17.46826

Visualisasi Hasil Forecast:

# === Plot Hasil Peramalan VAR ===

# Ambil hasil forecast dari objek fcst
fcst_income  <- as.data.frame(fcst$fcst$Income)
fcst_savings <- as.data.frame(fcst$fcst$Savings)

# Tambahkan periode waktu ke depan sesuai dengan data terakhir
last_date <- max(as.Date(data_var$Quarter))
future_quarters <- seq(last_date, by = "quarter", length.out = nrow(fcst_income) + 1)[-1]

# Plot Income Forecast
plot(future_quarters, fcst_income$fcst, type = "l", lwd = 2,
     ylim = range(c(fcst_income$lower, fcst_income$upper)),
     main = "Forecast Income (8 Kuartal ke Depan)",
     xlab = "Tahun", ylab = "Income (Forecast)")
lines(future_quarters, fcst_income$lower, lty = 2, col = "red")
lines(future_quarters, fcst_income$upper, lty = 2, col = "red")
points(future_quarters, fcst_income$fcst, pch = 16, col = "blue")
legend("topleft", legend = c("Prediksi", "Interval 95%"),
       col = c("blue", "red"), lty = c(1, 2), bty = "n")

# Plot Savings Forecast
plot(future_quarters, fcst_savings$fcst, type = "l", lwd = 2,
     ylim = range(c(fcst_savings$lower, fcst_savings$upper)),
     main = "Forecast Savings (8 Kuartal ke Depan)",
     xlab = "Tahun", ylab = "Savings (Forecast)")
lines(future_quarters, fcst_savings$lower, lty = 2, col = "red")
lines(future_quarters, fcst_savings$upper, lty = 2, col = "red")
points(future_quarters, fcst_savings$fcst, pch = 16, col = "blue")
legend("topleft", legend = c("Prediksi", "Interval 95%"),
       col = c("blue", "red"), lty = c(1, 2), bty = "n")

1.27 12. Evaluasi Model (MAPE In-Sample)

Mengukur keakuratan model dalam mencocokkan data historis.

# === MAPE In-Sample ===
# Ambil fitted value dari model VAR
fitted_values <- fitted(model1)

# Sesuaikan jumlah observasi aktual dengan fitted values
actual_values <- data_var[(nrow(data_var) - nrow(fitted_values) + 1):nrow(data_var), ]

# Hitung MAPE per variabel
mape_income  <- mean(abs((actual_values$Income  - fitted_values[, "Income"])  / actual_values$Income))  * 100
mape_savings <- mean(abs((actual_values$Savings - fitted_values[, "Savings"]) / actual_values$Savings)) * 100

cat("MAPE In-Sample Income :", round(mape_income, 3), "%\n")
## MAPE In-Sample Income : 109.003 %
cat("MAPE In-Sample Savings:", round(mape_savings, 3), "%\n")
## MAPE In-Sample Savings: 259.261 %
# === Plot Aktual vs Fitted (In-Sample) ===
# Ambil lagi fitted value dari model VAR
fitted_values <- fitted(model1)

# Sesuaikan indeks supaya panjangnya sama
idx_start <- nrow(data_var) - nrow(fitted_values) + 1

fit_df <- data.frame(
Date           = as.Date(data_var$Quarter[idx_start:nrow(data_var)]),
Income_Actual  = data_var$Income[idx_start:nrow(data_var)],
Income_Fitted  = fitted_values[, "Income"],
Savings_Actual = data_var$Savings[idx_start:nrow(data_var)],
Savings_Fitted = fitted_values[, "Savings"]
)

# Ubah ke format long untuk ggplot
fit_long <- fit_df %>%
pivot_longer(cols = -Date,
names_to = c("Variable", "Series"),
names_sep = "_",
values_to = "Value")

# Plot: Actual vs Fitted, difacet per variabel
ggplot(fit_long, aes(x = Date, y = Value, linetype = Series)) +
geom_line() +
facet_wrap(~ Variable, ncol = 1, scales = "free_y") +
labs(title    = "Aktual vs Fitted (In-Sample) - VAR",
x        = "Tahun",
y        = "Nilai",
linetype = NULL) +
scale_x_date(date_breaks = "2 years", date_labels = "%Y") +
theme_minimal(base_size = 11)

BAGIAN 3: ANALISIS SINGULAR SPECTRUM ANALYSIS (SSA)

Konsep Dasar: SSA adalah metode non-parametrik yang kuat untuk membedah data deret waktu menjadi komponen: Tren, Osilasi (Musiman), dan Noise. Metode ini tidak memerlukan asumsi distribusi data (seperti normalitas) dan bekerja dengan mengubah deret waktu menjadi matriks trajektori.

1.28 1. Persiapan Data

Data dummy digunakan untuk mendemonstrasikan algoritma SSA secara manual.

rm(list=ls())   # remove variables
graphics.off()  # close figures

# Input Data Contoh (1 s.d 15)
# Menggantikan read_excel sesuai permintaan
DataSSA <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)

Y <- as.matrix(DataSSA)
N <- length(Y)
L = 7           # Trial And error (L =< N/2)
K = N - L + 1

# Menampilkan data awal
print(t(Y))
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,]    1    2    3    4    5    6    7    8    9    10    11    12    13    14
##      [,15]
## [1,]    15

1.29 2. Tahap Embedding

Konsep: Embedding mengubah data deret waktu 1D menjadi matriks 2D yang disebut Matriks Hankel (Matriks Trajektori). Elemen diagonal dari matriks ini bernilai konstan. * \(L\) (Window Length): Menentukan seberapa besar “jendela” pengamatan pola.

# 1. Embedding
# To produce the Hankel matrix
X <- outer((1:L), (1:K), function(x,y) Y[(x+y-1)]) 

# Menampilkan sebagian matriks Hankel
print(head(X))
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    1    2    3    4    5    6    7    8    9
## [2,]    2    3    4    5    6    7    8    9   10
## [3,]    3    4    5    6    7    8    9   10   11
## [4,]    4    5    6    7    8    9   10   11   12
## [5,]    5    6    7    8    9   10   11   12   13
## [6,]    6    7    8    9   10   11   12   13   14

1.30 3. Tahap Dekomposisi (SVD)

Konsep: Matriks Trajektori dipecah menggunakan Singular Value Decomposition (SVD) menjadi Eigenvalues (\(\lambda\)) dan Eigenvectors. * Eigenvalues: Menunjukkan “energi” sinyal. Nilai yang besar biasanya merepresentasikan tren utama, sedangkan nilai kecil merepresentasikan noise.

# 2. SVD
# Hankelization R code function definition
UniHankel <- function(Y,L){
  k <- length(Y)-L+1
  outer((1:L), (1:k), function(x,y) Y[(x+y-1)])
}

# SVD R code function definition
SVD <- function(Y,L){
  X <- UniHankel(Y,L)
  svd(X)
}

# Eksekusi SVD pada Matriks X
SVD_result <- svd(X)

lambda <- SVD_result$d    # to see the eigenvalue
sing.value <- sqrt(lambda)
Lambda <- diag(lambda)    # Make singular matrix
U <- SVD_result$u
V <- SVD_result$v

# Menampilkan Eigenvalues
print(lambda)
## [1] 6.842069e+01 4.754857e+00 6.924323e-15 1.160704e-15 7.087244e-16
## [6] 4.326742e-16 1.169336e-16

1.31 4. Tahap Rekonstruksi

Tahap ini mengembalikan pecahan matriks tadi menjadi bentuk deret waktu melalui proses Grouping dan Diagonal Averaging.

1.31.1 Grouping

Konsep: Kita memilih komponen-komponen tertentu (berdasarkan besaran Eigenvalue) untuk direkonstruksi. Misalnya, memilih komponen ke-1 saja untuk mengambil tren, dan membuang sisanya sebagai noise.

# 1. Grouping
# Grouping R Code
Group <- function(Y, L, groups){
  I <- groups; p <- length(I)
  # Memanggil fungsi SVD yang didefinisikan sebelumnya
  SVD_res <- SVD(Y, L) 
  LambdaI <- matrix(diag(SVD_res$d)[I,I], p, p)
  SVD_res$u[,I] %*% LambdaI %*% t(SVD_res$v[,I])
}

# Diasumsikan menggunakan komponen pertama -> I <- c(1)
I <- c(1)
p <- length(I)

# Perhitungan manual komponen matriks
U[,I] %*% matrix(Lambda[I,I], p, p)
##          [,1]
## [1,] 16.60829
## [2,] 19.47749
## [3,] 22.34669
## [4,] 25.21589
## [5,] 28.08509
## [6,] 30.95429
## [7,] 33.82349
XI <- U[,I] %*% matrix(Lambda[I,I], p, p) %*% t(V[,I])   # Compute the SVD using component 1

print(XI[1:5, 1:5]) # Menampilkan sebagian hasil matriks rekonstruksi
##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 2.789869 3.416082 4.042295 4.668508 5.294721
## [2,] 3.271838 4.006234 4.740630 5.475026 6.209422
## [3,] 3.753808 4.596386 5.438965 6.281544 7.124123
## [4,] 4.235777 5.186539 6.137300 7.088062 8.038823
## [5,] 4.717747 5.776691 6.835635 7.894580 8.953524

1.31.2 Diagonal Averaging

Konsep: Matriks hasil rekonstruksi (\(XI\)) ditransformasikan kembali menjadi deret waktu satu dimensi dengan merata-ratakan nilai pada arah anti-diagonal (\(i+j=\text{konstan}\)).

# 2. Diagonal Averaging
# Diagonal Averaging R Code
DiagAver <- function(X){
  L <- nrow(X); k <- ncol(X); N <- k + L - 1
  D <- NULL
  for(j in 1:N){
    s1 <- max(1, (j - N + L))
    s2 <- min(L, j)
    place <- (s1:s2) + L * (((j + 1 - s1):(j + 1 - s2)) - 1)
    D[j] <- mean(X[place])
  }
  D
}

# SSA Reconstruction R Code (Fungsi gabungan)
SSA.Rec <- function(Y, L, groups){
  N <- length(Y)
  I <- groups; p <- length(I)
  XI <- Group(Y, L, groups)
  Approx <- DiagAver(XI)
  Resid <- Y - Approx
  list(Approximation = Approx, Residual = Resid)
}

# Implementasi manual loop Diagonal Averaging untuk XI
D <- NULL
N <- length(Y)
for(t in 1:N){
  s1 <- max(1, (t - N + L))
  s2 <- min(L, t)
  place <- (s1:s2) + L * (((t + 1 - s1):(t + 1 - s2)) - 1)
  D[t] <- mean(XI[place])
}

# Menampilkan hasil rekonstruksi (dibulatkan)
round(D, 2)
##  [1]  2.79  3.34  3.93  4.56  5.22  5.92  6.66  7.61  8.56  9.69 10.86 12.06
## [13] 13.30 14.57 15.88

1.32 5. Forecasting (SSA)

Konsep: Peramalan dilakukan menggunakan Linear Recurrent Formula (LRF). Metode ini memproyeksikan titik data masa depan sebagai kombinasi linear dari titik data sebelumnya, berdasarkan struktur eigen yang telah diekstraksi.

# Forecasting
# SSA LRF Forecasting R Code:
RSSA.Forecasting <- function(L, groups, h, Y){
  N <- length(Y)
  L <- min(L, (N - L + 1))
  X <- UniHankel(Y, L)
  U <- matrix(svd(X)$u, L, L)
  pi <- array(U[L, groups], dim = length(groups))
  V2 <- sum(pi^2)
  m <- length(groups)
  Udelta <- array(U[1:(L-1), groups], dim = c((L-1), m))
  A <- pi %*% t(Udelta) / (1 - V2)
  yhat <- array(0, dim = (N + h))
  yhat[1:N] <- SSA.Rec(Y, L, groups)$Approximation
  for(l in (N+1):(N+h))
    yhat[l] <- A %*% yhat[(l - L + 1):(l - 1)]
  yhat[(N+1):(N+h)]
}

# Menjalankan fungsi forecasting
# Parameter: L=7, Group=1, h=3 (forecast 3 periode ke depan)
Forecast_Result <- RSSA.Forecasting(L=7, groups=c(1), h=3, Y=Y)

print("Hasil Peramalan 3 periode ke depan:")
## [1] "Hasil Peramalan 3 periode ke depan:"
print(Forecast_Result)
## [1] 17.95759 19.90515 22.04894

2 KESIMPULAN AKHIR

Berdasarkan ketiga metode analisis yang telah dilakukan, dapat ditarik kesimpulan sebagai berikut:

  1. Analisis Fungsi Transfer: Metode ini sangat efektif ketika hubungan sebab-akibat (kausalitas) satu arah sudah diketahui dengan jelas, misalnya dampak Impor terhadap Harga. Kunci keberhasilan metode ini terletak pada proses pre-whitening untuk membersihkan autokorelasi input sehingga hubungan murni dapat teridentifikasi.
  2. Analisis VAR: Metode ini memberikan pandangan sistemik di mana variabel-variabel (seperti Income dan Savings) diperlakukan setara dan saling mempengaruhi. Keunggulannya terletak pada Impulse Response Function (IRF) yang memungkinkan kita melihat dampak guncangan (shock) dari satu variabel ke variabel lain dalam sistem.
  3. Analisis SSA: Metode ini menawarkan pendekatan non-parametrik yang kuat untuk membedah sinyal (tren/musiman) dari noise tanpa asumsi distribusi data yang ketat. Ini sangat berguna untuk data yang memiliki pola kompleks atau pendek, di mana metode parametrik klasik mungkin kurang optimal.

Pemilihan metode yang tepat sangat bergantung pada tujuan analisis: apakah untuk kontrol input-output (Fungsi Transfer), analisis kebijakan makroekonomi (VAR), atau dekomposisi sinyal data (SSA).

BAGIAN TAMBAHAN: KIAT & STRATEGI PENGERJAAN MANUAL (UAS)

Bagian ini disusun khusus untuk persiapan ujian tulis (tanpa laptop). Fokus utama adalah pemahaman logika, penulisan rumus manual, dan interpretasi output yang mungkin dicetak di soal.

2.1 1. STRATEGI FUNGSI TRANSFER (TF)

Tantangan Ujian: Biasanya diberikan plot CCF dan diminta menentukan orde \((b, r, s)\), atau diberikan persamaan dan diminta menghitung nilai ramalan manual.

2.1.1 A. Cara Membaca Plot CCF (Tanpa Komputer)

Jika soal memberikan grafik Cross Correlation Function (CCF) setelah pre-whitening: 1. Cari \(b\) (Delay): Lihat pada lag negatif (kiri). Di lag berapa batang pertama kali keluar garis putus-putus (signifikan)? * Contoh: Batang signifikan pertama di Lag -3 \(\rightarrow\) \(b = 3\). 2. Cari \(r\) (Numerator): Lihat berapa banyak batang signifikan setelah lag \(b\). * Rumus: Jumlah batang signifikan berurutan setelah \(b\). 3. Cari \(s\) (Denominator): Lihat pola batang setelah \(b+r\). * Jika mengecil secara eksponensial (melengkung turun) atau sinus (gelombang) \(\rightarrow\) \(s = 1\) atau \(s = 2\). * Jika langsung putus (tidak ada pola) \(\rightarrow\) \(s = 0\).

2.1.2 B. Menulis Persamaan Manual

Jangan terkecoh notasi. Ingat bentuk umum: \[Y_t = \text{Bagian Input} + \text{Bagian Noise}\] \[Y_t = \left( \frac{\omega_0 - \omega_1 B - \dots}{1 - \delta_1 B - \dots} \right) X_{t-b} + \frac{\theta(B)}{\phi(B)} a_t\]

Tips Menjawab: * Jika diminta forecast manual \(Y_{t+1}\): Masukkan nilai \(X\) dari \(b\) periode lalu. * Contoh (\(b=1\)): Harga bulan depan dipengaruhi Impor bulan ini.

2.2 2. STRATEGI VECTOR AUTOREGRESSION (VAR)

Tantangan Ujian: Menuliskan sistem persamaan matriks dan interpretasi hasil (Granger/IRF) dari tabel/gambar. Perhitungan manual biasanya terbatas pada VAR(1) dengan 2 variabel.

2.2.1 A. Menulis Persamaan VAR(1) Manual

Jika ada 2 variabel, misal \(Y_1\) (Income) dan \(Y_2\) (Savings), model VAR(1) bukan 1 persamaan, tapi 2 persamaan:

  1. Persamaan Income: \[Y_{1,t} = c_1 + \phi_{11}Y_{1,t-1} + \phi_{12}Y_{2,t-1} + e_{1,t}\] (Income dipengaruhi Income lalu DAN Savings lalu)

  2. Persamaan Savings: \[Y_{2,t} = c_2 + \phi_{21}Y_{1,t-1} + \phi_{22}Y_{2,t-1} + e_{2,t}\] (Savings dipengaruhi Income lalu DAN Savings lalu)

Tips Menjawab: * Perhatikan koefisien \(\phi_{12}\) dan \(\phi_{21}\). Jika \(\phi_{12}\) signifikan (tidak nol), berarti Savings mempengaruhi Income. Ini dasar Granger Causality.

2.2.2 B. Interpretasi Output (Tanpa Running Code)

  • Akar (Roots): Jika soal memberi angka modulus (misal: 0.8, 0.5), jawab “STABIL” karena < 1. Jika ada angka 1.05, jawab “TIDAK STABIL”.
  • IRF Graph:
    • Garis menyentuh 0 = Pengaruh hilang.
    • Garis jauh di atas/bawah 0 = Pengaruh kuat.
    • Lihat periode (sumbu X) kapan garis kembali ke 0. Itu durasi dampak guncangan.

2.3 3. STRATEGI SINGULAR SPECTRUM ANALYSIS (SSA)

Tantangan Ujian: Ini yang paling sering keluar hitungan manualnya. Biasanya diminta membuat Matriks Trajektori (Hankel) dan Diagonal Averaging untuk data kecil (N=5 atau N=10).

2.3.1 A. Langkah Hitung Manual (Wajib Hafal)

Misal Data: \(Y = \{1, 2, 3, 4, 5\}\), Panjang Jendela \(L=3\).

Langkah 1: Tentukan K \[N=5, L=3 \rightarrow K = N - L + 1 = 5 - 3 + 1 = 3\] Matriks akan berukuran \(L \times K\) (\(3 \times 3\)).

Langkah 2: Buat Matriks Trajektori (Embedding) Isi kolom pertama dengan data 1 sampai \(L\). Geser satu langkah untuk kolom kedua, dst. \[ X = \begin{bmatrix} 1 & 2 & 3 \\ 2 & 3 & 4 \\ 3 & 4 & 5 \end{bmatrix} \] Ciri Benar: Diagonal miring (kiri bawah ke kanan atas) angkanya sama (1, lalu 2-2, lalu 3-3-3…).

Langkah 3: Diagonal Averaging (Jika diberikan matriks hasil rekonstruksi) Misal dikasih matriks rekonstruksi \(Y_{rec}\): \[ Y_{rec} = \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \end{bmatrix} \] Cara mengembalikan ke deret angka: * \(y_1 = a\) * \(y_2 = (d+b)/2\) * \(y_3 = (g+e+c)/3\) * \(y_4 = (h+f)/2\) * \(y_5 = i\)

Tips Kalkulator: * Latihan merata-ratakan diagonal matriks dengan cepat. * Ingat pola indeks: \(y_k\) adalah rata-rata elemen matriks dimana indeks baris + kolom = \(k+1\).

2.3.2 B. Konsep SVD (Teori)

Jika ditanya konsep SVD tanpa hitungan: \[X = \sum \sqrt{\lambda_i} U_i V_i^T\] * \(\lambda\) (Eigenvalue): “Kekuatan” sinyal. Urutkan dari terbesar. * Grouping: Ambil \(\lambda\) besar untuk Tren. \(\lambda\) yang berpasangan biasanya Musiman. \(\lambda\) kecil adalah Noise.

3 TIPS UMUM MENJAWAB SOAL ESAI

  1. Stasioneritas adalah Kunci: Setiap kali ditanya langkah awal, jawaban wajibnya adalah “Cek Stasioneritas (Mean & Varians)”. Jika tidak stasioner \(\rightarrow\) Differencing/Box-Cox.
  2. Alasan Pemilihan Model:
    • Pilih TF jika tahu mana penyebab (X) dan akibat (Y).
    • Pilih VAR jika X dan Y saling mempengaruhi (ayam atau telur duluan?).
    • Pilih SSA jika datanya aneh, pendek, atau pola musiman tidak beraturan (non-parametrik).
  3. Interpretasi MAPE: Jangan cuma tulis angka. Tulis: “Model memiliki rata-rata kesalahan sebesar X%, yang berarti akurasinya (Sangat Baik/Baik/Cukup)”.
    • < 10% : Sangat Baik
    • 10-20% : Baik
    • 20-50% : Cukup

Semangat UAS Analisis Data Deret Waktu 2! Jangan lupa bawa kalkulator dan penggaris untuk menggambar plot manual jika diminta.