基于灰色–马尔可夫模型对中国结婚人数的预测

Prediction of the Number of Marriages in China Based on Gray-Markov Model

摘要: 中国的老龄化越来越严重，而中国的结婚人数又持续走低，这对中国的年龄结构有很大的影响，劳动力市场也会不完善。因此预测中国的结婚人数有利于制定相关的政策以及采取相应的措施去应对当前形势。本文主要采用灰色–马尔可夫链进行预测，并将其与直接由GM(1,1)模型得到的预测值以及时间序列预测所得的结果进行对比。结果表明利用灰色–马尔可夫链得到的预测值比直接用GM(1,1)模型以及时间序列预测得到的值要好。

Abstract: China’s aging is getting more and more serious, and the number of marriages in China continues to decline, which has a great impact on China’s age structure and the labor market will be imperfect. Therefore, predicting the number of marriages in China is conducive to making relevant policies and taking corresponding measures for the current situation. In this paper, gray-Markov chain is mainly used for forecasting, and the results are compared with those obtained from time series forecasting. The results show that the predicted values obtained by using gray-Markov chain are better than those obtained by using GM(1,1) model and time series directly.

1. 引言

中国的结婚人数持续走低，更是在2021年创下了历史新低。文献 [1] 研究了17~34岁年龄段青年的婚恋观；文献 [2] 调查了当代女大学生的婚恋观；文献 [3] 研究了农村70后、80后、90后的婚恋观状况及其影响因素；文献 [4] 研究了新生代农民工婚恋观现状及其影响因素。这些文献都是建立在影响因素上去探究的，没有从数量上去探究未来结婚人数的变化。

文献 [5] 利用BP神经网络，灰色预测模型，回归分析法对石油的消费情况进行了预测。文献 [6] 在原有的灰色–马尔可夫模型的基础上提出了新的状态的划分方式。文献 [7] 利用多元线性回归，灰色马尔可夫模型以及指数平滑法对对石油进行了预测。

中国的结婚人口不是单纯递增递减的形式，由于原始数据集年度时间数据集，所以排除利用回归模型对数据进行预测。由于当前的未来的结婚数据只会与当前是否结婚与离婚有关，并且从数据的趋势看，中国的结婚人口数据波动比较大。灰色–马尔可夫模型可以克服较大的随机性，并且只需要少量数据就可以对未来进行预测，因此选择利用灰色–马尔可夫模型进行预测。同时，由于数据是时间序列数据，因此也选用时间序列的方法对数据进行预测，比较两者的效果。

本文是从以往数据出发来预测未来结婚人数的高低。这有利于国家制定相应的政策来缓解当前结婚率低的压力。结婚人数创新低有很多主观和客观原因。主观原因是21世纪人们的受教育程度普遍提高，择偶标准发生变化。不再是古代的父母之命，媒妁之言。客观原因是社交圈狭小，没有机会去认识新的人。不婚主义者是少数，人们普遍不婚主要是受到物质和社交圈的限制而导致不婚。因此本文的研究具有一定的现实意义。

2. 模型介绍

2.1. GM(1,1)

第一：将所有的原始数据设置为与时间相关的序列

$X^{\left(0\right)}=\left(x^{\left(0\right)}\left(1\right),x^{\left(0\right)}\left(2\right),\cdots ,x^{\left(0\right)}\left(n\right)\right)$ (1)

第二：对原始序列进行累加随时间递增的单调序列，让 $X^{\left(1\right)}$ 作为 $X^{\left(0\right)}$ 的一个累加序列

$X^{\left(1\right)}=\left(x^{\left(1\right)}\left(1\right),x^{\left(1\right)}\left(2\right),\cdots ,x^{\left(1\right)}\left(n\right)\right)$ (2)

$X^{\left(1\right)}\left(k\right)={\displaystyle {\sum }_{i=0}^k{X}^{\left(0\right)}\left(i\right)}$ (3)

第三：令 $Z^{\left(1\right)}$ 是由 $X^{\left(1\right)}$ 生成的

