R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(FinTS)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(rugarch) 
## Loading required package: parallel
library(tseries) 
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(dynlm)
library(vars) 
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library(nlWaldTest) 
library(lmtest) 
library(broom) 
library(car) 
## Loading required package: carData
library(sandwich)
library(knitr) 
library(forecast)
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:FinTS':
## 
##     Acf
library(fGarch) 
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")
library(aTSA) 
## 
## Attaching package: 'aTSA'
## The following object is masked from 'package:forecast':
## 
##     forecast
## The following object is masked from 'package:vars':
## 
##     arch.test
## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test
## The following object is masked from 'package:graphics':
## 
##     identify
data <- read.csv("C:/Users/Me/OneDrive/Dokumen/data_inflansi_fiks.csv", sep = ";")
str(data)
## 'data.frame':    165 obs. of  3 variables:
##  $ No            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Period        : chr  "Januari 2012" "Februari 2012" "Maret 2012" "Apr-12" ...
##  $ Inflation.Data: num  3.65 3.56 3.97 4.5 4.45 4.53 4.56 4.58 4.31 4.61 ...
head(data)
##   No        Period Inflation.Data
## 1  1  Januari 2012           3.65
## 2  2 Februari 2012           3.56
## 3  3    Maret 2012           3.97
## 4  4        Apr-12           4.50
## 5  5      Mei 2012           4.45
## 6  6     Juni 2012           4.53
attach(data)
Inflation.TS <- ts(data$Inflation.Data, start = c(2012, 1), frequency = 12)
n <- length(Inflation.TS)      # total observasi,  165
train_size <- floor(0.8 * n)   # 132 data untuk training
test_size <- n - train_size    # 33 data untuk testing

inflation.train <- Inflation.TS[1:train_size]          # data training (1–132)
inflation.test  <- Inflation.TS[(train_size + 1):n]    # data testing (133–165)
length(inflation.train)  # harus 132
## [1] 132
length(inflation.test)   # harus 33
## [1] 33
str(inflation.train)
##  num [1:132] 3.65 3.56 3.97 4.5 4.45 4.53 4.56 4.58 4.31 4.61 ...
plot.ts(inflation.train)

str(data)
## 'data.frame':    165 obs. of  3 variables:
##  $ No            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Period        : chr  "Januari 2012" "Februari 2012" "Maret 2012" "Apr-12" ...
##  $ Inflation.Data: num  3.65 3.56 3.97 4.5 4.45 4.53 4.56 4.58 4.31 4.61 ...
hist(inflation.train, breaks=20,
     freq=FALSE, col="grey")

library(tseries)
library(forecast)
#UJI STASIONER (data asli)
adf.test(inflation.train)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag    ADF p.value
## [1,]   0 -0.388   0.532
## [2,]   1 -0.615   0.459
## [3,]   2 -0.472   0.508
## [4,]   3 -0.559   0.479
## [5,]   4 -0.514   0.495
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -1.56   0.500
## [2,]   1 -2.02   0.320
## [3,]   2 -1.65   0.465
## [4,]   3 -1.71   0.442
## [5,]   4 -1.66   0.462
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -1.82   0.646
## [2,]   1 -2.57   0.338
## [3,]   2 -1.95   0.592
## [4,]   3 -1.99   0.578
## [5,]   4 -1.99   0.575
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01
#UJI STATIONERITAS [data differencing ke-1]
adf.test(diff(inflation.train))
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -8.99    0.01
## [2,]   1 -8.42    0.01
## [3,]   2 -6.54    0.01
## [4,]   3 -5.76    0.01
## [5,]   4 -4.91    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -8.96    0.01
## [2,]   1 -8.39    0.01
## [3,]   2 -6.52    0.01
## [4,]   3 -5.74    0.01
## [5,]   4 -4.89    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -8.93    0.01
## [2,]   1 -8.36    0.01
## [3,]   2 -6.51    0.01
## [4,]   3 -5.74    0.01
## [5,]   4 -4.89    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01
acf(diff(inflation.train))

ndiffs(diff(inflation.train))
## [1] 0
pacf(diff(inflation.train))

