# =========================================================
# 1. PACKAGE
# =========================================================

library(readxl)
## Warning: package 'readxl' was built under R version 4.5.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(urca)
## Warning: package 'urca' was built under R version 4.5.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.5.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)
## Warning: package 'sandwich' was built under R version 4.5.3
library(strucchange)
## Warning: package 'strucchange' was built under R version 4.5.3
# =========================================================
# 2. IMPORT DATA
# =========================================================

data <- read_excel(
  "D:/Folder Farras/SEMESTER 4/ANALISIS RUNTUN WAKTU/uas/data kasus 3.xlsx",
  sheet = "Time_Series"
)

head(data)
## # A tibble: 6 × 10
##   Periode Tanggal             Tahun Triwulan Konsumsi_RT      PDB ln_Konsumsi_RT
##   <chr>   <dttm>              <dbl>    <dbl>       <dbl>    <dbl>          <dbl>
## 1 2010Q1  2010-03-31 00:00:00  2010        1     926098. 1642356.           13.7
## 2 2010Q2  2010-06-30 00:00:00  2010        2     937227. 1709132            13.8
## 3 2010Q3  2010-09-30 00:00:00  2010        3     959650. 1775110.           13.8
## 4 2010Q4  2010-12-31 00:00:00  2010        4     963089. 1737535.           13.8
## 5 2011Q1  2011-03-31 00:00:00  2011        1     964262. 1748731.           13.8
## 6 2011Q2  2011-06-30 00:00:00  2011        2     980157. 1816268.           13.8
## # ℹ 3 more variables: ln_PDB <dbl>, d_ln_Konsumsi_RT <dbl>, d_ln_PDB <dbl>
str(data)
## tibble [64 × 10] (S3: tbl_df/tbl/data.frame)
##  $ Periode         : chr [1:64] "2010Q1" "2010Q2" "2010Q3" "2010Q4" ...
##  $ Tanggal         : POSIXct[1:64], format: "2010-03-31" "2010-06-30" ...
##  $ Tahun           : num [1:64] 2010 2010 2010 2010 2011 ...
##  $ Triwulan        : num [1:64] 1 2 3 4 1 2 3 4 1 2 ...
##  $ Konsumsi_RT     : num [1:64] 926098 937227 959650 963089 964262 ...
##  $ PDB             : num [1:64] 1642356 1709132 1775110 1737535 1748731 ...
##  $ ln_Konsumsi_RT  : num [1:64] 13.7 13.8 13.8 13.8 13.8 ...
##  $ ln_PDB          : num [1:64] 14.3 14.4 14.4 14.4 14.4 ...
##  $ d_ln_Konsumsi_RT: num [1:64] NA 0.01195 0.02364 0.00358 0.00122 ...
##  $ d_ln_PDB        : num [1:64] NA 0.03985 0.03788 -0.02139 0.00642 ...
summary(data)
##    Periode             Tanggal                        Tahun         Triwulan   
##  Length:64          Min.   :2010-03-31 00:00:00   Min.   :2010   Min.   :1.00  
##  Class :character   1st Qu.:2014-03-08 12:00:00   1st Qu.:2014   1st Qu.:1.75  
##  Mode  :character   Median :2018-02-14 00:00:00   Median :2018   Median :2.50  
##                     Mean   :2018-02-13 18:00:00   Mean   :2018   Mean   :2.50  
##                     3rd Qu.:2022-01-22 12:00:00   3rd Qu.:2021   3rd Qu.:3.25  
##                     Max.   :2025-12-31 00:00:00   Max.   :2025   Max.   :4.00  
##                                                                                
##   Konsumsi_RT           PDB          ln_Konsumsi_RT      ln_PDB     
##  Min.   : 926098   Min.   :1642356   Min.   :13.74   Min.   :14.31  
##  1st Qu.:1131220   1st Qu.:2092345   1st Qu.:13.94   1st Qu.:14.55  
##  Median :1372887   Median :2530634   Median :14.13   Median :14.74  
##  Mean   :1348644   Mean   :2510329   Mean   :14.10   Mean   :14.72  
##  3rd Qu.:1513499   3rd Qu.:2826013   3rd Qu.:14.23   3rd Qu.:14.85  
##  Max.   :1819892   Max.   :3474461   Max.   :14.41   Max.   :15.06  
##                                                                     
##  d_ln_Konsumsi_RT        d_ln_PDB        
##  Min.   :-0.0674847   Min.   :-0.042804  
##  1st Qu.: 0.0005484   1st Qu.:-0.008712  
##  Median : 0.0064346   Median : 0.010442  
##  Mean   : 0.0107231   Mean   : 0.011894  
##  3rd Qu.: 0.0223945   3rd Qu.: 0.036038  
##  Max.   : 0.0458432   Max.   : 0.049240  
##  NA's   :1            NA's   :1
# =========================================================
# 3. MEMBUAT LOG / VARIABEL
# =========================================================

data <- data %>%
  mutate(
    ln_Konsumsi_RT = log(Konsumsi_RT),
    ln_PDB = log(PDB)
  )


