Proyek Ekonometrika - VECM

A. Library

library(readxl)
library(lubridate)
library(dplyr)
library(zoo)
library(psych)
library(patchwork)
library(ggplot2)
library(urca)
library(vars)
library(Metrics)
library(tseries)
library(sandwich)
library(lmtest)
library(car)

B. Data

2.1. Import Data

data <- read_xlsx("D:/UNY/MySta/SEM 5/Ekonometrika/Proyek/Dataset_Ekonom.xlsx")
head(data)

## # A tibble: 6 × 6
##   Date          Inflasi BI_rate  Kurs Money_supply  IHSG
##   <chr>           <dbl>   <dbl> <dbl>        <dbl> <dbl>
## 1 Januari 2013     4.57    5.75   NA      3268789. 4454.
## 2 Februari 2013    5.31    5.75   NA      3280420. 4796.
## 3 Maret 2013       5.9     5.75   NA      3322529. 4941.
## 4 April 2013       5.57    5.75   NA      3360928. 5034.
## 5 Mei 2013         5.47    5.75 9786.     3426305. 5069.
## 6 Juni 2013        5.9     6    9882.     3413379. 4819.

2.2. Data Structure, Character Variables, and Size

# Structure
str(data)

## tibble [144 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Date        : chr [1:144] "Januari 2013" "Februari 2013" "Maret 2013" "April 2013" ...
##  $ Inflasi     : num [1:144] 4.57 5.31 5.9 5.57 5.47 5.9 8.61 8.79 8.4 8.32 ...
##  $ BI_rate     : num [1:144] 5.75 5.75 5.75 5.75 5.75 6 6.5 7 7.25 7.25 ...
##  $ Kurs        : num [1:144] NA NA NA NA 9786 ...
##  $ Money_supply: num [1:144] 3268789 3280420 3322529 3360928 3426305 ...
##  $ IHSG        : num [1:144] 4454 4796 4941 5034 5069 ...

# Character variable
sapply(data, class)

##         Date      Inflasi      BI_rate         Kurs Money_supply         IHSG 
##  "character"    "numeric"    "numeric"    "numeric"    "numeric"    "numeric"

# Size
cat("Number of rows:", nrow(data), "\n")

## Number of rows: 144

cat("Number of columns:", ncol(data), "\n")

## Number of columns: 6

2.3. Data Summary

summary(data)

##      Date              Inflasi         BI_rate           Kurs      
##  Length:144         Min.   :1.320   Min.   :3.500   Min.   : 9786  
##  Class :character   1st Qu.:2.663   1st Qu.:4.438   1st Qu.:13313  
##  Mode  :character   Median :3.340   Median :5.750   Median :14114  
##                     Mean   :3.927   Mean   :5.543   Mean   :13920  
##                     3rd Qu.:4.947   3rd Qu.:6.500   3rd Qu.:14831  
##                     Max.   :8.790   Max.   :7.750   Max.   :16337  
##                                                     NA's   :4      
##   Money_supply          IHSG     
##  Min.   :3268789   Min.   :4195  
##  1st Qu.:4518614   1st Qu.:5082  
##  Median :5670877   Median :5944  
##  Mean   :6002823   Mean   :5853  
##  3rd Qu.:7657626   3rd Qu.:6612  
##  Max.   :9210816   Max.   :7671  
##

C. Data Preprocessing

3.1. Konversi kolom Date

bulan_id <- c("Januari"="01","Februari"="02","Maret"="03","April"="04",
              "Mei"="05","Juni"="06","Juli"="07","Agustus"="08",
              "September"="09","Oktober"="10","November"="11","Desember"="12")

data <- data %>%
  mutate(
    Month = sapply(strsplit(Date, " "), `[`, 1),
    Year  = sapply(strsplit(Date, " "), `[`, 2),
    Month_num = bulan_id[Month],
    Date2 = ymd(paste(Year, Month_num, "01", sep="-"))
  ) %>%
  dplyr::select(Date2, Inflasi, BI_rate, Kurs, Money_supply, IHSG)

# short by time
data <- data %>% arrange(Date2)

3.2. Missing and Duplicate Values

# Check missing values
colSums(is.na(data))

##        Date2      Inflasi      BI_rate         Kurs Money_supply         IHSG 
##            0            0            0            4            0            0

# Check duplicate values
sum(duplicated(data))

## [1] 0

Handling Missing Value: dengan na.fill(“extend”) karena metode ini stabil untuk time-series ekonomi, dengan mengisi NA awal menggunakan nilai valid terdekat tanpa menciptakan pola data yang tidak nyata.

# Handling missing value
data$Kurs <- na.fill(data$Kurs, "extend")
# Cek NA kembali
colSums(is.na(data))

