3.2 时间序列的指数平滑预测法
指数平滑法(Expinential smoothing method)的思想也是对时间序列进行修匀以消除不规则和随机的扰动。该方法是建立在如下基础上的加权平均法:即认为时间序列中的近期数据对未来值的影响比早期数据对未来值得影响更大。于是通过对时间序列的数据进行加权处理,越是近期的数据,其权数越大;反之,权数就越小。这样就将数据修匀了,并反映出时间序列中对预测时点值的影响程度。根据修匀的要求,可以有一次、二次甚至三次指数平滑。
3.3.1 一次指数平滑法
1.一次指数平滑法的计算公式及平滑系数a的讨论
设时间序列为x1,x2x3,,xN,一次指数平滑数列的递推公式为:
11Staxt(1a)St1,1S0x1,
0a1,1tN
(3-6)
式中,St1表示第t时点的一次指数平滑值,a称为平滑系数。递推公式(3-6)中,初
1
始值S0常用时间序列的首项x1(适用于历史数据个数较多,如50个历史数据及以上),如
果历史数据个数较少,如在15或20个数据及以下时,可以选用最初几期历-史数据的平均
1值作为初始值S0,这些选择都有一定的经验性和主观性。
下面讨论平滑系数a。将递推公式(3-6)展开可得:
St1axt(1a)St11axt(1a)axt1(1a)St12
axta(1a)xt1(1a)2St12
1
axta(1a)xt1a(1a)2xt2a(1a)t1x1(1a)tS0
容易看出,由于0a1,xi的系数a(1a)随着i的增加而递减。注意到这些系数之和为1,即:
i
a(1a)
i1
t
i1
1(1a)t
(1a)a(1a)t1
1(1a)
t
于是,递推公式(3-6)中的St1就是样本值x1,x2,,xt的一个加权平均。当用递推公式(3-6)来进行预测时,我们将用St1作为第t1时点的预测值。从上面的讨论可以看到,离预测时点t1最近的时点t的值xt,其权为a,最大,其次为xt1的权a(1a),,x1的权最小。
可见,公式(3-6)是在认为新近数据对未来影响大,远期数据对未来影响小的情况下对原时间序列的加权平均(修匀)。
11
若平滑系数a0,此时有St1St11S0(x1),这表明确定S0(x1)后,每个
1
时点的平滑值皆等于S0,此时,每个时点i的观察值xi不起任何作用。
若平滑系数a1,此时有St1xt。说明平滑后数列St1就是原时间序列,即没有对原时间序列进行任何处理和修匀。
对于平滑系数a来说,除了上述两种极端情况外,不可能出现xi与xj(ij)系数相等的情况。
综上所述,不难看到,0a1且比较接近1时,计算得到的一次指数平滑值对原历史数据的修匀程度将较小,平滑后的数列值St1能够比较快地反映出原时间序列的实际变化,因此,对于变化较大或趋势性较强的时间序列适合选择比较靠近1的数据作为平滑系数,比如取a0.95,0.90等。
若0a1且取值比较靠近0时,计算得到的一次指数平滑值对原历史数据的修匀程度将较大,平滑后的数列对原时间序列的变化反应较迟钝。因此,对于变化较小的、或接近平稳的时间序列,应选择比较靠近0的平滑系数使得平滑过程中的各数据的权数比较接近。
【例3-7】某电器销售企业1992-2003年某种电器销售额(万元)及其一次指数平滑数列计算列表如表3-2所示。
表3-2 一次指数平滑值计算表 (单位:万元)
年份 销售额
1992 50 51 51 51
52 50.8 50.50 50.20
47 51.04 51.25 51.64
51 50.23 49.13 47.93
49 50.38 50.07 50.39
48 50.10 49.54 49.28
51 49.68 48.77 48.26
40 49.94 49.89 50.45
48 47.95 44.95 42.09
51 47.96 46.48 46.82
51 48.57 48.74 50.16
59 49.06 49.87 50.83
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
0.2
0.5
0.8
1
本例中选择第一、第二期的历史数据平均值作为一次指数平滑的初始值S0。从表3-2
中可以看到,原时间序列变动较大,稍具周期性。若想指数平滑后数列敏感地反映出最新的一些观察值的变动,则应取较大的平滑系数a,如a0.8。若想消除其周期性变动,施加较大程度的修匀,以反映其长期趋势,则可取较小的a,如a0.2。
2.检验及预测
一次指数平滑法适用于对变化比较平稳的时间序列作短期的预测。第t1时点的预测值yt1等于按递推公式(3-6)计算的第t时点的一次指数平滑值St1,即:
yt1St1 (3-7) 类似于移动平均法,同一个问题随着平滑系数的不同可以有若干个一次指数平滑预测值。因此,应该在一个合适的评价标准基础上选择一个合理的平滑系数a。其方法是:首先计算与在原则上合理的多个平滑系数a相对应的平滑数列,然后分别计算其均方差MSE(见公式(3-7))或计算其平均绝对误差MAD(见公式(3-8)),以MSE或MAD最小者对应的平滑系数及其预测值为最合理的。
1N
MADet,而etxtSt11xtyt (3-8)
Nt1
公式(3-8)中的etxtSt11xtyt反映了各个时点的平滑值St11与实际值xt之间的误差。根据公式(3-2),计算误差数列et的自相关系数rt,若统计量
2
Q(N1)rt2(m1)
t1m
则说明误差数列具有随机性,可以认为此时的预测是有效的(可以利用前面的Matlab
程序funcoef.m判断)。一次指数平滑法的计算、检验及预测的Matlab程序如下(文件名funesm1.m)。
【例3-8】某商品2004年各月的销售量xt如表3-3所示,试预测2005年1月的销售量。
1利用funesm1.m计算的部分结果如下,其他部分数据如表3-3所示。程序中令S0x1423,
取L00.10,L10.05,L30.95,m7,0.05,对应的误差值et,也列于表3-3中。
表3-3 历史数据、一次指数平滑值及误差值
月 序号
0 0
1 1 423
2 2 358
3 3 434
4 4 445
5 5 527
6 6 429
7 7 426
8 8 502
9 9 480
10 10 385
11 11 427
12 12 446
t
xt
St1(0.1)
423
423 0
416.5 -65
418.8 17.5
420.9 26.7
431.5 106.1
431.3 -2.5
430.8 -5.3
437.9 71.2
442.1 42.1
436.4 -57.1
435.5 -9.4
436.5 10.5
et(0.1)
St1(0.2)
423
423 0
410 414.8
420.8 30.2
442 439.4
436.7 -13.4
449.8 65.3
455.8 30.2
441.6 -70.8
438.7 -14.6 436.9 -14.1
420.2 7.3
et(0.2)
-65 24 106.2
-13
St1(0.3)
423
423 0
403.5 -65
412.7 30.5
422.4 32.3
453.8 104.6
446.4 -24.8
440.3 -20.4
458.8 61.7
465.2 21.2
441.1 -80.2
439.8 9.1
et(0.3)
这就是说,取平滑系数0.10是合适的,此时预测值为2005年1月的销售量将是y = 436.4976。
funesm1.m
%一次指数平滑预测,文件名为funesm1.m
%输入时间序列想,平滑系数初值L0,步长L1,终值L2
%输入判断误差是否为随机误差时需要计算的自相关系数个数m,显著性水平alpha function ESM1=funesm1(x,L0,L1,L2,m,alpha)
s=zeros(round((L2-L0)/L1),length(x)); e=zeros(1,length(x));
MAD=zeros(1,round((L2-L0)/L1)); k=0;
for a=L0:L1:L2 k=k+1
s(k,1)=x(1); for i=2:length(x)
s(k,i)=a*x(i)+(1-a)*s(k,i-1); e(i)=x(i)-s(k,i-1); end
s(k,:),e
funcoef(e,m,alpha); MAD(k)=mean(abs(e)); end
disp ('The smallest MAD and the corresponding a and forecast') MAD,[MAD,k]=min(MAD);
a=L0+L1*(k-1),y=s(k,length(x)) end
利用上述程序计算得到的一次指数平滑预测值是在多个平滑系数(可自己确定)中计算、比较得到的。
3.3.2 二次指数平滑法
对于呈现出线性趋势的时间序列,在一次指数平滑数列的基础上用同一个平滑系数a再进行一次指数平滑,就是二次指数平滑。构成二次指数平滑数列St2的递推公式如下:
21令初始值S0,则有: S0x1(也可取其他值作为初始值)
11
Staxt(1a)St1
2 (3-9) 12
StaSt(1a)St1
二次指数平滑的目的是对原时间序列进行两次修匀,使得其不规则变动或周期变动尽量
消除掉,让时间序列的长期趋势性更能显示出来。对于平滑系数,同样有一个合理选区的问题。其方法与一次指数平滑法一样,先选取原则上较合理的多个a值分别计算,得到不同的数列St1和St2,再根据均方差MSE或MAD最小原则确定较为合理的a值,并得到相应的二次指数平滑值。
由于二次指数平滑的目的,二次指数平滑法较适用于具有线性趋势的时间序列,线性趋势预测模型为:
ytTatbtT(3-10)
上式中,T表示自T时点起向前预测的时点数;at,bt满足:
atSt1(St1St2)2St1St2,bt
a
(St1St2)(3-11) 1a
预测模型的有效性检验方法与一次指数平滑法一样。即通过自相关系数或2检验方法进行检验。利用二次指数平滑法预测具有线性趋势性的时间序列的基本步骤为:
1. 根据历史数据(时间序列)x按照公式(3-6)计算一次指数平滑值; 2.根据公式(3-9)计算二次指数平滑值;
3. 由公式(3-11)计算at,bt,并由(3-10)计算自t时点起先前T时期的预测值ytT。据此,编写的Matlab程序如下(文件名funesm2.m)。其中,程序中的s1,s2分别表示一次、二次指数平滑值,a2,b2表示相对于每个平滑系数的at,bt(公式(3-11)),y表示线性趋势的预测值。平滑系数从L2(步长L1),然后,从众多平滑系数计算的结果中挑选最小的MAD 所对应的平滑系数并最后得到预测值,自t时点起先前的时期数T由预测值确定并输入。
funesm2.m
%二次指数平滑,线性趋势预测模型
%输入时间序列x,平滑系数初值L0,步长L1,终值L2
%输入判断误差是否为随机误差时必须计算的自相关系数个数m,显著性水平alpha function ESM2=funesm2(x,L0,L1,L2,m,alpha) T=input('T=')
s1=zeros(round((L2-L0)/L1),length(x)); s2=zeros(round((L2-L0)/L1),length(x)); a2=zeros(round((L2-L0)/L1),length(x)); b2=zeros(round((L2-L0)/L1),length(x)); y=zeros(round((L2-L0)/L1),length(x)+T); e2=zeros(round((L2-L0)/L1),length(x)); MAD2=zeros(1,round((L2-L0)/L1));k=0; for a=L0:L1:L2 k=k+1;
s1(k,1)=x(1); s2(k,1)=x(1); for i=2:length(x)
s1(k,i)=a*x(i)+(1-a)*s1(k,i-1); s2(k,i)=a*s1(k,i)+(1-a)*s2(k,i-1); a2(k,i)=2*s1(k,i)-s2(k,i);
b2(k,i)=(s1(k,i)-s2(k,i))*(a/(1-a)); y(k,i+T)=a2(k,i)+b2(k,i)*T; if i+T
e2(k,i+T)=x(i+T)-y(k,i+T); end end
's1:',s1(k,:),'s2:',s2(k,:),
'a2:',a2(k,:),'b2:',b2(k,:),'y:',y(k,:),'e2:',e2(k,:), funcoef(e2(k,:),m,alpha);
MAD2(k)=mean(abs(e2(k,:))); end
MAD2;[MAD2,k]=min(MAD2);
a=L0+L1*(k-1),y(k,length(x)+T) end
【例3-9】某市1990-1999年工业产值为xt(t1,2,,12),如表(3-4)所示,由于原时间序列大体上呈现增长的趋势,取平滑系数a0.9利用二次指数平滑法计算的结果如表(包括at,bt值以及误差值et)3-4所示。
表3-4中的at,bt值按公式(3-11)计算。于是2000年的预测值为:
y13.700.51114.21。
而2005年的预测值为:y13.700.51616.76
表3-4 历史数据、二次指数平滑值及误差值
年份
0 10.1 10.1
1990 1 10.1 10.1 10.1
10.64 10.59 10.69 0.49
11.14 11.09 11.20 0.50 11.18 0.02
-0.002
-0.101
-0.220
-0/343
0.50 11.70
0.42 12.20
0.24 12.52
-0.040 12.54
12.17 0.434
12.91 0.296
-0.046 13.75
14.21
11.70
12.10
12.30
12.20
11.59
12.0
12.25
12.21
11.64
12.05
12.28
12.21
12.56 12.53 12.60 0.31
13.14 13.08 13.20 0.55
0.51
13.70
13.56
13.64
2 10.7
3 11.2
4 11.7
5 12.1
6 12.3
7 12.2
8 12.6
9 13.2
10 13.7
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
t xt
St1 St2
at bt yt et
按照表3-4中的误差值et容易看到,对所有的rt都有rk具有随机性(利用检验法也一样),预测有效。
2
1.96
0.69。说明误差数列利用Matlab程序funesm2.m计算结果如下:
>> x=[10.1,10.7,11.2,11.7,12.1,12.3,12.2,12.6,13.2,13.7]; >> funesm2(x,0.90,0.05,0.90,4,0.05)
T=1
ans = s1: ans =
Columns 1 through 7
10.1000 10.6400 11.1440 11.6444 12.0544 12.2754 12.2075 Columns 8 through 10 12.5608 13.1361 13.6436 ans = s2: ans =
Columns 1 through 7
10.1000 10.5860 11.0882 11.5888 12.0079 12.2487 12.2117 Columns 8 through 10 12.5258 13.0751 13.5868 ans = a2: ans =
Columns 1 through 7
0 10.6940 11.1998 11.7000 12.1010 12.3022 12.2034 Columns 8 through 10 12.5957 13.1971 13.7005 ans = b2: ans =
Columns 1 through 7
0 0.4860 0.5022 0.5006 0.4191 0.2408 -0.0370 Columns 8 through 10 0.3142 0.5492 0.5117 ans = y: ans =
Columns 1 through 7
0 0 11.1800 11.7020 12.2006 12.5201 12.5430 Columns 8 through 11
12.1664 12.9099 13.7463 14.2122 ans = e2: ans =
Columns 1 through 7
0 0 0.0200 -0.0020 -0.1006 -0.2201 -0.3430 Columns 8 through 10 0.4336 0.2901 -0.0463 These data are stochastic a = 0.9000
ans = 14.2122
3.3.3 布朗(Brown)非线性指数平滑法
当时间序列呈现非线性趋势时,可以采用布朗(Brown)非线性指数平滑法进行预测。其基本原理是在一次指数平滑、二次指数平滑数列的基础山,再进行一次指数平滑,即三次指数平滑,然后,以此对非线性模型(这里就是二次曲线模型)的参数进行估计,从而达到对非线性时间序列进行预测的目的。设时间序列为x1,x2,,xN,那么,第t时点的三次指数平滑数列St3(t1,2,N)按如下递推公式计算:
St1axt(1a)St11212
StaSt(1a)St1 (3-12)
S3aS2(1a)S2
tt1t
123
一般来说,取初始值S0S0S0x1,且取同一个平滑系数a。根据选定的较合理
的平滑系数a值的计算结果St1,St2和St3,按下式(3-13)计算系数。
123a3S3SStttta
(65a)St1(108a)St2(43a)St3 (3-13) bt2
2(1a)a2
(St12St2St3)ct2
2(1a)
合理的平滑系数a选定与前面的方法相同。即按均方差MSE或平均绝对误差MAD最小所对应的平滑系数为原则。用T表示自t时点起向前预测的时点期数,则此时的(布朗)
二次抛物线型预测模型为:
ytTatbtTctT2 (3-14)
预测模型的有效性检验方法同前。即先计算预测偏差et数列及其自相关系数rk,再按
2检验法来检验。或通过,rk1.N是否成立来检验。编写Matlab程序funesm3.m
如下:
funesm3.m
%布朗非线性指数平滑
%输入时间序列x,平滑系数初值L0,步长L1,终值L2
%输入判断误差是否为随机误差时必须计算的自相关系数个数m,显著性水平alpha function ESM3=funesm3(x,L0,L1,L2,m,alpha) T=input('T=')
s1=zeros(round((L2-L0)/L1),length(x)); s2=zeros(round((L2-L0)/L1),length(x)); s3=zeros(round((L2-L0)/L1),length(x)); a3=zeros(round((L2-L0)/L1),length(x));
b3=zeros(round((L2-L0)/L1),length(x)); c3=zeros(round((L2-L0)/L1),length(x)); y=zeros(round((L2-L0)/L1),length(x)+T); e3=zeros(round((L2-L0)/L1),length(x)); MAD3=zeros(1,round((L2-L0)/L1));k=0; for a=L0:L1:L2 k=k+1;
s1(k,1)=x(1); s2(k,1)=x(1); s3(k,1)=x(1); for i=2:length(x)
s1(k,i)=a*x(i)+(1-a)*s1(k,i-1); s2(k,i)=a*s1(k,i)+(1-a)*s2(k,i-1); s3(k,i)=a*s2(k,i)+(1-a)*s3(k,i-1); a3(k,i)=3*s1(k,i)-3*s2(k,i)+s3(k,i);
b3(k,i)=(a/(2*(1-a)^2))*((6-5*a)*s1(k,i)-(10-8*a)*s2(k,i)+(4-3*a)*s3(k,i)); c3(k,i)=(a^2/(2*(1-a)^2))*(s1(k,i)-2*s2(k,i)+s3(k,i)); y(k,i+T)=a3(k,i)+b3(k,i)*T+c3(k,i)*T^2; e3(k,i+T)=x(i+T)-y(k,i+T); end end
's1:',s1(k,:),'s2:',s2(k,:),'s3:',s3(k,:)
'a3:',a3(k,:),'b3:',b3(k,:),'c3:',c3(k,:),'y:',y(k,:),'e3:',e3(k,:), funcoef(e3(k,:),m,alpha); MAD3(k)=mean(abs(e3(k,:))); end
MAD3;[MAD3,k]=min(MAD3);
a=L0+L1*(k-1),y(k,length(x)+T) end
程序中的s1,s2,s3分别表示一次、二次以及三次指数平滑值,a3,b3,c3表示相对于每个平滑系数的at,bt,ct(公式(3-13)),y表示布朗非线性趋势的预测值。平滑系数从L0到L2(步长L1),然后,从众多平滑系数计算的结果中挑选最小的MAD所对应的平滑系数并最后得到预测值,自t时点起先前的时期数T由预测者确定并输入。
【例3-10】某国企1988-2000年的利润(万元)历史数据如表3-5所示。今取平滑系数0.5并利用三次指数平滑数列,用二次抛物线型预测模型(3-14)预测2001、2002年的利润。
表3-5 历史数据、三次指数平滑值0.50
年份 1988 1989 1990 1991 1992 1993 1994 t 0 1 2 3 4 5 6 7
xt
10.6
10.6 10.6
15.1 12.85
17.6 15.23
21.6 18.42
24.8 21.61
19.5 25.56
30.4 27.89
St1
St2 St31
t
10.6 10.6 8 33.0 30.49 27.78 24.99
10.6 10.6 9 34.5 32.49 30.14 27.56
11.73 11.17 10 52.4 42.45 36.29 31.93
13.48 12.32 11 67.9 55.17 45.73 38.98
15.95 14.14 12 79.3 67.24 56.49 47.66
18.78 16.46 13 89.8 78.52 67.50 57.58
22.17 19.31
25.07 22.19
xt
St1
St2 St31
利用手工计算如下。根据公式(3-13),有(系数的检验略):
a2000378.52367.557.5890.64
b2000c2000
0.5
(62.5)78.52(104)67.5(41.5)57.5813.772
20.250.5(78.52267.557.58)0.55 20.25
这样,从2000年起,二次抛物线预测模型为:
y2000T90.6413.77T0.55T2
即2001年塑料产量预测值为(T1):y200190.6413.770.55104.89 2002年塑料产量预测值为(T2):y200190.6413.7720.554120.38 实际预测时,应多选择平滑系数,并从中挑选合理的平滑系数对应的预测值。实际上,,对于这个例题,利用Matlab程序funesm3.m计算表明,在预测2001年的利润时,合理的平滑系数a0.6,而预测2002年的利润时,合理的平滑系数a0.6。参见计算结果(具体程序结果参见Matlab程序中运行结果)。
>> x=[10.6,15.1,17.6,21.6,24.8,19.5,30.4,33.0,34.5,52.4,67.9,79.3,89.8]; >> funesm3(x,0.10,0.05,0.90,7,0.05)
T=1 a=0.6000 ans=102.7434
>> funesm3(x,0.10,0.05,0.90,7,0.05) T=2 a=0.3500 ans=120.8212
由计算结果可知,2001年利润的预测值是102.74(万元),而2002年利润预测值是120.82(万元)。
3.2 时间序列的指数平滑预测法
指数平滑法(Expinential smoothing method)的思想也是对时间序列进行修匀以消除不规则和随机的扰动。该方法是建立在如下基础上的加权平均法:即认为时间序列中的近期数据对未来值的影响比早期数据对未来值得影响更大。于是通过对时间序列的数据进行加权处理,越是近期的数据,其权数越大;反之,权数就越小。这样就将数据修匀了,并反映出时间序列中对预测时点值的影响程度。根据修匀的要求,可以有一次、二次甚至三次指数平滑。
3.3.1 一次指数平滑法
1.一次指数平滑法的计算公式及平滑系数a的讨论
设时间序列为x1,x2x3,,xN,一次指数平滑数列的递推公式为:
11Staxt(1a)St1,1S0x1,
0a1,1tN
(3-6)
式中,St1表示第t时点的一次指数平滑值,a称为平滑系数。递推公式(3-6)中,初
1
始值S0常用时间序列的首项x1(适用于历史数据个数较多,如50个历史数据及以上),如
果历史数据个数较少,如在15或20个数据及以下时,可以选用最初几期历-史数据的平均
1值作为初始值S0,这些选择都有一定的经验性和主观性。
下面讨论平滑系数a。将递推公式(3-6)展开可得:
St1axt(1a)St11axt(1a)axt1(1a)St12
axta(1a)xt1(1a)2St12
1
axta(1a)xt1a(1a)2xt2a(1a)t1x1(1a)tS0
容易看出,由于0a1,xi的系数a(1a)随着i的增加而递减。注意到这些系数之和为1,即:
i
a(1a)
i1
t
i1
1(1a)t
(1a)a(1a)t1
1(1a)
t
于是,递推公式(3-6)中的St1就是样本值x1,x2,,xt的一个加权平均。当用递推公式(3-6)来进行预测时,我们将用St1作为第t1时点的预测值。从上面的讨论可以看到,离预测时点t1最近的时点t的值xt,其权为a,最大,其次为xt1的权a(1a),,x1的权最小。
可见,公式(3-6)是在认为新近数据对未来影响大,远期数据对未来影响小的情况下对原时间序列的加权平均(修匀)。
11
若平滑系数a0,此时有St1St11S0(x1),这表明确定S0(x1)后,每个
1
时点的平滑值皆等于S0,此时,每个时点i的观察值xi不起任何作用。
若平滑系数a1,此时有St1xt。说明平滑后数列St1就是原时间序列,即没有对原时间序列进行任何处理和修匀。
对于平滑系数a来说,除了上述两种极端情况外,不可能出现xi与xj(ij)系数相等的情况。
综上所述,不难看到,0a1且比较接近1时,计算得到的一次指数平滑值对原历史数据的修匀程度将较小,平滑后的数列值St1能够比较快地反映出原时间序列的实际变化,因此,对于变化较大或趋势性较强的时间序列适合选择比较靠近1的数据作为平滑系数,比如取a0.95,0.90等。
若0a1且取值比较靠近0时,计算得到的一次指数平滑值对原历史数据的修匀程度将较大,平滑后的数列对原时间序列的变化反应较迟钝。因此,对于变化较小的、或接近平稳的时间序列,应选择比较靠近0的平滑系数使得平滑过程中的各数据的权数比较接近。
【例3-7】某电器销售企业1992-2003年某种电器销售额(万元)及其一次指数平滑数列计算列表如表3-2所示。
表3-2 一次指数平滑值计算表 (单位:万元)
年份 销售额
1992 50 51 51 51
52 50.8 50.50 50.20
47 51.04 51.25 51.64
51 50.23 49.13 47.93
49 50.38 50.07 50.39
48 50.10 49.54 49.28
51 49.68 48.77 48.26
40 49.94 49.89 50.45
48 47.95 44.95 42.09
51 47.96 46.48 46.82
51 48.57 48.74 50.16
59 49.06 49.87 50.83
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
0.2
0.5
0.8
1
本例中选择第一、第二期的历史数据平均值作为一次指数平滑的初始值S0。从表3-2
中可以看到,原时间序列变动较大,稍具周期性。若想指数平滑后数列敏感地反映出最新的一些观察值的变动,则应取较大的平滑系数a,如a0.8。若想消除其周期性变动,施加较大程度的修匀,以反映其长期趋势,则可取较小的a,如a0.2。
2.检验及预测
一次指数平滑法适用于对变化比较平稳的时间序列作短期的预测。第t1时点的预测值yt1等于按递推公式(3-6)计算的第t时点的一次指数平滑值St1,即:
yt1St1 (3-7) 类似于移动平均法,同一个问题随着平滑系数的不同可以有若干个一次指数平滑预测值。因此,应该在一个合适的评价标准基础上选择一个合理的平滑系数a。其方法是:首先计算与在原则上合理的多个平滑系数a相对应的平滑数列,然后分别计算其均方差MSE(见公式(3-7))或计算其平均绝对误差MAD(见公式(3-8)),以MSE或MAD最小者对应的平滑系数及其预测值为最合理的。
1N
MADet,而etxtSt11xtyt (3-8)
Nt1
公式(3-8)中的etxtSt11xtyt反映了各个时点的平滑值St11与实际值xt之间的误差。根据公式(3-2),计算误差数列et的自相关系数rt,若统计量
2
Q(N1)rt2(m1)
t1m
则说明误差数列具有随机性,可以认为此时的预测是有效的(可以利用前面的Matlab
程序funcoef.m判断)。一次指数平滑法的计算、检验及预测的Matlab程序如下(文件名funesm1.m)。
【例3-8】某商品2004年各月的销售量xt如表3-3所示,试预测2005年1月的销售量。
1利用funesm1.m计算的部分结果如下,其他部分数据如表3-3所示。程序中令S0x1423,
取L00.10,L10.05,L30.95,m7,0.05,对应的误差值et,也列于表3-3中。
表3-3 历史数据、一次指数平滑值及误差值
月 序号
0 0
1 1 423
2 2 358
3 3 434
4 4 445
5 5 527
6 6 429
7 7 426
8 8 502
9 9 480
10 10 385
11 11 427
12 12 446
t
xt
St1(0.1)
423
423 0
416.5 -65
418.8 17.5
420.9 26.7
431.5 106.1
431.3 -2.5
430.8 -5.3
437.9 71.2
442.1 42.1
436.4 -57.1
435.5 -9.4
436.5 10.5
et(0.1)
St1(0.2)
423
423 0
410 414.8
420.8 30.2
442 439.4
436.7 -13.4
449.8 65.3
455.8 30.2
441.6 -70.8
438.7 -14.6 436.9 -14.1
420.2 7.3
et(0.2)
-65 24 106.2
-13
St1(0.3)
423
423 0
403.5 -65
412.7 30.5
422.4 32.3
453.8 104.6
446.4 -24.8
440.3 -20.4
458.8 61.7
465.2 21.2
441.1 -80.2
439.8 9.1
et(0.3)
这就是说,取平滑系数0.10是合适的,此时预测值为2005年1月的销售量将是y = 436.4976。
funesm1.m
%一次指数平滑预测,文件名为funesm1.m
%输入时间序列想,平滑系数初值L0,步长L1,终值L2
%输入判断误差是否为随机误差时需要计算的自相关系数个数m,显著性水平alpha function ESM1=funesm1(x,L0,L1,L2,m,alpha)
s=zeros(round((L2-L0)/L1),length(x)); e=zeros(1,length(x));
MAD=zeros(1,round((L2-L0)/L1)); k=0;
for a=L0:L1:L2 k=k+1
s(k,1)=x(1); for i=2:length(x)
s(k,i)=a*x(i)+(1-a)*s(k,i-1); e(i)=x(i)-s(k,i-1); end
s(k,:),e
funcoef(e,m,alpha); MAD(k)=mean(abs(e)); end
disp ('The smallest MAD and the corresponding a and forecast') MAD,[MAD,k]=min(MAD);
a=L0+L1*(k-1),y=s(k,length(x)) end
利用上述程序计算得到的一次指数平滑预测值是在多个平滑系数(可自己确定)中计算、比较得到的。
3.3.2 二次指数平滑法
对于呈现出线性趋势的时间序列,在一次指数平滑数列的基础上用同一个平滑系数a再进行一次指数平滑,就是二次指数平滑。构成二次指数平滑数列St2的递推公式如下:
21令初始值S0,则有: S0x1(也可取其他值作为初始值)
11
Staxt(1a)St1
2 (3-9) 12
StaSt(1a)St1
二次指数平滑的目的是对原时间序列进行两次修匀,使得其不规则变动或周期变动尽量
消除掉,让时间序列的长期趋势性更能显示出来。对于平滑系数,同样有一个合理选区的问题。其方法与一次指数平滑法一样,先选取原则上较合理的多个a值分别计算,得到不同的数列St1和St2,再根据均方差MSE或MAD最小原则确定较为合理的a值,并得到相应的二次指数平滑值。
由于二次指数平滑的目的,二次指数平滑法较适用于具有线性趋势的时间序列,线性趋势预测模型为:
ytTatbtT(3-10)
上式中,T表示自T时点起向前预测的时点数;at,bt满足:
atSt1(St1St2)2St1St2,bt
a
(St1St2)(3-11) 1a
预测模型的有效性检验方法与一次指数平滑法一样。即通过自相关系数或2检验方法进行检验。利用二次指数平滑法预测具有线性趋势性的时间序列的基本步骤为:
1. 根据历史数据(时间序列)x按照公式(3-6)计算一次指数平滑值; 2.根据公式(3-9)计算二次指数平滑值;
3. 由公式(3-11)计算at,bt,并由(3-10)计算自t时点起先前T时期的预测值ytT。据此,编写的Matlab程序如下(文件名funesm2.m)。其中,程序中的s1,s2分别表示一次、二次指数平滑值,a2,b2表示相对于每个平滑系数的at,bt(公式(3-11)),y表示线性趋势的预测值。平滑系数从L2(步长L1),然后,从众多平滑系数计算的结果中挑选最小的MAD 所对应的平滑系数并最后得到预测值,自t时点起先前的时期数T由预测值确定并输入。
funesm2.m
%二次指数平滑,线性趋势预测模型
%输入时间序列x,平滑系数初值L0,步长L1,终值L2
%输入判断误差是否为随机误差时必须计算的自相关系数个数m,显著性水平alpha function ESM2=funesm2(x,L0,L1,L2,m,alpha) T=input('T=')
s1=zeros(round((L2-L0)/L1),length(x)); s2=zeros(round((L2-L0)/L1),length(x)); a2=zeros(round((L2-L0)/L1),length(x)); b2=zeros(round((L2-L0)/L1),length(x)); y=zeros(round((L2-L0)/L1),length(x)+T); e2=zeros(round((L2-L0)/L1),length(x)); MAD2=zeros(1,round((L2-L0)/L1));k=0; for a=L0:L1:L2 k=k+1;
s1(k,1)=x(1); s2(k,1)=x(1); for i=2:length(x)
s1(k,i)=a*x(i)+(1-a)*s1(k,i-1); s2(k,i)=a*s1(k,i)+(1-a)*s2(k,i-1); a2(k,i)=2*s1(k,i)-s2(k,i);
b2(k,i)=(s1(k,i)-s2(k,i))*(a/(1-a)); y(k,i+T)=a2(k,i)+b2(k,i)*T; if i+T
e2(k,i+T)=x(i+T)-y(k,i+T); end end
's1:',s1(k,:),'s2:',s2(k,:),
'a2:',a2(k,:),'b2:',b2(k,:),'y:',y(k,:),'e2:',e2(k,:), funcoef(e2(k,:),m,alpha);
MAD2(k)=mean(abs(e2(k,:))); end
MAD2;[MAD2,k]=min(MAD2);
a=L0+L1*(k-1),y(k,length(x)+T) end
【例3-9】某市1990-1999年工业产值为xt(t1,2,,12),如表(3-4)所示,由于原时间序列大体上呈现增长的趋势,取平滑系数a0.9利用二次指数平滑法计算的结果如表(包括at,bt值以及误差值et)3-4所示。
表3-4中的at,bt值按公式(3-11)计算。于是2000年的预测值为:
y13.700.51114.21。
而2005年的预测值为:y13.700.51616.76
表3-4 历史数据、二次指数平滑值及误差值
年份
0 10.1 10.1
1990 1 10.1 10.1 10.1
10.64 10.59 10.69 0.49
11.14 11.09 11.20 0.50 11.18 0.02
-0.002
-0.101
-0.220
-0/343
0.50 11.70
0.42 12.20
0.24 12.52
-0.040 12.54
12.17 0.434
12.91 0.296
-0.046 13.75
14.21
11.70
12.10
12.30
12.20
11.59
12.0
12.25
12.21
11.64
12.05
12.28
12.21
12.56 12.53 12.60 0.31
13.14 13.08 13.20 0.55
0.51
13.70
13.56
13.64
2 10.7
3 11.2
4 11.7
5 12.1
6 12.3
7 12.2
8 12.6
9 13.2
10 13.7
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
t xt
St1 St2
at bt yt et
按照表3-4中的误差值et容易看到,对所有的rt都有rk具有随机性(利用检验法也一样),预测有效。
2
1.96
0.69。说明误差数列利用Matlab程序funesm2.m计算结果如下:
>> x=[10.1,10.7,11.2,11.7,12.1,12.3,12.2,12.6,13.2,13.7]; >> funesm2(x,0.90,0.05,0.90,4,0.05)
T=1
ans = s1: ans =
Columns 1 through 7
10.1000 10.6400 11.1440 11.6444 12.0544 12.2754 12.2075 Columns 8 through 10 12.5608 13.1361 13.6436 ans = s2: ans =
Columns 1 through 7
10.1000 10.5860 11.0882 11.5888 12.0079 12.2487 12.2117 Columns 8 through 10 12.5258 13.0751 13.5868 ans = a2: ans =
Columns 1 through 7
0 10.6940 11.1998 11.7000 12.1010 12.3022 12.2034 Columns 8 through 10 12.5957 13.1971 13.7005 ans = b2: ans =
Columns 1 through 7
0 0.4860 0.5022 0.5006 0.4191 0.2408 -0.0370 Columns 8 through 10 0.3142 0.5492 0.5117 ans = y: ans =
Columns 1 through 7
0 0 11.1800 11.7020 12.2006 12.5201 12.5430 Columns 8 through 11
12.1664 12.9099 13.7463 14.2122 ans = e2: ans =
Columns 1 through 7
0 0 0.0200 -0.0020 -0.1006 -0.2201 -0.3430 Columns 8 through 10 0.4336 0.2901 -0.0463 These data are stochastic a = 0.9000
ans = 14.2122
3.3.3 布朗(Brown)非线性指数平滑法
当时间序列呈现非线性趋势时,可以采用布朗(Brown)非线性指数平滑法进行预测。其基本原理是在一次指数平滑、二次指数平滑数列的基础山,再进行一次指数平滑,即三次指数平滑,然后,以此对非线性模型(这里就是二次曲线模型)的参数进行估计,从而达到对非线性时间序列进行预测的目的。设时间序列为x1,x2,,xN,那么,第t时点的三次指数平滑数列St3(t1,2,N)按如下递推公式计算:
St1axt(1a)St11212
StaSt(1a)St1 (3-12)
S3aS2(1a)S2
tt1t
123
一般来说,取初始值S0S0S0x1,且取同一个平滑系数a。根据选定的较合理
的平滑系数a值的计算结果St1,St2和St3,按下式(3-13)计算系数。
123a3S3SStttta
(65a)St1(108a)St2(43a)St3 (3-13) bt2
2(1a)a2
(St12St2St3)ct2
2(1a)
合理的平滑系数a选定与前面的方法相同。即按均方差MSE或平均绝对误差MAD最小所对应的平滑系数为原则。用T表示自t时点起向前预测的时点期数,则此时的(布朗)
二次抛物线型预测模型为:
ytTatbtTctT2 (3-14)
预测模型的有效性检验方法同前。即先计算预测偏差et数列及其自相关系数rk,再按
2检验法来检验。或通过,rk1.N是否成立来检验。编写Matlab程序funesm3.m
如下:
funesm3.m
%布朗非线性指数平滑
%输入时间序列x,平滑系数初值L0,步长L1,终值L2
%输入判断误差是否为随机误差时必须计算的自相关系数个数m,显著性水平alpha function ESM3=funesm3(x,L0,L1,L2,m,alpha) T=input('T=')
s1=zeros(round((L2-L0)/L1),length(x)); s2=zeros(round((L2-L0)/L1),length(x)); s3=zeros(round((L2-L0)/L1),length(x)); a3=zeros(round((L2-L0)/L1),length(x));
b3=zeros(round((L2-L0)/L1),length(x)); c3=zeros(round((L2-L0)/L1),length(x)); y=zeros(round((L2-L0)/L1),length(x)+T); e3=zeros(round((L2-L0)/L1),length(x)); MAD3=zeros(1,round((L2-L0)/L1));k=0; for a=L0:L1:L2 k=k+1;
s1(k,1)=x(1); s2(k,1)=x(1); s3(k,1)=x(1); for i=2:length(x)
s1(k,i)=a*x(i)+(1-a)*s1(k,i-1); s2(k,i)=a*s1(k,i)+(1-a)*s2(k,i-1); s3(k,i)=a*s2(k,i)+(1-a)*s3(k,i-1); a3(k,i)=3*s1(k,i)-3*s2(k,i)+s3(k,i);
b3(k,i)=(a/(2*(1-a)^2))*((6-5*a)*s1(k,i)-(10-8*a)*s2(k,i)+(4-3*a)*s3(k,i)); c3(k,i)=(a^2/(2*(1-a)^2))*(s1(k,i)-2*s2(k,i)+s3(k,i)); y(k,i+T)=a3(k,i)+b3(k,i)*T+c3(k,i)*T^2; e3(k,i+T)=x(i+T)-y(k,i+T); end end
's1:',s1(k,:),'s2:',s2(k,:),'s3:',s3(k,:)
'a3:',a3(k,:),'b3:',b3(k,:),'c3:',c3(k,:),'y:',y(k,:),'e3:',e3(k,:), funcoef(e3(k,:),m,alpha); MAD3(k)=mean(abs(e3(k,:))); end
MAD3;[MAD3,k]=min(MAD3);
a=L0+L1*(k-1),y(k,length(x)+T) end
程序中的s1,s2,s3分别表示一次、二次以及三次指数平滑值,a3,b3,c3表示相对于每个平滑系数的at,bt,ct(公式(3-13)),y表示布朗非线性趋势的预测值。平滑系数从L0到L2(步长L1),然后,从众多平滑系数计算的结果中挑选最小的MAD所对应的平滑系数并最后得到预测值,自t时点起先前的时期数T由预测者确定并输入。
【例3-10】某国企1988-2000年的利润(万元)历史数据如表3-5所示。今取平滑系数0.5并利用三次指数平滑数列,用二次抛物线型预测模型(3-14)预测2001、2002年的利润。
表3-5 历史数据、三次指数平滑值0.50
年份 1988 1989 1990 1991 1992 1993 1994 t 0 1 2 3 4 5 6 7
xt
10.6
10.6 10.6
15.1 12.85
17.6 15.23
21.6 18.42
24.8 21.61
19.5 25.56
30.4 27.89
St1
St2 St31
t
10.6 10.6 8 33.0 30.49 27.78 24.99
10.6 10.6 9 34.5 32.49 30.14 27.56
11.73 11.17 10 52.4 42.45 36.29 31.93
13.48 12.32 11 67.9 55.17 45.73 38.98
15.95 14.14 12 79.3 67.24 56.49 47.66
18.78 16.46 13 89.8 78.52 67.50 57.58
22.17 19.31
25.07 22.19
xt
St1
St2 St31
利用手工计算如下。根据公式(3-13),有(系数的检验略):
a2000378.52367.557.5890.64
b2000c2000
0.5
(62.5)78.52(104)67.5(41.5)57.5813.772
20.250.5(78.52267.557.58)0.55 20.25
这样,从2000年起,二次抛物线预测模型为:
y2000T90.6413.77T0.55T2
即2001年塑料产量预测值为(T1):y200190.6413.770.55104.89 2002年塑料产量预测值为(T2):y200190.6413.7720.554120.38 实际预测时,应多选择平滑系数,并从中挑选合理的平滑系数对应的预测值。实际上,,对于这个例题,利用Matlab程序funesm3.m计算表明,在预测2001年的利润时,合理的平滑系数a0.6,而预测2002年的利润时,合理的平滑系数a0.6。参见计算结果(具体程序结果参见Matlab程序中运行结果)。
>> x=[10.6,15.1,17.6,21.6,24.8,19.5,30.4,33.0,34.5,52.4,67.9,79.3,89.8]; >> funesm3(x,0.10,0.05,0.90,7,0.05)
T=1 a=0.6000 ans=102.7434
>> funesm3(x,0.10,0.05,0.90,7,0.05) T=2 a=0.3500 ans=120.8212
由计算结果可知,2001年利润的预测值是102.74(万元),而2002年利润预测值是120.82(万元)。