MATLAB时间序列----2019/8/10

时间序列

时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列

• 时间序列分析方法
    当时间序列的数值由于受周期变动和不规则变动的影响，起伏较大，不易显示出发展趋势时，可用移动平均法，消除这些因素的影响，分析、预测序列的长期趋势。
   1. 移动平均法：
     1）简单移动平均法
  
  注：简单移动平均法只适合做近期预测，而且是预测目标的发展趋势变化不大的情况。如果目标的发展趋势存在其它的变化，采用简单移动平均法就会产生较大的预测偏差和滞后。
  
  (此程序有问题，用细胞数组不行，还未找到调整办法，仅将程序打了一遍，如果有人知道此题正确程序，能否告诉我，谢谢。)

y=[533.8 574.6 606.9 649.8 705.1 772.0 816.4 892.7 963.9 1015.1 1102.7];
m=length(y);
n=[4,5];
%n为移动平均的项数

for i=1:length(n)
    %由于n的取值不同，yhat的长度不一致，下面使用了细胞数组
    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);
    s(i)=sqrt(mean((y(n(i)+1:m)-yhat{i}(1:end-1)).^2));
end
y12,s

2）加权移动平均法

       注意此题中的预测值修正的步骤
解：

y=[6.35 6.20 6.22 6.66 7.15 7.89 8.72 8.94 9.28 9.8];
w=[1/6;2/6;3/6];
m=length(y);
n=3;
for i=1:m-n+1
    yhat(i)=y(i:i+n-1)*w
end
%yhat为预测值，yhat(1)的值为预测的i+n-1年的值

err=abs(y(n+1:m)-yhat(1:end-1))./y(n+1:m)
%每年的预测值的相对误差

T_err=1-sum(yhat(1:end-1))/sum(y(n+1:m))
y1989=yhat(end)/(1-T_err)

  注：在加权移动平均法中， wt 的选择，同样具有一定的经验性。一般的原则是：近期数据的权数大，远期数据的权数小。至于大到什么程度和小到什么程度，则需要按照预测者对序列的了解和分析来确定。

3）趋势移动平均法

设时间序列 { y t } \{y_{t}\} {yt​} 从某时期开始具有直线趋势，且认为未来时期也按此直线趋势变
化，则可设此直线趋势预测模型为
$${\hat{y}}_{t+m}=a_t+b_{t}m,m=1,2\cdots$$
其中 $t$ 为当前时期数；$m$ 为由$t$ 至预测期的时期数；$a_t$ 为截距；$b_t$ 为斜率。两者又称为平滑系数。
平滑系数的计算公式为:
$$\{\begin{array}{l}{a}_t=2M_t^{(1)}-M_t^{(2)}\\ b_t=\frac{2}{N-1}(M_t^{(1)}-M_t^{(2)})\end{array}$$
例 3 我国 1965～1985 年的发电总量如表 3 所示，试预测 1986 年和 1987 年的发
电总量。

y=[676 825 774 716 940 1159 1384 1524 1668 1688 1958 2031 2234 2566 2820 3006 3093 3277 3514 3770 4107];
m1=length(y);
n=6;
%n为移动平均的项数

for i=1:m1-n+1
    yhat1(i)=sum(y(i:i+n-1))/n;
end
yhat1
m2=length(yhat1);
for i=1:m2-n+1
    yhat2(i)=sum(yhat1(i:i+n-1))/n;
end
yhat2
plot(1:21,y,'*')
a21=2*yhat1(end)-yhat2(end)
b21=2*(yhat1(end)-yhat2(end))/(n-1)
y1986=a21+b21;
y1987=a21+2*b21

2. 指数平滑法
   1).一次指数平滑法:
   一次指数平滑法公式为：
   $$S_t^{(1)}=\alpha y_t+(1-\alpha )S_{t-1}^{(1)}=S_{t-1}^{(1)}+\alpha (y_t-S_{t-1}^{(1)})$$
   预测模型为以第 $t$ 期指数平滑值作为 $t+1$期预测值：
   $${\hat{y}}_{t+1}=\alpha y_t+(1-\alpha )\widehat{y_t}$$
   用一次指数平滑法进行预测，除了选择合适的 $\alpha$ 外，还要确定初始值 $S_0^{(1)}$初始值
   是由预测者估计或指定的。当时间序列的数据较多，比如在 20 个以上时，初始值对以后的预测值影响很少，可选用第一期数据为初始值。如果时间序列的数据较少，在 20个以下时，初始值对以后的预测值影响很大，这时，就必须认真研究如何正确确定初始值。一般以最初几期实际值的平均值作为初始值
   

