Paket

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)

Data

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")

Transformasi BoxCox

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)

Statistik Hurst

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

Estimasi d dengan beberapa metode

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

cek semua nilai d

Jika d hasil GPH > 1, lakukan differencing manual

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

Definisikan nilai lambda yang sudah Anda hitung

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"

1. Ekstraksi Data dari Objek Bersarang

Mengakses elemen list pertama ([[1]]) dan kemudian sub-elemen yang benar

point_forecast <- as.vector(forecast_model[[1]]$Forecast)
se <- as.vector(forecast_model[[1]]$exactSD) # Menggunakan exactSD untuk Standard Error

2. Hitung Batas Interval Kepercayaan (CI)

(Z-score 80% = 1.282, Z-score 95% = 1.960)

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

3. Fungsi Transformasi Balik Box-Cox

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))
  }
}

4. Terapkan Transformasi Balik secara Berurutan

# 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)

5. Gabungkan Hasil ke dalam Tabel Final

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

PLOT BELUM

#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