setwd("C:/Users/favou/OneDrive/Documents")
require(tseries)
## Loading required package: tseries
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
require(forecast)
## Loading required package: forecast
library(readxl)
Data<-read_excel("icammda_pricemeaslesdata22.xls")
Data
## # A tibble: 84 × 4
## `Time(Month)` Year Month `Jere LGA, Borno State`
## <dbl> <dbl> <chr> <dbl>
## 1 1 2016 Jan 12
## 2 2 2016 Feb 7
## 3 3 2016 Mar 1
## 4 4 2016 Apr 17
## 5 5 2016 May 50
## 6 6 2016 Jun 2
## 7 7 2016 Jul 8
## 8 8 2016 Aug 166
## 9 9 2016 Sept 131
## 10 10 2016 Oct 59
## # ℹ 74 more rows
head(Data)
## # A tibble: 6 × 4
## `Time(Month)` Year Month `Jere LGA, Borno State`
## <dbl> <dbl> <chr> <dbl>
## 1 1 2016 Jan 12
## 2 2 2016 Feb 7
## 3 3 2016 Mar 1
## 4 4 2016 Apr 17
## 5 5 2016 May 50
## 6 6 2016 Jun 2
tail(Data)
## # A tibble: 6 × 4
## `Time(Month)` Year Month `Jere LGA, Borno State`
## <dbl> <dbl> <chr> <dbl>
## 1 79 2022 Jul 138
## 2 80 2022 Aug 158
## 3 81 2022 Sept 1
## 4 82 2022 Oct 0
## 5 83 2022 Nov 2
## 6 84 2022 Dec 0
Data$`Jere LGA, Borno State`
## [1] 12 7 1 17 50 2 8 166 131 59 39 9 0 4 0 1 1 0 0
## [20] 0 1 0 2 6 0 2 2 1 2 3 0 0 2 0 11 244 343 474
## [39] 658 468 56 34 0 3 0 1 1 1 1 0 41 30 546 275 310 58 262
## [58] 87 53 178 26 205 123 138 174 352 242 271 46 79 90 73 157 84 233 134
## [77] 162 83 138 158 1 0 2 0
Data_2<-Data$`Jere LGA, Borno State`
#date
dates_data<- seq(as.Date("2016/1/1"), as.Date("2022/12/1"), by = "month")
dates_data
## [1] "2016-01-01" "2016-02-01" "2016-03-01" "2016-04-01" "2016-05-01"
## [6] "2016-06-01" "2016-07-01" "2016-08-01" "2016-09-01" "2016-10-01"
## [11] "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01"
## [16] "2017-04-01" "2017-05-01" "2017-06-01" "2017-07-01" "2017-08-01"
## [21] "2017-09-01" "2017-10-01" "2017-11-01" "2017-12-01" "2018-01-01"
## [26] "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01" "2018-06-01"
## [31] "2018-07-01" "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01"
## [36] "2018-12-01" "2019-01-01" "2019-02-01" "2019-03-01" "2019-04-01"
## [41] "2019-05-01" "2019-06-01" "2019-07-01" "2019-08-01" "2019-09-01"
## [46] "2019-10-01" "2019-11-01" "2019-12-01" "2020-01-01" "2020-02-01"
## [51] "2020-03-01" "2020-04-01" "2020-05-01" "2020-06-01" "2020-07-01"
## [56] "2020-08-01" "2020-09-01" "2020-10-01" "2020-11-01" "2020-12-01"
## [61] "2021-01-01" "2021-02-01" "2021-03-01" "2021-04-01" "2021-05-01"
## [66] "2021-06-01" "2021-07-01" "2021-08-01" "2021-09-01" "2021-10-01"
## [71] "2021-11-01" "2021-12-01" "2022-01-01" "2022-02-01" "2022-03-01"
## [76] "2022-04-01" "2022-05-01" "2022-06-01" "2022-07-01" "2022-08-01"
## [81] "2022-09-01" "2022-10-01" "2022-11-01" "2022-12-01"
plot(`Jere LGA, Borno State`~dates_data, data=Data,col='red', xlab = "Time (months)", type='l', lwd=3, main="Time Series plot of Measles cases in Nigeria")
acf(Data$`Jere LGA, Borno State`
, main="ACF of Measles cases in Nigeria", lag.max=12)
pacf(Data$`Jere LGA, Borno State`, main="PACF of Measles cases in Nigeria", lag.max=12)
#Test
adf.test(Data$`Jere LGA, Borno State`)
##
## Augmented Dickey-Fuller Test
##
## data: Data$`Jere LGA, Borno State`
## Dickey-Fuller = -3.8246, Lag order = 4, p-value = 0.02172
## alternative hypothesis: stationary
#Measles difference 1
Measles_case_d1<-diff(Data$`Jere LGA, Borno State`)
adf.test(Measles_case_d1)
## Warning in adf.test(Measles_case_d1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: Measles_case_d1
## Dickey-Fuller = -5.0428, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
#plot
par(mfrow=c(1,1))
plot(Measles_case_d1, main="Differenced Measles Cases", type='l')
acf(Measles_case_d1, main="ACF of Measles cases in Jere LGA, Borno State", lag.max=12, ylim=c(-1,1))
pacf(Measles_case_d1, main="PACF of Measles cases in Jere LGA, Borno State", lag.max=12, ylim=c(-1,1))
#ARIMA fit
ts_measles<-ts(c(Data_2), start = c(2016), frequency = 12)
ts_measles
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2016 12 7 1 17 50 2 8 166 131 59 39 9
## 2017 0 4 0 1 1 0 0 0 1 0 2 6
## 2018 0 2 2 1 2 3 0 0 2 0 11 244
## 2019 343 474 658 468 56 34 0 3 0 1 1 1
## 2020 1 0 41 30 546 275 310 58 262 87 53 178
## 2021 26 205 123 138 174 352 242 271 46 79 90 73
## 2022 157 84 233 134 162 83 138 158 1 0 2 0
# to know various components
decomposed <- decompose(ts_measles)
plot(decomposed)
# view trend component
MeaslesTrend<-decomposed$trend
plot(MeaslesTrend, ylab = "Trend Component", xlab = "Time")
arima_model <- auto.arima(sqrt(ts_measles), ic = "aic", trace =T)
##
## ARIMA(2,1,2)(1,0,1)[12] with drift : 499.4077
## ARIMA(0,1,0) with drift : 499.3606
## ARIMA(1,1,0)(1,0,0)[12] with drift : 496.8347
## ARIMA(0,1,1)(0,0,1)[12] with drift : 497.1494
## ARIMA(0,1,0) : 497.3669
## ARIMA(1,1,0) with drift : 496.6547
## ARIMA(1,1,0)(0,0,1)[12] with drift : 496.5915
## ARIMA(1,1,0)(1,0,1)[12] with drift : 498.4702
## ARIMA(1,1,0)(0,0,2)[12] with drift : 498.4692
## ARIMA(1,1,0)(1,0,2)[12] with drift : Inf
## ARIMA(0,1,0)(0,0,1)[12] with drift : 497.3204
## ARIMA(2,1,0)(0,0,1)[12] with drift : 497.5704
## ARIMA(1,1,1)(0,0,1)[12] with drift : 497.8492
## ARIMA(2,1,1)(0,0,1)[12] with drift : 499.5704
## ARIMA(1,1,0)(0,0,1)[12] : 494.5929
## ARIMA(1,1,0) : 494.6626
## ARIMA(1,1,0)(1,0,1)[12] : 496.4714
## ARIMA(1,1,0)(0,0,2)[12] : 496.4702
## ARIMA(1,1,0)(1,0,0)[12] : 494.8373
## ARIMA(1,1,0)(1,0,2)[12] : 498.464
## ARIMA(0,1,0)(0,0,1)[12] : 495.3209
## ARIMA(2,1,0)(0,0,1)[12] : 495.5724
## ARIMA(1,1,1)(0,0,1)[12] : 495.8514
## ARIMA(0,1,1)(0,0,1)[12] : 495.1504
## ARIMA(2,1,1)(0,0,1)[12] : 497.5724
##
## Best model: ARIMA(1,1,0)(0,0,1)[12]
# Print the ARIMA model with the lowest AIC
summary(arima_model)
## Series: sqrt(ts_measles)
## ARIMA(1,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 sma1
## -0.1856 -0.1842
## s.e. 0.1115 0.1256
##
## sigma^2 = 21.49: log likelihood = -244.3
## AIC=494.59 AICc=494.9 BIC=501.85
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.02615164 4.5526 3.078198 NaN Inf 0.4495227 0.01862956
cases_fit<-arima(x=sqrt(ts_measles), order=c(1,1,0),
seasonal=list(order=c(0,0,1)))
forecast<-forecast(cases_fit,h=24)
forecast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2023 -0.20326985 -6.072703 5.666163 -9.179793 8.773253
## Feb 2023 0.05873514 -7.510793 7.628263 -11.517859 11.635329
## Mar 2023 -0.99065818 -10.052644 8.071328 -14.849770 12.868454
## Apr 2023 -0.32484884 -10.648064 9.998366 -16.112844 15.463146
## May 2023 -0.68124201 -12.130561 10.768077 -18.191464 16.828980
## Jun 2023 -0.15975139 -12.633413 12.313910 -19.236570 18.917067
## Jul 2023 -0.53597673 -13.956108 12.884155 -21.060297 19.988343
## Aug 2023 -0.65782193 -14.961919 13.646275 -22.534051 21.218407
## Sep 2023 1.74930633 -13.387224 16.885836 -21.400018 24.898631
## Oct 2023 1.90410117 -14.021409 17.829611 -22.451864 26.260067
## Nov 2023 1.63539480 -15.041811 18.312601 -23.870191 27.140980
## Dec 2023 1.88733016 -15.509121 19.283782 -24.718247 28.492907
## Jan 2024 1.84056413 -15.981054 19.662183 -25.415249 29.096378
## Feb 2024 1.84924517 -16.431257 20.129747 -26.108370 29.806860
## Mar 2024 1.84763373 -16.872452 20.567719 -26.782267 30.477534
## Apr 2024 1.84793286 -17.303106 20.998971 -27.441053 31.136919
## May 2024 1.84787733 -17.724363 21.420117 -28.085280 31.781035
## Jun 2024 1.84788764 -18.136726 21.832502 -28.715942 32.411717
## Jul 2024 1.84788573 -18.540755 22.236526 -29.333849 33.029621
## Aug 2024 1.84788608 -18.936930 22.632703 -29.939748 33.635520
## Sep 2024 1.84788602 -19.325695 23.021467 -30.534311 34.230083
## Oct 2024 1.84788603 -19.707448 23.403220 -31.118153 34.813925
## Nov 2024 1.84788603 -20.082558 23.778330 -31.691833 35.387606
## Dec 2024 1.84788603 -20.451358 24.147130 -32.255865 35.951637
plot(forecast(cases_fit,h=24), ylab="sqrt of cases", xlab="Time(years)",lwd=1)
legend("topleft", c("Measles cases"), cex=0.8, col=1,lty=1)
checkresiduals(cases_fit)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(0,0,1)[12]
## Q* = 27.214, df = 15, p-value = 0.02704
##
## Model df: 2. Total lags used: 17