$Z^{\left(1\right)}=\left(Z^{\left(1\right)}\left(1\right),Z^{\left(1\right)}\left(2\right),\cdots ,Z^{\left(1\right)}\left(n\right)\right)$ (4)

$Z^{\left(1\right)}\left(k\right)=0.5X^{\left(1\right)}\left(k\right)+X^{\left(1\right)}\left(k-1\right)$ (5)

第四：确定GM(1,1)的精确方程，即GM(1,1)的灰色微分方程模型为：

$x^{\left(1\right)}\left(k\right)+a{Z}^{\left(1\right)}\left(k\right)=b$ (6)

其中a称为发展系数，b为灰色作用量，假设 $\hat{a}$ 是要估计的系数，并且

$\hat{a}={\left(a,b\right)}^{\text{T}}$, $\hat{a}={\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}{Y}_n$ (7)

其中

$B=\left(\begin{array}{cc}-z^{\left(1\right)}\left(2\right)& 1\\ -z^{\left(1\right)}\left(3\right)& 1\\ ⋮& ⋮\\ -z^{\left(1\right)}\left(n\right)& 1\end{array}\right)$, $y_{\left(n\right)}=\left(\begin{array}{c}{x}^{\left(0\right)}\left(2\right)\\ x^{\left(0\right)}\left(3\right)\\ ⋮\\ x^{\left(0\right)}(\; n\; )\end{array}\right)$

$\frac{\text{d}{x}^{\left(1\right)}}{\text{d}t}+a{x}^{\left(1\right)}=b$，叫白化方程，也叫影子方程。

第五：根据以上步骤得到最终的预测方程

1)

$\widehat{x^{\left(1\right)}}\left(k\right)=\left(x^{\left(1\right)}\left(0\right)-\frac{b}{a}\right){\text{e}}^{-ak}+\frac{b}{a}$ (8)

2)

$\widehat{x^{\left(1\right)}}\left(k\right)=\left(x^{\left(1\right)}\left(0\right)-\frac{b}{a}\right){\text{e}}^{-ak}+\frac{b}{a},k=1,2,\cdots ,n$ (9)

3)

$\widehat{x^{\left(1\right)}}\left(k\right)=\left(x^{\left(1\right)}\left(1\right)-\frac{b}{a}\right){\text{e}}^{-ak}+\frac{b}{a},k=1,2,\cdots ,n$ (10)

4)

$\widehat{x^{\left(0\right)}}\left(k+1\right)=\widehat{x^{\left(1\right)}}\left(k+1\right)-\widehat{x^{\left(1\right)}}\left(k\right)$ (11)

2.2. 灰色–马尔可夫模型

GM(1,1)预测模型在灰色系统理论里占有重要地位，该模型研究的出发点是在自身时间序列中探索有价值的信息，发掘研究内容的规律，不需要考虑研究内容带来的影响。灰色系统模型以少量数据信息为研究对象，该模型对少量信息建模有较好的预测精度 [6]。马尔可夫链是指未来的状态只与当前的状态相关联，而与过去的状态无关。马尔可夫预测模型的原理是基于当前时间的状态和状态间的转移概率来预测未来时间的状态。马尔可夫链预测与灰色系统模型不同，它弥补了灰色预测模型的不足，可以对波动性大的数据进行研究 [6]。

在建立马尔可夫模型之前，需要对数据进行检验分别为级比检验和光滑度检验，级比的计算公式为：

$\rho \left(k\right)=\frac{X^{\left(0\right)}\left(k-1\right)}{X^{\left(0\right)}\left(k\right)}\left(k=2,3,\cdots ,n\right)$ (12)

如果所有级比都在 $\left({\text{e}}^{-\frac{2}{n+1}},{\text{e}}^{\frac{2}{n+1}}\right)$ 范围之内，则可以对原始数据直接建立灰色模型。光滑比的计算公式：

$\delta \left(k\right)=\frac{X^{\left(0\right)}\left(k\right)}{X^{\left(1\right)}\left(k-1\right)}\left(k=2,3,\cdots ,n\right)$ (13)

