移动加权平均修匀法端值数据外推的方法及其应用①

孙佳美 霍晶昌 赵星

（南开大学，天津 300071）

    ［摘要］随着人口老龄化问题的加剧，研究高龄人口的死亡率变化对保险公司正确计算退休后的养老金给付显得尤为重要。然而，高年龄段的投保人数一般较少，因此高年龄段死亡率的初始估计是不确切的，必须利用初始估计，结合先验观点对初始估计值进行修正，这种方法即为死亡率修匀。移动加权平均修匀方法是最早的非参数修匀方法。它的缺陷是，除非能够获得初始观察值以外的数据，否则前m个和后m个观察值的修匀值无法直接由基本公式得出。本文针对MWA修匀的上述缺陷，给出端值数据修匀的两种方法——矩阵向量法和外推法，着重给出这两种方法的计算机实现过程。最后以日本1996年生命表的编制为例，说明其中MWA修匀的过程。
    ［关键词］MWA修匀；端值数据；修匀矩阵；外推系数

    一、移动加权平均修匀法的基本公式及缺陷
    移动加权平均修匀 (MovingWeightedAverage Smoothing，也称MWA修匀) 法是一种经典的数据修匀方法，是DeForest于1870年提出的。该方法的基本公式为：ux=∑mj=-mcjyx+j（1）
其中cj=c-j，∑mj=-mcj=1，yx为初始观察值，ux为修匀值，m为给定的正整数，每一项的修匀值为原观察值的前后连续共2m+1项的加权和。
    公式（1）中的修匀系数cj一般由下述两个条件确定：
    1．ux再生三次多项式；
    2．极小化R2Z=Var(Δzux)Var(Δzyx)。
    特别地，当z=3时，如果令k=m+2，则得到2k-3项对称MWA修匀的修匀系数cj，结果如下：
cj=315［(k-1)2-j2］(k2-j2)［(k+1)2-j2］(3k2-16-11j2)8k(k2-1)(4k2-1)(4k2-9)(4k2-25)（2）
从上面的公式（1）看出，除非能够获得初始观察值以外的数据，否则前m个和后m个观察值的修匀值无法直接利用该公式得出。
由于MWA修匀公式的上述缺陷，很多学者在相关著作中对端值数据修匀的方法做出论述，但都没有取得一般共识。
Greville（1981）在MWA修匀过程中对端值数据处理方面进行了深入的研究，提出了一种矩阵向量法，这种方法可以记为u=Gy，其中y为观察值向量，u为相应的修匀值向量，G为n阶修匀矩阵。对于给定的MWA修匀公式，G是唯一确定的。
Greville还提出了另一种处理端值数据修匀的方法，即在观察值数据的两端各添加m个添加值，然后在
［作者简介］孙佳美，南开大学风险管理与保险系副教授；霍晶昌，南开大学风险管理与保险系2007级硕士研究生；赵星，南开大学风险管理与保险系2008级博士研究生。
所有数据范围内应用给定的MWA修匀公式，得到观察数据的全部修匀值。矩阵法和外推法并非不相关，实际上外推法的最终修匀结果也可以用矩阵表示。Greville（1981）已证明，对于同一组观察死亡率数据，矩阵法和外推法得到的最终修匀结果是一致的。
本文主要介绍上述两种端值数据修匀的方法，即矩阵向量法和外推法，最后以日本1996年生命表编制过程中的MWA修匀过程为例，说明如何用外推法对端值数据进行修匀。本文的重点是对端值数据修匀的两种方法，即矩阵向量法和外推法的具体算法给出说明，并结合具体例子给出两种方法的应用。
二、端值数据处理的矩阵向量法
T.N.Greville（1981）根据Whittaker修匀过程进行类推，并结合相关的数学概念，在MWA修匀的基础上提出一种用矩阵向量法对端值数据进行修匀的理论。矩阵向量法的关键是修匀矩阵G的确定，然后由u=Gy即可得到各个观察值的修匀值。
在介绍端值数据处理的矩阵向量法之前，首先注意到，如果对称MWA修匀公式（1）的精度是2s-1，则公式（1）可以改写为如下形式：ux=［1-(-1)sδ2sq(E)］yx（3）
其中yx为初始观察值，ux为修匀值，δ和E分别为中心差分算子和位移算子，q(E)=∑m-sj=-m+sqjEj，q-j=qj。
为了描述唯一的修匀矩阵G，我们进行如下定义：
定义1：对于一个N阶方阵M=(mij)Ni, j=0，如果存在非负整数h和k，当j-i>h并且i-j>k时，有mij=0，则矩阵M称为带状矩阵（Band Matrix）。注意，这里矩阵M的行和列都是从0开始标记的。如果再加上约束条件h+k≤N，则M称为严格的带状矩阵（Strictly Band Matrix）。如果mi+1, j+1=mij，(i, j=0,1,L,N-1)，M是个托伯利兹（Toeplitz）矩阵。
定义2：一个严格的带状矩阵如果满足下面条件时，就称作特伦特矩阵（Trench Matrix）。
令Mi(x)表示矩阵M第i行元素的母函数，Mi(x)=∑Nj=0mijxj，则当
Mi(x)=xiA(x)∑jj=0βjx-j,0≤i<k
xiA(x)B(1/x),k≤i≤N-h
xiB(1/x)∑Nj=N-iα jxj,N-h<i≤N
时，M是一个特伦特矩阵。其中A(x)=∑hj=0αjxj，B(x)=∑kj=0βjxj分别是次数为h和 k的多项式，并且α0β0≠0。除非特别说明，本文中所用到的特伦特矩阵都满足h=k并且αj和βj均为实数。
当且仅当A(x)和B(1/x)在复平面上没有公共的根时，特伦特矩阵（Trench Matrix）为非奇异的。而非奇异的特伦特矩阵（Trench Matrix）存在托伯利兹（Toeplitz）逆。
定义3：对于一个(N-s)×N阶的矩阵K，如果能把一个列向量转化为它的元素的s阶有限差分的形式，那么称矩阵K为差分矩阵。显然，K的每一行的非零元素是s阶差分的系数并改变符号，K的所有非零元素组成从左上到右下的对角带。
定理1：对于给定的2m+1项MWA修匀公式及任意的 N≥2m+1（N为初始观察值的个数），存在唯一的N阶矩阵G，且G具有如下五个性质：
（1）如果y是由一系列观察值组成的向量，并且u=Gy是相应的修匀向量，那么向量u的元素（除了前m项和后m项）可直接由修匀公式（1）计算得到；
（2）G是严格带状矩阵，并且h=k=m；
（3）G有如下形式：G=I-KTDK，D为(N-s)阶方阵，上标T表示矩阵转置；
（4）D为非奇异矩阵，并且存在托伯利兹（Toeplitz）逆矩阵；