##        Date2      Inflasi      BI_rate         Kurs Money_supply         IHSG 
##            0            0            0            0            0            0

3.3. Deteksi Outlier

# Outlier tiap variabel (IQR method)
vars <- c("Inflasi", "BI_rate", "Kurs", "Money_supply", "IHSG")

outlier_count <- sapply(vars, function(v){
  x <- data[[v]]
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR_val <- Q3 - Q1
  sum(x < (Q1 - 1.5*IQR_val) | x > (Q3 + 1.5*IQR_val), na.rm = TRUE)
})

# Tabel distribusi outlier
outlier_table <- data.frame(Variabel = vars, Jumlah_Outlier = outlier_count)
print(outlier_table)

##                  Variabel Jumlah_Outlier
## Inflasi           Inflasi              4
## BI_rate           BI_rate              0
## Kurs                 Kurs              8
## Money_supply Money_supply              0
## IHSG                 IHSG              0

# Boxplot
par(mfrow=c(2,3))
boxplot(data$Inflasi, main="Inflasi")
boxplot(data$BI_rate, main="BI Rate")
boxplot(data$Kurs, main="Kurs")
boxplot(data$Money_supply, main="Money Supply")
boxplot(data$IHSG, main="IHSG")

Interpretasi distribusi outlier: Outlier tidak perlu di-handling karena seluruh nilai ekstrem tersebut merupakan shock ekonomi yang nyata bukan kesalahan data, sehingga dipertahankan agar dinamika hubungan antarvariabel dalam model VECM tetap mencerminkan kondisi ekonomi sebenarnya.

3.4. Transformasi log

Pada time-series ekonomi, variabel Kurs, Money Supply, dan IHSG dilog karena ketiganya bersifat trending dan berskala besar, sehingga log transform membuat varians lebih stabil, hubungan lebih linear, dan hasil VECM lebih valid serta ekonomis untuk diinterpretasikan

data <- data %>%
  mutate(
    ln_Kurs  = log(Kurs),
    ln_M2    = log(Money_supply),
    ln_IHSG  = log(IHSG)
  )

# Dataset Final setelah Preprocessing
final_data <- data %>%
  dplyr::select(Date2, Inflasi, BI_rate, ln_Kurs, ln_M2, ln_IHSG)

D. Exploratory Data Analysis

4.1. Deskripsi Statistik

describe(final_data[, -1])

##         vars   n  mean   sd median trimmed  mad   min   max range  skew
## Inflasi    1 144  3.93 1.94   3.34    3.74 1.62  1.32  8.79  7.47  0.83
## BI_rate    2 144  5.54 1.34   5.75    5.55 1.48  3.50  7.75  4.25  0.03
## ln_Kurs    3 144  9.53 0.11   9.55    9.54 0.08  9.19  9.70  0.51 -1.23
## ln_M2      4 144 15.56 0.30  15.55   15.57 0.37 15.00 16.04  1.04 -0.10
## ln_IHSG    5 144  8.66 0.16   8.69    8.67 0.21  8.34  8.95  0.60 -0.18
##         kurtosis   se
## Inflasi    -0.20 0.16
## BI_rate    -1.15 0.11
## ln_Kurs     1.46 0.01
## ln_M2      -1.15 0.03
## ln_IHSG    -1.14 0.01

4.2. Plot Time Series untuk Setiap Variabel

# line plot time-series
plot_ts <- function(df, var, label){
ggplot(df, aes(x = Date2, y = .data[[var]])) +
geom_line(color = "steelblue", linewidth = 1) +
labs(title = paste("", label),
x = "Tahun",
y = label) +
theme_minimal()
}

p1 <- plot_ts(final_data, "Inflasi", "Inflasi (%)")
p2 <- plot_ts(final_data, "BI_rate", "BI Rate (%)")
p3 <- plot_ts(final_data, "ln_Kurs", "Kurs (Log)")
p4 <- plot_ts(final_data, "ln_M2", "Money Supply (Log)")
p5 <- plot_ts(final_data, "ln_IHSG", "IHSG (Log)")

(p1 + p2 + p3) /
(p4 + p5 + plot_spacer())

Plot di atas digunakan untuk menentukan deterministic terms yakni none, drift, atau tren. yang nantinya digunakan pada uji ADF. Dari plot di atas diperoleh term ADF drif untuk variabel Inflasi,ln_Kurs, dan BI_rate, sedangkan trend untuk variabel ln_M2, dan ln_IHSG

4.3. Plot IHSG per Tahun

final_data$Year  <- format(final_data$Date2, "%Y")
final_data$Bulan <- format(final_data$Date2, "%b")
final_data$Bulan <- factor(
  final_data$Bulan,
  levels = c("Jan","Feb","Mar","Apr","May","Jun",
             "Jul","Aug","Sep","Oct","Nov","Dec")
)