其中 $X^{\left(1\right)}\left(k-1\right)={\displaystyle {\sum }_{i=1}^{k-1}{X}^{\left(0\right)}\left(i\right)}\left(k=2,3,\cdots ,n\right)$，当原始数据序列满足 $\delta \left(k\right)$ 是k的递减函数时，才可以建模预测，只有当原始数据序列通过级比检验和光滑度检验时，才可以进行建模预测，否则需要对原始数据序列进行预处理。同时还要对发展系数计算发展系数−a，其关系见表1：

Table 1. Application range of development coefficient −a and GM(1,1) model

表1. 发展系数−a与GM(1,1)模型的适用范围

马尔科夫链预测的步骤如下：

第一、计算 $\gamma$ 值

$\gamma$ 值为实际值比上预测值，计算公式为：

$\gamma =\frac{X^{\left(0\right)}\left(k\right)}{X^{\left(1\right)}\left(k\right)}$ (14)

根据 $\gamma$ 值进行状态的划分，把 $\gamma$ 值划分为m个区间，即m个状态，任一状态区间可以表示为

$E_i=\left({\gamma }_i,{\gamma }_j\right)$

$\gamma$ 的所有状态区间集合可以表示为 $E=\left(E_1,E_2,\cdots ,E_m\right)$。

第二、计算状态转移概率

$P_{ij}\left(s\right)=\frac{M_{ij}\left(s\right)}{M_i}$ (15)

第三、计算状态转移概率矩阵

$P\left(s\right)=\left[\begin{array}{cccc}{p}_{11}\left(s\right)& p_{12}\left(s\right)& \cdots & p_{1m}\left(s\right)\\ p_{21}\left(s\right)& p_{22}\left(s\right)& \cdots & p_{2m}\left(s\right)\\ ⋮& ⋮& \ddots & ⋮\\ p_{m1}\left(s\right)& p_{m2}\left(s\right)& \cdots & p_{mm}\left(s\right)\end{array}\right]$

其中 $P_{i1}\left(1\right)+P_{i2}\left(s\right)+\cdots +P_{im}\left(s\right)=1$。

第四、确定对象转移状态

马尔科夫链的下一步状态和预测对象过去的状态没有关系，只和当前的状态有关。当预测对象的状态在 $E_i$ 时，只需观察状态转移概率矩阵中最大概率转入下一状态即可。如果有多个状态的转移概率相同，则可以具体分析应该转入哪一个状态即可。

第五、确定修正后的预测值

根据文献 [6] 提出的状态划分方式，每个状态区间的对应的数据不少于7个，且进行划分时，将两端的数据各划分为一个状态中间的数据平均分配。因此修正后的预测值的计算公式为：

$\hat{Y}\left(k\right)=\widehat{X^{\left(0\right)}}\left(k\right)\ast \left(\frac{{\gamma }_i+{\gamma }_j}{2}\right)$ (16)

3. 实证

本次的数据来源于中华人民共和国民政部。在建立GM(1,1)模型之前，已对数据进行级比检验和光滑度检验，且检验都通过，可以对数据进行灰色–马尔可夫链模型进行预测。同时计算出发展系数−a = 0.004164524，即原始数据可以进行中长期预测。表2给出了根据GM(1,1)模型得出的相关数据，并根据 $\gamma$ 值进行了状态的划分，结果见表2：

Table 2. Simulation data of GM(1,1) model

表2. GM(1,1)模型的模拟数据

对数据进行状态的划分是灰色–马尔可夫的重点，一般用残差的相对值进行状态的划分，而本篇文章根据 $\gamma$ 进行状态的划分，见文献 [6]，由以上数据我们可把数据分为5个状态。状态区间分别为：

$M_1:\left[0.75,0.92\right]$， $M_2:\left[0.92,1.007\right]$， $M_3:\left[1.007,1.05\right]$， $M_4:\left[1.05,1.14\right]$， $M_5:\left[1.14,1.36\right]$

一步状态转移概率矩阵为

$P=\left[\begin{array}{ccccc}\frac{6}{7}& \frac{1}{7}& 0& 0& 0\\ \frac{2}{7}& \frac{2}{7}& \frac{3}{7}& 0& 0\\ 0& \frac{2}{7}& \frac{3}{7}& \frac{2}{7}& 0\\ 0& \frac{1}{7}& \frac{1}{7}& \frac{4}{7}& \frac{1}{7}\\ 0& 0& 0& \frac{1}{7}& \frac{6}{7}\end{array}\right]$