① 本论文受到“南开大学2008年度人文社会科学校内文科青年项目（项目编号：NKQ08008）”的资助。（5）如果D-1=(ti j)，ti j=ti-j ，那么数列∑∞v=-∞tvzv在复平面的某个部分收敛（根据性质1-4可得，tv仅与v相关，与N相互独立）。
根据上述5个性质唯一确定的G具有如下性质：
（6）D是唯一的，并且G和D是对称矩阵；
（7）矩阵G的前m行和后m行非零元素仅由给定的MWA公式决定，与N无关；
（8）D为B(x)=A(x)条件下的Trench矩阵。并且有q(x)=A(x)A(1/x)，A(x)=0的m-s个解正好是方程q(x)=0的位于复平面单位圆外的解；
（9）F=I-G=KTDK是一个奇异的特伦特矩阵，可以用A＾(x)=B＾(x)=(x-1)sA(x)表示；
（10）∑∞v=-∞tvzv收敛于［q(z)］-1，在复平面中表示为包含单位圆的环形区域。
本文将略去定理1的具体证明，而只给出应用定理1及性质得到修匀矩阵G的方法。
假设m=6，对应的修匀系数cj由公式（2）得出，即：
可以将ux改写为公式（3）的形式（ 其中s=2），并且：
q(x)=0.0193498x-4+0.1052631x-3+0.30495356x-2+0.60014289x-1+0.8251964x0
+0.60014289x1+0.30495356x2+0.1052631x3+0.0193498x4（4）
利用计算机软件将q(x)进行如下分解：
q(x)=A(x)A(1/x)，其中A(x)=0的解为q(x)=0的位负平面单位圆之外的解。
则有A(x)=0.58617+0.57661x+0.35552x2+0.14711x3+0.03301x4
进一步地，A＾(x)=(x-1)2A(x)=0.58617-0.59572x-0.21154x2+0.012676x3+0.09432x4
+0.081085x5+0.033011x6
根据F=I-G=KTDK是一个奇异的特伦特矩阵（其中A＾(x)=B＾(x)=(x-1)sA(x)）及定义2得到特伦特矩阵F的前6行元素为：
0.343 6-0.349-0.1240.007 40.055 30.047 50.019 4000000…0
-0.3490.698 5-0.223-0.132-0.0490.0070.027 90.019 400000…0
-0.124-0.2230.743 2-0.226-0.151-0.06600.027 90.019 40000…0
0.007 4-0.132-0.2260.743 4-0.225-0.15-0.06500.027 90.019 4000…0
0.055 3-0.049-0.151-0.2250.752 3-0.217-0.147-0.06500.027 90.019 400…0
0.047 50.007-0.066-0.15-0.2170.758 9-0.214-0.147-0.06500.027 90.019 40…0G=I-F再利用，得到修匀矩阵G的前6行元素为：
0.656 40.349 20.124-0.007-0.055-0.048-0.019000000…0
0.349 20.301 50.223 20.131 50.048 8-0.007-0.028-0.01900000…0
0.1240.223 20.256 80.225 90.151 50.065 90-0.028-0.0190000…0
-0.0070.131 50.225 90.256 60.224 70.150 50.065 50-0.028-0.019000…0
-0.0550.048 80.151 50.224 70.247 70.2170.147 40.065 50-0.028-0.01900…0
-0.048-0.0070.065 90.150 50.2170.241 10.214 30.147 40.065 50-0.028-0.0190…0用同样的方法可得到修匀矩阵G的后6行元素。修匀矩阵G的其他行的非零元素为公式（2）对应的修匀系数c j，并使得G为一个N×N阶带状矩阵。
只要当N≥13，就可以根据u=Gy得到所有的修匀值。
三、端值数据处理的外推法
外推法是在观察值数据的两端各添加m个添加值，然后在所有数据范围内应用给定的MWA修匀公式，得到全部观察数据的MWA修匀值。外推法的关键是得到外推系数，由外推系数可以得到第一个观察值前面的m个添加值及最后一个观察值后面的m个添加值，再利用公式（1）就得到全部观察数据的修匀值。
下面从计算的角度简单说明端值数据修匀的外推法。假设给定一组从x=P到x=Q的观察死亡率数据，N=Q-P+1。同矩阵向量法一样，对于MWA修匀公式（1），首先把修匀公式（1）改写为公式（3）的形式，得到q(z)。定义p(z)为m-s次多项式，最高次幂的系数为1，p(z)=0的m-s个解为q(z)=0的所有单位圆内的解。再定义m次多项式a(z)，其系数a j由以下公式确定：
a(z)=(z-1)sp(z)=zm-∑mj=1ajzm-j(5)
连续利用yx=∑mj=1ajyx+j(x=P-1,P-2,…,P-m)可以得到第一个观察数据前面的m个添加值。同样利用yx=∑mj=1ajyx+j(x=Q+1,Q+2，…,Q+m)可以得到最后一个观察数据后面的m个添加值。最后，在所有数据范围（包括2m个添加值）内对观察数据应用MWA修匀公式（1），即得到全部观察数据的MWA修匀值。
例如，取m=6，采用前面的方法得到q(z)如公式（4）所示，其中s=2。取q(z)=0的所有在单位圆内的根组成多项式p(z)，则有：p(z)=0.056 316 54+0.250 963z+0.606 519z2+0.983 699z1+z4
再根据公式（5）及s=2，得到：
a(z)=(z-1)sp(z)=z6-1.016 300 692z5-0.360 879 5z4+0.021 624 5z3+0.160 908 7 z2+0.138 330 4z+
0.056 316 5(6)
比较（5）式和（6）式，得到m=6时13项MWA修匀公式对应的外推系数如下：
j123456aj1.016 3010.360 88-0.021 624 6-0.160 908 7-0.138 330 4-0.056 316 5连续利用yx=∑6j-1ajyx+j(x=P-1,P-2,…,P-6)可以得到第一个观察值前面的6个添加值。同样利用yx=∑6j-1ajyx-j(x=Q+1,Q+2，…,Q+6)可以得到观察数据后面的6个添加值。利用这12个添加值及原始观察数据和修匀公式（1）即得到所有观察数据的修匀值。
表1是极小化R3时移动加权平均修匀的修匀系数和外推系数。
四、应用举例
由于数据收集的限制，本文以日本1996年保险业生命表的编制过程为例说明其中MWA修匀的过程。在日本1996年生命表的编制过程中，对死亡率一共进行了三次修匀，其中第二次修匀就是采用的MWA修匀，对端值数据的修匀则采用了前面所述的外推法。
极小化R3时移动加权平均修匀（精度为3）的修匀系数cj和外推系数aj
表1
2m+15791113jcjajcjajcjajcjajcjaj00.5594410.4125870.3311390.2779450.24005710.2937062.0000000.2937061.6180420.2665571.3526130.2386931.1608110.2143371.0163012-0.073427-1.0000000.058741-0.2360730.1184700.1146970.1412670.2810790.1473570.3608803-0.058741-0.381969-0.009872-0.2872310.035723-0.1409680.065492-0.0216254-0.040724-0.180078-0.026792-0.2045460-0.1609095-0.027864-0.096377-0.027864-0.1383306-0.019350-0.0563172m+11517192123jcjajcjajcjajcjajcjaj00.2115410.1892310.1712660.1564690.14406010.1937420.9036650.1763900.8134420.1616910.7395860.1491360.6780000.1383180.62587920.1459040.3972960.1411120.4108850.1349650.4120920.1284230.4064950.1219490.39720630.0829180.0647500.0922930.1249330.0966580.1661610.0979560.1940250.0973950.21250140.024027-0.1007120.042093-0.0434560.0546850.0050950.0630380.0443140.0683030.0752375-0.014134-0.1354460.002467-0.1106440.017475-0.0782570.029628-0.0454380.038933-0.0153126-0.024499-0.094424-0.018639-0.106212-0.008155-0.0999740.003119-0.0840200.013430-0.0639277-0.013730-0.035128-0.020370-0.065896-0.018972-0.081844-0.012896-0.084711-0.004948-0.0787378-0.009960-0.023052-0.016601-0.047103-0.017614-0.063086-0.014527-0.0700639-0.007378-0.015756-0.013455-0.034444-0.015687-0.04897710-0.005570-0.011134-0.010918-0.02571411-0.004278-0.008092以1996年日本生命表男表的编制过程为例，第二次修匀采用了m=6的MWA修匀。由于粗死亡率的数据从x=0到x=78，因此应用外推法前需要将第一次修匀的结果向前外推至-6岁，向后外推至84岁，然后直接利用13项MWA修匀公式即得0岁—78岁粗死亡率的第二次修匀值。限于篇幅限制，本文只列出部分数据，主要说明外推的过程，如表2所示。直接利用前面给出的修匀矩阵G，并利用u=Gy同样可以得到与表2最后一列相同的修匀结果。
日本1996年生命表男表的前两次修匀过程
表2
 年龄粗死亡率第一次补整死亡率外推所需全部数据第二次修匀死亡率-60.003 426-50.003 055-40.002 681-30.002 298-20.001 9-10.001 49300.001 1450.001 490.001 490.001 110.000 1590.000 210.000 210.000 75820.000 5510.000 720.000 720.000 49630.000 0910.000 120.000 120.000 32640.000 2590.000 340.000 340.000 24250.000 1710.000 220.000 220.000 21760.000 2050.000 270.000 270.000 21670.000 1870.000 240.000 240.000 21180.000 1340.000 170.000 170.000 1990.000 1370.000 180.000 180.000 171100.000 1110.000 140.000 140.000 148110.000 1230.000 160.000 160.000 139120.000 1290.000 170.000 170.000 155130.000 1550.000 200.000 200.000 217140.000 1950.000 250.000 250.000 338650.013 9530.015 300.015 300.015 447660.0152 430.016 720.016 720.016 739670.016 6280.018 260.018 260.018 223680.018 2040.020 020.020 020.019 988690.020 2240.022 250.022 250.022 09700.022 0180.024 270.024 270.024 57710.024 6240.027 160.027 160.027 404720.027 6380.030 510.030 510.030 529730.031 6620.034 960.034 960.034 004740.034 1010.037 770.037 770.037 945750.038 2290.042 410.042 410.042 467760.041 5620.046 260.046 260.047 632770.047 4510.052 860.052 860.053 406780.054 8290.061 120.061 120.059 66790.066 173800.072 727810.079 178820.085 532830.091 848840.098 201［参考文献］
