数学建模

县域人口结构及劳动力供需状况分析与预测

楼方舟，数信学院，统计122，644465

陆怡君，数信学院，统计121，650453

殷飚，数信学院，统计122，631030

摘要

本文使用EXCEl对2010年该县第六次人口普查的相关数据进行处理，分析当年该县人口的总体规模和结构状况，包括：人口总数，性别比，未成年人数（0-14岁），劳动龄人口（15-60岁），老年人口（60岁以上），老龄化状况以及社会抚养比等。

从该县的实际情况和人口总数，生育、死亡和迁移情况的特点出发，针对人口规模，提出了Logistic、动态模拟等方法进行建模预测。

首先，在最简单的假设下，依照该县人口的历史数据，运用线形最小二乘法对其进行拟合，对2016至2020年的人口规模进行了预测。在此模型中，由于并没有考虑人口的年龄、出生人数男女比例等因素，只是粗略的进行了预测，所以只对中短期人口做了预测，理论上很好，实用性不强，有一定的局限性。

将人口增长的预测问题转化为对出生率、死亡率和城镇乡转移率的预测。通过原题数据的分析研究，发现影响人口增长的主要因素可以归结为出生率、死亡率和城镇乡转移率，并依此建立了不同参数随时间变化的递推数学模型，讨论了各个参数对人口增长的影响。

通过该县2003-2013各年实际从业人数以及2016-2020年预测劳动力需求量，使用灰色模型，估计2016-2020年各年劳动力供需缺口。

通过Eviews作图初步判断2003年到2013年从业人员的发展趋势，采用线性二次移动平均模型和GM(1，1)模型对2016年到2020年的从业人员进行预测。

最后，本文对模型进行了评价，特别指出了各个模型的优缺点，同时也对模型进行了合理性分析，针对该县的人口情况给政府提出了建议。

关键词: Logistic模型二次移动平均法微分方程 Eviews 灰色系统 EXCEL

1 问题重述

1.1问题背景

某县政府拟制定该县“十三五”（2016-2020年）规划总体发展战略，为此需要分析和预测该五年计划时段中该县劳动力的总体供需状况和结构性矛盾。劳动力供给的基础是人口，人口规模和结构决定了劳动力的规模和结构；劳动力需求的主要因素是经济发展所带动的产业发展。

1.2问题提出

通过前期大量调研获得相关数据;

1、通过2010年该县第六次人口普查的相关数据，基于这些数据分析当年该县人口的总体规模和结构状况，包括：人口总数，性别比，未成年人数（0-14岁），劳动龄人口（15-60岁），老年人口（60岁以上），老龄化状况以及社会抚养比等；

2、通过该县2000-2013年人口总数，生育、死亡和迁移情况，试建立合适的数学模型对2016-2020各年的人口规模进行估算，并对人口结构进行分析；

3、通过该县2010-2013年各学历段教育状况以及2010年六普公报公布的各学历阶段人数及比例，这里某一学历段的人数包括完成该学历段教育但未进入高一级学历段学习的人数和该学历段当前在校学生人数之和，如初中学历人数为当年初中毕业但未进入普通高中或职业高中的人数和当年初中在校学生的总和。对2011-2013以及2016-2020各年各学历人数作出合理的估计；

4、通过该县2003-2013各年实际从业人数以及2016-2020年预测劳动力需求量，估计2016-2020年各年劳动力供需缺口。

2 问题分析

对于问题（1），只需要通过excel将数据分析，将各个年龄段的人相加，再从网上寻找那些指标的意义，通过计算即可求的性别比，老龄化程度，以及社会抚养比。

对于问题（2），后者也可以通过excel将数据分析出来；最关键是对人口的预测，我组采用logistic模型，引入滞后变量，将各个指标的趋势进行建模，最后又使用二次移动平均的方式人口规模预测出来。

对于问题（3），我们建立的灰色预测模型，作为原始序列，在进行级比检验后，构造数据矩阵和数据向量，再通过计算得到发展系数和灰作用量的相应值，利用2010-202013年的数据预测出2014-2020年的数据。

对于问题（4），首先通过Eviews 作图初步判断2003年到2013年从业人员的发展趋势，采用线性二次移动平均模型和GM(1，1)模型对2016年到2020年的从业人员进行预测，然后将得到的结果与附件中给出的市场上的总需求人数进行比较，从而估计2016-2020年各年的劳动力供需缺口。

