Jawaharlal Nehru Krishi Vishwavidyalaya, India
# Imports
library(fpp2)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## -- Attaching packages ---------------------------------------------- fpp2 2.4 --
## v ggplot2 3.3.2 v fma 2.4
## v forecast 8.13 v expsmooth 2.3
##
library(forecast)
library(ggplot2)
library("readxl")
library(moments)
library(forecast)
library(here)
## here() starts at F:/Phd/Rain Water in Egypt
require(forecast)
require(tseries)
## Loading required package: tseries
require(markovchain)
## Loading required package: markovchain
## Package: markovchain
## Version: 0.8.5-3
## Date: 2020-12-03
## BugReport: https://github.com/spedygiorgio/markovchain/issues
require(data.table)
## Loading required package: data.table
#Global variable
Full_original_data <-read.csv("F:/Phd/Rain Water in Egypt/data/rain water.csv")
y_lab<- "Average Rainfall - (MM) In Egypt" # input name of data
Actual_date_interval <- c("1901/12/31","2016/12/31")
Forecast_date_interval <- c("2017/01/31","2017/12/31")
validation_data_days <-279
frequency<-"month"
# Data Preparation & calculate some of statistics measures
original_data<-(Full_original_data$Rainfall....MM.)
summary(original_data)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09002 1.15922 2.10026 2.77203 3.48025 23.72680
sd(original_data) # calculate standard deviation
## [1] 2.466202
skewness(original_data) # calculate Cofficient of skewness
## [1] 2.463856
kurtosis(original_data) # calculate Cofficient of kurtosis
## [1] 12.74796
rows <- NROW(original_data)
training_data<-original_data[1:(rows-validation_data_days)]
testing_data<-original_data[(rows-validation_data_days+1):rows]
AD<-fulldate<-seq(as.Date(Actual_date_interval[1]),as.Date(Actual_date_interval[2]), frequency) #input range for actual date
FD<-seq(as.Date(Forecast_date_interval[1]),as.Date(Forecast_date_interval[2]), frequency) #input range forecasting date
N_forecasting_days<-nrow(data.frame(FD))
validation_dates<-tail(AD,validation_data_days)
validation_data_by_name<-weekdays(validation_dates)
forecasting_data_by_name<-weekdays(FD)
data_series<-ts(training_data)
#NNAR Model
model_NNAR<-nnetar(data_series, size = 2)
saveRDS(model_NNAR, file = "model_NNAR.RDS")
my_model <- readRDS("model_NNAR.RDS")
accuracy(model_NNAR) # accuracy on training data
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0008149085 1.990764 1.339084 -73.45535 94.55947 0.6105124
## ACF1
## Training set 0.02006429
# Print Model Parameters
model_NNAR
## Series: data_series
## Model: NNAR(25,2)
## Call: nnetar(y = data_series, size = 2)
##
## Average of 20 networks, each of which is
## a 25-2-1 network with 55 weights
## options were - linear output units
##
## sigma^2 estimated as 3.963
# Testing Data Evaluation
forecasting_NNAR <- predict(model_NNAR, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_NNAR$mean,validation_data_days)
MAPE_Per_Day<-round( abs(((testing_data-validation_forecast)/testing_data)*100) ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using NNAR Model for ==> ",y_lab, sep=" ")
## [1] "MAPE % For 279 month by using NNAR Model for ==> Average Rainfall - (MM) In Egypt"
MAPE_Mean_All<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_NNAR<-paste(round(MAPE_Per_Day,3),"%")
MAPE_NNAR_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in NNAR Model for ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for 279 days in NNAR Model for ==> Average Rainfall - (MM) In Egypt"
paste(MAPE_Mean_All,"%")
## [1] "146.097 % MAPE 279 month Average Rainfall - (MM) In Egypt %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in NNAR Model for ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for 279 days in NNAR Model for ==> Average Rainfall - (MM) In Egypt"
plot(forecasting_NNAR)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)
graph1<-autoplot(forecasting_NNAR,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="black", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4, ylab=y_lab)
graph1
saveRDS(model_NNAR, file = "model_NNAR.RDS")
## Error of forecasting
Error_NNAR<-abs(testing_data-validation_forecast) # Absolute error of forecast (AEOF)
REOF_A_NNAR<-abs(((testing_data-validation_forecast)/testing_data)*100) #Relative error of forecast (divided by actual)(REOF_A)
REOF_F_NNAR<-abs(((testing_data-validation_forecast)/validation_forecast)*100) #Relative error of forecast (divided by forecast)(REOF_F)
correlation_NNAR<-cor(testing_data,validation_forecast, method = c("pearson")) # correlation coefficient between predicted and actual values
RMSE_NNAR<-sqrt(sum((Error_NNAR^2))/validation_data_days) # Root mean square forecast error
MSE_NNAR<-(sum((Error_NNAR^2))/validation_data_days) # Root mean square forecast error
MAD_NNAR<-abs((sum(testing_data-validation_forecast))/validation_data_days) # average forecast accuracy
AEOF_NNAR<-c(Error_NNAR)
REOF_ANNAR<-c(paste(round(REOF_A_NNAR,3),"%"))
REOF_FNNAR<-c(paste(round(REOF_F_NNAR,3),"%"))
data.frame(correlation_NNAR,MSE_NNAR,RMSE_NNAR,MAPE_Mean_All,MAD_NNAR) # analysis of Error by using NNAR Model shows result of correlation ,MSE ,MPER
## correlation_NNAR MSE_NNAR RMSE_NNAR
## 1 0.3288594 4.07024 2.017484
## MAPE_Mean_All MAD_NNAR
## 1 146.097 % MAPE 279 month Average Rainfall - (MM) In Egypt 0.5869489
##bats model
# Data Modeling
data_series<-ts(training_data)
autoplot(ts(original_data) ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
model_bats<-bats(data_series)
accuracy(model_bats) # accuracy on training data
## ME RMSE MAE MPE MAPE MASE
## Training set 0.6969259 2.585768 1.596551 -41.22989 81.82567 0.7278963
## ACF1
## Training set 0.008805843
# Print Model Parameters
model_bats
## BATS(0.001, {3,1}, -, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0.001118
## Alpha: 0.0216578
## AR coefficients: 0.819997 -0.166157 -0.066572
## MA coefficients: -0.667452
##
## Seed States:
## [,1]
## [1,] 0.9941995
## [2,] 0.0000000
## [3,] 0.0000000
## [4,] 0.0000000
## [5,] 0.0000000
## attr(,"lambda")
## [1] 0.001117534
##
## Sigma: 0.8014665
## AIC: 8949.284
plot(model_bats,xlab = paste ("Time in",frequency ,y_lab , sep=" "), col.main="black", col.lab="black", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4)
# Testing Data Evaluation
forecasting_bats <- predict(model_bats, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_bats$mean,validation_data_days)
MAPE_Per_Day<-round( abs(((testing_data-validation_forecast)/testing_data)*100) ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using bats Model for ==> ",y_lab, sep=" ")
## [1] "MAPE % For 279 month by using bats Model for ==> Average Rainfall - (MM) In Egypt"
MAPE_Mean_All<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_bats<-paste(round(MAPE_Per_Day,3),"%")
MAPE_bats_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in bats Model for ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for 279 days in bats Model for ==> Average Rainfall - (MM) In Egypt"
paste(MAPE_Mean_All,"%")
## [1] "104.867 % MAPE 279 month Average Rainfall - (MM) In Egypt %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in bats Model for ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for 279 days in bats Model for ==> Average Rainfall - (MM) In Egypt"
par(bg = "#F5F5F5")
plot(forecasting_bats)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)
graph1<-autoplot(forecasting_bats,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="blue", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4, ylab=y_lab)
graph1
## Error of forecasting
Error_bats<-abs(testing_data-validation_forecast) # Absolute error of forecast (AEOF)
MSE_Bats<-((sum(Error_bats^2)/validation_data_days))
MSE_Bats
## [1] 4.219481
REOF_A_bats<-abs(((testing_data-validation_forecast)/testing_data)*100) #Relative error of forecast (divided by actual)(REOF_A)
REOF_F_bats<-abs(((testing_data-validation_forecast)/validation_forecast)*100) #Relative error of forecast (divided by forecast)(REOF_F)
correlation_bats<-cor(testing_data,validation_forecast, method = c("pearson")) # correlation coefficient between predicted and actual values
RMSE_bats<-sqrt(sum((Error_bats^2))/validation_data_days) # Root mean square forecast error
MAD_bats<-abs((sum(testing_data-validation_forecast))/validation_data_days) # average forecast accuracy
AEOF_bats<-c(Error_bats)
REOF_Abats<-c(paste(round(REOF_A_bats,3),"%"))
REOF_Fbats<-c(paste(round(REOF_F_bats,3),"%"))
data.frame(correlation_bats,RMSE_bats,MAPE_Mean_All,MAD_bats) # analysis of Error by using Bats Model shows result of correlation ,MSE ,MPER
## correlation_bats RMSE_bats
## 1 0.06133969 2.054138
## MAPE_Mean_All MAD_bats
## 1 104.867 % MAPE 279 month Average Rainfall - (MM) In Egypt 0.3497803
## TBATS Model
# Data Modeling
data_series<-ts(training_data)
model_TBATS<-forecast:::fitSpecificTBATS(data_series,use.box.cox=FALSE, use.beta=TRUE, seasonal.periods=c(6),use.damping=FALSE,k.vector=c(2))
accuracy(model_TBATS) # accuracy on training data
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.07732442 2.560841 1.786943 -78.1541 109.2773 0.8146995 0.1089734
# Print Model Parameters
model_TBATS
## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 0.08391584
## Beta: 0.004256946
## Damping Parameter: 1
## Gamma-1 Values: 0.0007324185
## Gamma-2 Values: 0.002158919
##
## Seed States:
## [,1]
## [1,] 5.13130984
## [2,] -0.35706247
## [3,] 0.83815101
## [4,] 0.05626826
## [5,] -0.15465285
## [6,] 0.20866529
##
## Sigma: 2.560841
## AIC: 9920.676
plot(model_TBATS,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="black", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4, ylab=y_lab)
# Testing Data Evaluation
forecasting_tbats <- predict(model_TBATS, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_tbats$mean,validation_data_days)
MAPE_Per_Day<-round( abs(((testing_data-validation_forecast)/testing_data)*100) ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using TBATS Model for ==> ",y_lab, sep=" ")
## [1] "MAPE % For 279 month by using TBATS Model for ==> Average Rainfall - (MM) In Egypt"
MAPE_Mean_All<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_TBATS<-paste(round(MAPE_Per_Day,3),"%")
MAPE_TBATS_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in TBATS Model for ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for 279 days in TBATS Model for ==> Average Rainfall - (MM) In Egypt"
paste(MAPE_Mean_All,"%")
## [1] "266.145 % MAPE 279 month Average Rainfall - (MM) In Egypt %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in TBATS Model for ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for 279 days in TBATS Model for ==> Average Rainfall - (MM) In Egypt"
plot(forecasting_tbats)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)
graph2<-autoplot(forecasting_tbats,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="blue", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4, ylab=y_lab)
graph2
## Error of forecasting TBATS Model
Error_tbats<-abs(testing_data-validation_forecast) # Absolute error of forecast (AEOF)
sqrt(sum(Error_tbats^2)/validation_data_days)
## [1] 2.993632
REOF_A_tbats1<-abs(((testing_data-validation_forecast)/testing_data)*100) #Relative error of forecast (divided by actual)(REOF_A)
REOF_F_tbats<-abs(((testing_data-validation_forecast)/validation_forecast)*100) #Relative error of forecast (divided by forecast)(REOF_F)
correlation_tbats<-cor(testing_data,validation_forecast, method = c("pearson")) # correlation coefficient between predicted and actual values
RMSE_tbats<-sqrt(sum((Error_tbats^2))/validation_data_days) # Root mean square forecast error
MAD_tbats<-abs((sum(testing_data-validation_forecast))/validation_data_days) # average forecast accuracy
AEOF_tbats<-c(Error_tbats)
REOF_A_tbats<-c(paste(round(REOF_A_tbats1,3),"%"))
REOF_F_tbats<-c(paste(round(REOF_F_tbats,3),"%"))
data.frame(correlation_tbats,RMSE_tbats,MAPE_Mean_All,MAD_tbats) # analysis of Error by using Holt's linear model shows result of correlation ,MSE ,MPER
## correlation_tbats RMSE_tbats
## 1 0.1784506 2.993632
## MAPE_Mean_All MAD_tbats
## 1 266.145 % MAPE 279 month Average Rainfall - (MM) In Egypt 2.150012
## Holt's linear trend
# Data Modeling
data_series<-ts(training_data)
model_holt<-holt(data_series,h=N_forecasting_days+validation_data_days,lambda = "auto")
accuracy(model_holt) # accuracy on training data
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.841245 2.670937 1.645376 -36.80145 81.65181 0.7501562 0.1232587
# Print Model Parameters
summary(model_holt$model)
## Holt's method
##
## Call:
## holt(y = data_series, h = N_forecasting_days + validation_data_days,
##
## Call:
## lambda = "auto")
##
## Box-Cox transformation: lambda= -0.1044
##
## Smoothing parameters:
## alpha = 0.0153
## beta = 1e-04
##
## Initial states:
## l = 0.9983
## b = -5e-04
##
## sigma: 0.7772
##
## AIC AICc BIC
## 7252.354 7252.408 7277.428
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.841245 2.670937 1.645376 -36.80145 81.65181 0.7501562 0.1232587
# Testing Data Evaluation
forecasting_holt <- predict(model_holt, h=N_forecasting_days+validation_data_days,lambda = "auto")
validation_forecast<-head(forecasting_holt$mean,validation_data_days)
MAPE_Per_Day<-round( abs(((testing_data-validation_forecast)/testing_data)*100) ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using holt Model for ==> ",y_lab, sep=" ")
## [1] "MAPE % For 279 month by using holt Model for ==> Average Rainfall - (MM) In Egypt"
MAPE_Mean_All<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_holt<-paste(round(MAPE_Per_Day,3),"%")
MAPE_holt_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in holt Model for ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for 279 days in holt Model for ==> Average Rainfall - (MM) In Egypt"
paste(MAPE_Mean_All,"%")
## [1] "108.079 % MAPE 279 month Average Rainfall - (MM) In Egypt %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in holt Model for ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for 279 days in holt Model for ==> Average Rainfall - (MM) In Egypt"
plot(forecasting_holt,ylim=c(0,20))
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)
graph3<-autoplot(forecasting_holt,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="blue", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4, ylab=y_lab, ylim=c(0,20))
graph3
## Error of forecasting by using Holt's linear model
Error_Holt<-abs(testing_data-validation_forecast) # Absolute error of forecast (AEOF)
sqrt(sum(Error_Holt^2)/validation_data_days)
## [1] 2.048031
REOF_A_Holt1<-abs(((testing_data-validation_forecast)/testing_data)*100) #Relative error of forecast (divided by actual)(REOF_A)
REOF_F_Holt<-abs(((testing_data-validation_forecast)/validation_forecast)*100) #Relative error of forecast (divided by forecast)(REOF_F)
correlation_Holt<-cor(testing_data,validation_forecast, method = c("pearson")) # correlation coefficient between predicted and actual values
RMSE_Holt<-sqrt(sum((Error_Holt^2))/validation_data_days) # Root mean square forecast error
MAD_Holt<-abs((sum(testing_data-validation_forecast))/validation_data_days) # average forecast accuracy
AEOF_Holt<-c(Error_Holt)
REOF_A_Holt<-c(paste(round(REOF_A_Holt1,3),"%"))
REOF_F_Holt<-c(paste(round(REOF_F_Holt,3),"%"))
REOF_A_Holt11<-mean(abs(((testing_data-validation_forecast)/testing_data)*100))
data.frame(correlation_Holt,RMSE_Holt,MAPE_Mean_All,MAD_Holt) # analysis of Error by using Holt's linear model shows result of correlation ,MSE ,MPER
## correlation_Holt RMSE_Holt
## 1 0.001042441 2.048031
## MAPE_Mean_All MAD_Holt
## 1 108.079 % MAPE 279 month Average Rainfall - (MM) In Egypt 0.2968722
##################
#ARIMA Model
require(tseries) # need to install tseries tj test Stationarity in time series
paste ("tests For Check Stationarity in series ==> ",y_lab, sep=" ")
## [1] "tests For Check Stationarity in series ==> Average Rainfall - (MM) In Egypt"
kpss.test(data_series) # applay kpss test
## Warning in kpss.test(data_series): p-value smaller than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: data_series
## KPSS Level = 1.0065, Truncation lag parameter = 7, p-value = 0.01
pp.test(data_series) # applay pp test
## Warning in pp.test(data_series): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: data_series
## Dickey-Fuller Z(alpha) = -871.96, Truncation lag parameter = 7, p-value
## = 0.01
## alternative hypothesis: stationary
adf.test(data_series) # applay adf test
## Warning in adf.test(data_series): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: data_series
## Dickey-Fuller = -9.5506, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
ndiffs(data_series) # Doing first diffrencing on data
## [1] 1
#Taking the first difference
diff1_x1<-diff(data_series)
autoplot(diff1_x1, xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="black", col.sub="black", ylab=y_lab,main = "1nd differenced series")
## Warning: Ignoring unknown parameters: col.main, col.lab, col.sub
##Testing the stationary of the first differenced series
paste ("tests For Check Stationarity in series after taking first differences in ==> ",y_lab, sep=" ")
## [1] "tests For Check Stationarity in series after taking first differences in ==> Average Rainfall - (MM) In Egypt"
kpss.test(diff1_x1) # applay kpss test after taking first differences
## Warning in kpss.test(diff1_x1): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff1_x1
## KPSS Level = 0.013801, Truncation lag parameter = 7, p-value = 0.1
pp.test(diff1_x1) # applay pp test after taking first differences
## Warning in pp.test(diff1_x1): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: diff1_x1
## Dickey-Fuller Z(alpha) = -1302.3, Truncation lag parameter = 7, p-value
## = 0.01
## alternative hypothesis: stationary
adf.test(diff1_x1) # applay adf test after taking first differences
## Warning in adf.test(diff1_x1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff1_x1
## Dickey-Fuller = -24.701, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
#Taking the second difference
diff2_x1=diff(diff1_x1)
autoplot(diff2_x1, xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="black", col.sub="black", ylab=y_lab ,main = "2nd differenced series")
## Warning: Ignoring unknown parameters: col.main, col.lab, col.sub
##Testing the stationary of the first differenced series
paste ("tests For Check Stationarity in series after taking Second differences in",y_lab, sep=" ")
## [1] "tests For Check Stationarity in series after taking Second differences in Average Rainfall - (MM) In Egypt"
kpss.test(diff2_x1) # applay kpss test after taking Second differences
## Warning in kpss.test(diff2_x1): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff2_x1
## KPSS Level = 0.014869, Truncation lag parameter = 7, p-value = 0.1
pp.test(diff2_x1) # applay pp test after taking Second differences
## Warning in pp.test(diff2_x1): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: diff2_x1
## Dickey-Fuller Z(alpha) = -1593.8, Truncation lag parameter = 7, p-value
## = 0.01
## alternative hypothesis: stationary
adf.test(diff2_x1) # applay adf test after taking Second differences
## Warning in adf.test(diff2_x1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff2_x1
## Dickey-Fuller = -23.165, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
####Fitting an ARIMA Model
#1. Using auto arima function
model1 <- auto.arima(data_series,stepwise=FALSE, approximation=FALSE, trace=T, test = c("kpss", "adf", "pp")) #applaying auto arima
##
## ARIMA(0,1,0) : 5836.636
## ARIMA(0,1,0) with drift : 5838.633
## ARIMA(0,1,1) : 5237.635
## ARIMA(0,1,1) with drift : 5239.34
## ARIMA(0,1,2) : 5223.895
## ARIMA(0,1,2) with drift : Inf
## ARIMA(0,1,3) : Inf
## ARIMA(0,1,3) with drift : Inf
## ARIMA(0,1,4) : 5215.82
## ARIMA(0,1,4) with drift : Inf
## ARIMA(0,1,5) : 5202.892
## ARIMA(0,1,5) with drift : 5204.662
## ARIMA(1,1,0) : 5575.081
## ARIMA(1,1,0) with drift : 5577.079
## ARIMA(1,1,1) : Inf
## ARIMA(1,1,1) with drift : Inf
## ARIMA(1,1,2) : Inf
## ARIMA(1,1,2) with drift : Inf
## ARIMA(1,1,3) : Inf
## ARIMA(1,1,3) with drift : Inf
## ARIMA(1,1,4) : 5197.98
## ARIMA(1,1,4) with drift : 5199.776
## ARIMA(2,1,0) : 5525.909
## ARIMA(2,1,0) with drift : 5527.911
## ARIMA(2,1,1) : Inf
## ARIMA(2,1,1) with drift : Inf
## ARIMA(2,1,2) : 5212.798
## ARIMA(2,1,2) with drift : 5214.588
## ARIMA(2,1,3) : Inf
## ARIMA(2,1,3) with drift : Inf
## ARIMA(3,1,0) : 5488.79
## ARIMA(3,1,0) with drift : 5490.796
## ARIMA(3,1,1) : 5212.569
## ARIMA(3,1,1) with drift : 5214.244
## ARIMA(3,1,2) : 5191.777
## ARIMA(3,1,2) with drift : 5193.616
## ARIMA(4,1,0) : 5437.232
## ARIMA(4,1,0) with drift : 5439.243
## ARIMA(4,1,1) : 5203.876
## ARIMA(4,1,1) with drift : 5205.602
## ARIMA(5,1,0) : 5420.249
## ARIMA(5,1,0) with drift : 5422.264
##
##
##
## Best model: ARIMA(3,1,2)
model1 # show the result of autoarima
## Series: data_series
## ARIMA(3,1,2)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2
## 0.8226 -0.0269 -0.1564 -1.7079 0.7128
## s.e. 0.0615 0.0395 0.0308 0.0560 0.0555
##
## sigma^2 estimated as 6.178: log likelihood=-2589.85
## AIC=5191.7 AICc=5191.78 BIC=5221.78
#Make changes in the source of auto arima to run the best model
arima.string <- function (object, padding = FALSE)
{
order <- object$arma[c(1, 6, 2, 3, 7, 4, 5)]
m <- order[7]
result <- paste("ARIMA(", order[1], ",", order[2], ",",
order[3], ")", sep = "")
if (m > 1 && sum(order[4:6]) > 0) {
result <- paste(result, "(", order[4], ",", order[5],
",", order[6], ")[", m, "]", sep = "")
}
if (padding && m > 1 && sum(order[4:6]) == 0) {
result <- paste(result, " ", sep = "")
if (m <= 9) {
result <- paste(result, " ", sep = "")
}
else if (m <= 99) {
result <- paste(result, " ", sep = "")
}
else {
result <- paste(result, " ", sep = "")
}
}
if (!is.null(object$xreg)) {
if (NCOL(object$xreg) == 1 && is.element("drift", names(object$coef))) {
result <- paste(result, "with drift ")
}
else {
result <- paste("Regression with", result, "errors")
}
}
else {
if (is.element("constant", names(object$coef)) || is.element("intercept",
names(object$coef))) {
result <- paste(result, "with non-zero mean")
}
else if (order[2] == 0 && order[5] == 0) {
result <- paste(result, "with zero mean ")
}
else {
result <- paste(result, " ")
}
}
if (!padding) {
result <- gsub("[ ]*$", "", result)
}
return(result)
}
bestmodel <- arima.string(model1, padding = TRUE)
bestmodel <- substring(bestmodel,7,11)
bestmodel <- gsub(" ", "", bestmodel)
bestmodel <- gsub(")", "", bestmodel)
bestmodel <- strsplit(bestmodel, ",")[[1]]
bestmodel <- c(strtoi(bestmodel[1]),strtoi(bestmodel[2]),strtoi(bestmodel[3]))
bestmodel
## [1] 3 1 2
strtoi(bestmodel[3])
## [1] 2
#2. Using ACF and PACF Function
#par(mfrow=c(1,2)) # Code for making two plot in one graph
acf(diff2_x1,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), col.main="black", col.lab="black", col.sub="black", cex.main=1, cex.lab=1, cex.sub=1,font.main=4, font.lab=4, ylab=y_lab, main=paste("ACF-2nd differenced series ",y_lab, sep=" ",lag.max=20)) # plot ACF "auto correlation function after taking second diffrences
pacf(diff2_x1,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), ylab=y_lab,main=paste("PACF-2nd differenced series ",y_lab, sep=" ",lag.max=20)) # plot PACF " Partial auto correlation function after taking second diffrences
library(forecast) # install library forecast
x1_model1= arima(data_series, order=c(bestmodel)) # Run Best model of auto arima for forecasting
x1_model1 # Show result of best model of auto arima
##
## Call:
## arima(x = data_series, order = c(bestmodel))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2
## 0.8226 -0.0269 -0.1564 -1.7079 0.7128
## s.e. 0.0615 0.0395 0.0308 0.0560 0.0555
##
## sigma^2 estimated as 6.151: log likelihood = -2589.85, aic = 5191.7
paste ("accuracy of autoarima Model For ==> ",y_lab, sep=" ")
## [1] "accuracy of autoarima Model For ==> Average Rainfall - (MM) In Egypt"
accuracy(x1_model1) # aacuracy of best model from auto arima
## ME RMSE MAE MPE MAPE MASE
## Training set -0.05932172 2.478917 1.720522 -94.03555 119.6028 0.7844166
## ACF1
## Training set 0.001905031
x1_model1$x # show result of best model from auto arima
## NULL
checkresiduals(x1_model1,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), ylab=y_lab) # checkresiduals from best model from using auto arima
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,1,2)
## Q* = 4.4631, df = 5, p-value = 0.4848
##
## Model df: 5. Total lags used: 10
paste("Box-Ljung test , Ljung-Box test For Modelling for ==> ",y_lab, sep=" ")
## [1] "Box-Ljung test , Ljung-Box test For Modelling for ==> Average Rainfall - (MM) In Egypt"
Box.test(x1_model1$residuals^2, lag=20, type="Ljung-Box") # Do test for resdulas by using Box-Ljung test , Ljung-Box test For Modelling
##
## Box-Ljung test
##
## data: x1_model1$residuals^2
## X-squared = 38.613, df = 20, p-value = 0.007446
library(tseries)
jarque.bera.test(x1_model1$residuals) # Do test jarque.bera.test
##
## Jarque Bera Test
##
## data: x1_model1$residuals
## X-squared = 5742.9, df = 2, p-value < 2.2e-16
#Actual Vs Fitted
plot(data_series, col='red',lwd=2, main="Actual vs Fitted Plot", xlab='Time in (days)', ylab=y_lab) # plot actual and Fitted model
lines(fitted(x1_model1), col='black')
#Test data
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) ) # make testing data in time series and start from rows-6
forecasting_auto_arima <- forecast(x1_model1, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_auto_arima$mean,validation_data_days)
MAPE_Per_Day<-round(abs(((testing_data-validation_forecast)/testing_data)*100) ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using bats Model for ==> ",y_lab, sep=" ")
## [1] "MAPE % For 279 month by using bats Model for ==> Average Rainfall - (MM) In Egypt"
MAPE_Mean_All.ARIMA_Model<-round(mean(MAPE_Per_Day),3)
MAPE_Mean_All.ARIMA<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_auto_arima<-paste(round(MAPE_Per_Day,3),"%")
MAPE_auto.arima_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in bats Model for ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for 279 days in bats Model for ==> Average Rainfall - (MM) In Egypt"
paste(MAPE_Mean_All.ARIMA,"%")
## [1] "170.734 % MAPE 279 month Average Rainfall - (MM) In Egypt %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in bats Model for ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for 279 days in bats Model for ==> Average Rainfall - (MM) In Egypt"
data.frame(date_auto.arima=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_auto.arima=validation_forecast,MAPE_auto.arima_Model)
## date_auto.arima validation_data_by_name actual_data forecasting_auto.arima
## 1 1993-10-31 Sunday 4.74420 2.052765
## 2 1993-12-01 Wednesday 0.50180 2.551296
## 3 1993-12-31 Friday 4.47208 3.023406
## 4 1994-01-31 Monday 8.10726 3.200987
## 5 1994-03-03 Thursday 1.12499 3.256404
## 6 1994-03-31 Thursday 3.53859 3.223390
## 7 1994-05-01 Sunday 0.50400 3.166974
## 8 1994-05-31 Tuesday 0.95792 3.112790
## 9 1994-07-01 Friday 1.29018 3.074900
## 10 1994-07-31 Sunday 3.86952 3.054011
## 11 1994-08-31 Wednesday 1.46727 3.046321
## 12 1994-10-01 Saturday 0.96778 3.046482
## 13 1994-10-31 Monday 5.00519 3.050087
## 14 1994-12-01 Thursday 10.49740 3.054251
## 15 1994-12-31 Saturday 4.13700 3.057554
## 16 1995-01-31 Tuesday 1.00520 3.059595
## 17 1995-03-03 Friday 4.64740 3.060534
## 18 1995-03-31 Friday 1.20847 3.060736
## 19 1995-05-01 Monday 1.63924 3.060556
## 20 1995-05-31 Wednesday 0.95959 3.060257
## 21 1995-07-01 Saturday 1.30764 3.059984
## 22 1995-07-31 Monday 2.75944 3.059795
## 23 1995-08-31 Thursday 1.55201 3.059694
## 24 1995-10-01 Sunday 1.68352 3.059659
## 25 1995-10-31 Tuesday 0.94816 3.059662
## 26 1995-12-01 Friday 2.06694 3.059682
## 27 1995-12-31 Sunday 0.85773 3.059703
## 28 1996-01-31 Wednesday 3.92626 3.059720
## 29 1996-03-02 Saturday 1.06618 3.059730
## 30 1996-03-31 Sunday 3.27342 3.059734
## 31 1996-05-01 Wednesday 2.25442 3.059735
## 32 1996-05-31 Friday 0.58514 3.059734
## 33 1996-07-01 Monday 1.29002 3.059732
## 34 1996-07-31 Wednesday 2.41815 3.059731
## 35 1996-08-31 Saturday 1.37331 3.059730
## 36 1996-10-01 Tuesday 0.79516 3.059730
## 37 1996-10-31 Thursday 1.20975 3.059729
## 38 1996-12-01 Sunday 2.60112 3.059729
## 39 1996-12-31 Tuesday 1.12006 3.059730
## 40 1997-01-31 Friday 3.86200 3.059730
## 41 1997-03-03 Monday 2.84933 3.059730
## 42 1997-03-31 Monday 4.21092 3.059730
## 43 1997-05-01 Thursday 0.77283 3.059730
## 44 1997-05-31 Saturday 1.93075 3.059730
## 45 1997-07-01 Tuesday 1.82400 3.059730
## 46 1997-07-31 Thursday 2.53188 3.059730
## 47 1997-08-31 Sunday 1.37305 3.059730
## 48 1997-10-01 Wednesday 1.31306 3.059730
## 49 1997-10-31 Friday 6.88818 3.059730
## 50 1997-12-01 Monday 1.07951 3.059730
## 51 1997-12-31 Wednesday 2.48236 3.059730
## 52 1998-01-31 Saturday 2.35077 3.059730
## 53 1998-03-03 Tuesday 3.23532 3.059730
## 54 1998-03-31 Tuesday 2.65382 3.059730
## 55 1998-05-01 Friday 3.64956 3.059730
## 56 1998-05-31 Sunday 11.39700 3.059730
## 57 1998-07-01 Wednesday 1.28964 3.059730
## 58 1998-07-31 Friday 2.89250 3.059730
## 59 1998-08-31 Monday 3.33016 3.059730
## 60 1998-10-01 Thursday 2.18281 3.059730
## 61 1998-10-31 Saturday 0.27669 3.059730
## 62 1998-12-01 Tuesday 1.36996 3.059730
## 63 1998-12-31 Thursday 3.96099 3.059730
## 64 1999-01-31 Sunday 1.45933 3.059730
## 65 1999-03-03 Wednesday 6.14986 3.059730
## 66 1999-03-31 Wednesday 0.34525 3.059730
## 67 1999-05-01 Saturday 0.87614 3.059730
## 68 1999-05-31 Monday 3.92148 3.059730
## 69 1999-07-01 Thursday 1.28941 3.059730
## 70 1999-07-31 Saturday 2.75452 3.059730
## 71 1999-08-31 Tuesday 1.56583 3.059730
## 72 1999-10-01 Friday 1.30709 3.059730
## 73 1999-10-31 Sunday 0.40461 3.059730
## 74 1999-12-01 Wednesday 0.79104 3.059730
## 75 1999-12-31 Friday 1.06074 3.059730
## 76 2000-01-31 Monday 8.79867 3.059730
## 77 2000-03-02 Thursday 1.96720 3.059730
## 78 2000-03-31 Friday 0.94208 3.059730
## 79 2000-05-01 Monday 0.46715 3.059730
## 80 2000-05-31 Wednesday 1.39287 3.059730
## 81 2000-07-01 Saturday 1.28941 3.059730
## 82 2000-07-31 Monday 2.51911 3.059730
## 83 2000-08-31 Thursday 1.41707 3.059730
## 84 2000-10-01 Sunday 2.09922 3.059730
## 85 2000-10-31 Tuesday 3.18534 3.059730
## 86 2000-12-01 Friday 4.34451 3.059730
## 87 2000-12-31 Sunday 6.13581 3.059730
## 88 2001-01-31 Wednesday 2.30453 3.059730
## 89 2001-03-03 Saturday 1.98704 3.059730
## 90 2001-03-31 Saturday 1.51554 3.059730
## 91 2001-05-01 Tuesday 1.02323 3.059730
## 92 2001-05-31 Thursday 0.56850 3.059730
## 93 2001-07-01 Sunday 1.28941 3.059730
## 94 2001-07-31 Tuesday 2.51285 3.059730
## 95 2001-08-31 Friday 1.41755 3.059730
## 96 2001-10-01 Monday 0.77372 3.059730
## 97 2001-10-31 Wednesday 0.86521 3.059730
## 98 2001-12-01 Saturday 1.40814 3.059730
## 99 2001-12-31 Monday 3.16058 3.059730
## 100 2002-01-31 Thursday 6.39406 3.059730
## 101 2002-03-03 Sunday 5.24313 3.059730
## 102 2002-03-31 Sunday 1.00755 3.059730
## 103 2002-05-01 Wednesday 1.13408 3.059730
## 104 2002-05-31 Friday 0.50675 3.059730
## 105 2002-07-01 Monday 1.29804 3.059730
## 106 2002-07-31 Wednesday 2.43717 3.059730
## 107 2002-08-31 Saturday 1.54944 3.059730
## 108 2002-10-01 Tuesday 0.81303 3.059730
## 109 2002-10-31 Thursday 6.49295 3.059730
## 110 2002-12-01 Sunday 0.79751 3.059730
## 111 2002-12-31 Tuesday 3.44896 3.059730
## 112 2003-01-31 Friday 1.60621 3.059730
## 113 2003-03-03 Monday 3.74258 3.059730
## 114 2003-03-31 Monday 3.40988 3.059730
## 115 2003-05-01 Thursday 0.44270 3.059730
## 116 2003-05-31 Saturday 0.56668 3.059730
## 117 2003-07-01 Tuesday 1.41158 3.059730
## 118 2003-07-31 Thursday 3.09183 3.059730
## 119 2003-08-31 Sunday 2.91996 3.059730
## 120 2003-10-01 Wednesday 0.78733 3.059730
## 121 2003-10-31 Friday 0.23788 3.059730
## 122 2003-12-01 Monday 0.98593 3.059730
## 123 2003-12-31 Wednesday 4.39334 3.059730
## 124 2004-01-31 Saturday 7.90062 3.059730
## 125 2004-03-02 Tuesday 3.54160 3.059730
## 126 2004-03-31 Wednesday 2.82262 3.059730
## 127 2004-05-01 Saturday 0.37460 3.059730
## 128 2004-05-31 Monday 0.59284 3.059730
## 129 2004-07-01 Thursday 1.28963 3.059730
## 130 2004-07-31 Saturday 2.51388 3.059730
## 131 2004-08-31 Tuesday 1.61934 3.059730
## 132 2004-10-01 Friday 0.90617 3.059730
## 133 2004-10-31 Sunday 0.32351 3.059730
## 134 2004-12-01 Wednesday 1.60803 3.059730
## 135 2004-12-31 Friday 2.82753 3.059730
## 136 2005-01-31 Monday 3.79083 3.059730
## 137 2005-03-03 Thursday 4.93790 3.059730
## 138 2005-03-31 Thursday 5.90618 3.059730
## 139 2005-05-01 Sunday 5.32605 3.059730
## 140 2005-05-31 Tuesday 0.37812 3.059730
## 141 2005-07-01 Friday 1.28955 3.059730
## 142 2005-07-31 Sunday 3.32937 3.059730
## 143 2005-08-31 Wednesday 1.55622 3.059730
## 144 2005-10-01 Saturday 1.85197 3.059730
## 145 2005-10-31 Monday 0.54443 3.059730
## 146 2005-12-01 Thursday 1.06465 3.059730
## 147 2005-12-31 Saturday 2.56311 3.059730
## 148 2006-01-31 Tuesday 3.09579 3.059730
## 149 2006-03-03 Friday 2.55506 3.059730
## 150 2006-03-31 Friday 2.62328 3.059730
## 151 2006-05-01 Monday 3.43886 3.059730
## 152 2006-05-31 Wednesday 0.47531 3.059730
## 153 2006-07-01 Saturday 1.32460 3.059730
## 154 2006-07-31 Monday 2.44671 3.059730
## 155 2006-08-31 Thursday 1.76639 3.059730
## 156 2006-10-01 Sunday 0.92156 3.059730
## 157 2006-10-31 Tuesday 0.74955 3.059730
## 158 2006-12-01 Friday 0.85919 3.059730
## 159 2006-12-31 Sunday 2.70497 3.059730
## 160 2007-01-31 Wednesday 3.51390 3.059730
## 161 2007-03-03 Saturday 4.75467 3.059730
## 162 2007-03-31 Saturday 0.86146 3.059730
## 163 2007-05-01 Tuesday 1.80090 3.059730
## 164 2007-05-31 Thursday 1.45439 3.059730
## 165 2007-07-01 Sunday 1.50526 3.059730
## 166 2007-07-31 Tuesday 3.32094 3.059730
## 167 2007-08-31 Friday 1.92822 3.059730
## 168 2007-10-01 Monday 1.12446 3.059730
## 169 2007-10-31 Wednesday 1.18441 3.059730
## 170 2007-12-01 Saturday 1.28654 3.059730
## 171 2007-12-31 Monday 1.60943 3.059730
## 172 2008-01-31 Thursday 7.55999 3.059730
## 173 2008-03-02 Sunday 3.12595 3.059730
## 174 2008-03-31 Monday 0.26817 3.059730
## 175 2008-05-01 Thursday 3.13994 3.059730
## 176 2008-05-31 Saturday 0.56320 3.059730
## 177 2008-07-01 Tuesday 1.28941 3.059730
## 178 2008-07-31 Thursday 2.44671 3.059730
## 179 2008-08-31 Sunday 1.41721 3.059730
## 180 2008-10-01 Wednesday 2.98817 3.059730
## 181 2008-10-31 Friday 3.01147 3.059730
## 182 2008-12-01 Monday 0.30073 3.059730
## 183 2008-12-31 Wednesday 1.89878 3.059730
## 184 2009-01-31 Saturday 4.13442 3.059730
## 185 2009-03-03 Tuesday 5.61797 3.059730
## 186 2009-03-31 Tuesday 0.23356 3.059730
## 187 2009-05-01 Friday 0.14308 3.059730
## 188 2009-05-31 Sunday 0.54860 3.059730
## 189 2009-07-01 Wednesday 1.28941 3.059730
## 190 2009-07-31 Friday 2.44671 3.059730
## 191 2009-08-31 Monday 1.51752 3.059730
## 192 2009-10-01 Thursday 0.79499 3.059730
## 193 2009-10-31 Saturday 0.38449 3.059730
## 194 2009-12-01 Tuesday 0.60021 3.059730
## 195 2009-12-31 Thursday 5.60057 3.059730
## 196 2010-01-31 Sunday 2.24212 3.059730
## 197 2010-03-03 Wednesday 2.84198 3.059730
## 198 2010-03-31 Wednesday 0.24605 3.059730
## 199 2010-05-01 Saturday 0.27875 3.059730
## 200 2010-05-31 Monday 8.83731 3.059730
## 201 2010-07-01 Thursday 1.29149 3.059730
## 202 2010-07-31 Saturday 3.84991 3.059730
## 203 2010-08-31 Tuesday 8.03098 3.059730
## 204 2010-10-01 Friday 1.38485 3.059730
## 205 2010-10-31 Sunday 0.52419 3.059730
## 206 2010-12-01 Wednesday 0.26019 3.059730
## 207 2010-12-31 Friday 2.52660 3.059730
## 208 2011-01-31 Monday 1.79052 3.059730
## 209 2011-03-03 Thursday 2.35943 3.059730
## 210 2011-03-31 Thursday 2.76970 3.059730
## 211 2011-05-01 Sunday 2.99102 3.059730
## 212 2011-05-31 Tuesday 0.64547 3.059730
## 213 2011-07-01 Friday 1.64623 3.059730
## 214 2011-07-31 Sunday 2.50175 3.059730
## 215 2011-08-31 Wednesday 1.44966 3.059730
## 216 2011-10-01 Saturday 1.22527 3.059730
## 217 2011-10-31 Monday 0.40954 3.059730
## 218 2011-12-01 Thursday 4.53735 3.059730
## 219 2011-12-31 Saturday 1.27885 3.059730
## 220 2012-01-31 Tuesday 4.50333 3.059730
## 221 2012-03-02 Friday 2.71883 3.059730
## 222 2012-03-31 Saturday 1.84284 3.059730
## 223 2012-05-01 Tuesday 0.79004 3.059730
## 224 2012-05-31 Thursday 0.45443 3.059730
## 225 2012-07-01 Sunday 1.54115 3.059730
## 226 2012-07-31 Tuesday 2.48006 3.059730
## 227 2012-08-31 Friday 1.45472 3.059730
## 228 2012-10-01 Monday 1.08196 3.059730
## 229 2012-10-31 Wednesday 1.25241 3.059730
## 230 2012-12-01 Saturday 2.16875 3.059730
## 231 2012-12-31 Monday 5.50551 3.059730
## 232 2013-01-31 Thursday 9.83019 3.059730
## 233 2013-03-03 Sunday 2.08693 3.059730
## 234 2013-03-31 Sunday 0.75090 3.059730
## 235 2013-05-01 Wednesday 0.97940 3.059730
## 236 2013-05-31 Friday 0.84011 3.059730
## 237 2013-07-01 Monday 1.55332 3.059730
## 238 2013-07-31 Wednesday 2.47999 3.059730
## 239 2013-08-31 Saturday 9.08493 3.059730
## 240 2013-10-01 Tuesday 1.02467 3.059730
## 241 2013-10-31 Thursday 0.31919 3.059730
## 242 2013-12-01 Sunday 0.58748 3.059730
## 243 2013-12-31 Tuesday 6.73026 3.059730
## 244 2014-01-31 Friday 4.49754 3.059730
## 245 2014-03-03 Monday 11.76580 3.059730
## 246 2014-03-31 Monday 4.66164 3.059730
## 247 2014-05-01 Thursday 0.81749 3.059730
## 248 2014-05-31 Saturday 3.69818 3.059730
## 249 2014-07-01 Tuesday 1.55723 3.059730
## 250 2014-07-31 Thursday 3.29754 3.059730
## 251 2014-08-31 Sunday 1.49722 3.059730
## 252 2014-10-01 Wednesday 0.94467 3.059730
## 253 2014-10-31 Friday 0.87077 3.059730
## 254 2014-12-01 Monday 3.63473 3.059730
## 255 2014-12-31 Wednesday 0.71567 3.059730
## 256 2015-01-31 Saturday 3.95595 3.059730
## 257 2015-03-03 Tuesday 2.71103 3.059730
## 258 2015-03-31 Tuesday 2.63489 3.059730
## 259 2015-05-01 Friday 1.57308 3.059730
## 260 2015-05-31 Sunday 2.02107 3.059730
## 261 2015-07-01 Wednesday 1.94850 3.059730
## 262 2015-07-31 Friday 2.48587 3.059730
## 263 2015-08-31 Monday 1.90997 3.059730
## 264 2015-10-01 Thursday 1.23057 3.059730
## 265 2015-10-31 Saturday 6.97369 3.059730
## 266 2015-12-01 Tuesday 2.78345 3.059730
## 267 2015-12-31 Thursday 3.74498 3.059730
## 268 2016-01-31 Sunday 5.25677 3.059730
## 269 2016-03-02 Wednesday 1.14635 3.059730
## 270 2016-03-31 Thursday 2.81793 3.059730
## 271 2016-05-01 Sunday 2.13990 3.059730
## 272 2016-05-31 Tuesday 0.54040 3.059730
## 273 2016-07-01 Friday 1.55664 3.059730
## 274 2016-07-31 Sunday 2.77524 3.059730
## 275 2016-08-31 Wednesday 1.47609 3.059730
## 276 2016-10-01 Saturday 0.87034 3.059730
## 277 2016-10-31 Monday 3.08796 3.059730
## 278 2016-12-01 Thursday 1.11661 3.059730
## 279 2016-12-31 Saturday 2.59268 3.059730
## MAPE_auto.arima_Model
## 1 56.731 %
## 2 408.429 %
## 3 32.394 %
## 4 60.517 %
## 5 189.461 %
## 6 8.908 %
## 7 528.368 %
## 8 224.953 %
## 9 138.331 %
## 10 21.075 %
## 11 107.618 %
## 12 214.791 %
## 13 39.062 %
## 14 70.905 %
## 15 26.092 %
## 16 204.377 %
## 17 34.145 %
## 18 153.274 %
## 19 86.706 %
## 20 218.913 %
## 21 134.008 %
## 22 10.885 %
## 23 97.144 %
## 24 81.742 %
## 25 222.695 %
## 26 48.03 %
## 27 256.721 %
## 28 22.07 %
## 29 186.981 %
## 30 6.528 %
## 31 35.722 %
## 32 422.906 %
## 33 137.185 %
## 34 26.532 %
## 35 122.8 %
## 36 284.794 %
## 37 152.922 %
## 38 17.631 %
## 39 173.176 %
## 40 20.773 %
## 41 7.384 %
## 42 27.338 %
## 43 295.912 %
## 44 58.474 %
## 45 67.748 %
## 46 20.848 %
## 47 122.842 %
## 48 133.023 %
## 49 55.58 %
## 50 183.437 %
## 51 23.259 %
## 52 30.159 %
## 53 5.427 %
## 54 15.295 %
## 55 16.162 %
## 56 73.153 %
## 57 137.255 %
## 58 5.781 %
## 59 8.121 %
## 60 40.174 %
## 61 1005.833 %
## 62 123.344 %
## 63 22.753 %
## 64 109.667 %
## 65 50.247 %
## 66 786.236 %
## 67 249.228 %
## 68 21.975 %
## 69 137.297 %
## 70 11.08 %
## 71 95.406 %
## 72 134.087 %
## 73 656.217 %
## 74 286.798 %
## 75 188.452 %
## 76 65.225 %
## 77 55.537 %
## 78 224.785 %
## 79 554.978 %
## 80 119.671 %
## 81 137.297 %
## 82 21.461 %
## 83 115.919 %
## 84 45.756 %
## 85 3.943 %
## 86 29.572 %
## 87 50.133 %
## 88 32.77 %
## 89 53.984 %
## 90 101.89 %
## 91 199.027 %
## 92 438.211 %
## 93 137.297 %
## 94 21.763 %
## 95 115.846 %
## 96 295.457 %
## 97 253.64 %
## 98 117.289 %
## 99 3.191 %
## 100 52.147 %
## 101 41.643 %
## 102 203.68 %
## 103 169.798 %
## 104 503.795 %
## 105 135.719 %
## 106 25.544 %
## 107 97.473 %
## 108 276.337 %
## 109 52.876 %
## 110 283.66 %
## 111 11.285 %
## 112 90.494 %
## 113 18.245 %
## 114 10.269 %
## 115 591.152 %
## 116 439.94 %
## 117 116.759 %
## 118 1.038 %
## 119 4.787 %
## 120 288.621 %
## 121 1186.249 %
## 122 210.339 %
## 123 30.355 %
## 124 61.272 %
## 125 13.606 %
## 126 8.4 %
## 127 716.799 %
## 128 416.114 %
## 129 137.256 %
## 130 21.713 %
## 131 88.949 %
## 132 237.655 %
## 133 845.791 %
## 134 90.278 %
## 135 8.212 %
## 136 19.286 %
## 137 38.036 %
## 138 48.194 %
## 139 42.552 %
## 140 709.195 %
## 141 137.271 %
## 142 8.099 %
## 143 96.613 %
## 144 65.215 %
## 145 462.006 %
## 146 187.393 %
## 147 19.376 %
## 148 1.165 %
## 149 19.752 %
## 150 16.638 %
## 151 11.025 %
## 152 543.734 %
## 153 130.993 %
## 154 25.055 %
## 155 73.219 %
## 156 232.016 %
## 157 308.209 %
## 158 256.118 %
## 159 13.115 %
## 160 12.925 %
## 161 35.648 %
## 162 255.18 %
## 163 69.9 %
## 164 110.379 %
## 165 103.269 %
## 166 7.866 %
## 167 58.682 %
## 168 172.107 %
## 169 158.334 %
## 170 137.826 %
## 171 90.113 %
## 172 59.527 %
## 173 2.118 %
## 174 1040.966 %
## 175 2.555 %
## 176 443.276 %
## 177 137.297 %
## 178 25.055 %
## 179 115.898 %
## 180 2.395 %
## 181 1.603 %
## 182 917.434 %
## 183 61.142 %
## 184 25.994 %
## 185 45.537 %
## 186 1210.04 %
## 187 2038.475 %
## 188 457.734 %
## 189 137.297 %
## 190 25.055 %
## 191 101.627 %
## 192 284.877 %
## 193 695.789 %
## 194 409.777 %
## 195 45.368 %
## 196 36.466 %
## 197 7.662 %
## 198 1143.54 %
## 199 997.661 %
## 200 65.377 %
## 201 136.915 %
## 202 20.525 %
## 203 61.901 %
## 204 120.943 %
## 205 483.706 %
## 206 1075.96 %
## 207 21.101 %
## 208 70.885 %
## 209 29.681 %
## 210 10.472 %
## 211 2.297 %
## 212 374.031 %
## 213 85.863 %
## 214 22.304 %
## 215 111.065 %
## 216 149.719 %
## 217 647.114 %
## 218 32.566 %
## 219 139.256 %
## 220 32.056 %
## 221 12.538 %
## 222 66.033 %
## 223 287.288 %
## 224 573.312 %
## 225 98.535 %
## 226 23.373 %
## 227 110.331 %
## 228 182.795 %
## 229 144.307 %
## 230 41.083 %
## 231 44.424 %
## 232 68.874 %
## 233 46.614 %
## 234 307.475 %
## 235 212.409 %
## 236 264.206 %
## 237 96.98 %
## 238 23.377 %
## 239 66.321 %
## 240 198.606 %
## 241 858.592 %
## 242 420.823 %
## 243 54.538 %
## 244 31.969 %
## 245 73.995 %
## 246 34.364 %
## 247 274.283 %
## 248 17.264 %
## 249 96.485 %
## 250 7.212 %
## 251 104.361 %
## 252 223.894 %
## 253 251.382 %
## 254 15.82 %
## 255 327.534 %
## 256 22.655 %
## 257 12.862 %
## 258 16.124 %
## 259 94.506 %
## 260 51.392 %
## 261 57.03 %
## 262 23.085 %
## 263 60.198 %
## 264 148.643 %
## 265 56.125 %
## 266 9.926 %
## 267 18.298 %
## 268 41.794 %
## 269 166.911 %
## 270 8.581 %
## 271 42.985 %
## 272 466.197 %
## 273 96.56 %
## 274 10.251 %
## 275 107.286 %
## 276 251.556 %
## 277 0.914 %
## 278 174.02 %
## 279 18.014 %
data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_auto.arima=tail(forecasting_auto_arima$mean,N_forecasting_days))
## FD forecating_date forecasting_by_auto.arima
## 1 2017-01-31 Tuesday 3.05973
## 2 2017-03-03 Friday 3.05973
## 3 2017-03-31 Friday 3.05973
## 4 2017-05-01 Monday 3.05973
## 5 2017-05-31 Wednesday 3.05973
## 6 2017-07-01 Saturday 3.05973
## 7 2017-07-31 Monday 3.05973
## 8 2017-08-31 Thursday 3.05973
## 9 2017-10-01 Sunday 3.05973
## 10 2017-10-31 Tuesday 3.05973
## 11 2017-12-01 Friday 3.05973
## 12 2017-12-31 Sunday 3.05973
plot(forecasting_auto_arima)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)
graph4<-autoplot(forecasting_auto_arima,xlab = paste ("Time in ", frequency ,y_lab , sep=" "), ylab=y_lab)
graph4
## Error of forecasting
Error_auto.arima<-abs(testing_data-validation_forecast) # Absolute error of forecast (AEOF)
sqrt(sum(Error_auto.arima^2/validation_data_days))
## [1] 2.136216
REOF_A_auto.arima<-abs(((testing_data-validation_forecast)/testing_data)*100) #Relative error of forecast (divided by actual)(REOF_A)
REOF_F_auto.arima<-abs(((testing_data-validation_forecast)/validation_forecast)*100) #Relative error of forecast (divided by forecast)(REOF_F)
correlation_auto.arima<-cor(testing_data,validation_forecast, method = c("pearson")) # correlation coefficient between predicted and actual values
RMSE_auto.arima<-sqrt(sum((Error_auto.arima^2))/validation_data_days) # Root mean square forecast error
MAD_auto.arima<-abs((sum(testing_data-validation_forecast))/validation_data_days) # average forecast accuracy
AEOF_auto.arima<-c(Error_auto.arima)
REOF_auto.arima1<-c(paste(round(REOF_A_auto.arima,3),"%"))
REOF_auto.arima2<-c(paste(round(REOF_F_auto.arima,3),"%"))
data.frame(correlation_auto.arima,RMSE_auto.arima,MAPE_Mean_All,MAD_auto.arima) # analysis of Error by using Holt's linear model shows result of correlation ,MSE ,MPER
## correlation_auto.arima RMSE_auto.arima
## 1 -0.02744551 2.136216
## MAPE_Mean_All MAD_auto.arima
## 1 108.079 % MAPE 279 month Average Rainfall - (MM) In Egypt 0.6684562