##
===========================================
=================================
## TIME SERIES ANALYSIS OF NIGERIA'S
POPULATION DENSITY (1961-2023)
## Chapter Four - Complete R Script
##
## Author : [Your Name]
## Data :
Time_series_of_Nigeria_population_density.x
lsx
## (World Bank / FAO - Population
density, persons per sq km, 1961-2023)
##
## Required packages: readxl, tseries,
forecast, ggplot2, gridExtra
## Install once with:
##
[Link](c("readxl","tseries","fore
cast","ggplot2","gridExtra"))
##
===========================================
=================================
## ---- 0. SETUP --------------------------
-----------------------------------
library(readxl)
library(tseries)
library(forecast)
library(ggplot2)
library(gridExtra)
[Link](42)
options(scipen = 999)
## ---- 1. IMPORT AND PREPARE DATA --------
-----------------------------------
# Update the path below to the location of
your file
df <-
read_excel("Time_series_of_Nigeria_populati
on_density.xlsx")
colnames(df) <- c("Year", "Density")
df <- df[order(df$Year), ]
# Check completeness
sum([Link](df$Density)) #
should be 0
str(df)
# Convert to a formal annual time series
object
y <- ts(df$Density, start = df$Year[1],
frequency = 1)
print(y)
## ---- 2. DESCRIPTIVE STATISTICS (Table
4.1) --------------------------------
descr <- [Link](
n = length(y),
Mean = mean(y),
Median = median(y),
SD = sd(y),
Min = min(y),
Max = max(y),
Range = max(y) - min(y),
CV_pct = 100 * sd(y) / mean(y),
Skewness = e1071::skewness(y),
Kurtosis = e1071::kurtosis(y) #
excess kurtosis
)
print(round(descr, 3))
## ---- 3. TREND AND GROWTH RATE ANALYSIS
(Table 4.2) ------------------------
growth <- 100 * diff(y) / head(y, -1)
cat("Mean annual growth rate:",
round(mean(growth), 3), "%\n")
cat("SD of annual growth rate:",
round(sd(growth), 3), "%\n")
cat("Min growth:", round(min(growth), 3),
"% in", df$Year[[Link](growth) + 1],
"\n")
cat("Max growth:", round(max(growth), 3),
"% in", df$Year[[Link](growth) + 1],
"\n")
# Compound Annual Growth Rate (CAGR)
n_years <- max(df$Year) - min(df$Year)
cagr <- ([Link](tail(y, 1)) /
[Link](head(y, 1)))^(1 / n_years) - 1
cat("CAGR 1961-2023:", round(100 * cagr,
3), "%\n")
# Decade-wise average growth rate
decade <- (df$Year[-1] %/% 10) * 10
tapply(growth, decade, mean)
## ---- 4. TIME PLOT (Figure 4.1) ---------
-----------------------------------
p1 <- ggplot(df, aes(x = Year, y =
Density)) +
geom_line(color = "#1f4e79", linewidth =
1) +
geom_point(color = "#1f4e79", size = 1) +
labs(title = "Figure 4.1: Nigeria
Population Density, 1961-2023",
x = "Year", y = "Population Density
(persons/km2)") +
theme_minimal(base_size = 12)
print(p1)
ggsave("fig41_level.png", p1, width = 7,
height = 4, dpi = 150)
## ---- 5. DIFFERENCED SERIES (Figure 4.2)
-----------------------------------
d1 <- diff(y, differences = 1)
d2 <- diff(y, differences = 2)
par(mfrow = c(2, 1))
plot(d1, main = "(a) First Difference of
Population Density",
ylab = "Delta Density", col =
"#c0504d"); abline(h = 0)
plot(d2, main = "(b) Second Difference of
Population Density",
ylab = "Delta^2 Density", col =
"#548235"); abline(h = 0)
par(mfrow = c(1, 1))
## ---- 6. TREND-RESIDUAL DECOMPOSITION
(Figure 4.3) -------------------------
# Annual (non-seasonal) data: classical
stats::decompose() requires frequency > 1,
# so a 5-year centred moving average is
used to extract the trend component.
trend <- stats::filter(y, filter = rep(1/5,
5), sides = 2)
resid_dec <- y - trend
par(mfrow = c(3, 1))
plot(y, main = "Observed", col = "#1f4e79")
plot(trend, main = "Trend (5-Year Centred
Moving Average)", col = "#e36c09")
plot(resid_dec, main = "Residual", col =
"grey40"); abline(h = 0)
par(mfrow = c(1, 1))
## ---- 7. STATIONARITY TESTING: AUGMENTED
DICKEY-FULLER (Table 4.3) ---------
# k = trunc((n-1)^(1/3)) lag order, as used
by tseries::[Link] by default
[Link](y, k = trunc((length(y) -
1)^(1/3))) # level series
[Link](d1, k = trunc((length(d1) -
1)^(1/3))) # first difference
[Link](d2, k = trunc((length(d2) -
1)^(1/3))) # second difference
# Complementary Phillips-Perron test
[Link](y)
[Link](d1)
[Link](d2)
# KPSS test (complementary; null =
stationary)
[Link](y, null = "Trend")
[Link](d2, null = "Level")
## ---- 8. MODEL IDENTIFICATION: ACF / PACF
(Figure 4.4) ---------------------
par(mfrow = c(1, 2))
acf(d2, [Link] = 15, main = "ACF (Delta^2
series)")
pacf(d2, [Link] = 15, main = "PACF
(Delta^2 series)")
par(mfrow = c(1, 1))
## ---- 9. ARIMA MODEL ESTIMATION AND
SELECTION (Table 4.4) ------------------
# 9a. Candidate models identified from the
correlogram, compared by AIC/BIC
candidates <- list(
c(1, 2, 0), c(0, 2, 1), c(1, 2, 1),
c(2, 2, 0), c(2, 2, 1), c(1, 2, 2), c(0,
2, 2)
)
model_compare <- [Link](rbind,
lapply(candidates, function(ord) {
fit <- tryCatch(Arima(y, order = ord,
method = "ML"),
error = function(e)
NULL)
if ([Link](fit)) return(NULL)
[Link](
Model = paste0("ARIMA(", ord[1], ",",
ord[2], ",", ord[3], ")"),
AIC = AIC(fit), BIC = BIC(fit)
)
}))
model_compare <-
model_compare[order(model_compare$AIC), ]
print(model_compare)
# 9b. Cross-check with automatic selection
auto_fit <- [Link](y, d = 2, stepwise =
FALSE, approximation = FALSE,
seasonal = FALSE,
trace = TRUE)
summary(auto_fit)
# 9c. Fit the selected model explicitly:
ARIMA(0,2,2)
fit <- Arima(y, order = c(0, 2, 2), method
= "ML")
summary(fit)
coeftest_result <- lmtest::coeftest(fit)
# requires lmtest package
print(coeftest_result)
## ---- 10. DIAGNOSTIC CHECKING (Table 4.5,
Figure 4.5) ----------------------
checkresiduals(fit)
# combined diagnostic plot
res <- residuals(fit)
# Ljung-Box test at multiple lags
for (h in c(6, 10, 12)) {
print([Link](res, lag = h, type =
"Ljung-Box", fitdf = 2))
}
# Normality tests
[Link]([Link](res))
tseries::[Link]([Link](res))
# Identify largest residuals (outlier
inspection)
res_df <- [Link](Year = df$Year[-
c(1,2)], Residual = [Link](res))
res_df[order(-abs(res_df$Residual)), ][1:5,
]
## ---- 11. OUT-OF-SAMPLE FORECAST
VALIDATION (Table 4.6) --------------------
train <- window(y, end = 2018)
test <- window(y, start = 2019)
fit_train <- Arima(train, order = c(0, 2,
2), method = "ML")
fc_test <- forecast(fit_train, h =
length(test))
print([Link](Year = df$Year[df$Year >=
2019],
Actual =
[Link](test),
Forecast =
round([Link](fc_test$mean), 2)))
accuracy(fc_test, test) # RMSE, MAE,
MAPE, etc.
## ---- 12. FORECAST 2024-2030 (Table 4.7,
Figure 4.6) -----------------------
fc <- forecast(fit, h = 7, level = 95)
print(fc)
p6 <- autoplot(fc) +
labs(title = "Figure 4.6: Observed Series
and ARIMA(0,2,2) Forecast, 2024-2030",
x = "Year", y = "Population Density
(persons/km2)") +
theme_minimal(base_size = 12)
print(p6)
ggsave("fig46_forecast.png", p6, width =
7.5, height = 4.2, dpi = 150)
## ---- END OF SCRIPT ---------------------
------------------------------------