ggplot(final_data, aes(x = Bulan, y = ln_IHSG, color = Year, group = Year)) +
  geom_line(linewidth = 1) +
  labs(
    title = "Pergerakan IHSG Berdasarkan Tahun",
    x = "Bulan",
    y = "ln(IHSG)",
    color = "Tahun"
  ) +
  theme_minimal()

E. Hasil Analisis

1. Pengujian Stasioneritas

uji stasioneritas data dalam penelitian ini adalah uji Augmented Dickey Fuller dengan menggunakan taraf signifikan 5%. Pengujian hasil stasioneritas pada tingkat level dan first difference.

df <- final_data %>% 
  dplyr::select(Inflasi, BI_rate, ln_Kurs, ln_M2, ln_IHSG)

adf_test <- function(x, det_type){
  
  # LEVEL
  test_level <- ur.df(x, type = det_type, selectlags = "AIC")
  col_level  <- colnames(test_level@teststat)[1]
  adf_level  <- round(test_level@teststat["statistic", col_level], 4)
  crit_level <- test_level@cval[col_level, "5pct"]
  stat_level <- ifelse(adf_level < crit_level, "Stasioner", "Tidak Stasioner")
  
  # FIRST DIFFERENCE
  test_diff <- ur.df(diff(x), type = "drift", selectlags = "AIC")
  col_diff  <- colnames(test_diff@teststat)[1]
  adf_diff  <- round(test_diff@teststat["statistic", col_diff], 4)
  crit_diff <- test_diff@cval[col_diff, "5pct"]
  stat_diff <- ifelse(adf_diff < crit_diff, "Stasioner", "Tidak Stasioner")
  
  return(data.frame(
    ADF_Level      = adf_level,
    Crit_Level     = crit_level,
    Status_Level   = stat_level,
    ADF_FirstDiff  = adf_diff,
    Crit_FirstDiff = crit_diff,
    Status_Diff    = stat_diff
  ))
}


#  deterministic term 
det_terms <- c(
  Inflasi = "drift",
  BI_rate = "drift",
  ln_Kurs = "drift",
  ln_M2   = "trend",
  ln_IHSG = "trend"
)

adf_results <- do.call(rbind, lapply(names(df), function(v){
  adf_test(df[[v]], det_terms[[v]])
}))
rownames(adf_results) <- names(df)
adf_results

##         ADF_Level Crit_Level    Status_Level ADF_FirstDiff Crit_FirstDiff
## Inflasi   -2.0059      -2.88 Tidak Stasioner       -8.6798          -2.88
## BI_rate   -1.4064      -2.88 Tidak Stasioner       -5.0179          -2.88
## ln_Kurs   -2.8051      -2.88 Tidak Stasioner      -10.0603          -2.88
## ln_M2     -2.8359      -3.43 Tidak Stasioner      -10.7190          -2.88
## ln_IHSG   -3.1812      -3.43 Tidak Stasioner       -7.5032          -2.88
##         Status_Diff
## Inflasi   Stasioner
## BI_rate   Stasioner
## ln_Kurs   Stasioner
## ln_M2     Stasioner
## ln_IHSG   Stasioner

Interpretasi: variabel Y, X1, X2, X3, X4 merupakan data yang tidak stasioner pada tingkat level (data asli) karena nilai ADF lebih besar dibandingkan dengan tabel MacKinnon dengan tingkat signifikansi 5%. Sementara itu, dari hasil diferensi pertama terlihat bahwa nilai ADF lebih kecil dari tingkat signifikansi artinya semua variabel telah stasioner.

2. Pengujian Lag Optimal

Tahap ini menentukan panjang lag optimal pada penelitian ini didasarkan pada Akaike Information Criterion (AIC)

lag_sel <- VARselect(df, lag.max = 10, type = "const")

# Nilai AIC dari semua lag
AIC_values <- lag_sel$criteria["AIC(n)", ]

# Tabel
Tabel_Lag <- data.frame(
  Lag = 1:length(AIC_values),
  AIC = round(AIC_values, 5)
)
print(Tabel_Lag)

##    Lag       AIC
## 1    1 -28.43777
## 2    2 -28.79574
## 3    3 -28.75053
## 4    4 -28.51073
## 5    5 -28.42499
## 6    6 -28.33043
## 7    7 -28.36593
## 8    8 -28.28466
## 9    9 -28.17697
## 10  10 -28.16016

# Lag optimal
lag_optimal <- lag_sel$selection["AIC(n)"]
AIC_min <- AIC_values[lag_optimal]
cat("Lag optimal =", lag_optimal, "\n")

## Lag optimal = 2

cat("Nilai AIC terkecil =", round(AIC_min, 5), "\n")

## Nilai AIC terkecil = -28.79574