3 模型假设与符号说明

3.1模型的假设

（1）该县的各个指标数据均正确无误。

（2）每一年的人口总数，人口结构及分布和其他有关各量仅在年末发生变化，变化顺序是：一部分人先死亡，然后一部分人生小孩，最后一部分人迁移 （3）国际迁入迁出对于人口的影响较小 3.2符号的说明 问题二：

()i P t ---第t 年初年龄为i 的总人数

()i t μ---第t 年末年龄为i 的人的死亡率

()t β---第t 年末平均每个育龄女性的生育数 i h ---生育加权因子（即生育模式）

()t σ---第t 年末生育性别比

()i t λ---第t 年末迁移人口比例

()O t ---老龄化程度 ()A t ---城镇化程度

()G t ---性别比

问题四：

(),

F r t--------人口分布函数；即t时刻年龄小于r的人口。

(),

p r t-------人口密度函数；

()

N t--------t时刻的人口总数

(),r t

μ------死亡率；

()t

β-----总和生育率;

(),

Q r t---------迁移函数

()

f t------婴儿绝对出生函数

()

u t-------相对出生率函数；

(),

h r t---------t年r岁妇女的生育模式；

(),

k r t-----在t年r岁人口中女性所占比例

()

p r-----初始年代按年龄人口密度函数；

[]

,r r------育龄区

12

4 模型的准备

在对模型进行预测前，先将数据进行图像绘制，对于需要拟合的函数进行尝试，选择一个拟合度最好的函数。对于问题一而言，需要查阅相关文献，得知相应指标的意义。对于问题二的预测而言，通过对各个指标的分别建模，归纳出最适合的模型，最后将各个指标进行汇总而对总体进行预测。

5 模型的建立与求解

5.1问题一的求解

通过数据，可以直接得出人口总数为385310人；性别比（男：女）为99.809；未成年人数（0-14岁）为43456人，劳动龄人口（15-60岁）为264920人，老年人口（60岁以上）为76934人。

图1

（人口老龄化是指总人口中因年轻人口数量减少、年长人口数量增加而导致的老年人口比例相应增长的动态。两个含义：一是指老年人口相对增多，在总人口中所占比例不断上升的过程；二是指社会人口结构呈现老年状态，进入老龄化社会。国际上通常看法是，当一个国家或地区60岁以上老年人口占人口总数的10%，或65岁以上老年人口占人口总数的7%，即意味着这个国家或地区的人口处于老龄化社会。）因此我们发现，该县的老年人占总人口的20%，该地区已经进入老龄化。

（社会抚养比是指人口中非劳动年龄人口数对劳动年龄人口数之比 ,是衡量社会人均劳动年龄人口的抚养负担的指标。）

（社会扶养比 = (17岁及以下人口数 +6 3岁以上人口数 ) /（ 18～ 6 2岁人口数）= 少年抚养比与老年抚养比之和。）

社会抚养比=（58026+65497）/261787=0.471845

5.2问题二的求解

以下建立适用于各个子对象的通用模型。男、女性分别用下标m，w表示；

将第1

t 年初的人口分为在第t年末出生的人口与第t年末未死亡未迁出的非新生

人口两类，即：

(1)()()born other P t P t P t +=+ 设有年龄结构向量：

01n ()[(),(),,()]T P t P t P t P t =

其中，()i P t 为第t 年初年龄为i 的该类地区人数 5.2.1求第t 年末出生的人口数()born P t

设()i t μ为第t 年末年龄为i 的人的死亡率，()i h t 表示各年龄生育的女性占总生育女性的概率分布（即生育模式），()t β为第t 年末平均每个育龄女性的生育数，可以得到：

第t 年末年龄为i 的妇女人数为,,()[1()]w i w i P t t μ-

第t 年末年龄为i 的妇女生下的婴儿数量为,,()[1()]()()w i w i P t t h i t μβ- 5.2.2求第t 年末未死亡未迁出的非新生人口数()other P t

设从农村到城镇迁移率、从农村到城市迁移率分别为

t i ()t i ti t λ=第年末从农村迁往城镇的年龄为的人数

第年末年龄为的农村人数

t i ()t i ci t λ=

第年末从农村迁往城市的年龄为的人数

