Dự Báo Phát Triển Dân Số Việt Nam 2017 - 2026

Time Series Analysis

Để dự đoán phát triển dân số Việt Nam đến 2026, tôi đã lấy dữ liệu dân số VN từ http://www.worldometers.info/ trong khoảng thời gian từ năm 1960 đến 2016 và lưu trữ dưới dạng excel rồi phân tích và dự đoán dân số trong 10 năm tới. Package được sử dụng là forecast cho phân tích dữ liệu dạng time series.

1 Chuẩn bị dữ liệu

library("readxl")
dat <- read_excel("D:/Softwares/R/Data/Vnpop.xlsx")
dat <- as.data.frame(dat)
dat1 <- as.data.frame(t(dat))
library(data.table)
setDT(dat1, keep.rownames = TRUE)[]

##       rn       V1      V2
##  1: 1960 34743000 34.7430
##  2: 1961 35428000 35.4280
##  3: 1962 36123000 36.1230
##  4: 1963 36836000 36.8360
##  5: 1964 37574000 37.5740
##  6: 1965 38341000 38.3410
##  7: 1966 39142000 39.1420
##  8: 1967 39980000 39.9800
##  9: 1968 40856000 40.8560
## 10: 1969 41773000 41.7730
## 11: 1970 42729000 42.7290
## 12: 1971 43725000 43.7250
## 13: 1972 44758000 44.7580
## 14: 1973 45825000 45.8250
## 15: 1974 46918000 46.9180
## 16: 1975 48030000 48.0300
## 17: 1976 49158000 49.1580
## 18: 1977 50295000 50.2950
## 19: 1978 51436000 51.4360
## 20: 1979 52574000 52.5740
## 21: 1980 53700000 53.7000
## 22: 1981 54722000 54.7220
## 23: 1982 55687000 55.6870
## 24: 1983 56655000 56.6550
## 25: 1984 57692000 57.6920
## 26: 1985 58868000 58.8680
## 27: 1986 60249000 60.2490
## 28: 1987 61750000 61.7500
## 29: 1988 63263000 63.2630
## 30: 1989 64774000 64.7740
## 31: 1990 66016700 66.0167
## 32: 1991 67242400 67.2424
## 33: 1992 68450100 68.4501
## 34: 1993 69644500 69.6445
## 35: 1994 70824500 70.8245
## 36: 1995 71995500 71.9955
## 37: 1996 73156700 73.1567
## 38: 1997 74306900 74.3069
## 39: 1998 75456300 75.4563
## 40: 1999 76596700 76.5967
## 41: 2000 77630900 77.6309
## 42: 2001 78620500 78.6205
## 43: 2002 79537700 79.5377
## 44: 2003 80467400 80.4674
## 45: 2004 81436400 81.4364
## 46: 2005 82392100 82.3921
## 47: 2006 83311200 83.3112
## 48: 2007 84218500 84.2185
## 49: 2008 85118700 85.1187
## 50: 2009 86025000 86.0250
## 51: 2010 86932500 86.9325
## 52: 2011 87860300 87.8603
## 53: 2012 88809200 88.8092
## 54: 2013 89759500 89.7595
## 55: 2014 90728900 90.7289
## 56: 2015 91713300 91.7133
## 57: 2016 92701100 92.7011
##       rn       V1      V2

colnames(dat1) <- c("year", "pop", "Population")
as.numeric(dat1$year)

##  [1] 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973
## [15] 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987
## [29] 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001
## [43] 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015
## [57] 2016

Vẽ biểu đồ dân số Việt Nam từ 1960 đến 2016

p <- ggplot(dat1, aes(x= year, y = Population))
p + geom_line(group=1,colour = "blue") +
  theme(axis.text.x = element_text(color="#993333", size=7, angle=85)) +
  scale_y_continuous("Population in millions") +
  ggtitle("Vietnam Population Since 1960 To 2016") 

Biểu đồ có khuynh hướng tăng dần, ổn định (non-stationary trend). Tuy nhiên, không có chu kì biến thiên (có giai đoạn tăng và giai đoạn giảm, lặp lại trong một khoảng thời gian, gọi là seasional) như một biểu đồ Time Series đặc trưng.