Interpretasi: panjang lag optimal terletak pada lag 2 dikarenakan memiliki nilai AIC terkecil yaitu sebesar -28.79574. Penggunaan lag 2 sebagai lag optimal menyatakan bahwa semua variabel saling mempengaruhi satu sama lain tidak hanya pada periode yang sama, namun variabel-variabel tersebut saling berkaitan hingga pada dua periode sebelumnya.

3. Pengujian Stabilitas Model VAR

# lag optimal = 2
lag_optimal <- which.min(VARselect(df, lag.max=10, type="const")$criteria["AIC(n)", ])

# Estimasi VAR sesuai lag optimal
var_model <- VAR(df, p = lag_optimal, type = "const")

# Uji stabilitas (roots)
VAR_roots <- roots(var_model)

# Tabel modulus & roots
Stabilitas_VAR <- data.frame(
  Root = VAR_roots,
  Modulus = Mod(VAR_roots)
)

print(Stabilitas_VAR)

##          Root    Modulus
## 1  0.99661107 0.99661107
## 2  0.93111825 0.93111825
## 3  0.91246959 0.91246959
## 4  0.91246959 0.91246959
## 5  0.78482872 0.78482872
## 6  0.40582520 0.40582520
## 7  0.29049035 0.29049035
## 8  0.29049035 0.29049035
## 9  0.25863738 0.25863738
## 10 0.03495285 0.03495285

Interpretasi: Berdasarkan hasil uji stabilitas model VAR, seluruh nilai roots memiliki modulus kurang dari satu dan seluruh titik roots berada di dalam lingkaran satuan.Hal ini menunjukkan bahwa model VAR yang digunakan adalah stabil (stationary) sehingga analisis lebih lanjut seperti pengujian kointegrasi dan estimasi VECM dapat dilakukan.

4. Pengujian Kointegrasi

# 4. Pengujian Kointegrasi Johansen

# Johansen Trace Test
joh_trace <- ca.jo(
  df,
  type = "trace",
  K = lag_optimal, 
  ecdet = "const",
  spec = "transitory"
)

summary(joh_trace)

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 3.986501e-01 1.976767e-01 1.336932e-01 4.147073e-02 3.117743e-02
## [6] 2.103063e-16
## 
## Values of teststatistic and critical values of test:
## 
##            test 10pct  5pct  1pct
## r <= 4 |   4.50  7.52  9.24 12.97
## r <= 3 |  10.51 17.85 19.96 24.60
## r <= 2 |  30.89 32.00 34.91 41.07
## r <= 1 |  62.17 49.65 53.12 60.16
## r = 0  | 134.38 71.86 76.07 84.45
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##              Inflasi.l1   BI_rate.l1  ln_Kurs.l1     ln_M2.l1 ln_IHSG.l1
## Inflasi.l1    1.0000000    1.0000000   1.0000000     1.000000   1.000000
## BI_rate.l1   -0.3737654    0.4431524  -0.5022928    -8.772424  -5.384737
## ln_Kurs.l1   12.8706749  133.4946434   0.6286595   298.238227 -12.557590
## ln_M2.l1     -1.9984540   -8.5746679   2.4293010  -184.285081  17.261337
## ln_IHSG.l1    7.8652405  -56.8708544   0.9264255   205.483827 -33.842287
## constant   -165.4213080 -651.5049093 -52.5830612 -1704.777980 169.634669
##             constant
## Inflasi.l1   1.00000
## BI_rate.l1  -3.24294
## ln_Kurs.l1  56.54306
## ln_M2.l1   -46.10507
## ln_IHSG.l1  11.11568
## constant    95.77644
## 
## Weights W:
## (This is the loading matrix)
## 
##              Inflasi.l1    BI_rate.l1   ln_Kurs.l1      ln_M2.l1    ln_IHSG.l1
## Inflasi.d -0.0101296331 -1.712241e-02 -0.104446942  1.493577e-03  3.290964e-03
## BI_rate.d  0.0056740797 -3.296148e-03  0.030836741  7.032698e-04  2.301551e-03
## ln_Kurs.d -0.0008393762 -6.613512e-04  0.002070204 -7.021812e-05 -1.756638e-04
## ln_M2.d   -0.0025118928 -5.884919e-05  0.001026363  2.236242e-05 -4.067006e-05
## ln_IHSG.d -0.0007172964  8.293432e-04 -0.001974261 -1.698191e-04  7.026912e-04
##                constant
## Inflasi.d -7.752730e-15
## BI_rate.d  1.068305e-15
## ln_Kurs.d  2.420425e-15
## ln_M2.d    9.255746e-16
## ln_IHSG.d  1.589752e-15

# Johansen Max-Eigen Test
joh_eigen <- ca.jo(
  df,
  type = "eigen",
  K = lag_optimal,      
  ecdet = "const",
  spec = "transitory"
)