第年末年龄为的农村人数

则第t 年末未死亡未迁出的非新生人口数可表示为：

111,222000001()()()000()01()()()0

01()()()0()

c t c c other c t c cn tn cn c t t t P t t t t t t t P t μλλμλλμλλ?

???---?

???=---?

??

???---?

?

? 对该县人口，未死亡未迁出的非新生人口包括该县原来的剩下人口加上从外地迁移至该县的人口：

11

,220

00000

00001()000()000()01()0

0()0()00()1()0()0t t t other t t t v tn tn t t P t t P t t P t t t μλμλμλ????

????-?

???

????=-+?

??

??

??

?????-?

???

5.2.3死亡率模型

记死亡率()i t i t t i μ=

第年末年龄为的死亡人数

第年末年龄为的总人数

（1）年龄对死亡率的影响

众所周知，在同一时期的人口中，少年儿童的死亡率随年龄增长而下降，中青年人的死亡率变化较平稳，老年人的死亡率随年龄增长快速升高；总体来说，死亡率-年龄曲线为一“U ”字型。因此，考虑利用Kannisto 模型（文献[11]）通过数据拟合对死亡率随年龄增长的变化进行描述：

1(1)

i

i i

e e γγφμφ=+- 其中，i μ为年龄为i 的人的死亡率，φ和γ为参数，通过数据拟合得到

长期内， 由于()i t μ不断降低，最终将会趋于一个理想值，我们利用SIS 模型对其进行描述，体现制约与控制因素。有：

11()

()[1()]()d t t t t dt

μλμμγμ=?--?，0(0)μμ= 其中，1λ，1γ为参数，()i t μ的终值为11

1

λγμλ∞-=

5.2.4出生人口模型

出生人口的影响因素如图所示：

图8

生育人数对年龄符合一定的概率分布，每个年龄段的生育率归结为总和生育率与生育概率的乘积。这样对于生育率的预测，归结为总和生育率的预测上。

以下分别讨论各个子模型。 5.2.5生育模式模型

记()i h t 为各年龄生育的女性占总生育女性的概率分布，即生育模式（生育加权因子）。

由于生育模式主要受生理因素影响，假设其不随时间变化。考虑到生育的女性中中青年人口占较大比例的特点，我们设不同年龄女性的生育加权因子满足Γ分布：

1

11

11111()()

i i i i i e h αθαθα--

--=Γ，1i i >

5.2.6平均生育数模型

记平均生育数()t β=

每年总出生婴儿数

每年总生育女性数

通过分析数据，我们发现从2001年到2005年平均生育数上下波动，但总体来说有所减小。

出生人口

男女性别比k

生育性别比σ

平均生育数

β

生育模式i h

生育数

（1）中短期内，采用如下线性回归方程：

3()(0)t d t ββ=-

根据各年平均生育数数据，利用软件SAS 得到如下结果： 对人口：()0.002298 2.053c t t β=-

（2）长期内，平均生育数下降必将趋近一个定值。我们采用SIS 模型：

22()

()[1()]()d t t t t dt

βλββγβ=?--?，0(0)ββ= 终值为22

2

λμβλ∞-= 5.2.7迁移模型

记迁移率t i ()t i i t λ=

第年末该类地区迁移的年龄为的人数

第年该类地区末年龄为的总人数

()t t λ、()c t λ分别是从村到镇、从村到城的迁移率。

(),(),()c t v A t A t A t 分别为城、镇、乡外来人口比重。

假设迁移率没有性别差异。从数据得出迁移率的性别区别不明显，而且随着迁出的人在城镇的工作多元化，迁移率的性别区别将会减小。 时间的影响：

迁移率长期发展情况将是先增大后减小。这是由于当前城镇发展速度高于农村，所以中短期内迁移率将增加；但是，经济发展最终将消灭城乡差别，外来人口比重最终将变为零。

中短期内，外来人口比重变化速率改变不大，假设其随时间线性增加，通过线性回归得到如下结果：

()0.002380.006694c A t t =+ ()0.00110.01113t A t t =+ ()0.00652-0.01782v A t t =+

三个线性回归方程均通过显著性检验。

长期内外来人口比重减小并最终趋于零，我们采用SIS 模型描述其变化：

33()

()[1()]()dA t A t A t A t dt

λγ=?--， 0(0)A A =且33λγ≤ 5.2.8 老龄化程度

