#### Package ####
library(MASS)
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
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
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
library(quadprog)
library(zoo)
## Warning: package 'zoo' was built under R version 4.3.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(fracdiff)
## Warning: package 'fracdiff' was built under R version 4.3.3
library(fUnitRoots)
## Warning: package 'fUnitRoots' was built under R version 4.3.3
library(lmtest)
#### Data ####
data=read.table(file.choose(),header=T,sep = ",")
head(data)
## Price.All.
## 1 51.76
## 2 49.13
## 3 43.45
## 4 48.20
## 5 51.75
## 6 55.40
DataAll<-ts(data)
DataInsample<-ts(data[1:216,])
DataOutsample<-ts(data[217:240,],start = 217,end = 240)
plot(DataAll)
plot(DataInsample)
plot(DataOutsample)
#### Buat Model Dari InSample ####
# ACF dan PACF dari data aktual
acf(DataInsample,24);grid()
pacf(DataInsample,24);grid()
# Cek Stasioner Dalam Variance
boxcox(DataInsample~1)
powerTransform(DataInsample)
## Estimated transformation parameter
## DataInsample
## 0.489028
lambda = powerTransform(DataInsample)$lambda
DataTransform = DataInsample^(lambda)
head(DataTransform)
## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 6.889556 6.716081 6.324452 6.653607 6.888905 7.122380
boxcox(DataTransform~1)
powerTransform(DataTransform)
## Estimated transformation parameter
## DataTransform
## 0.9999998
plot(DataTransform)
# Cek Stasioner Dalam Rata-rata
adfTest(DataTransform)
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -0.2529
## P VALUE:
## 0.535
##
## Description:
## Fri Nov 15 00:58:05 2024 by user: mrio1
DataDiff = diff(DataTransform,differences = 1)
head(DataDiff)
## Time Series:
## Start = 2
## End = 7
## Frequency = 1
## [1] -0.1734743 -0.3916295 0.3291552 0.2352976 0.2334749 -0.3669769
plot(DataDiff)
adfTest(DataDiff)
## Warning in adfTest(DataDiff): p-value smaller than printed p-value
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -9.3892
## P VALUE:
## 0.01
##
## Description:
## Fri Nov 15 00:58:05 2024 by user: mrio1
#ACF dan PACF dari data Baru
acf(DataDiff);grid() #Terdapat ACF signifikan pada lag 1
pacf(DataDiff);grid() #Terdapat PACF signifikan pada lag 1
#Penaksiran dan Pengujian Signifikansi Parameter
model1=arima(x=DataTransform,order=c(0,1,1));model1
##
## Call:
## arima(x = DataTransform, order = c(0, 1, 1))
##
## Coefficients:
## ma1
## 0.2329
## s.e. 0.0657
##
## sigma^2 estimated as 0.1519: log likelihood = -102.53, aic = 209.07
coeftest(model1)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.232901 0.065691 3.5454 0.000392 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model2=arima(x=DataTransform,order=c(1,1,0));model2
##
## Call:
## arima(x = DataTransform, order = c(1, 1, 0))
##
## Coefficients:
## ar1
## 0.2263
## s.e. 0.0665
##
## sigma^2 estimated as 0.1522: log likelihood = -102.7, aic = 209.41
coeftest(model2)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.226254 0.066539 3.4003 0.000673 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model3=arima(x=DataTransform,order=c(1,1,1));model3
##
## Call:
## arima(x = DataTransform, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.0526 0.1830
## s.e. 0.3018 0.2973
##
## sigma^2 estimated as 0.1519: log likelihood = -102.52, aic = 211.04
coeftest(model3)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.052567 0.301763 0.1742 0.8617
## ma1 0.183031 0.297314 0.6156 0.5381
#### Pemeriksaan Diagnostik ####
#Autokorelasi Residual (LJung-Box)
tsdiag(model1,length(DataDiff)) #ARIMA (0,1,1)
tsdiag(model2,length(DataDiff)) #ARIMA (1,1,0)
#Normalitas Residual
shapiro.test(resid(model1))
##
## Shapiro-Wilk normality test
##
## data: resid(model1)
## W = 0.9519, p-value = 1.256e-06
shapiro.test(resid(model2))
##
## Shapiro-Wilk normality test
##
## data: resid(model2)
## W = 0.95103, p-value = 1.024e-06
#Pemilihan model terbaik (Ukuran Kesalahan Prediksi In Sample)
summary(model1)
##
## Call:
## arima(x = DataTransform, order = c(0, 1, 1))
##
## Coefficients:
## ma1
## 0.2329
## s.e. 0.0657
##
## sigma^2 estimated as 0.1519: log likelihood = -102.53, aic = 209.07
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005542218 0.3888798 0.2945416 -0.04576632 3.884844 0.9699076
## ACF1
## Training set 0.003340876
summary(model2)
##
## Call:
## arima(x = DataTransform, order = c(1, 1, 0))
##
## Coefficients:
## ar1
## 0.2263
## s.e. 0.0665
##
## sigma^2 estimated as 0.1522: log likelihood = -102.7, aic = 209.41
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005119543 0.389188 0.2945661 -0.04290725 3.884154 0.9699883
## ACF1
## Training set 0.009967086
#Peramalan (Out Sample) 24 periode
#ARIMA (0,1,1)
set.seed(34)
RawFcast1 = seq(1,24)
a_t = rnorm(24,0,sd(resid(model1)))
for (i in 1:length(RawFcast1)){
theta1 = summary(model1)$coef
if (i==1) {
a_t1 = resid(model1)[216]
RawFcast1[i] = DataTransform[216] + a_t[i] - theta1*a_t1
}
else {
a_t1 = a_t[i-1]
RawFcast1[i] = RawFcast1[i-1] + a_t[i] - theta1*a_t1
}
}
InversedFcast1 = ts(RawFcast1^(1/lambda),start = 217,end = 240)
Error1 = DataOutsample - InversedFcast1
ME1 = mean(Error1)
RMSE1 = sqrt(mean(Error1^2))
MAE1 = mean(abs(Error1))
MPE1 = mean(Error1/DataOutsample)*100
MAPE1 = mean(abs(Error1/DataOutsample))*100
MASE1 = mean(abs(Error1^2))
kesalahan = data.frame(ME1,RMSE1,MAE1,MPE1,MAPE1,MASE1)
print(kesalahan)
## ME1 RMSE1 MAE1 MPE1 MAPE1 MASE1
## 1 -4.745695 10.05708 7.444316 -6.444256 9.684252 101.1449
#ARIMA (1,1,0)
set.seed(407)
RawFcast2 = seq(1,24)
a_t = rnorm(24,0,sd(resid(model2)))
for (i in 1:length(RawFcast2)){
phi1 = summary(model2)$coef
if (i==1) {
RawFcast2[i] = phi1*DataTransform[216] + DataTransform[216] - phi1*DataTransform[215] + a_t[i]
}
else if (i==2){
RawFcast2[i] = phi1*RawFcast1[i-1] + RawFcast1[i-1] - phi1*DataTransform[216] + a_t[i]
}
else {
RawFcast2[i] = phi1*RawFcast1[i-1] + RawFcast1[i-1] - phi1*RawFcast1[i-2] + a_t[i]
}
}
InversedFcast2 = ts(RawFcast2^(1/lambda),start = 217,end = 240)
Error2 = DataOutsample - InversedFcast2
ME2 = mean(Error2)
RMSE2 = sqrt(mean(Error2^2))
MAE2 = mean(abs(Error2))
MPE2 = mean(Error2/DataOutsample)*100
MAPE2 = mean(abs(Error2/DataOutsample))*100
MASE2 = mean(abs(Error2^2))
kesalahan2 = data.frame(ME2,RMSE2,MAE2,MPE2,MAPE2,MASE2)
print(kesalahan2)
## ME2 RMSE2 MAE2 MPE2 MAPE2 MASE2
## 1 -2.573292 8.829646 7.099529 -3.436986 9.262223 77.96264
### Plot Predict and Forecast
InversedPred1=fitted(model1)^(1/lambda)
InversedPred2=fitted(model2)^(1/lambda)
warna = c("black","red","green","black") #Warna Data Aktual, Prediksi, Peramalan, Garis
#Model ARIMA (0,1,1)
plot(DataAll,xlim=c(-1,250),main="ARIMA(0,1,1)",xlab= "Periode",ylab="Harga Minyak Mentah Dunia",pch=16,col=warna[1])
points(DataAll,col=warna[4],lwd=1,pch=16,cex=1.2)
lines(InversedPred1,col=warna[4],lwd=1)
points(InversedPred1,cex=1,col=warna[2],pch=16)
points(InversedFcast1,cex=1,col=warna[3],pch=19)
lines(InversedFcast1,col=warna[4],lwd=1)
legend("topright",legend=c("Aktual","Prediksi 1","Peramalan 1"),cex=1.25,bty="n",col=warna[1:3],pch=c(16,16,19),text.font = 2);grid()
#Model ARIMA (1,1,0)
plot(DataAll,xlim=c(-1,250),main="ARIMA(1,1,0)",xlab= "Periode",ylab="Harga Minyak Mentah Dunia",pch=16,col=warna[1])
points(DataAll,col=warna[4],lwd=1,pch=16,cex=1.2)
lines(InversedPred2,col=warna[4],lwd=0.5)
points(InversedPred2,cex=1,col=warna[2],pch=16)
points(InversedFcast2,cex=1,col=warna[3],pch=19)
lines(InversedFcast2,col=warna[4],lwd=0.5)
legend("topright",legend=c("Aktual","Prediksi 2","Peramalan 2"),cex=1.25,bty="n",col=warna[1:3],pch=c(16,16,19),text.font = 2);grid()