summary(joh_eigen)

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 3.986501e-01 1.976767e-01 1.336932e-01 4.147073e-02 3.117743e-02
## [6] 2.103063e-16
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 4 |  4.50  7.52  9.24 12.97
## r <= 3 |  6.01 13.75 15.67 20.20
## r <= 2 | 20.38 19.77 22.00 26.81
## r <= 1 | 31.27 25.56 28.14 33.24
## r = 0  | 72.22 31.66 34.40 39.79
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##              Inflasi.l1   BI_rate.l1  ln_Kurs.l1     ln_M2.l1 ln_IHSG.l1
## Inflasi.l1    1.0000000    1.0000000   1.0000000     1.000000   1.000000
## BI_rate.l1   -0.3737654    0.4431524  -0.5022928    -8.772424  -5.384737
## ln_Kurs.l1   12.8706749  133.4946434   0.6286595   298.238227 -12.557590
## ln_M2.l1     -1.9984540   -8.5746679   2.4293010  -184.285081  17.261337
## ln_IHSG.l1    7.8652405  -56.8708544   0.9264255   205.483827 -33.842287
## constant   -165.4213080 -651.5049093 -52.5830612 -1704.777980 169.634669
##             constant
## Inflasi.l1   1.00000
## BI_rate.l1  -3.24294
## ln_Kurs.l1  56.54306
## ln_M2.l1   -46.10507
## ln_IHSG.l1  11.11568
## constant    95.77644
## 
## Weights W:
## (This is the loading matrix)
## 
##              Inflasi.l1    BI_rate.l1   ln_Kurs.l1      ln_M2.l1    ln_IHSG.l1
## Inflasi.d -0.0101296331 -1.712241e-02 -0.104446942  1.493577e-03  3.290964e-03
## BI_rate.d  0.0056740797 -3.296148e-03  0.030836741  7.032698e-04  2.301551e-03
## ln_Kurs.d -0.0008393762 -6.613512e-04  0.002070204 -7.021812e-05 -1.756638e-04
## ln_M2.d   -0.0025118928 -5.884919e-05  0.001026363  2.236242e-05 -4.067006e-05
## ln_IHSG.d -0.0007172964  8.293432e-04 -0.001974261 -1.698191e-04  7.026912e-04
##                constant
## Inflasi.d -7.752730e-15
## BI_rate.d  1.068305e-15
## ln_Kurs.d  2.420425e-15
## ln_M2.d    9.255746e-16
## ln_IHSG.d  1.589752e-15

Berdasarkan hasil uji kointegrasi Johansen menggunakan trace statistic pada tingkat signifikansi 5%, diperoleh bahwa nilai trace statistic pada hipotesis r = 0 dan r ≤ 1 lebih besar daripada nilai critical value 5%, sehingga H0 ditolak. Namun, pada hipotesis r ≤ 2, nilai trace statistic lebih kecil daripada nilai critical value 5%, sehingga H0 gagal ditolak.

baik Trace Test maupun Max-Eigen Test, keduanya secara konsisten menunjukkan jumlah hubungan kointegrasi = 2 (r = 2)

Dengan demikian, model yang tepat digunakan adalah VECM (Vector Error Correction Model)

5. Estimasi Vector Error Correction Model

vecm_model <- cajorls(joh_trace, r = 2)
summary(vecm_model)

##      Length Class  Mode   
## rlm  12     mlm    list   
## beta 12     -none- numeric

a. Estimasi Jangka Panjang (Long-run VECM)

# VECM Jangka Panjang

beta_raw <- vecm_model$beta   

# Normalisasi ke IHSG untuk kedua vektor
beta_norm1 <- beta_raw[, 1] / beta_raw["ln_IHSG.l1", 1]
beta_norm2 <- beta_raw[, 2] / beta_raw["ln_IHSG.l1", 2]
beta_norm <- cbind(beta_norm1, beta_norm2)

# Tabel
longrun_table <- data.frame(
  Variabel = rownames(beta_norm),
  Beta_VECM_1 = round(beta_norm[, 1], 6),
  Beta_VECM_2 = round(beta_norm[, 2], 6),
  row.names = NULL
)
longrun_table

##     Variabel Beta_VECM_1 Beta_VECM_2
## 1 Inflasi.l1   -0.045970    0.000000
## 2 BI_rate.l1    0.000000   -0.012619
## 3 ln_Kurs.l1   -3.128681   -1.863319
## 4   ln_M2.l1    0.230182    0.101585
## 5 ln_IHSG.l1    1.000000    1.000000
## 6   constant   17.827886    7.508695

b. Estimasi Jangka Pendek (Short-run VECM)

# VECM JANGKA PENDEK 

