一维平稳时间序列预报MATLAB实现

一维平稳时间序列预报

资料处理

​ 本次实验使用了课本某站15年6-8月降水量资料和长江中下游近30年5月平均降水量资料。将原始资料处理为距平序列。

计算各阶$r(\tau )$

​ 根据公式：
$$r(\tau )=\frac{\frac{1}{n-\tau }\sum _{i=1}^{n-\tau }(x_i-\bar{x})(x_{i+\tau }-\bar{x})}{\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x}{)}^2}$$
​ 分别计算$r(\tau )$（其中，一般取$\tau ⩽\frac{n}{2}$），并选出$\left|r(\tau )\right|$较大者所对应前期时刻的要素变量作为因子。

由正规方程组确定预报方程偏回归系数

​ 假设发现有$n-k,n-l,n-m$三个时间间隔的$\left|r(\tau )\right|$较其他值大得多，则取这三个前期变量作为因子，建立自回归方程：
$${\hat{x}}_{d,n}=b_k{x}_{d,n-k}+b_l{x}_{d,n-l}+b_m{x}_{d,n-m}$$
​ 仿照多元线性回归求偏回归系数的方法，得到求解自回归系数的正规方程组为：
$$\{\begin{array}{l}{b}_{k}R(k-k)+b_{l}R(k-l)+b_{m}R(k-m)=R(k)\\ b_{k}R(l-k)+b_{l}R(l-l)+b_{m}R(l-m)=R(k)\\ b_{k}R(m-k)+b_{l}R(m-l)+b_{m}R(m-m)=R(k)\end{array}$$
​ 最终可解出系数，建立自回归预报方程。

结果分析

​ 从图中可发现，计算结果比实测值要更接近于总体平均，即计算振幅比实际要小。趋势上总体保持了一致，但振幅偏小，且在某些区间出现了不同的偏离特征。总的来说，趋势预报的效果较好，但对于量值定量预报效果较差。
这是由于自回归方程是以最小二乘法为基础的，在正态平稳假设下，方程的计算值是在各因子出现在某一数值条件下预报对象的条件期望值；另一方面，由于有平稳遍历性假设，因此该法不能预报出特大值和特小值，虽然是定量计算，但也只能作为趋势预报使用效果较好。

clear all
a = 'YRV_JUN.dat';
fid=fopen(a);
X=fread(fid,[30 1],'float32');
fclose(fid);
x = X-mean(X);
x=[-5,-10,-10,-2,9,3,3,-1,-17,-13,-9,19,19,3,11];
ave = mean(x);
Dx=mean(((x-ave).^2));
n=length(x);
for tao = 1:floor(n/2)
    sum = 0;
    for i = 1:(n-tao)
        sum = sum + (x(i)-ave)*(x(i+tao)-ave);
    end
    r(tao)=sum/((n-tao)*Dx);
end
[B,I]=sort(abs(r),'descend');
I=num2cell(sort(I(1:3)));
[k l m]=deal(I{1:3});

A=[1,r(abs(k-l)),r(abs(k-m));
    r(abs(l-k)),1,r(abs(l-m));
    r(abs(m-k)),r(abs(m-l)),1];
B=[r(k);r(l);r(m)];
b=A\B;

for i = m+1:n
    x_pie(i)=b(1)*x(i-k)+b(2)*x(i-l)+b(3)*x(i-m);
    t(i)=i;
end

plot(t(m+1:end),x(m+1:end),'-','LineWidth',2);
hold on
plot(t(m+1:end),x_pie(m+1:end),'--','LineWidth',2);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
ax.Box = 'off';
legend('实际值','计算值')
xlabel('序号');ylabel('x');

参考资料：
《统计天气预报原理与方法》孔玉寿 钱建明 臧增亮 编著