y=[50 52 47 51 49 48 51 40 48 52 51 59];
n=length(y);
a=[0.2 0.5 0.8];
%不知如何选择α，所以设3个值

m=length(a);
yhat(1,1:m)=(y(1)+y(2))/2;
for i=2:n
     yhat(i,:)=a*y(i-1)+(1-a).*yhat(i-1,:);
end
%求预测的值
y=y';
%由于y是一行数据，而yhat中一个α 值的各年值是列数据，所以将y转换为列数据。
err=sqrt(mean((repmat(y,1,m)-yhat).^2))
%计算α 的预测标准误差 err ，选取使 err 较小的那个α 值
%repmat（Matrix,M,N）%把矩阵Matrix当成一个元素，以Matrix的内容堆叠在（MxN）的矩阵B中，B矩阵的大小由MxN及Matrix矩阵的内容决定
%由err计算内容可得err(1)数字最小，故选取α = 0.2

yhat1988=a*y(n)+(1-a).*yhat(n,:)
%故预测 1988 年该电器销售额为 yˆ1988 = 51.1754 。

2）二次指数平滑法

yt=[676 825 774 716 940 1159 1384 1524 1668 1688 1958 2031 2234 2566 2820 3006 3093 3277 3514 3770 4107];
 n=length(yt); 
a=0.8;
%由于数据变化较大，所以令a=0.8。

st1(1)=yt(1);
st2(1)=yt(1);
%时间序列的数据较多>20个,所以选用第一期数据为初始值

for i=2:n
    st1(i)=a*yt(i)+(1-a)*st1(i-1);
    st2(i)=a*st1(i)+(1-a)*st2(i-1);
end
at=2*st1-st2;
b=a/(1-a)*(st1-st2);
yhat=at+b;
str=char(['C',int2str(n+2)]);
at(n)+2*b(n)

3). 三次指数平滑法
  当时间序列的变动表现为二次曲线趋势时，则需要用三次指数平滑法。三次指数平滑是在二次指数平滑的基础上，再进行一次平滑，其计算公式为:



（数据在代码中）

n=length(y);
alph=[0.1,0.2,0.3];
m=length(alph);
for i=1:m
    st1(1,i)=mean(y(1:3));
    st2(1,i)=mean(y(1:3));
    st3(1,i)=mean(y(1:3));
end
%求初始值

for i=2:n
    st1(i,:)=alph*y(i)+(1-alph).*st1(i-1,:);
    st2(i,:)=alph*st1(i)+(1-alph).*st2(i-1,:);
    st3(i,:)=alph*st2(i)+(1-alph).*st3(i-1,:);
