New techniques for forecasting SAARC Milk production by using BATS,and TBATS Models
By
Dr:Pradeep Mishra
Department of Mathematics & Statistics
Jawaharlal Nehru Krishi Vishwavidyalaya University, Jabalpur, India
Abotaleb mostafa
Department of Electrical Engineering and Computer Science
South ural state university, Chelyabinsk, Russian federation
References
Time Series Investigation of Milk Production in Major States of India Using ARIMA Modeling
Modeling and Forecasting of Milk Production in Chhattisgarh and India
# Imports
# Imports
library(fpp2)## Warning: package 'fpp2' was built under R version 4.0.3## 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## Warning: package 'ggplot2' was built under R version 4.0.3## Warning: package 'forecast' was built under R version 4.0.3## library(forecast)
library(ggplot2)
library("readxl")## Warning: package 'readxl' was built under R version 4.0.3library(moments)## Warning: package 'moments' was built under R version 4.0.3library(forecast)
require(forecast)
require(tseries)## Loading required package: tseries## Warning: package 'tseries' was built under R version 4.0.3require(markovchain)## Loading required package: markovchain## Warning: package 'markovchain' was built under R version 4.0.3## Package: markovchain
## Version: 0.8.5-3
## Date: 2020-12-03
## BugReport: https://github.com/spedygiorgio/markovchain/issuesrequire(data.table)## Loading required package: data.tableFull_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in Pakistan" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$Pakistan
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4209000 5773750 10380500 12627086 19084250 28109000sd(original_data) # calculate standard deviation## [1] 7540331skewness(original_data) # calculate Cofficient of skewness## [1] 0.5334464kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 1.880945rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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
## Training set 25338.39 179751.4 120493.9 0.2595034 1.455225 0.3178252
## ACF1
## Training set -0.01911074# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 0.9868946
## Beta: 0.3584389
## Damping Parameter: 1
## Gamma-1 Values: -0.005289493
## Gamma-2 Values: 0.009289828
##
## Seed States:
## [,1]
## [1,] 4730320.9197
## [2,] 144576.0120
## [3,] 1161.5291
## [4,] -23156.1313
## [5,] 725.8252
## [6,] -25038.9612
##
## Sigma: 179751.4
## AIC: 1454.655plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in Pakistan"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Pakistan"paste(MAPE_Mean_All,"%")## [1] "1.746 % MAPE 7 years forecasting SAARC Milk production in Pakistan %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Pakistan"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 23652000 23502822
## 2 2013-12-31 вторник 24370000 24121649
## 3 2014-12-31 среда 25001000 24779868
## 4 2015-12-31 четверг 25744000 25331920
## 5 2016-12-31 суббота 26510000 25948036
## 6 2017-12-31 воскресенье 27298000 26600471
## 7 2018-12-31 понедельник 28109000 27149451
## MAPE_TBATS_Model
## 1 0.631 %
## 2 1.019 %
## 3 0.884 %
## 4 1.601 %
## 5 2.12 %
## 6 2.555 %
## 7 3.414 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 27768278
## 2 2020-12-31 четверг 28426497
## 3 2021-12-31 пятница 28978550
## 4 2022-12-31 суббота 29594665
## 5 2023-12-31 воскресенье 30247100
## 6 2024-12-31 вторник 30796080
## 7 2025-12-31 среда 31414907plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 ACF1
## Training set 16437.32 166425 99196.88 0.2896213 1.073816 0.2616504 -0.03808877# Print Model Parameters
model_bats## BATS(0.544, {0,0}, 0.998, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0.544361
## Alpha: 0.9867143
## Beta: 0.4119543
## Damping Parameter: 0.998458
##
## Seed States:
## [,1]
## [1,] 7707.8541597
## [2,] 0.0208483
## attr(,"lambda")
## [1] 0.5443614
##
## Sigma: 102.5477
## AIC: 1430.339plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in Pakistan"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 7 days in bats Model for ==> forecasting SAARC Milk production in Pakistan"paste(MAPE_Mean_All,"%")## [1] "0.751 % MAPE 7 years forecasting SAARC Milk production in Pakistan %"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 7 days in bats Model for ==> forecasting SAARC Milk production in Pakistan"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 23652000 23604177
## 2 2013-12-31 вторник 24370000 24261379
## 3 2014-12-31 среда 25001000 24925762
## 4 2015-12-31 четверг 25744000 25597270
## 5 2016-12-31 суббота 26510000 26275848
## 6 2017-12-31 воскресенье 27298000 26961440
## 7 2018-12-31 понедельник 28109000 27653994
## MAPE_bats_Model
## 1 0.202 %
## 2 0.446 %
## 3 0.301 %
## 4 0.57 %
## 5 0.883 %
## 6 1.233 %
## 7 1.619 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 28353455
## 2 2020-12-31 четверг 29059771
## 3 2021-12-31 пятница 29772889
## 4 2022-12-31 суббота 30492759
## 5 2023-12-31 воскресенье 31219327
## 6 2024-12-31 вторник 31952545
## 7 2025-12-31 среда 32692361plot(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 , sep=" "),ylab = y_lab)
graph1
Full_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in `China, mainland`" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$`China, mainland`
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 917000 1166250 1885000 1957335 2737500 3100000sd(original_data) # calculate standard deviation## [1] 791574skewness(original_data) # calculate Cofficient of skewness## [1] 0.09956085kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 1.450104rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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
## Training set 4547.584 38491.59 25951.82 0.4323523 1.426527 0.5495366
## ACF1
## Training set -0.03622683# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 0.9312581
## Beta: 0.1720789
## Damping Parameter: 1
## Gamma-1 Values: -0.009681016
## Gamma-2 Values: 0.008120177
##
## Seed States:
## [,1]
## [1,] 885291.763
## [2,] 2352.081
## [3,] -1283.867
## [4,] 6395.803
## [5,] -1743.830
## [6,] -1847.356
##
## Sigma: 38491.59
## AIC: 1297.459plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in `China, mainland`"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in `China, mainland`"paste(MAPE_Mean_All,"%")## [1] "6.772 % MAPE 7 years forecasting SAARC Milk production in `China, mainland` %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in `China, mainland`"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 3080000 3109069
## 2 2013-12-31 вторник 3050000 3141529
## 3 2014-12-31 среда 3100000 3181382
## 4 2015-12-31 четверг 2990666 3230758
## 5 2016-12-31 суббота 3005201 3262918
## 6 2017-12-31 воскресенье 2946374 3307867
## 7 2018-12-31 понедельник 3003323 3362639
## MAPE_TBATS_Model
## 1 0.944 %
## 2 3.001 %
## 3 2.625 %
## 4 8.028 %
## 5 8.576 %
## 6 12.269 %
## 7 11.964 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 3395099
## 2 2020-12-31 четверг 3434953
## 3 2021-12-31 пятница 3484329
## 4 2022-12-31 суббота 3516488
## 5 2023-12-31 воскресенье 3561438
## 6 2024-12-31 вторник 3616210
## 7 2025-12-31 среда 3648670plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 -1935.931 40750.93 23809.33 0.06616459 1.197824 0.5041687
## ACF1
## Training set -0.132177# Print Model Parameters
model_bats## BATS(0, {0,0}, 1, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0
## Alpha: 0.9732246
## Beta: 0.2432013
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 13.686930648
## [2,] 0.003080326
## attr(,"lambda")
## [1] 3.369018e-07
##
## Sigma: 0.01756693
## AIC: 1259.558plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in `China, mainland`"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 7 days in bats Model for ==> forecasting SAARC Milk production in `China, mainland`"paste(MAPE_Mean_All,"%")## [1] "6.443 % MAPE 7 years forecasting SAARC Milk production in `China, mainland` %"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 7 days in bats Model for ==> forecasting SAARC Milk production in `China, mainland`"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 3080000 3091986
## 2 2013-12-31 вторник 3050000 3133081
## 3 2014-12-31 среда 3100000 3174723
## 4 2015-12-31 четверг 2990666 3216919
## 5 2016-12-31 суббота 3005201 3259675
## 6 2017-12-31 воскресенье 2946374 3303000
## 7 2018-12-31 понедельник 3003323 3346900
## MAPE_bats_Model
## 1 0.389 %
## 2 2.724 %
## 3 2.41 %
## 4 7.565 %
## 5 8.468 %
## 6 12.104 %
## 7 11.44 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 3391384
## 2 2020-12-31 четверг 3436459
## 3 2021-12-31 пятница 3482133
## 4 2022-12-31 суббота 3528414
## 5 2023-12-31 воскресенье 3575311
## 6 2024-12-31 вторник 3622830
## 7 2025-12-31 среда 3670982plot(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 , sep=" "),ylab = y_lab)
graph1
Full_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in India" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$India
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10929000 14313750 28547500 34725088 49628250 91817140sd(original_data) # calculate standard deviation## [1] 22691830skewness(original_data) # calculate Cofficient of skewness## [1] 0.8041301kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 2.578639rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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 163237.7 607563 493032.2 0.7495706 2.317169 0.4204822 -0.04750128# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 1.135479
## Beta: 0.3232109
## Damping Parameter: 1
## Gamma-1 Values: -0.0005794413
## Gamma-2 Values: 0.0004366506
##
## Seed States:
## [,1]
## [1,] 10636498.27
## [2,] -333350.95
## [3,] -13434.89
## [4,] -232143.56
## [5,] -154919.57
## [6,] 55215.47
##
## Sigma: 607563
## AIC: 1578.879plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in India"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in India"paste(MAPE_Mean_All,"%")## [1] "4.051 % MAPE 7 years forecasting SAARC Milk production in India %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in India"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 67675432 67682808
## 2 2013-12-31 вторник 70442617 70563893
## 3 2014-12-31 среда 74709900 72816364
## 4 2015-12-31 четверг 76459000 74730907
## 5 2016-12-31 суббота 81266300 77355807
## 6 2017-12-31 воскресенье 86261680 79632439
## 7 2018-12-31 понедельник 91817140 81827328
## MAPE_TBATS_Model
## 1 0.011 %
## 2 0.172 %
## 3 2.535 %
## 4 2.26 %
## 5 4.812 %
## 6 7.685 %
## 7 10.88 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 84708413
## 2 2020-12-31 четверг 86960884
## 3 2021-12-31 пятница 88875426
## 4 2022-12-31 суббота 91500327
## 5 2023-12-31 воскресенье 93776959
## 6 2024-12-31 вторник 95971847
## 7 2025-12-31 среда 98852933plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 ACF1
## Training set 169496.8 689136.3 504977.5 0.7031322 2.22926 0.4306697 -0.02727572# Print Model Parameters
model_bats## BATS(0.77, {0,0}, 1, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0.770463
## Alpha: 1.115592
## Beta: 0.2071654
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 334769.596180
## [2,] -1.312595
## attr(,"lambda")
## [1] 0.7704627
##
## Sigma: 13822.66
## AIC: 1581.164plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in India"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 7 days in bats Model for ==> forecasting SAARC Milk production in India"paste(MAPE_Mean_All,"%")## [1] "4.546 % MAPE 7 years forecasting SAARC Milk production in India %"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 7 days in bats Model for ==> forecasting SAARC Milk production in India"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 67675432 67687320
## 2 2013-12-31 вторник 70442617 69935932
## 3 2014-12-31 среда 74709900 72201265
## 4 2015-12-31 четверг 76459000 74483033
## 5 2016-12-31 суббота 81266300 76780960
## 6 2017-12-31 воскресенье 86261680 79094785
## 7 2018-12-31 понедельник 91817140 81424253
## MAPE_bats_Model
## 1 0.018 %
## 2 0.719 %
## 3 3.358 %
## 4 2.584 %
## 5 5.519 %
## 6 8.308 %
## 7 11.319 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 83769121
## 2 2020-12-31 четверг 86129154
## 3 2021-12-31 пятница 88504126
## 4 2022-12-31 суббота 90893819
## 5 2023-12-31 воскресенье 93298021
## 6 2024-12-31 вторник 95716530
## 7 2025-12-31 среда 98149149plot(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 , sep=" "),ylab = y_lab)
graph1
Full_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in Nepal" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$Nepal
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 340000 482250 596618 678285 856086 1338277sd(original_data) # calculate standard deviation## [1] 279431.2skewness(original_data) # calculate Cofficient of skewness## [1] 0.7663647kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 2.412837rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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 2140.325 12074.91 8539.213 0.2344735 1.60696 0.5074679 -0.04518789# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 1.157697
## Beta: 0.209391
## Damping Parameter: 1
## Gamma-1 Values: -0.0003674468
## Gamma-2 Values: -0.0004999581
##
## Seed States:
## [,1]
## [1,] 299689.1851
## [2,] 10651.6178
## [3,] -935.5607
## [4,] 779.2784
## [5,] 696.6762
## [6,] 1923.7177
##
## Sigma: 12074.91
## AIC: 1179.209plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in Nepal"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Nepal"paste(MAPE_Mean_All,"%")## [1] "3.369 % MAPE 7 years forecasting SAARC Milk production in Nepal %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Nepal"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 1153838 1147760
## 2 2013-12-31 вторник 1188433 1180643
## 3 2014-12-31 среда 1167773 1209898
## 4 2015-12-31 четверг 1167154 1246415
## 5 2016-12-31 суббота 1210441 1281470
## 6 2017-12-31 воскресенье 1245954 1312594
## 7 2018-12-31 понедельник 1338277 1348808
## MAPE_TBATS_Model
## 1 0.527 %
## 2 0.655 %
## 3 3.607 %
## 4 6.791 %
## 5 5.868 %
## 6 5.349 %
## 7 0.787 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 1381691
## 2 2020-12-31 четверг 1410947
## 3 2021-12-31 пятница 1447463
## 4 2022-12-31 суббота 1482518
## 5 2023-12-31 воскресенье 1513642
## 6 2024-12-31 вторник 1549856
## 7 2025-12-31 среда 1582740plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 357.1187 11267.57 8026.917 -0.06524403 1.438529 0.4770232
## ACF1
## Training set 0.04559581# Print Model Parameters
model_bats## BATS(0.012, {0,0}, 1, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0.012029
## Alpha: 1.19226
## Beta: -0.01190352
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 13.6679832
## [2,] 0.0291863
## attr(,"lambda")
## [1] 0.01202911
##
## Sigma: 0.02441833
## AIC: 1167.955plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in Nepal"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 7 days in bats Model for ==> forecasting SAARC Milk production in Nepal"paste(MAPE_Mean_All,"%")## [1] "2.787 % MAPE 7 years forecasting SAARC Milk production in Nepal %"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 7 days in bats Model for ==> forecasting SAARC Milk production in Nepal"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 1153838 1140165
## 2 2013-12-31 вторник 1188433 1168947
## 3 2014-12-31 среда 1167773 1198446
## 4 2015-12-31 четверг 1167154 1228681
## 5 2016-12-31 суббота 1210441 1259669
## 6 2017-12-31 воскресенье 1245954 1291429
## 7 2018-12-31 понедельник 1338277 1323980
## MAPE_bats_Model
## 1 1.185 %
## 2 1.64 %
## 3 2.627 %
## 4 5.272 %
## 5 4.067 %
## 6 3.65 %
## 7 1.068 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 1357341
## 2 2020-12-31 четверг 1391533
## 3 2021-12-31 пятница 1426575
## 4 2022-12-31 суббота 1462488
## 5 2023-12-31 воскресенье 1499295
## 6 2024-12-31 вторник 1537017
## 7 2025-12-31 среда 1575676plot(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 , sep=" "),ylab = y_lab)
graph1
Full_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in `Sri Lanka`" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$`Sri Lanka`
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18320 29411 42535 45636 62135 85914sd(original_data) # calculate standard deviation## [1] 17528.07skewness(original_data) # calculate Cofficient of skewness## [1] 0.3694499kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 1.882833rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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 1305.45 8719.842 5599.515 1.577394 14.92159 1.126319 -0.03532745# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 0.9971949
## Beta: 0.01798248
## Damping Parameter: 1
## Gamma-1 Values: -0.01634612
## Gamma-2 Values: 0.008369212
##
## Seed States:
## [,1]
## [1,] 27277.6818
## [2,] -1678.4463
## [3,] 1781.1336
## [4,] 220.3718
## [5,] -1730.0898
## [6,] -1766.6832
##
## Sigma: 8719.842
## AIC: 1146.005plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in `Sri Lanka`"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in `Sri Lanka`"paste(MAPE_Mean_All,"%")## [1] "30.458 % MAPE 7 years forecasting SAARC Milk production in `Sri Lanka` %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in `Sri Lanka`"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 61710 44335.45
## 2 2013-12-31 вторник 54060 44215.58
## 3 2014-12-31 среда 45854 50277.39
## 4 2015-12-31 четверг 36118 48155.16
## 5 2016-12-31 суббота 66128 42665.94
## 6 2017-12-31 воскресенье 68591 43464.42
## 7 2018-12-31 понедельник 85914 41448.19
## MAPE_TBATS_Model
## 1 28.155 %
## 2 18.21 %
## 3 9.647 %
## 4 33.327 %
## 5 35.48 %
## 6 36.632 %
## 7 51.756 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 41328.32
## 2 2020-12-31 четверг 47390.14
## 3 2021-12-31 пятница 45267.90
## 4 2022-12-31 суббота 39778.68
## 5 2023-12-31 воскресенье 40577.17
## 6 2024-12-31 вторник 38560.94
## 7 2025-12-31 среда 38441.06plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 ACF1
## Training set 214.8835 9087.791 5102.434 -2.24465 13.93984 1.026333 -0.03144049# Print Model Parameters
model_bats## BATS(0.216, {0,0}, -, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0.216358
## Alpha: 0.9829001
##
## Seed States:
## [,1]
## [1,] 40.15351
## attr(,"lambda")
## [1] 0.2163576
##
## Sigma: 2.104403
## AIC: 1130.562plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in `Sri Lanka`"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 7 days in bats Model for ==> forecasting SAARC Milk production in `Sri Lanka`"paste(MAPE_Mean_All,"%")## [1] "25.282 % MAPE 7 years forecasting SAARC Milk production in `Sri Lanka` %"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 7 days in bats Model for ==> forecasting SAARC Milk production in `Sri Lanka`"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 61710 46339.52
## 2 2013-12-31 вторник 54060 46339.52
## 3 2014-12-31 среда 45854 46339.52
## 4 2015-12-31 четверг 36118 46339.52
## 5 2016-12-31 суббота 66128 46339.52
## 6 2017-12-31 воскресенье 68591 46339.52
## 7 2018-12-31 понедельник 85914 46339.52
## MAPE_bats_Model
## 1 24.908 %
## 2 14.281 %
## 3 1.059 %
## 4 28.3 %
## 5 29.925 %
## 6 32.441 %
## 7 46.063 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 46339.52
## 2 2020-12-31 четверг 46339.52
## 3 2021-12-31 пятница 46339.52
## 4 2022-12-31 суббота 46339.52
## 5 2023-12-31 воскресенье 46339.52
## 6 2024-12-31 вторник 46339.52
## 7 2025-12-31 среда 46339.52plot(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 , sep=" "),ylab = y_lab)
graph1
Full_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in Bangladesh" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$Bangladesh
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13090 20100 23377 24770 27750 39000sd(original_data) # calculate standard deviation## [1] 6580.095skewness(original_data) # calculate Cofficient of skewness## [1] 0.6500564kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 2.441656rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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
## Training set -45.62149 1570.798 1115.194 -0.758213 5.649507 1.09059
## ACF1
## Training set -0.003429975# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 0.9407896
## Beta: 0.1002629
## Damping Parameter: 1
## Gamma-1 Values: -0.002555942
## Gamma-2 Values: 0.01155586
##
## Seed States:
## [,1]
## [1,] 13161.3401
## [2,] 1130.1226
## [3,] -280.0889
## [4,] 101.4288
## [5,] 334.2858
## [6,] 186.5834
##
## Sigma: 1570.798
## AIC: 971.1757plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in Bangladesh"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Bangladesh"paste(MAPE_Mean_All,"%")## [1] "12.293 % MAPE 7 years forecasting SAARC Milk production in Bangladesh %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Bangladesh"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 38000 38218.04
## 2 2013-12-31 вторник 39000 38693.99
## 3 2014-12-31 среда 35173 39012.80
## 4 2015-12-31 четверг 35303 40355.38
## 5 2016-12-31 суббота 35432 41683.17
## 6 2017-12-31 воскресенье 35562 42555.17
## 7 2018-12-31 понедельник 35691 43599.09
## MAPE_TBATS_Model
## 1 0.574 %
## 2 0.785 %
## 3 10.917 %
## 4 14.311 %
## 5 17.643 %
## 6 19.665 %
## 7 22.157 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 44075.04
## 2 2020-12-31 четверг 44393.85
## 3 2021-12-31 пятница 45736.42
## 4 2022-12-31 суббота 47064.22
## 5 2023-12-31 воскресенье 47936.22
## 6 2024-12-31 вторник 48980.14
## 7 2025-12-31 среда 49456.09plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 ACF1
## Training set 340.776 1674.575 1064.698 0.8069479 5.202474 1.041208 -0.04360055# Print Model Parameters
model_bats## BATS(1, {0,0}, -, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Alpha: 1.031937
##
## Seed States:
## [,1]
## [1,] 19302.71
##
## Sigma: 1674.575
## AIC: 961.7012plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in Bangladesh"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 7 days in bats Model for ==> forecasting SAARC Milk production in Bangladesh"paste(MAPE_Mean_All,"%")## [1] "4.573 % MAPE 7 years forecasting SAARC Milk production in Bangladesh %"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 7 days in bats Model for ==> forecasting SAARC Milk production in Bangladesh"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 38000 37237.33
## 2 2013-12-31 вторник 39000 37237.33
## 3 2014-12-31 среда 35173 37237.33
## 4 2015-12-31 четверг 35303 37237.33
## 5 2016-12-31 суббота 35432 37237.33
## 6 2017-12-31 воскресенье 35562 37237.33
## 7 2018-12-31 понедельник 35691 37237.33
## MAPE_bats_Model
## 1 2.007 %
## 2 4.52 %
## 3 5.869 %
## 4 5.479 %
## 5 5.095 %
## 6 4.711 %
## 7 4.333 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 37237.33
## 2 2020-12-31 четверг 37237.33
## 3 2021-12-31 пятница 37237.33
## 4 2022-12-31 суббота 37237.33
## 5 2023-12-31 воскресенье 37237.33
## 6 2024-12-31 вторник 37237.33
## 7 2025-12-31 среда 37237.33plot(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 , sep=" "),ylab = y_lab)
graph1
Full_original_data<-read_excel("F:/Phd/ALL Russia Analysis/SAARC milk production final.xlsx")
y_lab<- "forecasting SAARC Milk production in Myanmar" # input name of data
Actual_date_interval <- c("1961/12/31","2018/12/31")
Forecast_date_interval <- c("2019/12/31","2025/12/31")
validation_data_days <-7
frequency<-"years"# Data Preparation & calculate some of statistics measures
original_data<-Full_original_data$Myanmar
summary(original_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18180 49035 93736 103086 159709 305631sd(original_data) # calculate standard deviation## [1] 74419.61skewness(original_data) # calculate Cofficient of skewness## [1] 0.9504659kurtosis(original_data) # calculate Cofficient of kurtosis## [1] 3.201248rows <- 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)
#plot
autoplot(data_series ,xlab=paste ("Time in ", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))
## 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
## Training set 1682.294 10178.43 6496.924 1.690186 8.622555 0.9309968
## ACF1
## Training set -0.006797351# Print Model Parameters
model_TBATS## TBATS(1, {0,0}, 1, {<6,2>})
##
## Call: NULL
##
## Parameters
## Alpha: 0.6496769
## Beta: 0.2621229
## Damping Parameter: 1
## Gamma-1 Values: 0.003815186
## Gamma-2 Values: -0.004483384
##
## Seed States:
## [,1]
## [1,] 12559.6504
## [2,] -861.8135
## [3,] 653.3233
## [4,] 1416.4450
## [5,] 3195.4670
## [6,] 1277.1974
##
## Sigma: 10178.43
## AIC: 1161.782plot(model_TBATS,main =paste(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 7 years by using TBATS Model for ==> forecasting SAARC Milk production in Myanmar"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Myanmar"paste(MAPE_Mean_All,"%")## [1] "115.646 % MAPE 7 years forecasting SAARC Milk production in Myanmar %"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 7 days in TBATS Model for ==> forecasting SAARC Milk production in Myanmar"data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)## date_TBATS validation_data_by_name actual_data forecasting_TBATS
## 1 2012-12-31 понедельник 171184 333191.1
## 2 2013-12-31 вторник 175526 351438.5
## 3 2014-12-31 среда 179751 370668.5
## 4 2015-12-31 четверг 184142 398602.2
## 5 2016-12-31 суббота 188490 421902.8
## 6 2017-12-31 воскресенье 192134 440604.3
## 7 2018-12-31 понедельник 193841 462956.4
## MAPE_TBATS_Model
## 1 94.639 %
## 2 100.22 %
## 3 106.212 %
## 4 116.465 %
## 5 123.833 %
## 6 129.321 %
## 7 138.833 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days))## FD forecating_date forecasting_by_TBATS
## 1 2019-12-31 вторник 481203.8
## 2 2020-12-31 четверг 500433.8
## 3 2021-12-31 пятница 528367.5
## 4 2022-12-31 суббота 551668.0
## 5 2023-12-31 воскресенье 570369.6
## 6 2024-12-31 вторник 592721.7
## 7 2025-12-31 среда 610969.0plot(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 , sep=" "),ylab = y_lab)
graph2
##bats model
# Data Modeling
data_series<-ts(training_data)
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 ACF1
## Training set 547.6818 10763.15 6829.284 -1.85141 7.643323 0.9786233 0.1023617# Print Model Parameters
model_bats## BATS(0.339, {0,0}, 1, -)
##
## Call: bats(y = data_series)
##
## Parameters
## Lambda: 0.338569
## Alpha: 0.82674
## Beta: -0.008284143
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 72.08607
## [2,] 2.79474
## attr(,"lambda")
## [1] 0.338569
##
## Sigma: 5.03121
## AIC: 1127.109plot(model_bats,main =paste(y_lab))
# 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 7 years by using bats Model for ==> forecasting SAARC Milk production in Myanmar"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 7 days in bats Model for ==> forecasting SAARC Milk production in Myanmar"paste(MAPE_Mean_All,"%")## [1] "95.673 % MAPE 7 years forecasting SAARC Milk production in Myanmar %"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 7 days in bats Model for ==> forecasting SAARC Milk production in Myanmar"data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)## date_bats validation_data_by_name actual_data forecasting_bats
## 1 2012-12-31 понедельник 171184 319048.1
## 2 2013-12-31 вторник 175526 332029.9
## 3 2014-12-31 среда 179751 345356.3
## 4 2015-12-31 четверг 184142 359031.8
## 5 2016-12-31 суббота 188490 373060.7
## 6 2017-12-31 воскресенье 192134 387447.5
## 7 2018-12-31 понедельник 193841 402196.5
## MAPE_bats_Model
## 1 86.377 %
## 2 89.163 %
## 3 92.13 %
## 4 94.975 %
## 5 97.921 %
## 6 101.655 %
## 7 107.488 %data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days))## FD forecating_date forecasting_by_bats
## 1 2019-12-31 вторник 417312.3
## 2 2020-12-31 четверг 432799.0
## 3 2021-12-31 пятница 448661.3
## 4 2022-12-31 суббота 464903.3
## 5 2023-12-31 воскресенье 481529.6
## 6 2024-12-31 вторник 498544.4
## 7 2025-12-31 среда 515952.3plot(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 , sep=" "),ylab = y_lab)
graph1