文献指出，60岁以上的人属于老年人。设第t 年末老年人总数为()old P t ，则

90

90

90

i=60

i=60

i=60

()()()()old ci ti vi P t P t P t P t =++∑∑∑

老龄化程度描述老年人口在总人口中的比例，所以老龄化程度表示为

()

()()

old P t O t P t =

其中，()P t 为第t 年末该县的总人口数。 5.2.9性别比

第t 年末的公民性别比为男性公民总人数与女性公民总人数之比，可表示为

()

()()

m w P t G t P t =

其中()m P t 、()w P t 分别为第t 年末男性公民、女性公民总人数。 最后对此进行预测：

人口发展方程：

()(),,p p r t Q r t r t

μ??+=-+?? ()()0,0p r p r =

()()()()0,p t f t u t N t ==

()()()()2

1

(),,,r r f t t h r t k r t p r t dr β=?

在进行数值计算时，将连续的人口模型离散化：

()i t N -------在t 年龄满i 周岁但不满1i +周岁的人口总数

()()i+1

i i

t =,N P r t dr ?（i =0，1，…….m。m 为最高年龄）

第一步：

用Excel 做出2000年到2013年人口总数的散点图如图所示，发现此散点图趋近于皮尔曲线，(资料曲线初期运行平缓，后发展迅速，达到一定程度，增长速度降低，终至平坦，渐

接高限的渐近线。)所以选用皮尔曲线作为预测模型：

^

^

1

t

t

K ab Y =+

第二步：

利用分段总和法估计模型的参数，a,K ,b

分段总和法的基本原理

此法是将原资料均分为三部分(资料项数不能被 3 整除时，将早期的一期或一、二期 Xt 剔除)，分别求取各部分总和，再据此计算参数估计值、a,K ,b ，建立预测方程，进行外推预测。

下面给出分段总和法的一般计算公式：

设：时间序列首、中、尾三段各期观察值的和分别为1

1t Y ∑、21t Y ∑、31

t

Y ∑表示。则参数的计算公式为：

3221212

11111111()(1)111()1t t m

t t

m t t m t Y Y

b Y Y b a Y Y b b k a m Y b ?-??=?

-?

??-=-?-?

???-?=-??

-??????∑∑∑∑∑∑∑

式中：m 为每一段中所包含的资料的项数，即m=N/3

利用附件2中给出的数据，用excel 计算得出

b=1.472,a=-8.23638E-10;k=2.62714E-06

第三步：将参数的估计值代入预测模型。通过excel 得到预测结果如下图所示

从表中可看出；2002年到2013年的预测值与附件中所给的值非常接近；所以我们可认为用该方法预测的结果是有效的。从而估计2016年的人口规模为409399.6人，2017年的人口规模为424535.9人，2018年的人口规模为448968.6人，2019年的人口规模为490521人，2020年的人口规模为567881.9人

从数据分析中得知，2000年-2013年平均人口为382578人；平均出生人数为2536人，平均死亡人数为2729人，其中大部分年数下男性死亡人数都要比女性死亡人数多；净迁入人口平均数为546人，2000年-2013年净迁入人口为7467人。 5.3问题三的求解

假设文盲的人数不变，先分别求出2010-2013年个学历的人数及比例。 在excel 中，就2011年的数据求解为例。 小学学历人数I4=I3-F3+F4,

初中学历人数N4=N3-K3-L3*(1-M3/100)+K4+L4*(1-M4/100)

高中学历人数U4=U3-P3-Q3*(1-R3/100)-S3-T3+P4+Q4*(1-R4/100)+S4+T4 大学及以上人数W4=W3+Q4*R4/100

这样就能求得2010-2013这四年个学历的人数。

灰色系统预测模型以微分方程为表述形式，揭示人数变化的连续过程。基本原理是：摈弃直接在历史数据中寻求统计规律和概率分布的传统方法，将无规律的原始数据通过一定的处理方式(比如1次或多次累加)，使其成为较有规律的时间序列，建立预测模型。数据的检验与处理：级比检验

以小学在校学生为例，原始序列如下式表示：

()()()(){}0

0000

=

1234x x x x X ，，，

带入具体数据得到:

{}0

=17542169931647515804X ，，，

求级比()k λ：

()()()0

k 1k =