利用马尔可夫模型算出预测值以后，同时以1986~2016年为训练集，2017~2021作为测试集，根据时间序列中的移动平均法和两参数法，利用R语言对数据进行预测，相关代码见附录，代码均来源于文献 [8]。不同模型的预测值见表3；GM(1,1)、灰色–马尔可夫模型与原始数据的对比见图1；移动平均法的预测结果见图2；两参数模型的预测结果见图3；移动平均法的趋势拟合见图4；两参数模型的趋势拟合图见图5。

Table 3. Prediction data of related models

表3. 相关模型的预测数据

Figure 1. Comparison of raw data with GM(1,1) and Markov model predictions

图1. 原始数据与GM(1,1)和马尔可夫模型预测值的对比图

Figure 2. Forecast trend of moving average method

图2. 移动平均法的预测趋势

Figure 3. Predicted trends of the two parameters

图3. 两参数的预测趋势

Figure 4. Fitting effect of moving average method

图4. 移动平均法的拟合效果

Figure 5. Fitting effect of the two-parameter model

图5. 两参数模型的拟合效果

4. 结论

从以上结果可知，单纯利用灰色模型和时间序列模型只能预测数据变化的趋势，不能克服随机波动因素的影响。而改进的灰色–马尔可夫模型可以克服这一缺点，将预测精度提高。通过观察，在状态进行转化时，预测误差会增大，这可能与选取的状态有关，因此在使用灰色–马尔可夫模型时，要注意状态的选取。同时利用灰色–马尔可夫模型对2022年的数据进行了预测，得出2022年中国的结婚人口数据会有一定的增大。根据时间序列模型预测的2022年结婚人口数据没有可比性，误差很大。综上，选用灰色–马尔可夫模型可以较好地预测中国的结婚人口数据。

附录

y1<-read.csv(C:/Users/Administrator/Desktop/data1.csv,header=F)

y2<-read.csv(C:/Users/Administrator/Desktop/data2.csv,header=F)

y3<-read.csv(C:/Users/Administrator/Desktop/data3.csv,header=F)

y11<-ts(y1,start=c(1986,1),frequency =1)

y22<-ts(y2,start=c(1986,1),frequency =1)

y33<-ts(y3,start=c(1986,1),frequency =1)

plot(y11,main=黑色：原始值 绿色：GM(1,1) 蓝色：马尔科夫链)

lines(y22,col=3)

lines(y33,col=4)

yy<-read.csv(C:/Users/Administrator/Desktop/data.csv,header=F)###读取数据

X<-read.csv(C:/Users/Administrator/Desktop/测试集.csv,header=F)

yy1<-ts(yy,start=c(1986,1),frequency =1)

install.packages(tseries)

library(tseries)#---------------###单位根检验

adf.test(yy1)#-------------------####平稳的

plot(yy1,main=时序图)#------------------------###时序图

acf(yy1,main=自相关图)

pacf(yy1,main=偏自相关图)

Box.test(yy1,type=Ljung-Box,lag=6)#---------###LB统计量纯随机序列检验

Box.test(yy1,type=Box-Pierce,lag=6)#--------###Q统计量纯随机序列检验

####平滑法（移动平均法）

install.packages(TTR)

library(TTR)

install.packages(forecast)

library(forecast)

y.fit2<-SMA(yy1,n=5)

plot(yy1,main=移动平均法的模拟效果)

lines(y.fit2,col=2)

y.fore2<-forecast(y.fit2,h=6)

y.fore2

plot(y.fore2,main=移动平均法的预测结果)

lines(X,col=3)###测试集

###两参数

y.fit3<-HoltWinters(yy1,gamma=F)

plot(y.fit3,main=两参数模型的模拟效果)

y.fore3<-forecast(y.fit3,h=6)

plot(y.fore3,main=两参数模型的模拟效果)

lines(yy1,col=3)