## ============================================================
## UAS ANALISIS RUNTUN WAKTU TERAPAN — STUDI KASUS 2
## ============================================================
library(urca)
library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: lmtest
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
# Tanggal : 2011-09-01, ... (tanggal 1 tiap bulan)
# Inflasi : inflasi y-on-y (%), diolah dari inflasi m-to-m BPS
# SukuBunga : BI Rate s.d. Jul 2016, lalu BI-7DRR (%)
# Kurs : kurs tengah Rp/USD rata-rata bulanan (Bank Indonesia)
# D_7DRR : 1 mulai Agustus 2016 (pergantian BI Rate -> BI-7DRR)
df <- read.csv("data_makro.csv")
df$Tanggal <- as.Date(df$Tanggal)
df <- df[order(df$Tanggal), ]
nrow(df)
## [1] 180
thn <- as.numeric(format(df$Tanggal[1], "%Y"))
bln <- as.numeric(format(df$Tanggal[1], "%m"))
sukubunga <- ts(df$SukuBunga, start = c(thn, bln), frequency = 12)
ln_kurs <- ts(100 * log(df$Kurs), start = c(thn, bln), frequency = 12) # 100 x ln
inflasi <- ts(df$Inflasi, start = c(thn, bln), frequency = 12)
# URUTAN KOLOM = URUTAN CHOLESKY (sukubunga -> kurs -> inflasi)
data_ts <- cbind(sukubunga, ln_kurs, inflasi)
colnames(data_ts) <- c("sukubunga", "ln_kurs", "inflasi")
# Dummy perubahan instrumen kebijakan BI
dum_step <- matrix(df$D_7DRR, ncol = 1, dimnames = list(NULL, "D_7DRR")) # untuk VAR level
dum_impulse <- matrix(c(0, diff(df$D_7DRR)), ncol = 1, dimnames = list(NULL, "I_Aug2016")) # untuk VECM
plot(data_ts, main = "Suku Bunga BI, 100 x Ln Kurs Rp/USD, dan Inflasi y-on-y")
summary(df[, c("SukuBunga", "Kurs", "Inflasi")])
## SukuBunga Kurs Inflasi
## Min. :3.500 Min. : 8766 Min. :-0.100
## 1st Qu.:4.750 1st Qu.:13105 1st Qu.: 2.640
## Median :5.750 Median :14114 Median : 3.385
## Mean :5.542 Mean :13737 Mean : 3.797
## 3rd Qu.:6.000 3rd Qu.:15139 3rd Qu.: 4.562
## Max. :7.750 Max. :18008 Max. : 8.790
uji_akar_unit <- function(x, nama) {
data.frame(Variabel = nama,
ADF_level = round(adf.test(x)$p.value, 4),
ADF_diff = round(adf.test(diff(x))$p.value, 4),
PP_level = round(pp.test(x)$p.value, 4),
PP_diff = round(pp.test(diff(x))$p.value, 4))
}
tabel_akar_unit <- rbind(uji_akar_unit(sukubunga, "Suku bunga"),
uji_akar_unit(ln_kurs, "Ln kurs"),
uji_akar_unit(inflasi, "Inflasi"))
## Warning in pp.test(diff(x)): p-value smaller than printed p-value
## Warning in adf.test(diff(x)): p-value smaller than printed p-value
## Warning in pp.test(diff(x)): p-value smaller than printed p-value
## Warning in adf.test(diff(x)): p-value smaller than printed p-value
## Warning in pp.test(diff(x)): p-value smaller than printed p-value
tabel_akar_unit
## Variabel ADF_level ADF_diff PP_level PP_diff
## 1 Suku bunga 0.4209 0.018 0.7798 0.01
## 2 Ln kurs 0.2505 0.010 0.5913 0.01
## 3 Inflasi 0.1479 0.010 0.2670 0.01
lag_select <- VARselect(data_ts, lag.max = 12, type = "const", exogen = dum_step)
lag_select$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 3 2 3
# Untuk lag K = 2 dan K = 3: rank kointegrasi (trace & eigen, 5%)
# dan uji autokorelasi residual VECM (Portmanteau).
rank_johansen <- function(jo) {
cv <- jo@cval[, "5pct"]; st <- jo@teststat
n <- length(st); r <- 0
for (i in n:1) { if (st[i] > cv[i]) r <- r + 1 else break } # baca dari r = 0
r
}
bandingkan <- data.frame(K = 2:3, ecdet = "const", stringsAsFactors = FALSE)
hasil_spes <- do.call(rbind, lapply(seq_len(nrow(bandingkan)), function(i) {
K <- bandingkan$K[i]; ec <- bandingkan$ecdet[i]
jt <- ca.jo(data_ts, type = "trace", ecdet = ec, K = K, spec = "transitory", dumvar = dum_impulse)
je <- ca.jo(data_ts, type = "eigen", ecdet = ec, K = K, spec = "transitory", dumvar = dum_impulse)
rt <- rank_johansen(jt); re <- rank_johansen(je)
r <- max(1, min(rt, re))
vv <- vec2var(jt, r = r)
data.frame(K = K, ecdet = ec, r_trace = rt, r_eigen = re, r_dipakai = r,
Portmanteau_p = round(serial.test(vv, lags.pt = 12, type = "PT.asymptotic")$serial$p.value, 4),
ARCH_p = round(arch.test(vv, lags.multi = 5)$arch.mul$p.value, 4),
JB_p = round(normality.test(vv)$jb.mul$JB$p.value, 4))
}))
hasil_spes
## K ecdet r_trace r_eigen r_dipakai Portmanteau_p ARCH_p JB_p
## Chi-squared 2 const 1 1 1 0.1013 0.0002 0
## Chi-squared1 3 const 3 0 1 0.1859 0.0112 0
# Kriteria pilih: (1) rank trace & eigen konsisten, (2) Portmanteau_p > 0,05,
# (3) lag sesuai kriteria informasi. Hasil: K = 2
# (r trace = r eigen = 1, Portmanteau p = 0,10, lag pilihan SC).
# K = 3 ditolak: trace memberi r = 3 (full rank, bertentangan dengan I(1))
# sedangkan eigen memberi r = 0.
K_pilih <- 2
ecdet_pilih <- "const"
jo_trace <- ca.jo(data_ts, type = "trace", ecdet = ecdet_pilih, K = K_pilih,
spec = "transitory", dumvar = dum_impulse)
jo_eigen <- ca.jo(data_ts, type = "eigen", ecdet = ecdet_pilih, K = K_pilih,
spec = "transitory", dumvar = dum_impulse)
summary(jo_trace)
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 1.247961e-01 5.068796e-02 3.815853e-02 1.110223e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 2 | 6.93 7.52 9.24 12.97
## r <= 1 | 16.18 17.85 19.96 24.60
## r = 0 | 39.91 32.00 34.91 41.07
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## sukubunga.l1 ln_kurs.l1 inflasi.l1 constant
## sukubunga.l1 1.00000000 1.000000 1.00000000 1.0000000
## ln_kurs.l1 -0.06308535 -1.154511 0.02988741 -0.4341021
## inflasi.l1 -1.35082836 1.105225 -0.07857562 0.9165834
## constant 59.54308910 1113.744084 -33.58489848 396.5043886
##
## Weights W:
## (This is the loading matrix)
##
## sukubunga.l1 ln_kurs.l1 inflasi.l1 constant
## sukubunga.d -0.01820402 2.503685e-05 -0.01971871 -5.820781e-16
## ln_kurs.d -0.04570161 1.278014e-02 -0.03177282 5.013580e-16
## inflasi.d 0.08165053 4.671031e-04 -0.06863250 6.084860e-15
summary(jo_eigen)
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 1.247961e-01 5.068796e-02 3.815853e-02 1.110223e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 2 | 6.93 7.52 9.24 12.97
## r <= 1 | 9.26 13.75 15.67 20.20
## r = 0 | 23.73 19.77 22.00 26.81
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## sukubunga.l1 ln_kurs.l1 inflasi.l1 constant
## sukubunga.l1 1.00000000 1.000000 1.00000000 1.0000000
## ln_kurs.l1 -0.06308535 -1.154511 0.02988741 -0.4341021
## inflasi.l1 -1.35082836 1.105225 -0.07857562 0.9165834
## constant 59.54308910 1113.744084 -33.58489848 396.5043886
##
## Weights W:
## (This is the loading matrix)
##
## sukubunga.l1 ln_kurs.l1 inflasi.l1 constant
## sukubunga.d -0.01820402 2.503685e-05 -0.01971871 -5.820781e-16
## ln_kurs.d -0.04570161 1.278014e-02 -0.03177282 5.013580e-16
## inflasi.d 0.08165053 4.671031e-04 -0.06863250 6.084860e-15
r_rank <- 1 # trace & eigen (5%) sama-sama r = 1
vecm <- cajorls(jo_trace, r = r_rank)
cat("\n=== Vektor kointegrasi (dinormalisasi ke sukubunga) ===\n")
##
## === Vektor kointegrasi (dinormalisasi ke sukubunga) ===
print(round(vecm$beta, 4))
## ect1
## sukubunga.l1 1.0000
## ln_kurs.l1 -0.0631
## inflasi.l1 -1.3508
## constant 59.5431
cat("\n=== Persamaan jangka pendek & koefisien penyesuaian (alpha) ===\n")
##
## === Persamaan jangka pendek & koefisien penyesuaian (alpha) ===
print(summary(vecm$rlm))
## Response sukubunga.d :
##
## Call:
## lm(formula = sukubunga.d ~ ect1 + I_Aug2016 + sukubunga.dl1 +
## ln_kurs.dl1 + inflasi.dl1 - 1, data = data.mat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36867 -0.05763 -0.00214 0.03081 0.47991
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ect1 -0.018204 0.006337 -2.873 0.00458 **
## I_Aug2016 -1.198003 0.137464 -8.715 2.30e-15 ***
## sukubunga.dl1 0.382043 0.058029 6.584 5.28e-10 ***
## ln_kurs.dl1 0.009738 0.005616 1.734 0.08468 .
## inflasi.dl1 0.001858 0.019337 0.096 0.92356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1367 on 173 degrees of freedom
## Multiple R-squared: 0.4976, Adjusted R-squared: 0.4831
## F-statistic: 34.27 on 5 and 173 DF, p-value: < 2.2e-16
##
##
## Response ln_kurs.d :
##
## Call:
## lm(formula = ln_kurs.d ~ ect1 + I_Aug2016 + sukubunga.dl1 + ln_kurs.dl1 +
## inflasi.dl1 - 1, data = data.mat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1676 -0.4523 0.2929 1.3143 9.6686
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ect1 -0.04570 0.08508 -0.537 0.59185
## I_Aug2016 0.92886 1.84558 0.503 0.61540
## sukubunga.dl1 0.26556 0.77910 0.341 0.73363
## ln_kurs.dl1 0.23201 0.07539 3.077 0.00243 **
## inflasi.dl1 0.32337 0.25962 1.246 0.21461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.836 on 173 degrees of freedom
## Multiple R-squared: 0.08108, Adjusted R-squared: 0.05453
## F-statistic: 3.053 on 5 and 173 DF, p-value: 0.01146
##
##
## Response inflasi.d :
##
## Call:
## lm(formula = inflasi.d ~ ect1 + I_Aug2016 + sukubunga.dl1 + ln_kurs.dl1 +
## inflasi.dl1 - 1, data = data.mat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.15713 -0.27375 -0.00932 0.24757 2.57115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ect1 0.0816505 0.0245592 3.325 0.00108 **
## I_Aug2016 -0.5127978 0.5327359 -0.963 0.33710
## sukubunga.dl1 0.3473209 0.2248899 1.544 0.12432
## ln_kurs.dl1 0.0006877 0.0217630 0.032 0.97483
## inflasi.dl1 0.2468503 0.0749396 3.294 0.00120 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5299 on 173 degrees of freedom
## Multiple R-squared: 0.1033, Adjusted R-squared: 0.07736
## F-statistic: 3.985 on 5 and 173 DF, p-value: 0.001905
model_irf <- vec2var(jo_trace, r = r_rank)
# Error correction term (ECT) sepanjang waktu: jarak suku bunga dari keseimbangan
b <- vecm$beta[, 1]
ect <- ts(cbind(data_ts, 1) %*% b, start = c(thn, bln), frequency = 12)
plot(ect, ylab = "ECT (poin persen)", xlab = "",
main = "Error correction term: suku bunga aktual - keseimbangan")
abline(h = 0, lty = 2)
round(c(min = min(ect), max = max(ect), terakhir = tail(ect, 1)), 2)
## min max terakhir
## -4.92 4.23 -0.79
# 6a. Stabilitas: modulus akar matriks companion VECM.
# Harus ada (3 - r) akar = 1 (tren stokastik bersama); sisanya < 1.
akar_vecm <- function(v) {
K <- v$K; p <- v$p
A <- do.call(cbind, v$A)
comp <- rbind(A, cbind(diag(K * (p - 1)), matrix(0, K * (p - 1), K)))
sort(Mod(eigen(comp)$values), decreasing = TRUE)
}
round(akar_vecm(model_irf), 4)
## [1] 1.0000 1.0000 0.8385 0.4067 0.2546 0.2546
# 6b. Autokorelasi, normalitas, heteroskedastisitas residual
serial.test(model_irf, lags.pt = 12, type = "PT.asymptotic")
##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object model_irf
## Chi-squared = 110.73, df = 93, p-value = 0.1013
normality.test(model_irf)$jb.mul$JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object model_irf
## Chi-squared = 441.98, df = 6, p-value < 2.2e-16
arch.test(model_irf, lags.multi = 5)$arch.mul
##
## ARCH (multivariate)
##
## data: Residuals of VAR object model_irf
## Chi-squared = 256.23, df = 180, p-value = 0.0001632
v <- colnames(data_ts)
set.seed(123)
irf_all <- irf(model_irf, n.ahead = 24, ortho = TRUE, boot = TRUE, ci = 0.95, runs = 500)
plot(irf_all) # 3 grafik: respons terhadap shock sukubunga, ln_kurs, inflasi
# Angka IRF bulan ke-0, 1, 3, 6, 12, 24 beserta batas CI 95%
tabel_irf <- function(imp) {
h <- c(1, 2, 4, 7, 13, 25)
do.call(cbind, lapply(v, function(res) {
out <- cbind(irf_all$irf[[imp]][h, res], irf_all$Lower[[imp]][h, res], irf_all$Upper[[imp]][h, res])
colnames(out) <- paste0(res, c("", "_low", "_up")); round(out, 4)
})) -> m
rownames(m) <- paste0("h=", h - 1); m
}
lapply(setNames(v, v), tabel_irf)
## $sukubunga
## sukubunga sukubunga_low sukubunga_up ln_kurs ln_kurs_low ln_kurs_up
## h=0 0.1348 0.1150 0.1500 0.3042 0.0102 0.5493
## h=1 0.1901 0.1562 0.2169 0.4481 -0.0004 0.8108
## h=3 0.2263 0.1725 0.2706 0.5618 -0.0631 1.0413
## h=6 0.2412 0.1711 0.3013 0.5963 -0.0903 1.1780
## h=12 0.2515 0.1655 0.3335 0.6064 -0.0918 1.2162
## h=24 0.2562 0.1603 0.3445 0.6102 -0.0804 1.2442
## inflasi inflasi_low inflasi_up
## h=0 0.1114 0.0114 0.2331
## h=1 0.1830 0.0545 0.3385
## h=3 0.2202 0.0811 0.3816
## h=6 0.2044 0.0834 0.3335
## h=12 0.1769 0.0649 0.2884
## h=24 0.1635 0.0486 0.2721
##
## $ln_kurs
## sukubunga sukubunga_low sukubunga_up ln_kurs ln_kurs_low ln_kurs_up
## h=0 0.0000 0.0000 0.0000 1.7840 1.3787 2.1128
## h=1 0.0206 0.0020 0.0374 2.2209 1.6122 2.6534
## h=3 0.0460 0.0060 0.0821 2.3710 1.6497 2.8707
## h=6 0.0613 0.0079 0.1167 2.3968 1.6377 2.9769
## h=12 0.0734 0.0051 0.1460 2.4074 1.6370 3.0460
## h=24 0.0790 0.0044 0.1635 2.4119 1.6220 3.0849
## inflasi inflasi_low inflasi_up
## h=0 0.0463 -0.0043 0.1061
## h=1 0.0447 -0.0581 0.1460
## h=3 0.0274 -0.1146 0.1488
## h=6 -0.0022 -0.1454 0.1147
## h=12 -0.0356 -0.1737 0.0977
## h=24 -0.0515 -0.1906 0.0898
##
## $inflasi
## sukubunga sukubunga_low sukubunga_up ln_kurs ln_kurs_low ln_kurs_up
## h=0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## h=1 0.0134 -0.0042 0.0329 0.1957 -0.0480 0.4498
## h=3 0.0568 0.0105 0.1057 0.3522 -0.1322 0.8536
## h=6 0.1081 0.0232 0.1813 0.4060 -0.3640 1.2083
## h=12 0.1571 0.0290 0.2652 0.4451 -0.6110 1.6248
## h=24 0.1800 0.0292 0.3317 0.4633 -0.6911 1.8816
## inflasi inflasi_low inflasi_up
## h=0 0.5083 0.4081 0.5833
## h=1 0.5777 0.4337 0.6857
## h=3 0.4754 0.2807 0.5943
## h=6 0.3282 0.1139 0.5159
## h=12 0.1883 0.0299 0.4305
## h=24 0.1228 0.0009 0.3760
fevd_res <- fevd(model_irf, n.ahead = 24)
lapply(fevd_res, function(m) round(m[c(1, 6, 12, 24), ] * 100, 2))
## $sukubunga
## sukubunga ln_kurs inflasi
## [1,] 100.00 0.00 0.00
## [2,] 89.98 3.40 6.61
## [3,] 79.31 4.72 15.97
## [4,] 70.50 5.43 24.07
##
## $ln_kurs
## sukubunga ln_kurs inflasi
## [1,] 2.83 97.17 0.00
## [2,] 4.80 93.52 1.69
## [3,] 5.30 92.38 2.32
## [4,] 5.56 91.63 2.81
##
## $inflasi
## sukubunga ln_kurs inflasi
## [1,] 4.54 0.79 94.67
## [2,] 14.00 0.40 85.60
## [3,] 19.70 0.41 79.89
## [4,] 27.11 1.18 71.71
warna <- c("grey20", "grey55", "grey85")
par(mfrow = c(3, 1), mar = c(4, 4, 2.5, 9), xpd = TRUE)
for (nm in names(fevd_res)) {
barplot(t(fevd_res[[nm]] * 100), col = warna, border = NA, names.arg = 1:24,
ylab = "Persen (%)", xlab = "Horizon (bulan)",
main = paste("Dekomposisi varians", nm))
legend("right", inset = c(-0.2, 0), legend = colnames(fevd_res[[nm]]),
fill = warna, bty = "n", cex = 0.9)
}
par(mfrow = c(1, 1), xpd = FALSE)