k X X λ-

()()()()()=2345λλλλλ，，，

根据级比检验要求，只有当()k λ均落在可容覆盖2

21

2

(,)n n e

e

-

++内时，原始序列(0)X 可

以作为模型GM （1,1）进行预测。根据计算，所有的级比均满足()()0.9216451.128426k λ∈　，落在可容覆盖（0.8187,1.1993）内，则原始序列可以进行预测。 GM （1,1）模型的建立

（1）对原始序列(0)X 进行一次累加

(1)(1)(1)(1){(1),(2),,(4)}X x x x =

代入具体数据得到：

(1)

(17542249653588545791)X

= 其中，(1)

()()k

i x k x i ==∑。

（2）构造数据矩阵B 及数据向量Y

(1)(1)

(1)(1)(1)(1)

1((1)(2)) 21(()()) 2 1((4)(5)) 2x x x x B x x ??-+ ?

?

?-+ ?= ?

? ?-+ ???

１２３ １ １ (0)(0)(0)(2)(3)(5)x x Y x ?? ? ?

= ? ? ?

??

（3）计算参数向量 u

10.0446742[ ]()12225T T T u a b B B B Y -??=== ???

　

于是得到0.0446742a =，12225b =。其中a 为发展系数，b 为灰作用量。

（4）建立模型

根据(0)()x k 的GM （1,1）白化形式的微分方程为：

(1)

(1)dx ax b dt

+= 将,a b 值代入上式（5）得到：

(1)

(1)(0.0446742)12225dx x dt

+= 求解得：

(1)(0)0.0446742(1)((1))26098.927420.1ak k b b

x k x e e a a --+=-+=-+

（5）求生成数列值 (1)(1)x

k +及模型还原值 (0)(1)x k + 由于满足下式（7）：

(0)(1)(1)(1)(1)()x k x k x k +=+- 1,2,,4k =

我们令 (1)(0)(0)(1)(1)(1)13212x

x x ===，得： (0)(0)(0)(0)((1),(2),,(4))x

x x x =

（6）模型检验

我们通过设定残差值来评价原始值与模型值的精确度，下面给出残差的定义： |-|

=

原始值模型值残差原始值

预测结果：利用Matlab 编程计算得到2014年-2020年在校学生人数（其余的见附件二）

年份 小学在校学生 2014 15371 2015 14887 2016 14653 2017 14526 2018 14298 2019 13987 2020

13896

5.4问题四的求解

用eviews 做出从业人员随年份变化的散点图和线图，可看出该时间序列时间序列具有长期且呈直线型。因此可采用线性二次移动平均法对2016年到2020年的从业人员进行估计。预测的模型为:^

t T t t X a b T +=+?

计算步骤：

第一步：计算一次移动平均数和二次移动平均数

()111

t

t t n t X X X M n

--++++= ()

()()()111211t t t n t

M M M M n

--++++=

第二步：计算t a ，t b 的值

()(

)

1

22t t t a M M =-

()()()

122

1

t t t b M M n =

-- 第三步：根据计算出来的参数值进行估算。

为了避免由于使用不同时刻

t a ，t b 的值进行预算所导致的预测结果不同，我们取

n=6，通过使用Excel 得到的预测结果如下图所示。

通过线性二次移动平均法得到的结果为：2016年的从业总数为257312.2，2017年的从业总数为257655，2018年的从业总数为257997.8，2019年的从业总数为258340.6，2020年的从业总数为258683.3，

利用GM(1，1)模型预测2016-2020年的从业人员，附件4关于从业人员的时间序列

()0X 有n=11个观察值，()()()()()()()()(){}

000001,2,3.....11X X X X X =，通过累加生成新序

列。()()()()()()()()(){}

111111,2,3.....11X X X X X =，则GM(1，1)模型相应的微分方程为：

()

()1

1dX aX dt

μ+=；a 为发展灰度，μ为内生控制灰度。 设^

a 为待估参数向量，()^

,T

a a μ=。利用最小二乘法求解可得：()

^

1

T

T n a B B

B Y -=

其中：()()()()()()()()()()()()1111111

12121

2312

1101112

X X B X X X X ????

?-+?

? ?

???

=-+ ??? ? ? ???-+ ?????

;()()()()()

()()

0002,311T n Y X X X = 求解微分方程，即可得预测模型：