#Identifikasi Model ARIMA
auto.arima(inflation.train,trace = T,d = 1)
## 
##  ARIMA(2,1,2) with drift         : 214.778
##  ARIMA(0,1,0) with drift         : 217.7413
##  ARIMA(1,1,0) with drift         : 212.8452
##  ARIMA(0,1,1) with drift         : 209.9714
##  ARIMA(0,1,0)                    : 215.7668
##  ARIMA(1,1,1) with drift         : 211.1361
##  ARIMA(0,1,2) with drift         : 210.7404
##  ARIMA(1,1,2) with drift         : 212.6486
##  ARIMA(0,1,1)                    : 207.9304
##  ARIMA(1,1,1)                    : 209.0681
##  ARIMA(0,1,2)                    : 208.6793
##  ARIMA(1,1,0)                    : 210.8053
##  ARIMA(1,1,2)                    : 210.5587
## 
##  Best model: ARIMA(0,1,1)
## Series: inflation.train 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.3122
## s.e.  0.0909
## 
## sigma^2 = 0.2794:  log likelihood = -101.92
## AIC=207.84   AICc=207.93   BIC=213.59
#Pemilihan Model
modelarima012<-arima(inflation.train,
                     order=c(0,1,2), method="ML")
(modelarima012)
## 
## Call:
## arima(x = inflation.train, order = c(0, 1, 2), method = "ML")
## 
## Coefficients:
##          ma1      ma2
##       0.2704  -0.1019
## s.e.  0.0884   0.0870
## 
## sigma^2 estimated as 0.2744:  log likelihood = -101.25,  aic = 208.49
modelarima210<-arima(inflation.train,
                     order=c(2,1,0), method="ML")
(modelarima210)
## 
## Call:
## arima(x = inflation.train, order = c(2, 1, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2
##       0.2695  -0.1798
## s.e.  0.0857   0.0856
## 
## sigma^2 estimated as 0.2742:  log likelihood = -101.19,  aic = 208.37
modelarima011<-arima(inflation.train,
                     order=c(0,1,1), method="ML")
(modelarima011)
## 
## Call:
## arima(x = inflation.train, order = c(0, 1, 1), method = "ML")
## 
## Coefficients:
##          ma1
##       0.3122
## s.e.  0.0909
## 
## sigma^2 estimated as 0.2773:  log likelihood = -101.92,  aic = 207.84
#Analisis Residual
checkresiduals(modelarima210$residuals)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 1.6356, df = 10, p-value = 0.9984
## 
## Model df: 0.   Total lags used: 10
jarque.bera.test(residuals(modelarima210))
## 
##  Jarque Bera Test
## 
## data:  residuals(modelarima210)
## X-squared = 267.74, df = 2, p-value < 2.2e-16
library(aTSA)
arch.test(modelarima210)
## ARCH heteroscedasticity test for residuals 
## alternative: heteroscedastic 
## 
## Portmanteau-Q test: 
##      order    PQ p.value
## [1,]     4  5.64 0.22811
## [2,]     8  9.49 0.30283
## [3,]    12 27.01 0.00770
## [4,]    16 32.47 0.00869
## [5,]    20 48.35 0.00038
## [6,]    24 48.99 0.00190
## Lagrange-Multiplier test: 
##      order      LM  p.value
## [1,]     4 167.023 0.00e+00
## [2,]     8  65.457 1.22e-11
## [3,]    12  38.080 7.58e-05
## [4,]    16  12.383 6.50e-01
## [5,]    20   1.966 1.00e+00
## [6,]    24   0.946 1.00e+00

#Pemeriksaan efek
library(FinTS)
alpha <- 0.05
q <- length(coef(modelarima210))-1
Chicr <- qchisq(1-alpha, q)
Chicr
## [1] 3.841459
infArchTest <-
  ArchTest(modelarima210$residuals, lags = 1,
           demean = TRUE)
infArchTest
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  modelarima210$residuals
## Chi-squared = 5.174, df = 1, p-value = 0.02293
for (i in 1:10){
  ArchTest <- ArchTest(diff(inflation.train), lags=i, demean=TRUE)
  cat("P Value LM Test lag ke", i,"adalah" , ArchTest$p.value, "\n")}
