# Trend Analysis
# 1. descriptive Analysis
# data Pre-processing
library(readxl)
livestock<-read_xlsx("~/badeed.xlsx")
# checking missing
any(is.na(livestock))
## [1] FALSE
# Declaring Time series
is.ts(livestock)
## [1] FALSE
# converting to Time series
livestock_ts<-ts(livestock$export, start=c(2012,1), end= c(2023,12), frequency = 12)
is.ts(livestock_ts)
## [1] TRUE
# descriptive
library(psych)
describe(livestock_ts)
## vars n mean sd median trimmed mad min max range
## X1 1 144 210617.2 288985.3 127549.5 139258 102883.5 8545 1593849 1585304
## skew kurtosis se
## X1 2.9 8.39 24082.11
library(TSstudio)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
ts.plot(livestock_ts, col="purple")
ts.plot(livestock_ts,xlab="Year",ylab="export data",main="Number of livestocks exported from somaliland
2012 -2023", col="magenta")
# Splitting Data sets
# Test + Train
library(TSstudio)
splitdata<-ts_split(livestock_ts,sample.out = 24)
train<-splitdata$train
train
## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2012 141262 288139 148713 118107 88395 210183 248533 179825 274086
## 2013 217617 168020 118326 73415 140644 226719 253359 230338 879725
## 2014 222553 142834 172535 143341 82279 280423 242883 281287 1235587
## 2015 165686 164504 220457 129177 152708 290257 165012 649526 1316095
## 2016 183440 146106 177433 144725 45098 321055 165424 1211406 478958
## 2017 41223 41631 33356 34467 51572 53659 57906 1272478 14342
## 2018 25262 34785 17444 32587 32391 103578 241896 842898 8920
## 2019 50292 68750 62667 84853 61328 39953 893599 454654 29675
## 2020 48775 64343 68689 55810 60841 172894 398729 86914 95209
## 2021 128551 137701 144710 152199 80101 560929 410966 58405 70708
## Oct Nov Dec
## 2012 1593849 48583 202758
## 2013 694393 37108 122131
## 2014 186338 201707 197694
## 2015 30330 42383 135040
## 2016 125284 69851 66149
## 2017 37742 27659 35877
## 2018 8545 17486 34139
## 2019 37499 59947 55214
## 2020 59416 63729 91785
## 2021 58866 98707 106891
test<-splitdata$test
test
## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2022 104728 172380 132969 165145 150354 791462 267923 57046 75544
## 2023 168787 139831 157523 97143 346864 1407445 117624 111897 97692
## Oct Nov Dec
## 2022 79359 112012 126548
## 2023 152918 123840 205914
ts.plot(train,xlab="Year",ylab="export",main="Number of livestocks exported from somaliland
2012 -2023", col="darkblue")
# diagnostic
library(tseries)
adf.test(train)
## Warning in adf.test(train): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: train
## Dickey-Fuller = -6.3548, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
pp.test(train)
## Warning in pp.test(train): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: train
## Dickey-Fuller Z(alpha) = -96.337, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
# How to check number of difference
library(forecast)
ndiffs(train)
## [1] 1
diff1<-diff(train)
par(mfrow=c(1,3))
ts.plot(livestock_ts,xlab="Year",ylab="export data",main="number of livestocks exported from somaliland before spliting
2012 -2023", col="purple")
ts.plot(train,xlab="Year",ylab="export",main="number of livestocks exported from somaliland
2012 -2023", col="darkred")
ts.plot(diff1, xlab="Years", ylab=" Number of exported", main="First differencing of train data")
#ACF and PACF
par(mfrow=c(1,4))
acf(train)
pacf(train)
acf(diff1)
pacf(diff1)
# Decomposition
plot(decompose(train,type = "additive"))
plot(decompose(train,type = "multiplicative"))
# III. Predictive Analytics
# Fit the specified SARIMA model
sarima1 <- arima(train, order=c(4, 1, 0), seasonal=list(order=c(2, 0, 0), period=12))
sarima2 <- arima(train, order=c(3, 1, 0), seasonal=list(order=c(2, 0, 0), period=12))
sarima3 <- arima(train, order=c(2, 1, 0), seasonal=list(order=c(2, 0, 0), period=12))
sarima4 <- arima(train, order=c(1, 1, 0), seasonal=list(order=c(2, 0, 0), period=12))
sarima5 <- arima(train, order=c(0, 1, 0), seasonal=list(order=c(2, 0, 0), period=12))
auto.arima(train, seasonal=TRUE)
## Series: train
## ARIMA(4,1,0)(2,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 sar1 sar2
## -0.9507 -0.7179 -0.4506 -0.1618 1.1869 -0.4504
## s.e. 0.1071 0.1318 0.1222 0.0907 0.0891 0.0914
##
## sigma^2 = 3.57e+10: log likelihood = -1621.47
## AIC=3256.93 AICc=3257.94 BIC=3276.39
AIC(sarima1,sarima2,sarima3,sarima4,sarima5)
## df AIC
## sarima1 7 3256.933
## sarima2 6 3258.062
## sarima3 5 3267.226
## sarima4 4 3280.582
## sarima5 3 3308.935
BIC(sarima1,sarima2,sarima3,sarima4,sarima5)
## df BIC
## sarima1 7 3276.387
## sarima2 6 3274.737
## sarima3 5 3281.122
## sarima4 4 3291.698
## sarima5 3 3317.273
#Fit Individual Models on Training Data:
# SARIMA model
fit_sarima_t <- auto.arima(train, seasonal = TRUE)
# ETS model
fit_ets <- ets(train)
# TBATS model
fit_tbats <- tbats(train)
# ANN model (Neural Network for Forecasting)
fit_ann <- nnetar(train)
#Forecast on Test Data
# Forecast using ETS
fc_ets <- forecast(fit_ets, h =24)
# Forecast using TBATS
fc_tbats <- forecast(fit_tbats, h =24)
# Forecast using ANN
fc_ann <- forecast(fit_ann, h = 24)
# Forecast using SARIMA
fc_sarima <- forecast(fit_sarima_t, h = 24)
#Evaluate Accuracy
# Function to calculate accuracy
evaluate_forecast <- function(forecast, actual) {
return(accuracy(forecast, actual))
}
# Evaluate individual models
accuracy_ets <- evaluate_forecast(fc_ets$mean, test)
accuracy_tbats <- evaluate_forecast(fc_tbats$mean, test)
accuracy_ann <- evaluate_forecast(fc_ann$mean, test)
accuracy_sarima <- evaluate_forecast(fc_sarima$mean, test)
# Print accuracies
print(list(
accuracy_ets = accuracy_ets,
accuracy_tbats = accuracy_tbats,
accuracy_ann = accuracy_ann,
accuracy_sarima = accuracy_sarima
))
## $accuracy_ets
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 55663.6 231622.8 114159.8 -4.45238 45.55983 -0.1385541 0.7770662
##
## $accuracy_tbats
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 119799.1 309359.9 129088.1 21.72566 35.25137 0.07261207 1.075958
##
## $accuracy_ann
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 23317.62 227730.9 106305.3 -8.001652 40.54603 -0.2311541 0.8533414
##
## $accuracy_sarima
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 53665.81 200193.6 86549.68 3.785126 29.97349 0.02046479 0.6577297
#Combine Forecasts for All Possible Hybrids:
# Function to combine forecasts
combine_forecasts <- function(..., weights = NULL) {
forecasts <- list(...)
if (is.null(weights)) {
weights <- rep(1/length(forecasts), length(forecasts))
}
combined <- Reduce(`+`, Map(`*`, forecasts, weights))
return(combined)
}
# Combine pairs
combined_ets_tbats <- combine_forecasts(fc_ets$mean, fc_tbats$mean)
combined_ets_ann <- combine_forecasts(fc_ets$mean, fc_ann$mean)
combined_ets_sarima <- combine_forecasts(fc_ets$mean, fc_sarima$mean)
combined_tbats_ann <- combine_forecasts(fc_tbats$mean, fc_ann$mean)
combined_tbats_sarima <- combine_forecasts(fc_tbats$mean, fc_sarima$mean)
combined_ann_sarima <- combine_forecasts(fc_ann$mean, fc_sarima$mean)
# Combine triples
combined_ets_tbats_ann <- combine_forecasts(fc_ets$mean, fc_tbats$mean, fc_ann$mean)
combined_ets_tbats_sarima <- combine_forecasts(fc_ets$mean, fc_tbats$mean, fc_sarima$mean)
combined_ets_ann_sarima <- combine_forecasts(fc_ets$mean, fc_ann$mean, fc_sarima$mean)
combined_tbats_ann_sarima <- combine_forecasts(fc_tbats$mean, fc_ann$mean, fc_sarima$mean)
# Combine all
combined_all <- combine_forecasts(fc_ets$mean, fc_tbats$mean, fc_ann$mean, fc_sarima$mean)
#Evaluate Each Hybrid:
# Function to calculate accuracy
evaluate_forecast <- function(forecast, actual) {
return(accuracy(forecast, actual))
}
# Evaluate pairs
accuracy_ets_tbats <- evaluate_forecast(combined_ets_tbats, test)
accuracy_ets_ann <- evaluate_forecast(combined_ets_ann, test)
accuracy_ets_sarima <- evaluate_forecast(combined_ets_sarima, test)
accuracy_tbats_ann <- evaluate_forecast(combined_tbats_ann, test)
accuracy_tbats_sarima <- evaluate_forecast(combined_tbats_sarima, test)
accuracy_ann_sarima <- evaluate_forecast(combined_ann_sarima, test)
# Evaluate triples
accuracy_ets_tbats_ann <- evaluate_forecast(combined_ets_tbats_ann, test)