［1］Greville,T.N.E.(1981).Movingweightedaverage smoothing extended to the extremities of the data.Ⅰ.Theory.Scand.Actuarial J.1981.39-56
［2］Greville,T.N.E.(1981).Movingweightedaverage smoothing extended to the extremities of the data.Ⅱ.Methods .Scand.Actuarial J.1981.65-81
［3］Greville,T.N.E.(1981).Movingweightedaverage smoothing extended to the extremities of the data.Ⅲ.Stability and optimal properties .Journal of Approximation theory 33, 43-58
［4］贾利新 王爱兰. 常系数非齐次线性差分方程的算子解法, 数学的实践与认识. 2000.10.(30-4).485-490

Abstract： With the accelerated aging of the population, it becomes more important for  insurance companies to study the mortality rate of the elderly population, because the calculation of retirement pensions depend on the mortality rate of the elderly population. However, the number of elderly insureds are often small, so the original estimates of highage mortality rate is inaccurate, therefose they must be smoothed. The use    of a symmetrical moving weighted average of terms to smooth equally spaced observations of a function of one variable does not yield smoothed value of the first m and the last m observations, unless additional data beyond the range of the original observations are available. This article gives two methods for smoothing the extremities of the data which is the matrix approach and the  extrapolation approach, especially focusing on the process of realization by computer software. Finally, takeing Japan′s life insurance industry mortality 1996 for example, the MWA smoothing process is interpreted.
Key words：MWA smoothing; extremities of the data; graduation matrix; extension coefficients
［编辑：施敏］保险研究2008年增刊2