## P Value LM Test lag ke 1 adalah 0.03361739 
## P Value LM Test lag ke 2 adalah 0.1052701 
## P Value LM Test lag ke 3 adalah 0.2006683 
## P Value LM Test lag ke 4 adalah 0.2446249 
## P Value LM Test lag ke 5 adalah 0.1090221 
## P Value LM Test lag ke 6 adalah 0.1679345 
## P Value LM Test lag ke 7 adalah 0.2381361 
## P Value LM Test lag ke 8 adalah 0.3241266 
## P Value LM Test lag ke 9 adalah 0.3959315 
## P Value LM Test lag ke 10 adalah 0.494796
#Pemodelan ARCH/GARCH
library(fGarch)
#arch1
arch1<-
  garchFit(~arma(2,0)+garch(1,0),data = diff(inflation.train), trace = F)
library(texreg)
## Version:  1.39.4
## Date:     2024-07-23
## Author:   Philip Leifeld (University of Manchester)
## 
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").
extract(arch1,
        include.nobs = FALSE,
        include.aic = TRUE,
        include.loglik = FALSE)
## 
## 
##                coef.       s.e.            p
## mu     -0.0133485909 0.03002730 6.566461e-01
## ar1     0.3514101504 0.07583747 3.591406e-06
## ar2    -0.0003227572 0.10631846 9.975778e-01
## omega   0.0911342610 0.02531067 3.174522e-04
## alpha1  0.9999999900 0.30187505 9.242357e-04
## 
##          GOF dec. places
## AIC 1.303773        TRUE
#garch(1,1)
garch1<-
  garchFit(~arma(2,0)+garch(1,1),data = diff(inflation.train), trace = F)
library(texreg)
extract(garch1,include.nobs = FALSE,include.aic = TRUE,include.loglik = FALSE)
## 
## 
##              coef.       s.e.            p
## mu     -0.02599381 0.02327261 2.640255e-01
## ar1     0.34342135 0.10206110 7.658240e-04
## ar2     0.06604696 0.07689285 3.903692e-01
## omega   0.03214646 0.01305048 1.376875e-02
## alpha1  0.99999999 0.21355094 2.830885e-06
## beta1   0.22212530 0.08148402 6.410733e-03
## 
##          GOF dec. places
## AIC 1.248905        TRUE
#(garch1,2)
garch12<-
  garchFit(~arma(2,0)+garch(1,2),data = diff(inflation.train), trace = F)
library(texreg)
extract(garch12,include.nobs = FALSE,include.aic = TRUE,include.loglik = FALSE)
## 
## 
##              coef.       s.e.            p
## mu     -0.02060143 0.02492805 4.085570e-01
## ar1     0.32761137 0.10237823 1.374227e-03
## ar2     0.04283613 0.07590821 5.725399e-01
## omega   0.03044215 0.01270027 1.653132e-02
## alpha1  0.99999999 0.22022995 5.606551e-06
## beta1   0.10822442 0.06877296 1.155689e-01
## beta2   0.09372267 0.08389276 2.639207e-01
## 
##          GOF dec. places
## AIC 1.253429        TRUE
#garch(2,1)
garch21<-
  garchFit(~arma(2,0)+garch(2,1),data = diff(inflation.train), trace = F)
library(texreg)
extract(garch21,include.nobs = FALSE,include.aic = TRUE,include.loglik = FALSE)
## 
## 
##              coef.       s.e.            p
## mu     -0.02794780 0.02354105 2.351509e-01
## ar1     0.32809650 0.10273026 1.404294e-03
## ar2     0.07261092 0.07706788 3.461064e-01
## omega   0.03072632 0.01301493 1.823301e-02
## alpha1  0.99999999 0.21350902 2.818208e-06
## alpha2  0.00000001 0.22225759 1.000000e+00
## beta1   0.22237737 0.14329872 1.206994e-01
## 
##          GOF dec. places
## AIC 1.267143        TRUE
#Model ARIMA (0,1,1)-GARCH(1,2)
garch12<-
  garchFit(~arma(2,0)+garch(1,2),data = diff(inflation.train), trace = F)