vecm_eq <- summary(vecm_model$rlm)

# Tabel
extract_shortrun <- function(model_summary, varname) {
  coef_table <- model_summary$coefficients
  data.frame(
    Variabel   = rownames(coef_table),
    Estimate   = round(coef_table[, 1], 6),
    t_stat     = round(coef_table[, 3], 6),
    Signifikan = ifelse(abs(coef_table[,3]) > 1.96, "Signifikan (5%)", "Tidak Signifikan"),
    Persamaan  = varname,
    row.names  = NULL
  )
}

# Persamaan
short_tables <- rbind(
  extract_shortrun(vecm_eq[[1]], "Inflasi.d"),
  extract_shortrun(vecm_eq[[2]], "BI_rate.d"),
  extract_shortrun(vecm_eq[[3]], "ln_Kurs.d"),
  extract_shortrun(vecm_eq[[4]], "ln_M2.d"),
  extract_shortrun(vecm_eq[[5]], "ln_IHSG.d")
)

short_tables

##       Variabel  Estimate    t_stat       Signifikan Persamaan
## 1         ect1 -0.027252 -2.124444  Signifikan (5%) Inflasi.d
## 2         ect2 -0.003802 -0.772040 Tidak Signifikan Inflasi.d
## 3  Inflasi.dl1  0.147692  1.750209 Tidak Signifikan Inflasi.d
## 4  BI_rate.dl1  0.095043  0.411366 Tidak Signifikan Inflasi.d
## 5  ln_Kurs.dl1  0.545120  0.220536 Tidak Signifikan Inflasi.d
## 6    ln_M2.dl1 -2.973735 -0.876242 Tidak Signifikan Inflasi.d
## 7  ln_IHSG.dl1 -1.029607 -0.776155 Tidak Signifikan Inflasi.d
## 8         ect1  0.002378  0.570381 Tidak Signifikan BI_rate.d
## 9         ect2 -0.003581 -2.237892  Signifikan (5%) BI_rate.d
## 10 Inflasi.dl1  0.024776  0.903417 Tidak Signifikan BI_rate.d
## 11 BI_rate.dl1  0.392198  5.223134  Signifikan (5%) BI_rate.d
## 12 ln_Kurs.dl1  1.904001  2.370135  Signifikan (5%) BI_rate.d
## 13   ln_M2.dl1  1.107189  1.003836 Tidak Signifikan BI_rate.d
## 14 ln_IHSG.dl1  0.074194  0.172093 Tidak Signifikan BI_rate.d
## 15        ect1 -0.001501 -3.379210  Signifikan (5%) ln_Kurs.d
## 16        ect2  0.000021  0.121131 Tidak Signifikan ln_Kurs.d
## 17 Inflasi.dl1  0.000647  0.221535 Tidak Signifikan ln_Kurs.d
## 18 BI_rate.dl1  0.003212  0.401546 Tidak Signifikan ln_Kurs.d
## 19 ln_Kurs.dl1  0.062203  0.726888 Tidak Signifikan ln_Kurs.d
## 20   ln_M2.dl1  0.109983  0.936086 Tidak Signifikan ln_Kurs.d
## 21 ln_IHSG.dl1 -0.227865 -4.961588  Signifikan (5%) ln_Kurs.d
## 22        ect1 -0.002571 -8.371823  Signifikan (5%)   ln_M2.d
## 23        ect2  0.000913  7.743540  Signifikan (5%)   ln_M2.d
## 24 Inflasi.dl1  0.002996  1.483341 Tidak Signifikan   ln_M2.d
## 25 BI_rate.dl1  0.012203  2.206403  Signifikan (5%)   ln_M2.d
## 26 ln_Kurs.dl1  0.033109  0.559557 Tidak Signifikan   ln_M2.d
## 27   ln_M2.dl1 -0.399405 -4.916436  Signifikan (5%)   ln_M2.d
## 28 ln_IHSG.dl1 -0.016248 -0.511660 Tidak Signifikan   ln_M2.d
## 29        ect1  0.000112  0.122962 Tidak Signifikan ln_IHSG.d
## 30        ect2  0.000636  1.817127 Tidak Signifikan ln_IHSG.d
## 31 Inflasi.dl1 -0.002221 -0.370509 Tidak Signifikan ln_IHSG.d
## 32 BI_rate.dl1  0.003995  0.243388 Tidak Signifikan ln_IHSG.d
## 33 ln_Kurs.dl1  0.015442  0.087948 Tidak Signifikan ln_IHSG.d
## 34   ln_M2.dl1  0.026437  0.109662 Tidak Signifikan ln_IHSG.d
## 35 ln_IHSG.dl1  0.177753  1.886330 Tidak Signifikan ln_IHSG.d

6. Impulse Response Function (IRF)