Để phân tích và dự đoán dữ liệu dân số này ta set dữ liệu này thành Time Series data bằng các hàm của package forecast như sau:

df <- (dat[2,])
colnames(df) <- NULL
df <- unlist(c(df))
popseries <- ts(df, start = c(1960))
plot.ts(popseries)

Để tính toán các giá trị của dân số theo những thời điểm nhất định (dự báo) chúng ta dùng hàm HoltWinters của package forecast.

popseriesforecasts <- HoltWinters(popseries, beta=FALSE, gamma=FALSE)
plot(popseriesforecasts)

Đường biểu diễn màu đỏ cho thấy hàm HoltWinters tính toán và dự báo các giá trị dân số trong mẫu trên (in-sample prediction) khá tương đồng với giá trị thật quan sát được (Holt-Winters exponential smoothing function).

Trên cơ sở tính toán này, các giá trị dân số trong các thời điểm tương lai ngắn được dự báo. Hàm forecast Holt-Winters giúp ta dự báo trong tương lai gần của phát triển dân số.

popseriesforecasts2 <- forecast:::forecast.HoltWinters(popseriesforecasts, h=10)
popseriesforecasts2

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2017       92.70105 92.46256 92.93954 92.33631 93.06578
## 2018       92.70105 92.36379 93.03831 92.18525 93.21685
## 2019       92.70105 92.28799 93.11411 92.06933 93.33277
## 2020       92.70105 92.22409 93.17800 91.97161 93.43049
## 2021       92.70105 92.16780 93.23430 91.88551 93.51659
## 2022       92.70105 92.11690 93.28519 91.80767 93.59442
## 2023       92.70105 92.07010 93.33200 91.73610 93.66600
## 2024       92.70105 92.02654 93.37556 91.66947 93.73263
## 2025       92.70105 91.98562 93.41648 91.60690 93.79520
## 2026       92.70105 91.94692 93.45518 91.54771 93.85439

plot(popseriesforecasts2)

h = 10 cho biết yêu cầu muốn dự báo trong 10 time points sắp tới, từ 2017 đến 2026. Các giá trị của dân số Việt Nam và khoảng tin cậy 80%, 95% được trình bày trong bảng kết quả trên. Tuy nhiên đường biểu diễn đi không theo hướng phát triển của những năm trước 2016, mà có khuynh hướng giảm gia tăng. Liệu có cách nào khác để cải thiện mức độ chính xác của mô hình dự báo này không?

Hai mô hình thường được nhắc tới trong phân tích Time Series là AR (Automatic Regression) và MA (Moving Average). Mô hình AR muốn nói tới những dữ liệu có tương quan giữa các thời điểm. Có nghĩa dân số năm tới liên quan trực tiếp với dân số hiện có trong năm nay hoặc trong vài năm vừa qua.

Chúng ta cần một mô hình dự báo được cải thiện hơn tính chính xác nếu có thể. Để làm điều đó, chúng ta tiếp tục phân tích sâu hơn về residuals của mô hình prediction trên, để xem những sai số này có tương quan với nhau hay không. Kiểm định sự tương quan của khoảng từ 1 đến 20 giá trị về phía trước như sau:

popseriesforecasts2$residuals

## Time Series:
## Start = 1960 
## End = 2016 
## Frequency = 1 
##  [1]        NA 0.6850000 0.6950355 0.7130361 0.7380370 0.7670383 0.8010398
##  [8] 0.8380416 0.8760435 0.9170455 0.9560476 0.9960496 1.0330517 1.0670536
## [15] 1.0930554 1.1120567 1.1280577 1.1370585 1.1410590 1.1380592 1.1260591
## [22] 1.0220584 0.9650530 0.9680501 1.0370502 1.1760538 1.3810610 1.5010717
## [29] 1.5130779 1.5110785 1.2427784 1.2257645 1.2077636 1.1944627 1.1800620
## [36] 1.1710612 1.1612608 1.1502603 1.1494597 1.1404597 1.0342592 0.9896537
## [43] 0.9172514 0.9297476 0.9690482 0.9557503 0.9191496 0.9073477 0.9002471
## [50] 0.9063467 0.9075470 0.9278471 0.9489482 0.9503492 0.9694493 0.9844503
## [57] 0.9878511