()()()()^

1011ak X k X e a a μμ-?

?+=-+???

?

6 模型的优缺点

6.1模型的优点

1.可操作性

引入了影响死亡率，和生育率的子因素，将人口总量的变化转化到这些子因素的变化上。子因素的数据更易获取。

2.健壮性

将死亡率转化为死亡分布率和死亡减小率，这样将死亡变化转化为总和死亡率的变化上，仅仅预测这一个参数就可以得到死亡率，这样大大地加强了模型的强壮性。将生育率转化为总和生育率和生育率分布的乘积，这样把生育率这个随时间和年龄变化的量转化为仅随时间变化的量，并在短期内通过拟合得到变化规律，并在长期内通过一个经典的非线性模型来描述。这样增加了模型的健壮性。

3.可控性

在做长期参数预测时，将变化规律用参数来描述，这些参数体现了国家政策的控制，这样通过这些参数的变化就可以体现人口总题的变化趋势。这样，就可以根据国家对人口的要求制定不同的政策。

6. 2模型的缺点

除了上面的优点还有如下的缺点：

1.模型中是基于概率分布假设的，但对于概率分布的形式很难预测，我们只能对其进行大致的拟合，但可能并不是实际情形。死亡率的拟合函数参数值稳定性太差。出生模式的拟合曲线不能完全符合实际情况，迁移模式的拟合曲线太理想化。

2.总和生育率不能完全准确的描述，只能用近似形式的曲线来描述。可能实际总和生育率的变化并不是那个形式。

3.并没有对90岁以外的人进行预测研究。

参考文献

[1] 韩中庚，数学建模方法及其应用（第二版），北京：高等教育出版社，2009。

[2] 全国大学生数学建模竞赛组委会，高教社杯全国大学生数学建模竞赛论文格式规范，北京，2009.

[3] 姜启源,谢金星,叶俊.数学模型（第三版）[M].北京：高等教育出版社,2003

[4] 王学民.应用概率统计[M].上海：上海财经大学出版社，2005

[5] 梅长林,范金城.数据分析方法[M].北京：高等教育出版社，2006

[6] Chahnazarian :Determinants of the Sex Ratio at Birth :Review of Recent Literature, Social Ecology 112Win2ter ;35 (3 —4) :214 —235.

[7] Poston D L , Fertility and family planning in China : an analysis of provincial patterns, Journal of BiosocialScience, 1988 .20 .4 Oct ,379 —391.

[8] Arnold ,F. et al ,Sex preference ,fertility ,and family Planning in China ,Population and Dvelopment Review ,1986 ,12 ,2 ,J une ,221 —246.

附录1

clc,clear

x0=[17542 16993 16475 15804];

n=length(x0);

lamda=x0(1:n-1)./x0(2:n)

range=minmax(lamda)

x1=cumsum(x0)

for i=2:n

z(i)=0.5*(x1(i)+x1(i-1));

end

B=[-z(2:n)',ones(n-1,1)];

Y=x0(2:n)';

u=B\Y

x=dsolve('Dx+a*x=b','x(0)=x0');

x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});

yuce1=subs(x,'t',[0:n-1]);

digits(6),y=vpa(x) %为提高预测精度，先计算预测值，再显示微分方程的解yuce=[x0(1),diff(yuce1)]

epsilon=x0-yuce %计算残差

delta=abs(epsilon./x0) %计算相对误差

rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda %计算级比偏差值

附件二

教育状况文盲小学

年份人口总数人数

比例

（%）在校学生毕业生

升学率

（%）学历人数比例（%）

2010 385310 24739 6.42 17542 3395 100 104830 27.21 2011 385878 24739 6.41 16993 3062 100 104281 27.02 2012 386178 24739 6.41 16475 3067 100 103763 26.87 2013 386064 24739 6.41 15804 2806 100 103092 26.70 2014 393397 24739 6.29 15371 2822 100 102659 26.10 2015 399717 24739 6.19 14887 2794 100 102175 25.56 2016 409400 24739 6.04 14653 2666 100 101941 24.90 2017 424536 24739 5.83 14526 2538 100 101814 23.98 2018 448969 24739 5.51 14298 2510 100 101586 22.63 2019 490521 24739 5.04 13987 2482 100 101275 20.65 2020 567882 24739 4.36 13896 2354 100 101184 17.82