# Konversi Johansen menjadi VAR(p)
variables <- colnames(df)
var_model <- vec2var(joh_trace, r = 2)

irf <- irf(
  var_model,
  impulse = variables,
  response = variables,
  n.ahead = 10,
  boot = TRUE,
  ci = 0.95
)

plot(irf)

## 7. Forcast Eror Variance Decomposition (FEVD)

# FEVD
fevd <- fevd(var_model, n.ahead = 10)
plot(fevd)

fevd_IHSG <- fevd$ln_IHSG

fevd_table <- data.frame(
  Period = 1:10,
  Inflasi = fevd_IHSG[,"Inflasi"] * 100,
  BI_rate = fevd_IHSG[,"BI_rate"] * 100,
  ln_Kurs = fevd_IHSG[,"ln_Kurs"] * 100,
  ln_M2 = fevd_IHSG[,"ln_M2"] * 100,
  ln_IHSG = fevd_IHSG[,"ln_IHSG"] * 100
)

knitr::kable(fevd_table, caption="FEVD IHSG (%)")

FEVD IHSG (%)
Period	Inflasi	BI_rate	ln_Kurs	ln_M2	ln_IHSG
1	1.444728	1.2221457	16.794190	0.2725501	80.26639
2	1.862657	0.9878017	14.864125	0.2456270	82.03979
3	2.063668	0.8261104	13.173696	0.2283131	83.70821
4	2.184536	0.7089725	11.759897	0.2162220	85.13037
5	2.282448	0.6222441	10.600034	0.2063713	86.28890
6	2.381186	0.5574765	9.656398	0.1977790	87.20716
7	2.491702	0.5093647	8.890842	0.1899386	87.91815
8	2.619690	0.4744219	8.270380	0.1826362	88.45287
9	2.768580	0.4502681	7.768094	0.1757620	88.83730
10	2.940671	0.4352231	7.362589	0.1692648	89.09225

Interpretasi: IHSG pada periode awal dipengaruhi oleh shock internal sebesar 80 persen. Dalam jangka panjang kontribusi shock internal IHSG meningkat menjadi 89 persen. Hal ini menunjukkan bahwa dinamika IHSG lebih banyak dijelaskan oleh pergerakan IHSG itu sendiri dibandingkan variabel makro lainnya. Dengan kata lain, variabel makro seperti inflasi, BI rate, nilai tukar, dan jumlah uang beredar memiliki pengaruh yang relatif kecil dalam menjelaskan variasi IHSG, terutama dalam horizon jangka panjang.

8. Pengujian Asumsi Klasik

res <- vecm_model$rlm$residuals
colnames(res) <- c("Inflasi", "BI_rate", "ln_Kurs", "ln_M2", "ln_IHSG")

a. Multikolinearitas

# Model OLS untuk VIF
model <- lm(ln_IHSG ~ Inflasi + BI_rate + ln_Kurs + ln_M2, data = df)
vif(model)

##  Inflasi  BI_rate  ln_Kurs    ln_M2 
## 2.552916 2.131576 5.491630 6.212092

semua VIF < 10 maka tidak ada multikolinearitas.

b. Autokorelasi

# Autokorelasi untuk ln_IHSG
dwtest(res[, "ln_IHSG"] ~ 1)

## 
##  Durbin-Watson test
## 
## data:  res[, "ln_IHSG"] ~ 1
## DW = 2.069, p-value = 0.6605
## alternative hypothesis: true autocorrelation is greater than 0

DW ≈ 2 mengindikasikan tidak ada autokorelasi.

c. Homoskedastisitas

# Residual ln_IHSG
y_res <- res[, "ln_IHSG"]

white_test <- bptest(y_res ~ 1, ~ y_res + I(y_res^2))
white_test

## 
##  studentized Breusch-Pagan test
## 
## data:  y_res ~ 1
## BP = 142, df = 2, p-value < 2.2e-16

chi_crit <- qchisq(0.95, df=2)
chi_crit

## [1] 5.991465

P-value << 0.05 → terdapat heteroskedastisitas pada residual ln_IHSG.

# robust t dan p-value
robust_coeftest <- function(y_res) {
  model <- lm(y_res ~ 1)
  robust_vcov <- vcovHC(model, type = "HC1")
  coeftest(model, vcov = robust_vcov)
}

# Terapkan ke residual
res_robust <- lapply(as.data.frame(res), robust_coeftest)
res_robust

## $Inflasi
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.027489   0.041768 -0.6581   0.5115
## 
## 
## $BI_rate
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.011304   0.013562  0.8335    0.406
## 
## 
## $ln_Kurs
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00042953 0.00144781  0.2967   0.7671
## 
## 
## $ln_M2
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00041622 0.00100076  0.4159   0.6781
## 
## 
## $ln_IHSG
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0015245  0.0029688 -0.5135   0.6084

