时间序列分析

确定性时间序列分析方法

时间序列预测：将预测对象按时间顺序排列，构成一个“时间序列”，再从该时间序列过去的变化规律，推断今后变化的可能性和变化趋势、变化规律。但一般只适用于短期预测，也不能反映事物的内在联系，不能分析两个因素的相互关系。

一个时间序列往往使一下几种变化类型的叠加或耦合：
（1）长期趋势变动：时间序列朝一定方向持续上升或下降或停留在某一水          平，反映了客观事物的主要变化趋势。
（2）季节变动
（3）循环变动
（4）不规则变动
通常用 T t T_{t} Tt​表示长期趋势变动， s t s_{t} st​表示季节变动趋势， C t C_{t} Ct​表示循环变动趋势， R t R_{t} Rt​表示随机干扰。常见确定性时间序列模型有：
（1）加法模型
y t = T t + S t + C t + R t y_{t}=T_{t}+S_{t}+C_{t}+R_{t} yt​=Tt​+St​+Ct​+Rt​
（2）乘法模型
y t = T t ⋅ S t ⋅ C t ⋅ R t y_{t}=T_{t}\cdot S_{t}\cdot C_{t}\cdot R_{t} yt​=Tt​⋅St​⋅Ct​⋅Rt​
（3）混合模型
y t = T t ⋅ S t + R t y_{t}=T_{t}\cdot S_{t}+R_{t} yt​=Tt​⋅St​+Rt​ y t = S t + T t ⋅ C t ⋅ R t y_{t}=S_{t}+T_{t}\cdot C_{t}\cdot R_{t} yt​=St​+Tt​⋅Ct​⋅Rt​

一、移动平均法

即以 N N N项的平均值作为预测值：
设观测序列为： y 1 , ⋅ ⋅ ⋅ , y T y_{1},\cdot \cdot \cdot ,y_{T} y1​,⋅⋅⋅,yT​,取移动平均的项数 N < T N<T N<T.
一次移动平均：
M t ( 1 ) = 1 N ( y t + y t − 1 + ⋅ ⋅ ⋅ + y t − N + 1 ) M_{t}^{(1)}=\frac{1}{N}(y_{t}+y_{t-1}+\cdot \cdot \cdot +y_{t-N+1}) Mt(1)​=N1​(yt​+yt−1​+⋅⋅⋅+yt−N+1​)
= 1 N ( y t − 1 + ⋅ ⋅ ⋅ + y t − N ) + 1 N ( y t − y t − N ) =\frac{1}{N}(y_{t-1}+\cdot \cdot \cdot +y_{t-N})+\frac{1}{N}(y_{t}-y_{t-N}) =N1​(yt−1​+⋅⋅⋅+yt−N​)+N1​(yt​−yt−N​)
= M t − 1 ( 1 ) + 1 N ( y t − y t − N ) =M_{t-1}^{(1)}+\frac{1}{N}(y_{t}-y_{t-N}) =Mt−1(1)​+N1​(yt​−yt−N​)
二次移动平均：
M t ( 2 ) = 1 N ( M t ( 1 ) + ⋅ ⋅ ⋅ + M t − N + 1 ( 1 ) ) M_{t}^{(2)}=\frac{1}{N}(M_{t}^{(1)}+\cdot \cdot \cdot +M_{t-N+1}^{(1)}) Mt(2)​=N1​(Mt(1)​+⋅⋅⋅+Mt−N+1(1)​)
= M t − 1 ( 2 ) + 1 N ( M t ( 1 ) − M t − N ( 1 ) ) =M_{t-1}^{(2)}+\frac{1}{N}(M_{t}^{(1)}-M_{t-N}^{(1)}) =Mt−1(2)​+N1​(Mt(1)​−Mt−N(1)​)
当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均法建立预测模型：
y ^ t + 1 = M t ( 1 ) , t = N , N + 1 , ⋅ ⋅ ⋅ T \hat{y}_{t+1}=M_{t}^{(1)},t=N,N+1,\cdot \cdot \cdot T y^​t+1​=Mt(1)​,t=N,N+1,⋅⋅⋅T
预测标误差：
s = ∑ t = N + 1 T ( y ^ t − y i ) 2 T − N s=\sqrt{\frac{\sum_{t=N+1}^{T}(\hat{y}_{t}-y_{i})^2}{T-N}} s=T−N∑t=N+1T​(y^​t​−yi​)2​ ​
一般 N N N的取值范围是 [ 5 , 20 ] [5,20] [5,20],历史序列的基本变化趋势不大，随机变动成分较多时， N N N取值大一些。

