####-----------------PACKAGE-----------------####
library(MASS)
library(hms)
library(car)
## Loading required package: carData
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(tseries)
library(timeSeries)
## Loading required package: timeDate
##
## Attaching package: 'timeSeries'
## The following objects are masked from 'package:graphics':
##
## lines, points
library(quadprog)
library(zoo)
##
## Attaching package: 'zoo'
## The following object is masked from 'package:timeSeries':
##
## time<-
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(fracdiff)
library(fUnitRoots)
library(lmtest)
library(Metrics)
##
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
##
## accuracy
WorldCup<-read.csv("Data GoogleTrends (World Cup).csv", sep = ";")
WorldCup$Month <- as.Date(paste0(WorldCup$Month, "-1"))
head(WorldCup);tail(WorldCup)
## Month World.Cup D1 D2 D3 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 M12
## 1 2005-05-01 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## 2 2005-06-01 2 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## 3 2005-07-01 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 4 2005-08-01 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## 5 2005-09-01 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## 6 2005-10-01 2 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## Month World.Cup D1 D2 D3 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 M12
## 235 2024-11-01 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## 236 2024-12-01 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## 237 2025-01-01 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## 238 2025-02-01 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## 239 2025-03-01 2 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## 240 2025-04-01 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
Zt = as.numeric(WorldCup$World.Cup)
Zt_zoo <- zoo(Zt, order.by = WorldCup$Month)
{plot(Zt_zoo, lwd = 5, col = "black",ylim = c(0,110),
main="Plot Pencarian World Cup",ylab = "Keyword WorldCup", xlab = "Tahun")
text(x = WorldCup$Month[14],
y = WorldCup$World.Cup[14]+5,
labels = format(as.Date(WorldCup$Month[14]), "%Y-%b"),
col = "red")
text(x = WorldCup$Month[62],
y = WorldCup$World.Cup[62]+5,
labels = format(as.Date(WorldCup$Month[62]), "%Y-%b"),
col = "red")
text(x = WorldCup$Month[110],
y = WorldCup$World.Cup[110]+5,
labels = format(as.Date(WorldCup$Month[110]), "%Y-%b"),
col = "red")
text(x = WorldCup$Month[158],
y = WorldCup$World.Cup[158]+5,
labels = format(as.Date(WorldCup$Month[158]), "%Y-%b"),
col = "red")
text(x = WorldCup$Month[212],
y = WorldCup$World.Cup[212]+5,
labels = format(as.Date(WorldCup$Month[212]), "%Y-%b"),
col = "red")
}
#In Sample 85%(204) & Out Sample 15%(36)
Insample =WorldCup$World.Cup[1:204]
Outsample=WorldCup$World.Cup[205:240]
t=c(1:204)
Insample_zoo = zoo(Insample, order.by = WorldCup$Month[1:204])
plot(Insample_zoo, lwd = 5, col = "green", ylab = "Insample", xlab = "Tahun")
####-----------------Naive Trend Linear & Exponential-----------------####
### Naive Trend Linear
PredNaiveLinear = c(NA,NA,rep(0,202))
for (i in 2:203){
PredNaiveLinear[i+1] = (Insample[i] + (Insample[i]- Insample[i-1]))
#if (PredNaiveLinear[i+1] < 0) {PredNaiveLinear[i+1] = 0}
}
ForeNaiveLinear = rep(0,36)
for (i in 1:36){
if (i == 1){
ForeNaiveLinear[i] = (Insample[length(Insample)] + (Insample[length(Insample)] - Insample[length(Insample)-1]))
}
else if (i == 2){
ForeNaiveLinear[i] = (ForeNaiveLinear[1] + (ForeNaiveLinear[1] - Insample[length(Insample)]))
}
else {
ForeNaiveLinear[i] = (ForeNaiveLinear[i-1]+ (ForeNaiveLinear[i-1] - ForeNaiveLinear[i-2]))
}
}
NaiveLinear_zoo = zoo(c(PredNaiveLinear,ForeNaiveLinear), order.by = WorldCup$Month)
{plot(Zt_zoo, lwd = 5, col = "black",lty=1,
main = "Metode Naive Linear",ylab = "Keyword World Cup", xlab = "Tahun")
lines(NaiveLinear_zoo[1:204],lty=2,lwd = 4, col = "yellow2")
lines(NaiveLinear_zoo[205:240],lty=2,lwd = 4, col = "green")
abline(v = WorldCup$Month[204] , col = "red", lty = 3, lwd = 2)
text(x = WorldCup$Month[204], y = max(WorldCup$World.Cup), labels = "Batas Insample",srt = 90, pos = 2, col = "red", cex = 1.2)
text(x = WorldCup$Month[204],y = par("usr")[3]+1.5, pos = 2,labels = "2022 Apr",col = "red",cex = 1.2)
legend("topleft" , lwd = 3, bty = "n",
legend = c("Data Aktual (Insample)", "Pred Insample", "Pred Outsample"),
col = c("black","yellow2","green"), lty = c(1,2))
}
### Naive Trend Exponential
PredNaiveExp = c(NA,NA,rep(0,202))
for (i in 2:203){
PredNaiveExp[i+1] = abs(Insample[i] * (Insample[i]/Insample[i-1]))
#if (PredNaiveExp[i+1] < 0) {PredNaiveExp[i+1] = 0}
}
ForeNaiveExp = rep(0,36)
for (i in 1:36){
if (i == 1){
ForeNaiveExp[i] = abs(Insample[length(Insample)] * (Insample[length(Insample)] / Insample[length(Insample)-1]))
}
else if (i == 2){
ForeNaiveExp[i] = abs(ForeNaiveExp[1] * (ForeNaiveExp[1]/Insample[length(Insample)]))
}
else {
ForeNaiveExp[i] = abs(ForeNaiveExp[i-1] * (ForeNaiveExp[i-1]/ForeNaiveExp[i-2]))
}
}
NaiveExp_zoo = zoo(c(PredNaiveExp,ForeNaiveExp), order.by = WorldCup$Month)
{plot(Zt_zoo, lwd = 5, col = "black",lty=1,
main = "Metode Naive Exponential",ylab = "Keyword World Cup", xlab = "Tahun")
lines(NaiveExp_zoo[1:204],lty=2,lwd = 4, col = "yellow2")
lines(NaiveExp_zoo[205:240],lty=2,lwd = 4, col = "green")
abline(v = WorldCup$Month[204] , col = "red", lty = 3, lwd = 2)
text(x = WorldCup$Month[204], y = max(WorldCup$World.Cup), labels = "Batas Insample",srt = 90, pos = 2, col = "red", cex = 1.2)
text(x = WorldCup$Month[204],y = par("usr")[3]+1.5, pos = 2,labels = "2022 Apr",col = "red",cex = 1.2)
legend("topleft" , lwd = 3,bty = "n",
legend = c("Data Aktual (Insample)", "Pred Insample", "Pred Outsample"),
col = c("black","yellow2","green"), lty = c(1,2))
}
####-----------------Pemodelan SARIMA-----------------####
#Pemodelan SARIMA dari WorldCup Residual TSR
Insample_zoo = zoo(Insample, order.by = WorldCup$Month[1:204])
plot(Insample_zoo, lwd = 5, col = "green", ylab = "Insample", xlab = "Tahun")
acf(Insample,150,col="blue",lwd=1);grid() #Membuat ACF
###-----------Periodogram-----------###
s <- spectrum(Insample, log = "no") #fungsi spectrum() untuk mendapatkan periodogram
points(x = s$freq[which.max(s$spec)],
y = max(s$spec),
pch = 16, col = "red", cex = 1.5)
text(x = s$freq[which.max(s$spec)],
y = max(s$spec)-30,
labels = "Freq = 0.04166667",
col = "red")
frekuensi_dominan <- s$freq[which.max(s$spec)] #Ambil frekuensi dengan spektrum tertinggi
periode <- 1 / frekuensi_dominan #Hitung periode
# Tampilkan hasil
cat("Frekuensi dominan:", frekuensi_dominan, "\n")
## Frekuensi dominan: 0.04166667
cat("Periode yang terdeteksi:", periode, "\n")
## Periode yang terdeteksi: 24
#Stasioneritas dalam Variansi
boxcox(Insample~1) #belum stasioner dalam varians
p<-powerTransform(Insample);p
## Estimated transformation parameter
## Insample
## -2.182963
lambda = p$lambda
Yt = ts(Insample^lambda)
plot(zoo(Yt, order.by = WorldCup$Month[1:204]), lwd = 5, col = "green",main = "Data Transformasi",ylab = "Insample", xlab = "Tahun")
boxcox(Yt~1);powerTransform(Yt) #sudah stasioner dalam varians setelah pangkat lambda
## Estimated transformation parameter
## Yt
## 1.000001
#Stasioneritas dalam Rata-Rata ADF
library(urca)
##
## Attaching package: 'urca'
## The following objects are masked from 'package:fUnitRoots':
##
## punitroot, qunitroot, unitrootTable
adfTest(Yt)#belum stasioner dalam rata-rata
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -2.1571
## P VALUE:
## 0.03171
##
## Description:
## Fri Apr 25 09:27:27 2025 by user: mrio1
datadiff=diff(Yt) #NON MUSIMAN
adfTest(datadiff)#sudah stasioner dalam rata-rata setelah differencing 1 kali (d=1)
## 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: -12.9118
## P VALUE:
## 0.01
##
## Description:
## Fri Apr 25 09:27:27 2025 by user: mrio1
plot(zoo(datadiff, order.by = WorldCup$Month[1:204]), lwd = 3, col = "black",main = "Data Differencing",ylab = "Insample", xlab = "Tahun")
acf(datadiff,50,col="blue",main="Data Diff Non-Musiman",cex=2,lwd=2);grid() #Membuat ACF (cut-off setelah lag 1 -> q=1)
pacf(datadiff,50,col="red",main="Data Diff Non-Musiman",cex=2,lwd=2);grid() #Membuat PACF (cut-off setelah lag 1 -> p=1))
datadiffSeasonal=diff(datadiff,lag=periode) #MUSIMAN
adfTest(datadiffSeasonal)#Seasonal sudah stasioner dalam rata-rata setelah differencing 1 kali (D=1)
## Warning in adfTest(datadiffSeasonal): p-value smaller than printed p-value
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -11.5043
## P VALUE:
## 0.01
##
## Description:
## Fri Apr 25 09:27:27 2025 by user: mrio1
acf(datadiffSeasonal,50,col="blue",main="Data Diff Musiman",cex=2,lwd=2);grid() #Membuat ACF (cut-off setelah lag 1 -> Q=1)
pacf(datadiffSeasonal,50,col="red",main="Data Diff Musiman",cex=2,lwd=2);grid() #Membuat PACF (cut-off setelah lag 1 -> P=1))
#Model Sementara
fit1=arima(Yt,order=c(0,1,1),seasonal=list(order=c(0,1,0),period=24))
fit2=Arima(Yt,order=c(1,1,0),seasonal=list(order=c(0,1,0),period=24))
fit3=Arima(Yt,order=c(1,1,1),seasonal=list(order=c(0,1,0),period=24))
fit4=Arima(Yt,order=c(0,1,1),seasonal=list(order=c(1,1,0),period=24))
fit5=Arima(Yt,order=c(1,1,0),seasonal=list(order=c(1,1,0),period=24))
#fit6=Arima(Yt,order=c(1,1,1),seasonal=list(order=c(1,1,0),period=24))
fit7=Arima(Yt,order=c(0,1,1),seasonal=list(order=c(1,1,1),period=24))
fit8=Arima(Yt,order=c(1,1,0),seasonal=list(order=c(1,1,1),period=24))
fit9=Arima(Yt,order=c(1,1,1),seasonal=list(order=c(1,1,1),period=24))
coeftest(fit1)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.545240 0.076691 -7.1096 1.164e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit2)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.377384 0.069551 -5.426 5.762e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit3)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.225356 0.125397 1.7971 0.07231 .
## ma1 -0.708065 0.091004 -7.7806 7.216e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit4)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.77876 0.06480 -12.018 < 2.2e-16 ***
## sar1 -0.79366 0.04396 -18.054 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit5)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.445731 0.068178 -6.5377 6.246e-11 ***
## sar1 -0.730791 0.046441 -15.7358 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coeftest(fit6)
coeftest(fit7)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.748047 0.074173 -10.0852 < 2.2e-16 ***
## sar1 -0.686923 0.124395 -5.5221 3.349e-08 ***
## sma1 -0.270112 0.263324 -1.0258 0.305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit8)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.446934 0.067882 -6.5840 4.580e-11 ***
## sar1 -0.530424 0.091166 -5.8182 5.947e-09 ***
## sma1 -0.473182 0.140040 -3.3789 0.0007277 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit9)
## Warning in sqrt(diag(se)): NaNs produced
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.055773 NaN NaN NaN
## ma1 -0.750309 NaN NaN NaN
## sar1 -0.656952 NaN NaN NaN
## sma1 -0.312631 NaN NaN NaN
#Pemeriksaan Diagnostik
resid1=resid(fit1)
resid2=resid(fit2)
resid4=resid(fit4)
resid5=resid(fit5)
resid8=resid(fit8)
#WN
x<-resid1
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box Resid1",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
abline(v = seq(0, 200, by = 24), col = "blue", lty = 2)
Box.test(x,lag=12,type="Ljung")$p.value
## [1] 0.03126445
Box.test(x,lag=24,type="Ljung")$p.value
## [1] 0
Box.test(x,lag=36,type="Ljung")$p.value
## [1] 0
Box.test(x,lag=48,type="Ljung")$p.value
## [1] 0
Box.test(x,lag=60,type="Ljung")$p.value
## [1] 0
x<-resid2
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box Resid2",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
abline(v = seq(0, 200, by = 24), col = "blue", lty = 2)
Box.test(x,lag=12,type="Ljung")$p.value
## [1] 0.00132457
Box.test(x,lag=24,type="Ljung")$p.value
## [1] 0
Box.test(x,lag=36,type="Ljung")$p.value
## [1] 0
Box.test(x,lag=48,type="Ljung")$p.value
## [1] 0
Box.test(x,lag=60,type="Ljung")$p.value
## [1] 0
x<-resid4
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box Resid7",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
abline(v = seq(0, 200, by = 24), col = "blue", lty = 2) #Lolos
Box.test(x,lag=12,type="Ljung")$p.value
## [1] 0.8883301
Box.test(x,lag=24,type="Ljung")$p.value
## [1] 0.96117
Box.test(x,lag=36,type="Ljung")$p.value
## [1] 0.8284312
Box.test(x,lag=48,type="Ljung")$p.value
## [1] 0.6887752
Box.test(x,lag=60,type="Ljung")$p.value
## [1] 0.7302078
x<-resid5
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box Resid15",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
abline(v = seq(0, 200, by = 24), col = "blue", lty = 2)
Box.test(x,lag=12,type="Ljung")$p.value
## [1] 0.03745838
Box.test(x,lag=24,type="Ljung")$p.value
## [1] 0.09247981
Box.test(x,lag=36,type="Ljung")$p.value
## [1] 0.02735657
Box.test(x,lag=48,type="Ljung")$p.value
## [1] 0.0102357
Box.test(x,lag=60,type="Ljung")$p.value
## [1] 0.02223235
x<-resid8
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box Resid19",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
abline(v = seq(0, 200, by = 24), col = "blue", lty = 2) #Lolos
Box.test(x,lag=12,type="Ljung")$p.value
## [1] 0.09512407
Box.test(x,lag=24,type="Ljung")$p.value
## [1] 0.267634
Box.test(x,lag=36,type="Ljung")$p.value
## [1] 0.1946229
Box.test(x,lag=48,type="Ljung")$p.value
## [1] 0.1642025
Box.test(x,lag=60,type="Ljung")$p.value
## [1] 0.3713284
#Cek Normalitas (4)
ks.test(resid1,"pnorm",mean(resid1),sqrt(var(resid1)))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid1
## D = 0.13556, p-value = 0.001109
## alternative hypothesis: two-sided
ks.test(resid2,"pnorm",mean(resid2),sqrt(var(resid2)))
## Warning in ks.test.default(resid2, "pnorm", mean(resid2), sqrt(var(resid2))):
## ties should not be present for the one-sample Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid2
## D = 0.21893, p-value = 6.425e-09
## alternative hypothesis: two-sided
ks.test(resid4,"pnorm",mean(resid4),sqrt(var(resid4))) #Normal dengan alpha=0.01
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid4
## D = 0.11148, p-value = 0.01256
## alternative hypothesis: two-sided
ks.test(resid5,"pnorm",mean(resid5),sqrt(var(resid5)))
## Warning in ks.test.default(resid5, "pnorm", mean(resid5), sqrt(var(resid5))):
## ties should not be present for the one-sample Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid5
## D = 0.18229, p-value = 2.59e-06
## alternative hypothesis: two-sided
ks.test(resid8,"pnorm",mean(resid8),sqrt(var(resid8)))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid8
## D = 0.15373, p-value = 0.0001299
## alternative hypothesis: two-sided
#Predict & Forecast SARIMA
PredSARIMA = abs(fitted(fit4))^(1/lambda)
ForeSARIMA = abs(as.numeric(forecast(fit4, h = 36)$mean))^(1/lambda)
SARIMA_zoo = zoo(c(PredSARIMA,ForeSARIMA), order.by = WorldCup$Month[1:240])
{plot(Zt_zoo, lwd = 5, col = "black",lty=1,
main = "Metode SARIMA",ylab = "Keyword World Cup", xlab = "Tahun")
lines(SARIMA_zoo[1:204],lty=2,lwd = 4, col = "yellow2")
lines(SARIMA_zoo[205:240],lty=2,lwd = 4, col = "green")
abline(v = WorldCup$Month[204] , col = "red", lty = 3, lwd = 2)
text(x = WorldCup$Month[204], y = max(WorldCup$World.Cup), labels = "Batas Insample",srt = 90, pos = 2, col = "red", cex = 1.2)
text(x = WorldCup$Month[204],y = par("usr")[3]+1.5, pos = 2,labels = "2022 Apr",col = "red",cex = 1.2)
legend("topleft" , lwd = 3,
legend = c("Data Aktual (Insample)", "Pred Insample", "Pred Outsample"),
col = c("black","yellow2","green"), lty = c(1,2))
}
####-----------------Time Series Regression-----------------####
D1=WorldCup$D1[1:204]
D2=WorldCup$D2[1:204]
D3=WorldCup$D3[1:204]
M1=WorldCup$M1[1:204]
M2=WorldCup$M2[1:204]
M3=WorldCup$M3[1:204]
M4=WorldCup$M4[1:204]
M5=WorldCup$M5[1:204]
M6=WorldCup$M6[1:204]
M7=WorldCup$M7[1:204]
M8=WorldCup$M8[1:204]
M9=WorldCup$M9[1:204]
M10=WorldCup$M10[1:204]
M11=WorldCup$M11[1:204]
M12=WorldCup$M12[1:204]
model1=lm(Insample~t+D1+D2+D3+M1+M2+M3+M4+M5+M6+M7+M8+M9+M10+M11+M12)
step(model1)
## Start: AIC=610.84
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M3 + M4 + M5 + M6 + M7 +
## M8 + M9 + M10 + M11 + M12
##
##
## Step: AIC=610.84
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M3 + M4 + M5 + M6 + M7 +
## M8 + M9 + M10 + M11
##
## Df Sum of Sq RSS AIC
## - M3 1 0.0 3482.8 608.84
## - M10 1 0.0 3482.8 608.85
## - M11 1 0.2 3483.0 608.86
## - M5 1 0.3 3483.1 608.86
## - M8 1 0.8 3483.6 608.89
## - M1 1 0.8 3483.6 608.89
## - M9 1 0.8 3483.6 608.89
## - M2 1 0.8 3483.6 608.89
## - D3 1 3.0 3485.8 609.02
## - M4 1 5.3 3488.1 609.15
## <none> 3482.8 610.84
## - t 1 39.0 3521.8 611.12
## - D1 1 102.3 3585.1 614.75
## - M7 1 166.1 3648.9 618.35
## - M6 1 184.3 3667.1 619.36
## - D2 1 11767.7 15250.5 910.11
##
## Step: AIC=608.84
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M4 + M5 + M6 + M7 + M8 +
## M9 + M10 + M11
##
## Df Sum of Sq RSS AIC
## - M10 1 0.0 3482.8 606.85
## - M11 1 0.3 3483.1 606.86
## - M5 1 0.5 3483.3 606.87
## - M8 1 0.9 3483.7 606.90
## - M1 1 1.0 3483.8 606.90
## - M9 1 1.0 3483.8 606.90
## - M2 1 1.0 3483.8 606.90
## - D3 1 3.0 3485.8 607.02
## - M4 1 6.7 3489.5 607.24
## <none> 3482.8 608.84
## - t 1 39.0 3521.8 609.12
## - D1 1 102.3 3585.1 612.75
## - M7 1 215.2 3698.1 619.08
## - M6 1 241.0 3723.8 620.50
## - D2 1 11767.7 15250.5 908.11
##
## Step: AIC=606.85
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M4 + M5 + M6 + M7 + M8 +
## M9 + M11
##
## Df Sum of Sq RSS AIC
## - M11 1 0.3 3483.1 604.86
## - M5 1 0.6 3483.4 604.88
## - M8 1 1.0 3483.8 604.90
## - M1 1 1.0 3483.8 604.90
## - M9 1 1.0 3483.8 604.90
## - M2 1 1.0 3483.9 604.91
## - D3 1 3.0 3485.8 605.02
## - M4 1 7.2 3490.0 605.27
## <none> 3482.8 606.85
## - t 1 39.0 3521.9 607.12
## - D1 1 102.3 3585.1 610.75
## - M7 1 238.2 3721.1 618.34
## - M6 1 269.3 3752.1 620.04
## - D2 1 11767.7 15250.5 906.11
##
## Step: AIC=604.86
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M4 + M5 + M6 + M7 + M8 +
## M9
##
## Df Sum of Sq RSS AIC
## - M8 1 0.8 3483.9 602.91
## - M1 1 0.8 3483.9 602.91
## - M9 1 0.8 3483.9 602.91
## - M5 1 0.8 3483.9 602.91
## - M2 1 0.9 3483.9 602.91
## - D3 1 3.0 3486.1 603.04
## - M4 1 7.0 3490.1 603.27
## <none> 3483.1 604.86
## - t 1 39.1 3522.2 605.14
## - D1 1 102.3 3585.4 608.77
## - M7 1 248.3 3731.4 616.91
## - M6 1 289.6 3772.7 619.15
## - D2 1 11767.7 15250.8 904.11
##
## Step: AIC=602.91
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M4 + M5 + M6 + M7 + M9
##
## Df Sum of Sq RSS AIC
## - M9 1 0.5 3484.3 600.93
## - M1 1 0.6 3484.4 600.94
## - M2 1 0.6 3484.5 600.94
## - M5 1 1.2 3485.1 600.98
## - D3 1 2.4 3486.2 601.04
## - M4 1 6.4 3490.2 601.28
## <none> 3483.9 602.91
## - t 1 39.1 3523.0 603.18
## - D1 1 102.3 3586.2 606.81
## - M7 1 251.4 3735.3 615.12
## - M6 1 304.6 3788.5 618.01
## - D2 1 11767.8 15251.6 902.12
##
## Step: AIC=600.93
## Insample ~ t + D1 + D2 + D3 + M1 + M2 + M4 + M5 + M6 + M7
##
## Df Sum of Sq RSS AIC
## - M1 1 0.4 3484.8 598.96
## - M2 1 0.5 3484.8 598.96
## - M5 1 1.4 3485.8 599.02
## - D3 1 2.0 3486.3 599.05
## - M4 1 6.0 3490.4 599.29
## <none> 3484.3 600.93
## - t 1 39.1 3523.4 601.21
## - D1 1 102.3 3586.6 604.84
## - M7 1 253.1 3737.4 613.24
## - M6 1 313.5 3797.9 616.51
## - D2 1 11767.7 15252.1 900.13
##
## Step: AIC=598.96
## Insample ~ t + D1 + D2 + D3 + M2 + M4 + M5 + M6 + M7
##
## Df Sum of Sq RSS AIC
## - M2 1 0.3 3485.1 596.98
## - M5 1 1.7 3486.5 597.06
## - D3 1 2.2 3487.0 597.09
## - M4 1 5.7 3490.4 597.29
## <none> 3484.8 598.96
## - t 1 39.0 3523.8 599.23
## - D1 1 102.3 3587.1 602.86
## - M7 1 255.4 3740.2 611.39
## - M6 1 323.6 3808.3 615.07
## - D2 1 11767.7 15252.5 898.13
##
## Step: AIC=596.98
## Insample ~ t + D1 + D2 + D3 + M4 + M5 + M6 + M7
##
## Df Sum of Sq RSS AIC
## - M5 1 1.9 3487.0 595.09
## - D3 1 2.4 3487.5 595.12
## - M4 1 5.4 3490.5 595.30
## <none> 3485.1 596.98
## - t 1 39.0 3524.1 597.25
## - D1 1 102.3 3587.4 600.88
## - M7 1 257.1 3742.2 609.50
## - M6 1 331.4 3816.5 613.51
## - D2 1 11767.7 15252.8 896.14
##
## Step: AIC=595.09
## Insample ~ t + D1 + D2 + D3 + M4 + M6 + M7
##
## Df Sum of Sq RSS AIC
## - D3 1 2.1 3489.1 593.21
## - M4 1 7.3 3494.4 593.52
## <none> 3487.0 595.09
## - t 1 38.8 3525.9 595.35
## - D1 1 127.3 3614.3 600.41
## - M7 1 264.5 3751.5 608.01
## - M6 1 329.7 3816.7 611.52
## - D2 1 11767.5 15254.6 894.16
##
## Step: AIC=593.21
## Insample ~ t + D1 + D2 + M4 + M6 + M7
##
## Df Sum of Sq RSS AIC
## - M4 1 7.8 3496.9 591.67
## <none> 3489.1 593.21
## - t 1 38.0 3527.1 593.42
## - D1 1 126.4 3615.5 598.47
## - M7 1 269.2 3758.3 606.37
## - M6 1 327.6 3816.7 609.52
## - D2 1 11766.9 15256.0 892.18
##
## Step: AIC=591.67
## Insample ~ t + D1 + D2 + M6 + M7
##
## Df Sum of Sq RSS AIC
## <none> 3496.9 591.67
## - t 1 36.5 3533.5 591.79
## - D1 1 118.9 3615.8 596.49
## - M7 1 263.6 3760.6 604.50
## - M6 1 337.7 3834.6 608.47
## - D2 1 11765.6 15262.6 890.27
##
## Call:
## lm(formula = Insample ~ t + D1 + D2 + M6 + M7)
##
## Coefficients:
## (Intercept) t D1 D2 M6 M7
## 0.433533 0.007212 3.957088 43.901628 5.000525 -4.418451
#Seleksi Variabel
model2=lm(formula = Insample ~ t + D1 + D2 + M6 + M7)
summary(model2)
##
## Call:
## lm(formula = Insample ~ t + D1 + D2 + M6 + M7)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.025 -0.604 -0.043 0.424 32.525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.433533 0.618229 0.701 0.483970
## t 0.007212 0.005015 1.438 0.151987
## D1 3.957088 1.525278 2.594 0.010186 *
## D2 43.901628 1.700922 25.810 < 2e-16 ***
## M6 5.000525 1.143640 4.372 1.98e-05 ***
## M7 -4.418451 1.143596 -3.864 0.000151 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.203 on 198 degrees of freedom
## Multiple R-squared: 0.8178, Adjusted R-squared: 0.8132
## F-statistic: 177.7 on 5 and 198 DF, p-value: < 2.2e-16
res=resid(model2) #Residual dari model TSR
#Predict And Forecast
tnext = c(205:240)
PredTSR = abs(fitted(model2))
ForeTSR = abs(0.433533+0.007212*tnext+3.957088*WorldCup$D1[205:240]+43.901628*WorldCup$D2[205:240]+5.000525*WorldCup$M6[205:240]-4.418451*WorldCup$M7[205:240])
TSR_zoo = zoo(c(PredTSR,ForeTSR), order.by = WorldCup$Month[1:240])
{plot(Zt_zoo, lwd = 5, col = "black",lty=1,
main = "Metode TSR",ylab = "Keyword World Cup", xlab = "Tahun")
lines(TSR_zoo[1:204],lty=2,lwd = 4, col = "yellow2")
lines(TSR_zoo[205:240],lty=2,lwd = 4, col = "green")
abline(v = WorldCup$Month[204] , col = "red", lty = 3, lwd = 2)
text(x = WorldCup$Month[204], y = max(WorldCup$World.Cup), labels = "Batas Insample",srt = 90, pos = 2, col = "red", cex = 1.2)
text(x = WorldCup$Month[204],y = par("usr")[3]+1.5, pos = 2,labels = "2022 Apr",col = "red",cex = 1.2)
legend("topleft" , lwd = 3,
legend = c("Data Aktual (Insample)", "Pred Insample", "Pred Outsample"),
col = c("black","yellow2","green"), lty = c(1,2))
}
####-----------------Perhitungan Metrik-----------------####
###INSAMPLE
Metrik_Insample = data.frame(RMSE =c(rep(0,4)),
MAPE =c(rep(0,4)),
SMAPE=c(rep(0,4)))
row.names(Metrik_Insample) = c("Naive Linear",
"Naive Exponential",
"SARIMA",
"TSR")
Metrik_Insample$RMSE[1] = rmse(Insample[3:204],PredNaiveLinear[3:204])
Metrik_Insample$RMSE[2] = rmse(Insample[3:204],PredNaiveExp[3:204])
Metrik_Insample$RMSE[3] = rmse(Insample[1:204],PredSARIMA)
Metrik_Insample$RMSE[4] = rmse(Insample[1:204],PredTSR)
Metrik_Insample$MAPE[1] = mape(Insample[3:204],PredNaiveLinear[3:204])
Metrik_Insample$MAPE[2] = mape(Insample[3:204],PredNaiveExp[3:204])
Metrik_Insample$MAPE[3] = mape(Insample[1:204],PredSARIMA)
Metrik_Insample$MAPE[4] = mape(Insample[1:204],PredTSR)
Metrik_Insample$SMAPE[1] = smape(Insample[3:204],PredNaiveLinear[3:204])
Metrik_Insample$SMAPE[2] = smape(Insample[3:204],PredNaiveExp[3:204])
Metrik_Insample$SMAPE[3] = smape(Insample[1:204],PredSARIMA)
Metrik_Insample$SMAPE[4] = smape(Insample[1:204],PredTSR)
###OUTSAMPLE
Metrik_Outsample = data.frame(RMSE =c(rep(0,4)),
MAPE =c(rep(0,4)),
SMAPE=c(rep(0,4)))
row.names(Metrik_Outsample) = c("Naive Linear",
"Naive Exponential",
"SARIMA",
"TSR")
Metrik_Outsample$RMSE[1] = rmse(Outsample,ForeNaiveLinear)
Metrik_Outsample$RMSE[2] = rmse(Outsample,ForeNaiveExp)
Metrik_Outsample$RMSE[3] = rmse(Outsample,ForeSARIMA)
Metrik_Outsample$RMSE[4] = rmse(Outsample,ForeTSR)
Metrik_Outsample$MAPE[1] = mape(Outsample,ForeNaiveLinear)
Metrik_Outsample$MAPE[2] = mape(Outsample,ForeNaiveExp)
Metrik_Outsample$MAPE[3] = mape(Outsample,ForeSARIMA)
Metrik_Outsample$MAPE[4] = mape(Outsample,ForeTSR)
Metrik_Outsample$SMAPE[1] = smape(Outsample,ForeNaiveLinear)
Metrik_Outsample$SMAPE[2] = smape(Outsample,ForeNaiveExp)
Metrik_Outsample$SMAPE[3] = smape(Outsample,ForeSARIMA)
Metrik_Outsample$SMAPE[4] = smape(Outsample,ForeTSR)
Metrik_Insample;Metrik_Outsample
## RMSE MAPE SMAPE
## Naive Linear 16.398654 0.9233378 0.4707805
## Naive Exponential 110.806949 1.0041659 0.3776846
## SARIMA 8.989157 0.1604673 0.1875439
## TSR 4.017813 0.6821315 0.4739551
## RMSE MAPE SMAPE
## Naive Linear 21.71085 1.0861933 0.7222614
## Naive Exponential 21.71085 1.0861933 0.7222614
## SARIMA 22.18820 0.4589255 0.5969660
## TSR 11.78169 0.7491415 0.5683089