Penanganan robust dilakukan untuk mengoreksi standard errors agar t-statistik dan p-value menjadi valid meskipun residual mengalami heteroskedastisitas

d. Normalitas

# Residual ln_IHSG
y_res <- res[, "ln_IHSG"]

# Uji Jarque–Bera
jb_test <- jarque.bera.test(y_res)
jb_test

## 
##  Jarque Bera Test
## 
## data:  y_res
## X-squared = 91.861, df = 2, p-value < 2.2e-16

residual tidak normal (pvalue < 0.05), inference konvensional bisa kurang valid, namun pada sampel besar (n > 100) efeknya relatif kecil dan koefisien VECM tetap konsisten

F. Peramalan

a. Forecast 6 periode kedepan

# Training data: Januari 2013 – Juni 2024
train_df <- final_data[1:(nrow(final_data)-6), ]
# Testing data: Juli – Desember 2024
test_df  <- final_data[(nrow(final_data)-5):nrow(final_data), ]

train_model <- train_df %>% dplyr::select(Inflasi, BI_rate, ln_Kurs, ln_M2, ln_IHSG)
test_model  <- test_df %>% dplyr::select(Inflasi, BI_rate, ln_Kurs, ln_M2, ln_IHSG)

# Estimasi VECM di train 
# lag optimal di train set
lag_sel <- VARselect(train_model, lag.max = 10, type = "const")
lag_optimal <- lag_sel$selection["AIC(n)"]

# Estimasi VAR untuk uji kointegrasi
var_train <- VAR(train_model, p = lag_optimal, type = "const")

# Johansen test
joh_trace <- ca.jo(train_model, type="trace", K=lag_optimal, ecdet="const", spec="transitory")

# VECM
vecm_model <- cajorls(joh_trace, r = 2)

# Konversi VECM ke VAR untuk forecast
var_model <- vec2var(joh_trace, r = 2)
forecast_horizon <- 6

vecm_forecast <- predict(var_model, n.ahead = forecast_horizon, ci = 0.95)
ihsg_forecast <- vecm_forecast$fcst$ln_IHSG

# Tanggal forecast 
forecast_dates <- test_df$Date2

forecast_table <- data.frame(
  Periode = forecast_dates,
  Hasil_Ramalan = round(exp(ihsg_forecast[, "fcst"]), 2),
  Lower_CI = round(exp(ihsg_forecast[, "lower"]), 2),
  Upper_CI = round(exp(ihsg_forecast[, "upper"]), 2),
  Data_Aktual = exp(test_df$ln_IHSG)  # kembalikan ke level IHSG
)
forecast_table

##      Periode Hasil_Ramalan Lower_CI Upper_CI Data_Aktual
## 1 2024-07-01       7134.80  6657.75  7646.03     7255.76
## 2 2024-08-01       7190.17  6484.38  7972.79     7670.73
## 3 2024-09-01       7236.75  6383.04  8204.63     7527.93
## 4 2024-10-01       7278.07  6320.57  8380.64     7574.02
## 5 2024-11-01       7315.96  6280.10  8522.67     7114.27
## 6 2024-12-01       7351.41  6252.86  8642.97     7079.90

b. plot forecast

# data historis IHSG dan forecast
plot_data <- data.frame(
  Periode = c(final_data$Date2, forecast_dates),
  IHSG = c(exp(final_data$ln_IHSG), forecast_table$Hasil_Ramalan),
  Type = c(rep("Actual", nrow(final_data)), rep("Forecast", forecast_horizon))
)

plot_data$Lower_CI <- c(rep(NA, nrow(final_data)), forecast_table$Lower_CI)
plot_data$Upper_CI <- c(rep(NA, nrow(final_data)), forecast_table$Upper_CI)
forecast_start <- forecast_dates[1]

# Plot
ggplot(plot_data, aes(x=Periode)) +
  geom_line(aes(y=IHSG, color=Type, group=1), size=1.2) +
  geom_point(aes(y=IHSG, color=Type), size=2) +
  geom_ribbon(aes(ymin=Lower_CI, ymax=Upper_CI), fill="blue", alpha=0.1) +
  geom_vline(xintercept=as.numeric(forecast_start), linetype="dashed", color="black") +
  labs(title="IHSG: Data Historis dan Forecast 6 Bulan Terakhir",
       y="IHSG", x="Waktu", color="Legend") +
  scale_color_manual(values=c("Actual"="black", "Forecast"="blue")) +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

c. Evaluasi

# MAPE
mape <- mean(abs(forecast_table$Hasil_Ramalan - forecast_table$Data_Aktual) / forecast_table$Data_Aktual) * 100
cat("MAPE =", round(mape, 2), "%\n")

## MAPE = 3.73 %