当预测目标的基本趋势与某一线性模型相吻合时，常采用二次平均法。

例：预测12月销售收入。

-	-	-	-	-	-	-
月份	1	2	3	4	5	6
销售收入	533.8	574.6	606.9	649.8	705.1	772.0
月份	7	8	9	10	11	12
销售收入	816.4	892.7	963.9	1015.1	1102.7

分别取 N = 4 、 5 N=4、5 N=4、5带入上诉公式可以得到预测值，但计算标准误差 s s s后，发现 N = 4 N=4 N=4时的误差更小，所以取 N = 4 N=4 N=4。

clc,clear
y=[533.8  574.6  606.9  649.8   705.1  772.0  816.4  892.7  963.9  1015.1  1102.7];
x=1:11;
plot(x,y,'-');%画如直观观察变化趋势
m=length(y);   
n=[4,5];             %n为移动平均的项数
for i=1:length(n)     %由于n的取值不同，下面使用了细胞数组
    for j=1:m-n(i)+1
        yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i); 
    end
    y12(i)=yhat{i}(end);       %提出第12月份的预测值
    s(i)=sqrt(mean((y(n(i)+1:end)-yhat{i}(1:end-1)).^2)); %求预测的标准误差
end
y12, s   %分别显示两种方法的预测值和预测的标准误差

细胞数组：参考1、参考2

做出收入变化趋势图：

N = 4 、 5 N=4、5 N=4、5预测结果及误差：

y12 =

  993.6000  958.1600


s =

  150.5121  182.3851

二、指数平滑法

一次移动平均认为最近 N N N期的数据对未来值影响相同，但实际上，历史数据对未来的影响一般都是随时间的增长而递减的，所以要对各期观测值按时间顺序进行加权平均作为预测值。

1. 一次指数平滑法

1. 预测模型
设α为加权系数 0 < α < 1 0<α<1 0<α<1,一次指数平滑公式为：
S t ( 1 ) = α y t + ( 1 − α ) S t − 1 ( 1 ) = S t − 1 ( 1 ) + α ( y t − S t − 1 ( 1 ) ) S_{t}^{(1)}=αy_{t}+(1-α)S_{t-1}^{(1)}=S_{t-1}^{(1)}+α(y_{t}-S_{t-1}^{(1)}) St(1)​=αyt​+(1−α)St−1(1)​=St−1(1)​+α(yt​−St−1(1)​)
展开：
S t ( 1 ) = α y t + ( 1 − α ) [ α y t − 1 + ( 1 − α ) S t − 2 ( 1 ) ] = α ∑ j = 0 ∞ ( 1 − α ) j y t − i S_{t}^{(1)}=αy_{t}+(1-α)[αy_{t-1}+(1-α)S_{t-2}^{(1)}]=α\sum_{j=0}^{∞}(1-α)^jy_{t-i} St(1)​=αyt​+(1−α)[αyt−1​+(1−α)St−2(1)​]=αj=0∑∞​(1−α)jyt−i​
这说明 S t ( 1 ) S_{t}^{(1)} St(1)​是所有历史数据的加权平均，加权系数为 α , α ( 1 − α ) , α ( 1 − α ) 2 , ⋅ ⋅ ⋅ , α,α(1-α),α(1-α)^2,\cdot \cdot \cdot, α,α(1−α),α(1−α)2,⋅⋅⋅,它们的和为1。加权系数符合指数规律，故称之为指数平滑。预测模型为：
y ^ t + 1 = S t ( 1 ) = α y t + ( 1 − α ) y ^ t \hat{y}_{t+1}=S_{t}^{(1)}=αy_{t}+(1-α)\hat{y}_{t} y^​t+1​=St(1)​=αyt​+(1−α)y^​t​
即第 t t t期的数据平滑值作为 t + 1 t+1 t+1期的预测值。