# =========================================================
# 4. MEMBUAT TIME SERIES
# =========================================================

konsumsi <- ts(
  data$ln_Konsumsi_RT,
  start = c(2010, 1),
  frequency = 4
)

pdb <- ts(
  data$ln_PDB,
  start = c(2010, 1),
  frequency = 4
)


# =========================================================
# 5. UJI STASIONERITAS LEVEL
# =========================================================

adf.test(konsumsi)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  konsumsi
## Dickey-Fuller = -1.5984, Lag order = 3, p-value = 0.7377
## alternative hypothesis: stationary
adf.test(pdb)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  pdb
## Dickey-Fuller = -1.9887, Lag order = 3, p-value = 0.5795
## alternative hypothesis: stationary
# =========================================================
# 6. FIRST DIFFERENCE
# =========================================================

d_konsumsi <- diff(konsumsi)
d_pdb <- diff(pdb)

adf.test(d_konsumsi)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  d_konsumsi
## Dickey-Fuller = -3.3666, Lag order = 3, p-value = 0.06926
## alternative hypothesis: stationary
adf.test(d_pdb)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  d_pdb
## Dickey-Fuller = -3.2162, Lag order = 3, p-value = 0.09312
## alternative hypothesis: stationary
# =========================================================
# 7. ESTIMASI HUBUNGAN JANGKA PANJANG
# =========================================================

model_longrun <- lm(
  ln_Konsumsi_RT ~ ln_PDB,
  data = data
)

summary(model_longrun)
## 
## Call:
## lm(formula = ln_Konsumsi_RT ~ ln_PDB, data = data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0248602 -0.0084976  0.0000157  0.0073021  0.0249801 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.365758   0.116805   3.131  0.00265 ** 
## ln_PDB      0.933097   0.007937 117.570  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01278 on 62 degrees of freedom
## Multiple R-squared:  0.9955, Adjusted R-squared:  0.9955 
## F-statistic: 1.382e+04 on 1 and 62 DF,  p-value: < 2.2e-16
# =========================================================
# 8. MEMBENTUK RESIDUAL / ECT
# =========================================================

data$ECT <- residuals(model_longrun)


# =========================================================
# 9. UJI KOINTEGRASI
# =========================================================

adf.test(na.omit(data$ECT))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(data$ECT)
## Dickey-Fuller = -0.84948, Lag order = 3, p-value = 0.9525
## alternative hypothesis: stationary
# =========================================================
# 10. MEMBENTUK FIRST DIFFERENCE DAN ECT LAG
# =========================================================

data <- data %>%
  mutate(
    d_ln_Konsumsi = ln_Konsumsi_RT -
      lag(ln_Konsumsi_RT, 1),
    
    d_ln_PDB = ln_PDB -
      lag(ln_PDB, 1),
    
    ECT_lag = lag(ECT, 1)
  )


# =========================================================
# 11. ESTIMASI ECM
# =========================================================

model_ecm <- lm(
  d_ln_Konsumsi ~ d_ln_PDB + ECT_lag,
  data = data
)

summary(model_ecm)
## 
## Call:
## lm(formula = d_ln_Konsumsi ~ d_ln_PDB + ECT_lag, data = data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0306795 -0.0030449  0.0004639  0.0058199  0.0117619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.003865   0.001242   3.112  0.00285 ** 
## d_ln_PDB     0.580679   0.047422  12.245  < 2e-16 ***
## ECT_lag     -0.633102   0.087958  -7.198 1.15e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.008786 on 60 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7545, Adjusted R-squared:  0.7463 
## F-statistic: 92.18 on 2 and 60 DF,  p-value: < 2.2e-16
# =========================================================
# 12. UJI AUTOKORELASI
# =========================================================

dwtest(model_ecm)
## 
##  Durbin-Watson test
## 
## data:  model_ecm
## DW = 2.0506, p-value = 0.5637
## alternative hypothesis: true autocorrelation is greater than 0
bgtest(
  model_ecm,
  order = 4
)
## 
##  Breusch-Godfrey test for serial correlation of order up to 4
## 
## data:  model_ecm
## LM test = 25.247, df = 4, p-value = 4.487e-05
# =========================================================
# 13. UJI HETEROSKEDASTISITAS
# =========================================================

bptest(model_ecm)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ecm
## BP = 5.5302, df = 2, p-value = 0.06297
# =========================================================
# 14. UJI NORMALITAS
# =========================================================

jarque.bera.test(
  residuals(model_ecm)
)
## 
##  Jarque Bera Test
## 
## data:  residuals(model_ecm)
## X-squared = 24.24, df = 2, p-value = 5.45e-06
# =========================================================
# 15. STANDARD ERROR ROBUST
# =========================================================

coeftest(
  model_ecm,
  vcov = vcovHC(
    model_ecm,
    type = "HC1"
  )
)
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  0.0038648  0.0014785  2.6139    0.0113 *  
## d_ln_PDB     0.5806794  0.0645590  8.9946 1.008e-12 ***
## ECT_lag     -0.6331019  0.1026728 -6.1662 6.472e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1