Box.test(popseriesforecasts2$residuals, lag=20, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  popseriesforecasts2$residuals
## X-squared = 175.78, df = 20, p-value < 2.2e-16

p-value < 2.2e-16. Vậy đã có bằng chứng là các errors trong khoảng lùi 1 - 20 thời điểm có tương quan với nhau.

Ngoài ra, người ta còn đánh giá xem các errors trong 20 thời điểm cuối có phân bố chuẩn và biến thiên hằng định hay không bằng cách vẽ biểu đồ của những residuals. Trong những residual ở trên ở vị trí thứ nhất là NA, nên chúng ta loại giá trị này ra để hàm acf chạy không bị lỗi.

popseriesforecasts2$residuals <- popseriesforecasts2$residuals[-1]
acf(popseriesforecasts2$residuals,lag.max=20)

pacf(popseriesforecasts2$residuals, lag.max=20)

Biểu đồ ACF cho thấy sự tương quan giữa các giá trị dân số giảm dần khi đi xa dần về phía các thời điểm trước đó. Mô hình AR(2) phù hợp với biến thiên của dữ liệu này.

Để đánh giá xem errors có phân phối chuẩn, mean = 0 và biến thiên hằng định không chúng ta phân tích phân phối của những errors như sau:

# Function for forecasterror
plotForecastErrors <- function(forecasterrors)
{
  # make a histogram of the forecast errors:
  mybinsize <- IQR(forecasterrors)/4
  mysd <- sd(forecasterrors)
  mymin <- min(forecasterrors) - mysd*5
  mymax <- max(forecasterrors) + mysd*3
  # generate normally distributed data with mean 0 and standard deviation mysd
  mynorm <- rnorm(10000, mean=0, sd=mysd)
  mymin2 <- min(mynorm)
  mymax2 <- max(mynorm)
  if (mymin2 < mymin) { mymin <- mymin2 }
  if (mymax2 > mymax) { mymax <- mymax2 }
  # make a red histogram of the forecast errors, with the normally distributed
  #→data overlaid:
  mybins <- seq(mymin, mymax, mybinsize)
  hist(forecasterrors, col="red", freq=FALSE, breaks=mybins)
  # freq=FALSE ensures the area under the histogram = 1
  # generate normally distributed data with mean 0 and standard deviation mysd
  myhist <- hist(mynorm, plot=FALSE, breaks=mybins)
  # plot the normal curve as a blue line on top of the histogram of forecast
  #→errors:
  points(myhist$mids, myhist$density, type="l", col="blue", lwd=2)
}
plotForecastErrors(popseriesforecasts2$residuals)

Biểu đồ plotForecastErrors cho thấy các errors không phân phối chuẩn và mean khác 0.

2 Dựng Mô hình ARIMA

Đây là mô hình dự báo các giá trị tương lai gần của dữ liệu Time Series. Mô hình mới này cần ba thông số p, d và q. Hàm auto.arima() giúp ta xác định 3 thông số này.

auto.arima(popseries)

## Series: popseries 
## ARIMA(2,2,2)                    
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       0.8277  -0.5780  -0.3921  0.6941
## s.e.  0.1964   0.1677   0.1616  0.1583
## 
## sigma^2 estimated as 0.002612:  log likelihood=87.57
## AIC=-165.14   AICc=-163.92   BIC=-155.11

Ba thông số đó là ARIMA(2,2,2) với p = 2, d = 2 và q = 2. Đây là một mô hình hỗn hợp của AR và MA.

Từ đó viết code cho mô hình ARIMA để dự báo gần cho dân số như sau:

popseriesarima <- arima(popseries, order=c(2,2,2))
popseriesforecasts <- forecast(popseriesarima, h=10)
popseriesforecasts

##      Point Forecast     Lo 80     Hi 80    Lo 95     Hi 95
## 2017       93.68824  93.62570  93.75077 93.59260  93.78388
## 2018       94.67311  94.50846  94.83775 94.42131  94.92490
## 2019       95.65648  95.33863  95.97432 95.17037  96.14258
## 2020       96.63992  96.12871  97.15113 95.85809  97.42175
## 2021       97.62429  96.89576  98.35282 96.51010  98.73848
## 2022       98.60939  97.64988  99.56890 97.14195 100.07683
## 2023       99.59455  98.39226 100.79684 97.75580 101.43329
## 2024      100.57934  99.11964 102.03904 98.34692 102.81176
## 2025      101.56379  99.82929 103.29830 98.91110 104.21649
## 2026      102.54818 100.52118 104.57518 99.44815 105.64821

plot(popseriesforecasts)

Biểu đồ dân số Việt Nam giai đoạn 2017 - 2026 được vẽ tiếp theo kèm với khoàng tin cậy 80% và 95%. Biểu đồ này phần nào phản ánh được khuynh hướng gia tăng dân số của Việt Nam theo tương quan với những năm trước đó. Theo kết quả dự báo này, đến năm 2024 dân số nước ta vượt mốc 100 triệu người.

Để vẽ chi tiết hơn giai đoạn này ta xử lí dữ liệu như sau:

library(data.table)
pop2026 <- as.data.frame(popseriesforecasts)
setDT(pop2026, keep.rownames = TRUE)[]

##       rn Point Forecast     Lo 80     Hi 80    Lo 95     Hi 95
##  1: 2017       93.68824  93.62570  93.75077 93.59260  93.78388
##  2: 2018       94.67311  94.50846  94.83775 94.42131  94.92490
##  3: 2019       95.65648  95.33863  95.97432 95.17037  96.14258
##  4: 2020       96.63992  96.12871  97.15113 95.85809  97.42175
##  5: 2021       97.62429  96.89576  98.35282 96.51010  98.73848
##  6: 2022       98.60939  97.64988  99.56890 97.14195 100.07683
##  7: 2023       99.59455  98.39226 100.79684 97.75580 101.43329
##  8: 2024      100.57934  99.11964 102.03904 98.34692 102.81176
##  9: 2025      101.56379  99.82929 103.29830 98.91110 104.21649
## 10: 2026      102.54818 100.52118 104.57518 99.44815 105.64821

colnames(pop2026) <- c("year", "population", "Lo.80", "Hi.80", "Lo.95", "Hi.95")

Vẽ biểu đồ dân số giai đoạn này:

library(ggplot2)
g <- ggplot(pop2026, aes(x = year, group = 1))
g + geom_point(aes(y=population)) +
  geom_line(aes(y=population), colour = "blue", size=1) +
  geom_line(aes(y=Lo.80), colour="orange", linetype = "dotted", size=1) +
  geom_line(aes(y=Hi.80), colour="orange", linetype = "dotted", size=1) +
  geom_line(aes(y=Lo.95), colour="green", linetype = "dashed", size=1) +
  geom_line(aes(y=Hi.95), colour="green", linetype = "dashed", size=1) +
  ggtitle("Forecast for Vietnam Population in 2017 - 2026")+
  scale_y_continuous("Population in millions") +
  theme(axis.line = element_line(colour="black"))

Biểu đồ dân số Việt Nam giai đoạn 2017 - 2026 được vẽ tiếp theo kèm với khoàng tin cậy 80% và 95%. Đến năm 2024, dân số Việt Nam được dự báo vượt trên cột mốc 100 triệu người.

Bài phân tích này chỉ dựa vào một dataset cụ thể, chưa phản ánh được những đặc trưng đầy đủ của các loại time series data. Trong y học, một số dữ liệu như tần suất sanh theo tháng, mùa, năm hoặc tỉ suất bệnh mới theo tháng, năm là những ví dụ của loại dữ liệu time series. Theo dõi số lượng bệnh nhân đến khám ở một cơ sở y tế theo tháng, mùa, năm cũng là những dữ liệu thuộc loại Time Series.

Để biết thêm về loại phân tích này, Time Series Analysis, các bạn có thể tải xuống và đọc tài liệu sau đây: https://media.readthedocs.org/pdf/a-little-book-of-r-for-time-series/latest/a-little-book-of-r-for-time-series.pdf

Cám ơn các anh chị, các bạn đã dành thời gian để đọc bài viết này. Mong được góp ý, bàn luận để được hiểu biết rõ hơn.