2. 加权系数的选择

①时间序列波动小，比较平稳，α取小一点；
②时间序列具有迅速且明显的变动倾向，α取大一点。
实际运用中，多取结果α值进行比较，哪个误差小就用哪个。

3. 初始值的确定
初始值 S 0 ( 1 ) S_{0}^{(1)} S0(1)​由预测者估计或指定。时间序列比较多（大于20）时，一般取第一期作数据为初始值，较少时，一般取最初几期实际值的平均值作为初始值。
例：预测1998年销售额

年份	t	实际销售额 y i y_{i} yi​	预测值 y ^ i \hat{y}_{i} y^​i​ α=0.2	预测值 y ^ i \hat{y}_{i} y^​i​ α=0.5	预测值 y ^ i \hat{y}_{i} y^​i​ α=0.8
1976	1	50	51	51	51
1977	2	52	50.8	50.5	50.2
1978	3	47	51.04	51.25	51.64
1979	4	51	50.23	49.13	47.93
1980	5	49	50.39	50.06	50.39
1981	6	48	50.11	49.53	49.28
1982	7	51	49.69	48.77	28.26
1983	8	40	49.95	49.88	50.45
1984	9	48	47.96	44.94	42.09
1985	10	52	47.97	46.47	46.82
1986	11	51	48.77	49.24	50.96
1987	12	59	49.22	50.12	50.99

clc,clear
yt=load('dianqi.txt');   %实际销售额数据以列向量的方式存放在纯文本文件中
n=length(yt); alpha=[0.2 0.5 0.8]; m=length(alpha);
yhat(1,[1:m])=(yt(1)+yt(2))/2;
for i=2:n
    yhat(i,:)=alpha*yt(i-1)+(1-alpha).*yhat(i-1,:);
end
yhat
err=sqrt(mean((repmat(yt,1,m)-yhat).^2)) 
xlswrite('dianqi.xls',yhat) %把预测数据写到Excel文件，准备在word表格中使用
yhat1988=alpha*yt(n)+(1-alpha).*yhat(n,:)

输出：

hat =

   51.0000   51.0000   51.0000
   50.8000   50.5000   50.2000
   51.0400   51.2500   51.6400
   50.2320   49.1250   47.9280
   50.3856   50.0625   50.3856
   50.1085   49.5313   49.2771
   49.6868   48.7656   48.2554
   49.9494   49.8828   50.4511
   47.9595   44.9414   42.0902
   47.9676   46.4707   46.8180
   48.7741   49.2354   50.9636
   49.2193   50.1177   50.9927


err =

    4.5029    4.5908    4.8426


yhat1988 =

   51.1754   54.5588   57.3985

>> 

则选取α=0.2,1998年预测值为51.1754

2. 二次指数平滑法