end
err=sqrt(mean((repmat(y',1,m)-st3).^2));
%求α 的预测标准误差 err ，选取使 err 较小的那个α 值
%得：

alpha=0.3;

%由于alpha=0.3，在alph矩阵中是最后一列，所以以下各公式中的st1,st2,st3均取最后一列。

a=3*st1(:,m)-3*st2(:,m)+st3(:,m); 
b=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1(:,m)-2*(5-4*alpha)*st2(:,m)+(4-3*alpha)*st3(:,m)); 
c=0.5*alpha^2/(1-alpha)^2*(st1(:,m)-2*st2(:,m)+st3(:,m));
yt_1=(3-3*alpha+alpha^2)/(1-alpha)^2*st1(:,m)-(3-alpha)/(1-alpha)^2*st2(:,m)+1/(1-alpha)^2*st3(:,m)
%求y_t+1的值，1989年的预测值是yt_1矩阵的最后一个值

plot(1:n,y,'*',1:n,yt_1(1:n),'o')
%画图

yt_2=a+b*2+c*2^2
%求y_t+2的值，1990年的预测值是yt_2矩阵的最后一个值

在EXCEL中也可求指数平滑（较为便捷，不用编程序）：
具体步骤如下：
（1） 选择工具菜单中的数据—数据分析命令，弹出数据分析对话框。
（2） 在分析工具列表框中，选择指数平滑工具。
（3） 出现指数平滑对话框，输入 $\alpha$ 的值，得到结果

• 差分指数平滑法
  1）一阶差分指数平滑法
  
  
  （代码是我自己写的，书上没有）

y=[24 26 27 30 32 33 36 40 41 44];
n=length(y);
plot(1:n,y,'o')
%由图可知：数据呈直线增长，所以采用一阶差分指数平滑法

alpha=0.4;
yt(1)=0;
for i=2:n
    yt(i)=y(i)-y(i-1);
end
st1(2)=0;
for i=3:n+1
    st1(i)=alpha*yt(i-1)+(1-alpha)*st1(i-1);
end
for i=3:n+1
   y_t(i)=st1(i)+y(i-1);
end
y_t(n+1)

2)二阶差分指数平滑法

• 自适应滤波法
  
  
  例：。设有一个时间序列包括 10 个观测值，如表 9 所示。试用自适应滤波法，以两个权数来求第 11 期的预测值。
  

yt=0.1:0.1:1;
m=length(yt);
k=0.9;
N=2;
Terr=10000;
w=ones(1,N)/N;
while abs(Terr)>0.00001
    Terr=[];
    for j=N+1:m-1
        yhat(j)=w*yt(j-1:-1:j-N)';
        err=yt(j)-yhat(j);
        Terr=[Terr,abs(err)];
        w=w+2*k*err*yt(j-1:-1:j-N);
    end
    %j表示t的取值
    Terr=max(Terr);
end
w,yhat
% w：求出的“最佳”权数
%yhat:预测值
y_end=w(1)*yt(end)+w(2)*yt(end-1)

• 趋势外推预测方法
  1) 指数曲线法
  
  
  
  
  
  解：

y=[42.1 47.5 52.7 57.7 62.5 67.1 71.5 75.7 79.8 83.7 87.5 91.1 94.6 97.9 101.1];
n=length(y);
m=n/3;
for i=2:n-1
    b=(y(i+1)-y(i))/(y(i)-y(i-1));
end
range=minmax(b)
s1=sum(y(1:m));
s2=sum(y(m+1:2*m));
s3=sum(y(2*m+1:end));
b=((s3-s2)/(s2-s1))^(1/m)
a=(s2-s1)*(b-1)/(b*(b^m-1)^2)
k=(s1-a*b*(b^m-1)/(b-1))/m
t1=[1:18];
yt=k+a*b.^t1;

• Compertz 曲线
  
  
  

y=[42.1 47.5 52.7 57.7 62.5 67.1 71.5 75.7 79.8 83.7 87.5 91.1 94.6 97.9 101.1];
y=log(y);
%求导
n=length(y);
m=n/3;
s1=sum(y(1:m));
s2=sum(y(m+1:2*m));
s3=sum(y(2*m+1:end));
b=((s3-s2)/(s2-s1))^(1/m)
a=(s2-s1)*(b-1)/(b*(b^m-1)^2)
k=(s1-a*b*(b^m-1)/(b-1))/m
a=exp(a)
k=exp(k)
t=[1:18]
yt=k*a.^(b.^t);
yt(18)

• Logistic 曲线（生长曲线）

y=[42.1 47.5 52.7 57.7 62.5 67.1 71.5 75.7 79.8 83.7 87.5 91.1 94.6 97.9 101.1];
y=1./y;
n=length(y);
m=n/3;
s1=sum(y(1:m));
s2=sum(y(m+1:2*m));
s3=sum(y(2*m+1:end));
b=((s3-s2)/(s2-s1))^(1/m)
a=(s2-s1)*(b-1)/(b*(b^m-1)^2)
k=(s1-a*b*(b^m-1)/(b-1))/m
t=[1:18];
yt=1./(k+a*b.^t)
yt(18)