## ============================================================
## STUDI KASUS 2 — VAR/VECM: Suku Bunga, Kurs, dan Inflasi
## UAS Analisis Runtun Waktu Terapan
## ============================================================
## Spesifikasi akhir: VECM, lag K = 2 (VAR level), r = 1, konstanta
## terestriksi dalam kointegrasi, impulse dummy BI-7DRR Agustus 2016.
# install.packages(c("urca","vars","tseries"))
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
# Drift kurs: apakah rata-rata perubahan bulanan ln kurs signifikan?
t.test(diff(ln_kurs))   # jika signifikan -> ada tren naik (depresiasi rata-rata)
## 
##  One Sample t-test
## 
## data:  diff(ln_kurs)
## t = 2.8753, df = 178, p-value = 0.00453
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.1246155 0.6699564
## sample estimates:
## mean of x 
## 0.3972859
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 K = 2, 3 dan ecdet = "const"/"none": rank (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 <- expand.grid(K = 2:3, ecdet = c("const", "none"), 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
## Chi-squared2 2  none       1       1         1        0.0650 0.0002    0
## Chi-squared3 3  none       1       0         1        0.2149 0.0188    0
# Kriteria pilih: (1) rank trace & eigen konsisten, (2) Portmanteau_p > 0,05,
# (3) lag sesuai kriteria informasi. Hasil: K = 2 & ecdet = "const"
# (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"

# Uji formal: konstanta cukup di dalam hubungan kointegrasi (ecdet = "const")
# atau perlu drift/tren linear di level (ecdet = "none")?
# H0: tidak ada tren linear di level -> p > 0,05 berarti "const" tepat.
jo_lt <- ca.jo(data_ts, type = "trace", ecdet = "const", K = K_pilih,
               spec = "transitory", dumvar = dum_impulse)
lttest(jo_lt, r = 1)
## LR-test for no linear trend
## 
## H0: H*2(r<=1)
## H1: H2(r<=1)
## 
## Test statistic is distributed as chi-square
## with 2 degress of freedom
##         test statistic p-value
## LR test           5.26    0.07
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
# Kecepatan penyesuaian & half-life per variabel:
# koreksi efektif = alpha_i * beta_i (harus negatif agar kembali ke keseimbangan)
alpha <- coef(vecm$rlm)["ect1", ]
beta  <- vecm$beta[1:3, 1]
koreksi <- alpha * beta
data.frame(alpha = round(alpha, 4), beta = round(beta, 4),
           koreksi = round(koreksi, 4),
           half_life_bulan = round(ifelse(koreksi < 0, log(0.5) / log(1 + koreksi), NA), 1))
##               alpha    beta koreksi half_life_bulan
## sukubunga.d -0.0182  1.0000 -0.0182            37.7
## ln_kurs.d   -0.0457 -0.0631  0.0029              NA
## inflasi.d    0.0817 -1.3508 -0.1103             5.9
model_irf <- vec2var(jo_trace, r = r_rank)
# 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
# Variabel I(1): VAR level dengan lag K + 1, uji Wald hanya pada K lag pertama.
var_ty <- VAR(data_ts, p = K_pilih + 1, type = "const", exogen = dum_step)
granger_ty <- function(sebab, akibat) {
  m <- var_ty$varresult[[akibat]]
  idx <- paste0(sebab, ".l", 1:K_pilih)
  b <- coef(m)[idx]; V <- vcov(m)[idx, idx]
  W <- as.numeric(t(b) %*% solve(V) %*% b)
  data.frame(H0 = paste(sebab, "tidak Granger-cause", akibat),
             Wald = round(W, 3), df = K_pilih,
             p_value = round(pchisq(W, K_pilih, lower.tail = FALSE), 4))
}
v <- colnames(data_ts)
tabel_granger <- do.call(rbind, lapply(v, function(s)
  do.call(rbind, lapply(setdiff(v, s), function(a) granger_ty(s, a)))))
tabel_granger
##                                      H0  Wald df p_value
## 1 sukubunga tidak Granger-cause ln_kurs 2.602  2  0.2723
## 2 sukubunga tidak Granger-cause inflasi 1.492  2  0.4742
## 3 ln_kurs tidak Granger-cause sukubunga 7.174  2  0.0277
## 4   ln_kurs tidak Granger-cause inflasi 0.327  2  0.8492
## 5 inflasi tidak Granger-cause sukubunga 1.892  2  0.3884
## 6   inflasi tidak Granger-cause ln_kurs 2.940  2  0.2299
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
# Urutan inflasi -> kurs -> suku bunga (BI bereaksi terhadap inflasi & kurs
# dalam bulan yang sama). Jika arah respons tetap sama, kesimpulan IRF robust.
data_alt <- data_ts[, c("inflasi", "ln_kurs", "sukubunga")]
jo_alt   <- ca.jo(data_alt, type = "trace", ecdet = ecdet_pilih, K = K_pilih,
                  spec = "transitory", dumvar = dum_impulse)
irf_alt  <- irf(vec2var(jo_alt, r = r_rank), impulse = "sukubunga", n.ahead = 24,
                ortho = TRUE, boot = TRUE, ci = 0.95, runs = 500, seed = 123)
round(cbind(inflasi = irf_alt$irf$sukubunga[c(1, 2, 4, 7, 13, 25), "inflasi"],
            low     = irf_alt$Lower$sukubunga[c(1, 2, 4, 7, 13, 25), "inflasi"],
            up      = irf_alt$Upper$sukubunga[c(1, 2, 4, 7, 13, 25), "inflasi"],
            ln_kurs = irf_alt$irf$sukubunga[c(1, 2, 4, 7, 13, 25), "ln_kurs"]), 4)
##      inflasi     low     up ln_kurs
## [1,]  0.0000  0.0000 0.0000  0.0000
## [2,]  0.0559 -0.0006 0.1206  0.0286
## [3,]  0.1148  0.0131 0.2206  0.0829
## [4,]  0.1334  0.0288 0.2273  0.1015
## [5,]  0.1397  0.0375 0.2245  0.1018
## [6,]  0.1423  0.0426 0.2205  0.1011
plot(irf_alt)

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)