当时间序列的变动出现直线趋势时，用一次指数平滑法会出现明显的滞后偏差。所以采用二次指数平滑：
{ S 1 ( 1 ) = α y t + ( 1 − α ) S t − 1 ( 1 ) S t ( 2 ) = α S t ( 1 ) + ( 1 − α ) S t − 1 ( 2 ) \left\{\begin{matrix}S_{1}^{(1)}=αy_{t}+(1-α)S_{t-1}^{(1)} \\ S_{t}^{(2)}=αS_{t}^{(1)}+(1-α)S_{t-1}^{(2)} \end{matrix}\right. {S1(1)​=αyt​+(1−α)St−1(1)​St(2)​=αSt(1)​+(1−α)St−1(2)​​
当时间序列从某时刻开始具有直线趋势时，可用直线趋势模型：
y ^ t + 1 = a t + b t m , m = 1 , 2 , ⋅ ⋅ ⋅ \hat {y}_{t+1}=a_{t}+b_{t}m,m=1,2,\cdot \cdot \cdot y^​t+1​=at​+bt​m,m=1,2,⋅⋅⋅ { a t = s S t ( 1 ) − S t ( 2 ) b t = α 1 − α ( S t ( 1 ) − S t ( 2 ) ) \left\{\begin{matrix}a_{t}=sS_{t}^{(1)}-S_{t}^{(2)} \\ b_{t}=\frac{α}{1-α}(S_{t}^{(1)}-S_{t}^{(2)}) \end{matrix}\right. {at​=sSt(1)​−St(2)​bt​=1−αα​(St(1)​−St(2)​)​
进行预测。
例：预测1986和1987年的发电量

解：
取α=0.3，初始值 S 0 ( 1 ) = S 0 ( 2 ) = 676 S_{0}^{(1)}=S_{0}^{(2)}=676 S0(1)​=S0(2)​=676（取时间序列的第一项）。计算 S t ( 1 ) 、 S t ( 2 ) S_{t}^{(1)}、S_{t}^{(2)} St(1)​、St(2)​的值如上表。
将t=21时的值带入直线趋势模型，可得系数：
a 21 = 2 × 3523.1 − 3032.6 = 4013.7 a_{21}=2×3523.1-3032.6=4013.7 a21​=2×3523.1−3032.6=4013.7 b 21 = 0.3 1 − 0.3 ( 3523.1 − 3032.6 ) = 210.21 b_{21}=\frac{0.3}{1-0.3}(3523.1-3032.6)=210.21 b21​=1−0.30.3​(3523.1−3032.6)=210.21
则得到t=21时的直线趋势方程：
y ^ 21 + m = 4013.7 + 210.24 m \hat{y}_{21+m}=4013.7+210.24m y^​21+m​=4013.7+210.24m
m = 1 、 2 m=1、2 m=1、2时即对应 1986年的1987年的预测值。
程序：

clc,clear
yt=load('fadian.txt');  %原始发电总量数据以列向量的方式存放在纯文本文件中
n=length(yt), alpha=0.3; st1(1)=yt(1); st2(1)=yt(1);
for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
end
xlswrite('fadian.xls',[st1',st2']) %把数据写入表单Sheet1中的前两列
at=2*st1-st2;
bt=alpha/(1-alpha)*(st1-st2);
yhat=at+bt;  %最后的一个分量为1986年的预测值
xlswrite('fadian.xls',yhat','Sheet1','C2') %把预测值写入第3列
str=['C',int2str(n+2)]; %准备写1987年预测值位置的字符串
xlswrite('fadian.xls',at(n)+2*bt(n),'Sheet1',str)%把1987年预测值写到相应位置

其中：
fadian.txt

676
825
774
716
940
1159
1384
1524
1668
1688
1958
2031
2234
2566
2820
3006
3093
3277
3514
3770
4107

fadian.xls

3. 三次指数平滑法

当时间序列的变动表现为二次曲线趋势时，使用三次指数平滑。
{ S t ( 1 ) = α y t + ( 1 − α ) S t − 1 ( 1 ) S t ( 2 ) = α S t ( 1 ) + ( 1 − α ) S t − 1 ( 2 ) S t ( 3 ) = α S t ( 2 ) + ( 1 − α ) S t − 1 ( 3 ) \left\{\begin{matrix}S_{t}^{(1)}=αy_{t}+(1-α)S_{t-1}^{(1)} \\ S_{t}^{(2)}=αS_{t}^{(1)}+(1-α)S_{t-1}^{(2)} \\ S_{t}^{(3)}=αS_{t}^{(2)}+(1-α)S_{t-1}^{(3)} \end{matrix}\right. ⎩⎪⎨⎪⎧​St(1)​=αyt​+(1−α)St−1(1)​St(2)​=αSt(1)​+(1−α)St−1(2)​St(3)​=αSt(2)​+(1−α)St−1(3)​​
三次指数平滑的预测模型为：
y ^ t + m = a t + b t m + c t m 2 , m = 1 、 2 , ⋅ ⋅ ⋅ \hat{y}_{t+m}=a_{t}+b_{t}m+c_{t}m^2,m=1、2,\cdot \cdot \cdot y^​t+m​=at​+bt​m+ct​m2,m=1、2,⋅⋅⋅
其中：
{ a t = 3 S t ( 1 ) − 3 S t ( 2 ) + S t ( 3 ) b t = α 1 − α [ ( 6 − 5 α ) S t ( 1 ) − 2 ( 5 − 4 α ) S t ( 2 ) + ( 4 − 3 α ) S t ( 3 ) ] c t = α 2 2 ( 1 − α ) 2 [ S t ( 1 ) − 2 S t ( 2 ) + S t ( 3 ) ] \left\{\begin{matrix}a_{t}=3S_{t}^{(1)}-3S_{t}^{(2)}+S_{t}^{(3)} \\ b_{t}=\frac{α}{1-α}[(6-5α)S_{t}^{(1)}-2(5-4α)S_{t}^{(2)}+(4-3α)S_{t}^{(3)}] \\ c_{t}=\frac{α^2}{2(1-α)^2}[S_{t}^{(1)}-2S_{t}^{(2)}+S_{t}^{(3)}] \end{matrix}\right. ⎩⎪⎨⎪⎧​at​=3St(1)​−3St(2)​+St(3)​bt​=1−αα​[(6−5α)St(1)​−2(5−4α)St(2)​+(4−3α)St(3)​]ct​=2(1−α)2α2​[St(1)​−2St(2)​+St(3)​]​
例：预测1989和1990年固定资产投资总额。

解：
做投资总额趋势图：

可以看出投资总额呈二次曲线上升，故采用三次指数平滑。
取α=0.3，初始值 S 0 ( 1 ) = S 0 ( 2 ) = S 0 ( 3 ) = y 1 + y 2 + y 3 3 = 21.94 S_{0}^{(1)}=S_{0}^{(2)}=S_{0}^{(3)}=\frac{y_{1}+y_{2}+y_{3}}{3}=21.94 S0(1)​=S0(2)​=S0(3)​=3y1​+y2​+y3​​=21.94
再计算 S t ( 1 ) 、 S t ( 2 ) 、 S t ( 3 ) S_{t}^{(1)}、S_{t}^{(2)}、S_{t}^{(3)} St(1)​、St(2)​、St(3)​,如上表。
根据t=11时的值可以计算三次平滑模型的各项系数。得到预测模型为： y ^ t + m = 219.91 + 38.38 m + 1.62 m 2 \hat{y}_{t+m}=219.91+38.38m+1.62m^2 y^​t+m​=219.91+38.38m+1.62m2m取相应的值即可得到所要预测的投资总额。
程序：

clc,clear
yt=load('touzi.txt'); %原始投资总额数据以列向量的方式存放在纯文本文件中
n=length(yt); 
alpha=0.3;            %取α=0.3
st0=mean(yt(1:3));    %初值取前三期平均值

st1(1)=alpha*yt(1)+(1-alpha)*st0; %计算一二三次指数平滑的第一个平滑值
st2(1)=alpha*st1(1)+(1-alpha)*st0;
st3(1)=alpha*st2(1)+(1-alpha)*st0;

for i=2:n                        %计算所有一二三次指数平滑值
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
    st3(i)=alpha*st2(i)+(1-alpha)*st3(i-1);
end

xlswrite('touzi.xls',[st1',st2',st3'])  %把数据写在前三列

at=3*st1-3*st2+st3;                     %计算三次指数平滑预测模型的三个系数
bt=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1-2*(5-4*alpha)*st2+(4-3*alpha)*st3);
ct=0.5*alpha^2/(1-alpha)^2*(st1-2*st2+st3);
xishu=[ct(end),bt(end),at(end)]          %二次预测多项式的系数向量

yhat=at+bt+ct;   %预测模型

xlswrite('touzi.xls',yhat','Sheet1','D2')  %把数据写在第4列第2行开始的位置

plot(1:n,yt,'D',2:n,yhat(1:end-1),'*')
legend('实际值','预测值','Location','northwest')  %图注显示在左上角

yhat1989=polyval(xishu,1)  %求预测多项式m=1时的值
yhat1990=polyval(xishu,2)  %求预测多项式m=2时的值

输出:

xishu =

    1.6205   38.3849  219.9120


yhat1989 =

  259.9174


yhat1990 =

  303.1637

touzi.xls

根据历史经验：α取0.1~0.3为宜，α越大意味着采用的数据越少。一些基本准则：

1. 序列的基本趋势比较稳，预测偏差由随机因素造成，α取值小一些。此时预测模型包含更多的历史数据信息。
2. 基本趋势已发生系统的变化，α的取值应大一些。

三、差分指数平滑法

上面讲到，当时间序列的变化趋势有直线趋势时，我们采用二次指数平滑模型。但我们也可以从数据变换的角度来改进，即在使用指数平滑法前先对数据做一些技术上的处理，让处理后的数据能适合于一次指数平滑模型，然后再对输出结果做一定的数据处理，使其恢复为原变量的形态。（可以联系自动控制原理中的拉普拉斯变换的思想）差分方法是改变数据变动趋势的简易方法。差分概念的介绍

1.一阶差分指数平滑法

当时间序列呈直线增加时，可以采用一阶差分指数平滑模型来预测：
▽ y t = y t − y t − 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ( 1 ) \triangledown y_{t}=y_{t}-y_{t-1} ..................................(1) ▽yt​=yt​−yt−1​..................................(1) ▽ y ^ t + 1 = α ▽ y t + ( 1 − α ) ▽ y t ^ . . . . . . . . . . . . . . . . ( 2 ) \triangledown \hat{y}_{t+1}=α\triangledown y_{t}+(1-α)\triangledown \hat{y_{t}}................(2) ▽y^​t+1​=α▽yt​+(1−α)▽yt​^​................(2) y ^ t + 1 = ▽ y ^ t + 1 + y t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ( 3 ) \hat{y}_{t+1}=\triangledown \hat{y}_{t+1}+y_{t}................................(3) y^​t+1​=▽y^​t+1​+yt​................................(3)
(1)表示把呈直线增加的序列做一阶差分，构成一个平稳的新序列；
（3）表示把经过一阶差分后新序列的指数平滑预测值与变量的当前值做叠加，作为变量下一期的预测值。
例：预测1987年的燃料消耗量。

程序：

clc,clear
yt=load('ranliao.txt'); %实际燃料消耗量数据以列向量的方式存放在纯文本文件中
n=length(yt); alpha=0.4; 
dyt=diff(yt); %求yt的一阶向前差分
dyt=[0;dyt]; %这里使用的是一阶向后差分，加“0”补位
dyhat(2)=dyt(2); %指数平滑值的初始值
for i=2:n
    dyhat(i+1)=alpha*dyt(i)+(1-alpha)*dyhat(i);
end
for i=1:n
    yhat(i+1)=dyhat(i+1)+yt(i);  
end
yhat
xlswrite('ranliao.xls',[yt,dyt]) 
xlswrite('ranliao.xls',[dyhat',yhat'],'Sheet1','C1')

预测值：

yhat =

         0   26.0000   28.0000   28.6000   32.1600   34.0960   34.6576   38.1946   42.9167   43.1500   46.4900

1987年预测值为46.49.

2.二阶差分指数平滑模型

时间序列呈二次曲线增长时，采用二阶差分指数平滑模型来预测：
▽ y t = y t − y t − 1 \triangledown y_{t}=y_{t}-y_{t-1} ▽yt​=yt​−yt−1​ ▽ 2 y t = ▽ y t − ▽ y t − 1 \triangledown ^2y_{t}=\triangledown y_{t}-\triangledown y_{t-1} ▽2yt​=▽yt​−▽yt−1​ ▽ 2 y ^ t + 1 = α ▽ 2 y t + ( 1 − α ) ▽ 2 y t ^ \triangledown ^2\hat{y}_{t+1}=α\triangledown ^2y_{t}+(1-α)\triangledown ^2\hat{y_{t}} ▽2y^​t+1​=α▽2yt​+(1−α)▽2yt​^​ y ^ t + 1 = ▽ 2 y ^ t + 1 + ▽ y t + y t \hat{y}_{t+1}=\triangledown ^2\hat{y}_{t+1}+\triangledown y_{t}+y_{t} y^​t+1​=▽2y^​t+1​+▽yt​+yt​

四、具有季节特点的时间序列的预测

季节，可以是自然季节，也可以是某产品的销售季节等。可以采用季节系数法进行预测，步骤为：

1. 收集m年的每年各季度或各月份（每年n个季度）的时间序列样本数据 a i j a_{ij} aij​。 i i i为年份， j j j为季度或月份。
2. 计算每年所有季度或月份的算术平均值 a ˉ \bar{a} aˉ。
3. 计算所有年份的同季度或月份的算术平均值 a ˉ . j \bar{a}_{.j} aˉ.j​。
4. 计算季度系数或月份系数 b j = a ˉ . j a ˉ b_{j}=\frac{\bar{a}_{.j}}{\bar{a}} bj​=aˉaˉ.j​​。
5. 预测计算：
   先求预测年份的年加权平均：
   y m + 1 = ∑ i = 1 m w i y i ∑ i = 1 m w i y_{m+1}=\frac{\sum_{i=1}^{m}w_{i}y_{i}}{\sum_{i=1}^{m}w_{i}} ym+1​=∑i=1m​wi​∑i=1m​wi​yi​​其中： y i = ∑ i = 1 n a i j y_{i}=\sum_{i=1}^{n}a_{ij} yi​=∑i=1n​aij​为第 i i i年的年合计数， w i w_{i} wi​为第 i i i年的权值 （ w i = i ） （w_{i}=i） （wi​=i）
   
   再计算预测年份的季度平均值：
   y ˉ m + 1 = y m + 1 n \bar{y}_{m+1}=\frac{y_{m+1}}{n} yˉ​m+1​=nym+1​​最后预测该年底 j j j季度的预测值：
   y m + 1 , j = b j y ˉ m + 1 y_{m+1,j}=b_{j}\bar{y}_{m+1} ym+1,j​=bj​yˉ​m+1​
   例：预测2004年各季度的销售额。

程序:

clc, clear
format long g
a=load('jijie.txt');
[m,n]=size(a);
a_mean=mean(mean(a));  %计算所有数据的算术平均值
aj_mean=mean(a);  %计算同季节的算术平均值
bj=aj_mean/a_mean  %计算季节系数
w=1:m;
yhat=w*sum(a,2)/sum(w)  %预测下一年的年加权平均值,这里是求行和
yjmean=yhat/n  %计算预测年份的季节平均值
yjhat=yjmean*bj  %预测年份的季节预测值
format   %恢复默认的显示格式

预测值：

yjhat =

          145572.530322036            201169.1663137          272694.982172952          183904.321191313