library(texreg)
extract(garch12,include.nobs = FALSE,include.aic = TRUE,include.loglik = FALSE)
## 
## 
##              coef.       s.e.            p
## mu     -0.02060143 0.02492805 4.085570e-01
## ar1     0.32761137 0.10237823 1.374227e-03
## ar2     0.04283613 0.07590821 5.725399e-01
## omega   0.03044215 0.01270027 1.653132e-02
## alpha1  0.99999999 0.22022995 5.606551e-06
## beta1   0.10822442 0.06877296 1.155689e-01
## beta2   0.09372267 0.08389276 2.639207e-01
## 
##          GOF dec. places
## AIC 1.253429        TRUE
#Diagnostik GARCH (1,2)
plot(residuals(garch12), type='h', ylab="Standardized Residuals")

qqnorm(residuals(garch12)); qqline(residuals(garch12))

# Pendugaan variansi
fitted.volatility <- volatility(garch12)
plot(fitted.volatility, type="l", xlab="", ylab="", xaxt="n")

#PERAMALAN ARIMA (2,0,1)-GARCH(1,2)
library(rugarch)
garch012.spec <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 1), include.mean = TRUE),
  distribution.model = "norm"
)
garch012.fit <- ugarchfit(spec = garch012.spec, data = diff(inflation.train))
f.garch <- ugarchforecast(garch012.fit, n.ahead = 11)
f.garch
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 11
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=1970-05-12]:
##        Series  Sigma
## T+1   0.01184 0.3711
## T+2  -0.02137 0.4496
## T+3  -0.02137 0.4859
## T+4  -0.02137 0.5236
## T+5  -0.02137 0.5583
## T+6  -0.02137 0.5909
## T+7  -0.02137 0.6218
## T+8  -0.02137 0.6512
## T+9  -0.02137 0.6793
## T+10 -0.02137 0.7063
## T+11 -0.02137 0.7322
# PERAMALAN ARIMA (0,1,1)

# Data testing (33 observasi)
inflation.test <- Inflation.TS[(train_size + 1):n]

# Cek panjang data testing
length(inflation.test)
## [1] 33
# Forecast 11 periode ke depan
f.arima <- predict(modelarima011, n.ahead = 11)

# Ambil 11 observasi pertama dari data testing
actual <- inflation.test[1:11]

# Buat tabel perbandingan
hasil_perbandingan <- data.frame(
  Actual = as.numeric(actual),
  Forecast = as.numeric(f.arima$pred),
  Variance = as.numeric(f.arima$se^2)
)

# Tampilkan hasil
print(hasil_perbandingan)
##    Actual Forecast  Variance
## 1    5.28 5.546046 0.2773001
## 2    5.47 5.546046 0.7548095
## 3    4.97 5.546046 1.2323189
## 4    4.33 5.546046 1.7098282
## 5    4.00 5.546046 2.1873376
## 6    3.52 5.546046 2.6648470
## 7    3.08 5.546046 3.1423563
## 8    3.27 5.546046 3.6198657
## 9    2.28 5.546046 4.0973751
## 10   2.56 5.546046 4.5748845
## 11   2.86 5.546046 5.0523938
# HITUNG MAPE

MAPE <- mean(
  abs(
    (hasil_perbandingan$Actual -
       hasil_perbandingan$Forecast) /
      hasil_perbandingan$Actual
  )
) * 100

cat("Nilai MAPE =", round(MAPE, 4), "%\n")
## Nilai MAPE = 58.7083 %
# HITUNG RMSE

RMSE <- sqrt(
  mean(
    (hasil_perbandingan$Actual -
       hasil_perbandingan$Forecast)^2
  )
)

cat("Nilai RMSE =", round(RMSE, 4), "\n")
## Nilai RMSE = 2.0548
# HITUNG MAE

MAE <- mean(
  abs(
    hasil_perbandingan$Actual -
      hasil_perbandingan$Forecast
  )
)

cat("Nilai MAE =", round(MAE, 4), "\n")
## Nilai MAE = 1.7624
# PLOT HASIL FORECAST VS AKTUAL

plot(
  actual,
  type = "l",
  col = "blue",
  lwd = 2,
  ylab = "Inflasi",
  xlab = "Periode",
  main = "Perbandingan Aktual dan Forecast ARIMA(0,1,1)"
)

lines(
  f.arima$pred,
  col = "red",
  lwd = 2
)

legend(
  "topright",
  legend = c("Aktual", "Forecast"),
  col = c("blue", "red"),
  lwd = 2
)

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.