library(arfima)
## Warning: package 'arfima' was built under R version 4.3.3
## Loading required package: ltsa
## Warning: package 'ltsa' was built under R version 4.3.3
## Note that the arfima package has new defaults starting with
## 1.4-0: type arfimachanges() for a list, as well as some other notes.
## NOTE: some of these are quite important!
##
## Attaching package: 'arfima'
## The following object is masked from 'package:stats':
##
## BIC
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
##
## Attaching package: 'forecast'
## The following object is masked from 'package:arfima':
##
## arfima
library(fracdiff)
## Warning: package 'fracdiff' was built under R version 4.3.3
library(pracma)
## Warning: package 'pracma' was built under R version 4.3.3
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
library(TSA)
## Warning: package 'TSA' was built under R version 4.3.3
## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library(nortest)
df <- read.csv("D:/SEMESTER 5/Metopen/Data Jabodetabekk.csv", sep = ",")
head(df)
## X Jabodetabek
## 1 2019-1 27768
## 2 2019-2 25305
## 3 2019-3 28366
## 4 2019-4 28062
## 5 2019-5 28369
## 6 2019-6 25816
data <- df[,2]
plot(data, type = "l",xlab = "Periode",
ylab = "Jumlah Penumpang", main = "Plot Jumlah Penumpang Kereta Api Jabodetabek")
lamda <- BoxCox.lambda(data) ; lamda
## [1] 1.999924
ts_data <- BoxCox(data, lamda)
lamda1 <- BoxCox.lambda(ts_data) ; lamda1
## [1] 1.44628
ts_data2 <- BoxCox(ts_data, lamda1)
lamda2 <- BoxCox.lambda(ts_data2) ; lamda2
## [1] 0.9522903
datats <- ts_data2
acf(datats)
hurstexp(datats)
## Simple R/S Hurst estimation: 0.7870558
## Corrected R over S Hurst exponent: 1.037753
## Empirical Hurst exponent: 0.7493313
## Corrected empirical Hurst exponent: 0.7678741
## Theoretical Hurst exponent: 0.5197842
dgph <- fdGPH(datats)$d
dgphta <- fdGPH(spec.taper(datats))$d
dspr <- fdSperio(datats)$d
cat("Nilai d GPH :", dgph,"Nilai d GPH-TA :", dgphta, "Nilai d SPR :", dspr, "\n")
## Nilai d GPH : 1.178347 Nilai d GPH-TA : 0.9887196 Nilai d SPR : 0.8141534
if(dgph > 1){
datats <- diff(datats) # buat stasioner
d_GPH <- fdGPH(datats)$d # update nilai d
} else {
d_GPH <- dgph$d
}
Karena hasil d GPH > 1, dilakukan pemecahan differencing. Differencing angka bulat dan differencing fractional
datagph <- diff(datats)
d_GPH <- fdGPH(datagph)$d ; d_GPH
## [1] -0.6700006
cat("Nilai d GPH setelah dilakukan differencing :", d_GPH, "\n")
## Nilai d GPH setelah dilakukan differencing : -0.6700006
diff_gph <- diffseries(datagph, d_GPH)
acf(diff_gph); pacf(diff_gph)
#Differencing dari d GPH-TA
diff_gphta <- diffseries(datats, dgphta)
acf(diff_gphta) ; pacf(diff_gphta)
#Differencing dari d SPR
diff_spr <- diffseries(datats, dspr)
acf(diff_spr) ; pacf(diff_spr)
#Pemodelan ARFIMA #Model dengan d GPH
fit1 <- arfima::arfima(z=datats, order = c(2,0,2), dmean = FALSE,
fixed = list(theta = c(0,NA),
phi = c(0,NA),
d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit1)
##
## Call:
##
## arfima::arfima(z = datats, order = c(2, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), phi = c(0, NA), d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) 0.00000e+00 NA 3.90000e-01 NA NA
## phi(2) 2.51690e-01 2.84165e-01 3.07522e-01 0.88572 0.3757695
## theta(1) 0.00000e+00 NA 2.35053e-01 NA NA
## theta(2) -2.50780e-01 2.66646e-01 2.40375e-01 -0.94050 0.3469623
## d.f -3.14205e-01 1.10635e-01 2.65876e-01 -2.84003 0.0045109 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 7.0318e+22; Log-likelihood = -2075.76; AIC = 4165.51; BIC = 4182.1
##
## Numerical Correlations of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 0.00 NaN 0.00 NaN NaN
## phi(2) NaN 1.00 NaN 0.87 -0.41
## theta(1) 0.00 NaN 0.00 NaN NaN
## theta(2) NaN 0.87 NaN 1.00 -0.15
## d.f NaN -0.41 NaN -0.15 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 1.00 0.53 0.67 0.23 -0.81
## phi(2) 0.53 1.00 0.13 0.82 -0.66
## theta(1) 0.67 0.13 1.00 0.06 -0.20
## theta(2) 0.23 0.82 0.06 1.00 -0.28
## d.f -0.81 -0.66 -0.20 -0.28 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 1.06 0.00 -0.94 0.00 1.09
## phi(2) 0.00 1.07 0.00 -0.94 0.58
## theta(1) -0.94 0.00 1.06 0.00 -0.92
## theta(2) 0.00 -0.94 0.00 1.07 -0.45
## d.f 1.09 0.58 -0.92 -0.45 1.64
fit2 <- arfima::arfima(z=datats, order = c(2,0,0), dmean = FALSE,
fixed = list(phi = c(0,NA),
d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit2)
##
## Call:
##
## arfima::arfima(z = datats, order = c(2, 0, 0), dmean = FALSE, fixed = list(phi = c(0,
## NA), d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) 0.00000e+00 NA 3.24654e-01 NA NA
## phi(2) 4.51083e-01 1.29112e-01 1.96926e-01 3.49373 0.00047633 ***
## d.f -3.18301e-01 1.10363e-01 3.19844e-01 -2.88413 0.00392499 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.94133e+22; Log-likelihood = -2076.25; AIC = 4162.5; BIC = 4174.35
##
## Numerical Correlations of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 0.00 NaN NaN
## phi(2) NaN 1.00 -0.54
## d.f NaN -0.54 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 1.00 0.82 -0.95
## phi(2) 0.82 1.00 -0.86
## d.f -0.95 -0.86 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 1.26 0.00 1.22
## phi(2) 0.00 1.26 0.67
## d.f 1.22 0.67 1.65
fit3 <- arfima::arfima(z=datats, order = c(0,0,2), dmean = FALSE,
fixed = list(theta = c(0,NA),
d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit3)
##
## Call:
##
## arfima::arfima(z = datats, order = c(0, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## theta(1) 0.00000e+00 NA 1.30996e-01 NA NA
## theta(2) -4.48667e-01 1.13842e-01 1.07981e-01 -3.94114 8.1094e-05 ***
## d.f -2.77259e-01 9.45304e-02 1.19373e-01 -2.93301 0.0033569 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.9056e+22; Log-likelihood = -2076.07; AIC = 4162.14; BIC = 4173.99
##
## Numerical Correlations of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 0.00 NaN NaN
## theta(2) NaN 1.00 0.36
## d.f NaN 0.36 1.00
##
## Theoretical Correlations of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 1.00 0.23 0.64
## theta(2) 0.23 1.00 0.36
## d.f 0.64 0.36 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 1.25 0.00 -0.88
## theta(2) 0.00 1.25 -0.41
## d.f -0.88 -0.41 1.64
fit4 <- arfima::arfima(z=datats, order = c(1,0,1), dmean = FALSE,
fixed = list(d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit4)
##
## Call:
##
## arfima::arfima(z = datats, order = c(1, 0, 1), dmean = FALSE, fixed = list(d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) -6.04726e-01 1.67212e-01 1.96596e-01 -3.61651 0.0002986 ***
## theta(1) -1.55768e-01 2.33602e-01 3.05221e-01 -0.66681 0.5048939
## d.f 1.34192e-01 1.53264e-01 1.36602e-01 0.87556 0.3812716
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 7.19594e+22; Log-likelihood = -2077.53; AIC = 4165.05; BIC = 4176.9
##
## Numerical Correlations of Coefficients:
## phi(1) theta(1) d.f
## phi(1) 1.00 0.71 0.13
## theta(1) 0.71 1.00 0.64
## d.f 0.13 0.64 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) theta(1) d.f
## phi(1) 1.00 0.85 0.39
## theta(1) 0.85 1.00 0.68
## d.f 0.39 0.68 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) theta(1) d.f
## phi(1) 1.58 -1.11 0.78
## theta(1) -1.11 1.03 -0.93
## d.f 0.78 -0.93 1.65
fit5 <- arfima::arfima(z=datats, order = c(1,0,0), dmean = FALSE,
fixed = list(d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit5)
##
## Call:
##
## arfima::arfima(z = datats, order = c(1, 0, 0), dmean = FALSE, fixed = list(d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) -5.23589e-01 1.34201e-01 1.13402e-01 -3.90152 9.5592e-05 ***
## d.f 1.99860e-01 1.28413e-01 1.03779e-01 1.55638 0.11962
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 7.13145e+22; Log-likelihood = -2077.74; AIC = 4163.49; BIC = 4172.97
##
## Numerical Correlations of Coefficients:
## phi(1) d.f
## phi(1) 1.00 -0.66
## d.f -0.66 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) d.f
## phi(1) 1.00 -0.53
## d.f -0.53 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) d.f
## phi(1) 1.38 0.80
## d.f 0.80 1.64
fit6 <- arfima::arfima(z=datats, order = c(0,0,1), dmean = FALSE,
fixed = list(d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit6)
##
## Call:
##
## arfima::arfima(z = datats, order = c(0, 0, 1), dmean = FALSE, fixed = list(d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## theta(1) 2.76277e-01 1.64560e-01 2.25055e-01 1.67888 0.093175 .
## d.f 5.17299e-02 1.65891e-01 1.82592e-01 0.31183 0.755170
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 7.69614e+22; Log-likelihood = -2080.57; AIC = 4169.15; BIC = 4178.62
##
## Numerical Correlations of Coefficients:
## theta(1) d.f
## theta(1) 1.00 0.84
## d.f 0.84 1.00
##
## Theoretical Correlations of Coefficients:
## theta(1) d.f
## theta(1) 1.00 0.88
## d.f 0.88 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## theta(1) d.f
## theta(1) 1.08 -1.17
## d.f -1.17 1.64
fit7 <- arfima::arfima(z=datats, order = c(1,0,2), dmean = FALSE,
fixed = list(theta = c(0,NA),
d = dgph))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit7)
##
## Call:
##
## arfima::arfima(z = datats, order = c(1, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), d = dgph))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) -2.11191e-01 2.07629e-01 3.25269e-01 -1.01715 0.309081
## theta(1) 0.00000e+00 NA 2.45784e-01 NA NA
## theta(2) -3.84495e-01 1.48358e-01 1.45362e-01 -2.59167 0.009551 **
## d.f -1.20223e-01 1.87643e-01 1.55981e-01 -0.64070 0.521717
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.93163e+22; Log-likelihood = -2075.59; AIC = 4163.19; BIC = 4177.4
##
## Numerical Correlations of Coefficients:
## phi(1) theta(1) theta(2) d.f
## phi(1) 1.00 NaN -0.56 -0.84
## theta(1) NaN 0.00 NaN NaN
## theta(2) -0.56 NaN 1.00 0.65
## d.f -0.84 NaN 0.65 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) theta(1) theta(2) d.f
## phi(1) 1.00 0.82 -0.63 -0.60
## theta(1) 0.82 1.00 -0.39 -0.18
## theta(2) -0.63 -0.39 1.00 0.62
## d.f -0.60 -0.18 0.62 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) theta(1) theta(2) d.f
## phi(1) 1.04 -0.98 0.21 0.90
## theta(1) -0.98 1.17 0.00 -0.89
## theta(2) 0.21 0.00 1.17 -0.42
## d.f 0.90 -0.89 -0.42 1.64
#Model dengan d GPH TA
fit8 <- arfima::arfima(z=datats, order = c(2,0,2), dmean = FALSE,
fixed = list(theta = c(0, NA),
phi = c(0, NA),
d = dgphta))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit8)
##
## Call:
##
## arfima::arfima(z = datats, order = c(2, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), phi = c(0, NA), d = dgphta))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) 0.00000e+00 NA 3.90000e-01 NA NA
## phi(2) 2.51690e-01 2.84165e-01 3.07522e-01 0.88572 0.3757695
## theta(1) 0.00000e+00 NA 2.35053e-01 NA NA
## theta(2) -2.50780e-01 2.66646e-01 2.40375e-01 -0.94050 0.3469623
## d.f -3.14205e-01 1.10635e-01 2.65876e-01 -2.84003 0.0045109 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 7.0318e+22; Log-likelihood = -2075.76; AIC = 4165.51; BIC = 4182.1
##
## Numerical Correlations of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 0.00 NaN 0.00 NaN NaN
## phi(2) NaN 1.00 NaN 0.87 -0.41
## theta(1) 0.00 NaN 0.00 NaN NaN
## theta(2) NaN 0.87 NaN 1.00 -0.15
## d.f NaN -0.41 NaN -0.15 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 1.00 0.53 0.67 0.23 -0.81
## phi(2) 0.53 1.00 0.13 0.82 -0.66
## theta(1) 0.67 0.13 1.00 0.06 -0.20
## theta(2) 0.23 0.82 0.06 1.00 -0.28
## d.f -0.81 -0.66 -0.20 -0.28 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 1.06 0.00 -0.94 0.00 1.09
## phi(2) 0.00 1.07 0.00 -0.94 0.58
## theta(1) -0.94 0.00 1.06 0.00 -0.92
## theta(2) 0.00 -0.94 0.00 1.07 -0.45
## d.f 1.09 0.58 -0.92 -0.45 1.64
fit9 <- arfima::arfima(z=datats, order = c(2,0,0), dmean = FALSE,
fixed = list(phi = c(0, NA),
d = dgphta))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit9)
##
## Call:
##
## arfima::arfima(z = datats, order = c(2, 0, 0), dmean = FALSE, fixed = list(phi = c(0,
## NA), d = dgphta))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) 0.00000e+00 NA 3.24654e-01 NA NA
## phi(2) 4.51083e-01 1.29112e-01 1.96926e-01 3.49373 0.00047633 ***
## d.f -3.18301e-01 1.10363e-01 3.19844e-01 -2.88413 0.00392499 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.94133e+22; Log-likelihood = -2076.25; AIC = 4162.5; BIC = 4174.35
##
## Numerical Correlations of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 0.00 NaN NaN
## phi(2) NaN 1.00 -0.54
## d.f NaN -0.54 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 1.00 0.82 -0.95
## phi(2) 0.82 1.00 -0.86
## d.f -0.95 -0.86 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 1.26 0.00 1.22
## phi(2) 0.00 1.26 0.67
## d.f 1.22 0.67 1.65
fit10 <- arfima::arfima(z=datats, order = c(0,0,2), dmean = FALSE,
fixed = list(theta = c(0, NA),
d = dgphta))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit10)
##
## Call:
##
## arfima::arfima(z = datats, order = c(0, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), d = dgphta))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## theta(1) 0.00000e+00 NA 1.30996e-01 NA NA
## theta(2) -4.48667e-01 1.13842e-01 1.07981e-01 -3.94114 8.1094e-05 ***
## d.f -2.77259e-01 9.45304e-02 1.19373e-01 -2.93301 0.0033569 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.9056e+22; Log-likelihood = -2076.07; AIC = 4162.14; BIC = 4173.99
##
## Numerical Correlations of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 0.00 NaN NaN
## theta(2) NaN 1.00 0.36
## d.f NaN 0.36 1.00
##
## Theoretical Correlations of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 1.00 0.23 0.64
## theta(2) 0.23 1.00 0.36
## d.f 0.64 0.36 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 1.25 0.00 -0.88
## theta(2) 0.00 1.25 -0.41
## d.f -0.88 -0.41 1.64
#Model dengan d SPR
fit11 <- arfima::arfima(z=datats, order = c(2,0,2), dmean = FALSE,
fixed = list(theta = c(0, NA),
phi = c(0, NA),
d = dspr))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit11)
##
## Call:
##
## arfima::arfima(z = datats, order = c(2, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), phi = c(0, NA), d = dspr))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) 0.00000e+00 NA 3.90000e-01 NA NA
## phi(2) 2.51690e-01 2.84165e-01 3.07522e-01 0.88572 0.3757695
## theta(1) 0.00000e+00 NA 2.35053e-01 NA NA
## theta(2) -2.50780e-01 2.66646e-01 2.40375e-01 -0.94050 0.3469623
## d.f -3.14205e-01 1.10635e-01 2.65876e-01 -2.84003 0.0045109 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 7.0318e+22; Log-likelihood = -2075.76; AIC = 4165.51; BIC = 4182.1
##
## Numerical Correlations of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 0.00 NaN 0.00 NaN NaN
## phi(2) NaN 1.00 NaN 0.87 -0.41
## theta(1) 0.00 NaN 0.00 NaN NaN
## theta(2) NaN 0.87 NaN 1.00 -0.15
## d.f NaN -0.41 NaN -0.15 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 1.00 0.53 0.67 0.23 -0.81
## phi(2) 0.53 1.00 0.13 0.82 -0.66
## theta(1) 0.67 0.13 1.00 0.06 -0.20
## theta(2) 0.23 0.82 0.06 1.00 -0.28
## d.f -0.81 -0.66 -0.20 -0.28 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) phi(2) theta(1) theta(2) d.f
## phi(1) 1.06 0.00 -0.94 0.00 1.09
## phi(2) 0.00 1.07 0.00 -0.94 0.58
## theta(1) -0.94 0.00 1.06 0.00 -0.92
## theta(2) 0.00 -0.94 0.00 1.07 -0.45
## d.f 1.09 0.58 -0.92 -0.45 1.64
fit12 <- arfima::arfima(z=datats, order = c(2,0,0), dmean = FALSE,
fixed = list(phi = c(0, NA),
d = dspr))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit12)
##
## Call:
##
## arfima::arfima(z = datats, order = c(2, 0, 0), dmean = FALSE, fixed = list(phi = c(0,
## NA), d = dspr))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## phi(1) 0.00000e+00 NA 3.24654e-01 NA NA
## phi(2) 4.51083e-01 1.29112e-01 1.96926e-01 3.49373 0.00047633 ***
## d.f -3.18301e-01 1.10363e-01 3.19844e-01 -2.88413 0.00392499 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.94133e+22; Log-likelihood = -2076.25; AIC = 4162.5; BIC = 4174.35
##
## Numerical Correlations of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 0.00 NaN NaN
## phi(2) NaN 1.00 -0.54
## d.f NaN -0.54 1.00
##
## Theoretical Correlations of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 1.00 0.82 -0.95
## phi(2) 0.82 1.00 -0.86
## d.f -0.95 -0.86 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## phi(1) phi(2) d.f
## phi(1) 1.26 0.00 1.22
## phi(2) 0.00 1.26 0.67
## d.f 1.22 0.67 1.65
fit13 <- arfima::arfima(z=datats, order = c(0,0,2), dmean = FALSE,
fixed = list(theta = c(0, NA),
d = dspr))
## Note: only one starting point. Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.
summary(fit13)
##
## Call:
##
## arfima::arfima(z = datats, order = c(0, 0, 2), dmean = FALSE, fixed = list(theta = c(0,
## NA), d = dspr))
##
##
## Mode 1 Coefficients:
## Estimate Std. Error Th. Std. Err. z-value Pr(>|z|)
## theta(1) 0.00000e+00 NA 1.30996e-01 NA NA
## theta(2) -4.48667e-01 1.13842e-01 1.07981e-01 -3.94114 8.1094e-05 ***
## d.f -2.77259e-01 9.45304e-02 1.19373e-01 -2.93301 0.0033569 **
## zbar 2.92454e+09 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 6.9056e+22; Log-likelihood = -2076.07; AIC = 4162.14; BIC = 4173.99
##
## Numerical Correlations of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 0.00 NaN NaN
## theta(2) NaN 1.00 0.36
## d.f NaN 0.36 1.00
##
## Theoretical Correlations of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 1.00 0.23 0.64
## theta(2) 0.23 1.00 0.36
## d.f 0.64 0.36 1.00
##
## Expected Fisher Information Matrix of Coefficients:
## theta(1) theta(2) d.f
## theta(1) 1.25 0.00 -0.88
## theta(2) 0.00 1.25 -0.41
## d.f -0.88 -0.41 1.64
#Pengujian Asumsi #Uji Asumsi Model Fit2 #Uji Normalitas
r2 <- residuals(fit2)
str(r2)
## List of 1
## $ Mode1: num [1:79] -3.76e+11 3.48e+11 2.10e+11 -9.87e+10 -4.30e+11 ...
n <- length(r2$Mode1)
mean <- mean(r2$Mode1)
sd <- sd(r2$Mode1)
res <- rnorm(n,mean,sd)
ks.test(r2$Mode1,res)
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: r2$Mode1 and res
## D = 0.10127, p-value = 0.8161
## alternative hypothesis: two-sided
#Uji Autokorelasi
Box.test(r2$Mode1,type=c("Ljung-Box"))
##
## Box-Ljung test
##
## data: r2$Mode1
## X-squared = 0.020934, df = 1, p-value = 0.885
#Uji Asumsi Model Fit3 #Uji Normalitas
r3 <- residuals(fit3)
str(r3)
## List of 1
## $ Mode1: num [1:79] -3.84e+11 3.76e+11 1.68e+11 -7.22e+10 -4.94e+11 ...
n <- length(r3$Mode1)
mean <- mean(r3$Mode1)
sd <- sd(r3$Mode1)
res <- rnorm(n,mean,sd)
ks.test(r3$Mode1,res)
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: r3$Mode1 and res
## D = 0.27848, p-value = 0.00419
## alternative hypothesis: two-sided
#Uji Autokorelasi
Box.test(r3$Mode1,type=c("Ljung-Box"))
##
## Box-Ljung test
##
## data: r3$Mode1
## X-squared = 0.12172, df = 1, p-value = 0.7272
#Uji Asumsi Model Fit9 #Uji Normalitas
r9 <- residuals(fit9)
str(r9)
## List of 1
## $ Mode1: num [1:79] -3.76e+11 3.48e+11 2.10e+11 -9.87e+10 -4.30e+11 ...
n <- length(r9$Mode1)
mean <- mean(r9$Mode1)
sd <- sd(r9$Mode1)
res <- rnorm(n,mean,sd)
ks.test(r9$Mode1,res)
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: r9$Mode1 and res
## D = 0.20253, p-value = 0.07815
## alternative hypothesis: two-sided
#Uji Autokorelasi
Box.test(r9$Mode1,type=c("Ljung-Box"))
##
## Box-Ljung test
##
## data: r9$Mode1
## X-squared = 0.020934, df = 1, p-value = 0.885
#Uji Asumsi Model Fit10 #Uji Normalitas
r10 <- residuals(fit10)
str(r10)
## List of 1
## $ Mode1: num [1:79] -3.84e+11 3.76e+11 1.68e+11 -7.22e+10 -4.94e+11 ...
n <- length(r10$Mode1)
mean <- mean(r10$Mode1)
sd <- sd(r10$Mode1)
res <- rnorm(n,mean,sd)
ks.test(r10$Mode1,res)
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: r10$Mode1 and res
## D = 0.16456, p-value = 0.2361
## alternative hypothesis: two-sided
#Uji Autokorelasi
Box.test(r10$Mode1,type=c("Ljung-Box"))
##
## Box-Ljung test
##
## data: r10$Mode1
## X-squared = 0.12172, df = 1, p-value = 0.7272
#Uji Asumsi Model Fit12 #Uji Normalitas
r12 <- residuals(fit12)
str(r12)
## List of 1
## $ Mode1: num [1:79] -3.76e+11 3.48e+11 2.10e+11 -9.87e+10 -4.30e+11 ...
n <- length(r12$Mode1)
mean <- mean(r12$Mode1)
sd <- sd(r12$Mode1)
res <- rnorm(n,mean,sd)
ks.test(r12$Mode1,res)
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: r12$Mode1 and res
## D = 0.18987, p-value = 0.116
## alternative hypothesis: two-sided
#Uji Autokorelasi
Box.test(r12$Mode1,type=c("Ljung-Box"))
##
## Box-Ljung test
##
## data: r12$Mode1
## X-squared = 0.020934, df = 1, p-value = 0.885
#Uji Asumsi Model Fit13 #Uji Normalitas
r13 <- residuals(fit13)
str(r13)
## List of 1
## $ Mode1: num [1:79] -3.84e+11 3.76e+11 1.68e+11 -7.22e+10 -4.94e+11 ...
n <- length(r13$Mode1)
mean <- mean(r13$Mode1)
sd <- sd(r13$Mode1)
res <- rnorm(n,mean,sd)
ks.test(r13$Mode1,res)
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: r13$Mode1 and res
## D = 0.17722, p-value = 0.1677
## alternative hypothesis: two-sided
#Uji Autokorelasi
Box.test(r13$Mode1,type=c("Ljung-Box"))
##
## Box-Ljung test
##
## data: r13$Mode1
## X-squared = 0.12172, df = 1, p-value = 0.7272
#Nilai AIC
model <- list(fit2,fit3,fit9,fit10,fit12,fit13)
sapply(model, AIC)
## [1] 4162.500 4162.139 4162.500 4162.139 4162.500 4162.139
lamda <- 1.999924
lamda1 <- 1.44628
lamda2 <- 0.9522903
h <- 12
forecast_model <- predict(fit3, n.ahead = h)
str(forecast_model)
## List of 6
## $ :List of 7
## ..$ SDForecasts: num [1, 1:12] 2.63e+11 2.73e+11 2.88e+11 2.92e+11 2.93e+11 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr "79"
## .. .. ..$ : chr [1:12] "1" "2" "3" "4" ...
## ..$ Forecast : Named num [1:12] 2.84e+11 -2.28e+11 -5.13e+09 -1.52e+10 -1.39e+10 ...
## .. ..- attr(*, "names")= chr [1:12] "1" "2" "3" "4" ...
## ..$ exactVar : num [1:12] 6.91e+22 7.44e+22 8.28e+22 8.51e+22 8.56e+22 ...
## ..$ exactSD : num [1:12] 2.63e+11 2.73e+11 2.88e+11 2.92e+11 2.93e+11 ...
## ..$ limitVar : num [1:12] 6.91e+22 7.44e+22 8.28e+22 8.50e+22 8.55e+22 ...
## ..$ limitSD : num [1:12] 2.63e+11 2.73e+11 2.88e+11 2.92e+11 2.92e+11 ...
## ..$ sigma2 : num 6.91e+22
## $ z : num [1:79] -4.26e+11 5.41e+11 -5.90e+10 5.96e+10 -4.59e+11 ...
## $ limiting: logi TRUE
## $ predint : num 0.95
## $ name : chr "fit3"
## $ m : int 1
## - attr(*, "class")= chr "predarfima"
point_forecast <- as.vector(forecast_model[[1]]$Forecast)
se <- as.vector(forecast_model[[1]]$exactSD) # Menggunakan exactSD untuk Standard Error
forecast_lo80 <- point_forecast - 1.282 * se
forecast_hi80 <- point_forecast + 1.282 * se
forecast_lo95 <- point_forecast - 1.960 * se
forecast_hi95 <- point_forecast + 1.960 * se
inverse_boxcox <- function(x, lamda) {
if (abs(lamda) < 1e-6) {
return(exp(x))
} else {
result <- (lamda * x + 1)
# Penanganan nilai negatif/nol untuk menghindari error akar
result[result <= 0] <- 1e-10
return(result^(1 / lamda))
}
}
# Ramalan Titik
forecast_mean_final <- inverse_boxcox(inverse_boxcox(inverse_boxcox(point_forecast, lamda2), lamda1), lamda)
# Lo 80%
forecast_lo80_final <- inverse_boxcox(inverse_boxcox(inverse_boxcox(forecast_lo80, lamda2), lamda1), lamda)
# Hi 80%
forecast_hi80_final <- inverse_boxcox(inverse_boxcox(inverse_boxcox(forecast_hi80, lamda2), lamda1), lamda)
# Lo 95%
forecast_lo95_final <- inverse_boxcox(inverse_boxcox(inverse_boxcox(forecast_lo95, lamda2), lamda1), lamda)
# Hi 95%
forecast_hi95_final <- inverse_boxcox(inverse_boxcox(inverse_boxcox(forecast_hi95, lamda2), lamda1), lamda)
final_forecast_table <- data.frame(
"Point.Forecast" = forecast_mean_final,
"Lo.80" = forecast_lo80_final,
"Hi.80" = forecast_hi80_final,
"Lo.95" = forecast_lo95_final,
"Hi.95" = forecast_hi95_final
)
print(final_forecast_table)
## Point.Forecast Lo.80 Hi.80 Lo.95 Hi.95
## 1 22729.585824 1.732065 30186.38 1.732065 33080.18
## 2 1.732065 1.732065 16704.84 1.732065 23359.42
## 3 1.732065 1.732065 24852.60 1.732065 29045.14
## 4 1.732065 1.732065 24727.36 1.732065 29000.16
## 5 1.732065 1.732065 24786.09 1.732065 29055.08
## 6 1.732065 1.732065 24844.92 1.732065 29104.04
## 7 1.732065 1.732065 24892.91 1.732065 29142.80
## 8 1.732065 1.732065 24931.47 1.732065 29173.50
## 9 1.732065 1.732065 24962.80 1.732065 29198.25
## 10 1.732065 1.732065 24988.65 1.732065 29218.57
## 11 1.732065 1.732065 25010.32 1.732065 29235.54
## 12 1.732065 1.732065 25028.73 1.732065 29249.92
#Langkah 1: Memisahkan data
n_total <- length(datats)
n_test <- 8 # Asumsi data tes 8 periode (dari 80 total data, 80% train = 64, 20% test = 16. Jika n_test=8, train 72). Kita asumsikan n_total=80.
n_train <- n_total - n_test
train_data <- datats[1:n_train]
test_data <- datats[(n_train + 1):n_total]
# Langkah 2: Membuat model terbaik (ARFIMA(0,0,2)) dengan data pelatihan
# Menggunakan fracdiff untuk estimasi d, p=0, q=2
model_mape <- fracdiff(x = train_data, nar = 0, nma = 2)
## Warning: unable to compute correlation matrix; maybe change 'h'
# Langkah 3: Melakukan peramalan untuk periode pengujian (h = panjang test_data)
forecast_mape <- forecast(model_mape, h = length(test_data))
# Langkah 4: Menghitung metrik akurasi
# AKURASI DIHITUNG PADA SKALA YANG DITRANSFORMASI (datats)
accuracy_result <- accuracy(forecast_mape, test_data)
print(accuracy_result)
## ME RMSE MAE MPE MAPE
## Training set -1587423385 255243562277 181102907075 -376.32845 671.52715
## Test set 13441261688 325123925297 231240144800 96.81618 96.81618
## MASE ACF1
## Training set 0.508153 -0.0245641
## Test set 0.648832 NA