预测方法
15 预测方法
15.1 微分方程模型
以 15.1 美日硫磺岛战役模型为例
微分方程如下
代码如下
dxy=@(t,x)[-0.0544*x(2)+54000*(t>=0&t<1)+6000*(t>=2&t<3)+13000*(t>=5&t<6);-0.0106*x(1)];
[t,xy]=ode45(dxy,[0:36],[0,21500])
subplot(211);
plot(t,xy(:,1),'r*',t,xy(:,2),'gD');
subplot(212);
plot(xy(:,1),xy(:,2));
其中 dxy=@(t,x)[] 是构造了一个具有 \(t,x\) 的函数表达式。方括号内用分号隔开了两个式子，分别表示两个微分方程（注意微分形式一定要在左边，方程个数和微分方程个数相同）
ode45是用于求微分方程数值解的函数
\([t,x]=ode45(Fun,tspan,x0)\) ，其中 \(Fun\) 表示函数名， \(tspan\) 表示求解区间（例中就是 \(t\) 的取值）， \(x0\) 表示初始值（例中就是0和21500），它返回的 \(x\) 是一个 \(t\times2\) 的矩阵，用于存放不同 \(t\) 对应的两个变量。
由解出的 \(x\) ，就可以画出微分方程组的轨线。
15.2 灰色预测模型
灰色模型 \(Grey Model,GM\) ，根据较少的数据得到近似指数规律再进行建模的方法，但是只适用于中短期并且满足指数增长的预测。
累加生成：
令 \(x^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n))\) 为原始序列
其生成序列为\(x^{(1)}=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(n))\)
其中 \(x^{(1)}(k)=\sum_{i=1}^{k}x^{(0)}(i)\) ，也就是 \(x^{(1)}\) 是 \(x^{(0)}\) 的一个前缀和序列
\(r\) 次累加生成则有 \(x^{(r)}(k)=\sum_{i=1}^{k}x^{(r-1)}(i)\)
累加生成可以让非负数据变成连续单增的序列。如果有负数，可以进行移轴操作。
累减生成：
令 \(x^{(r)}\) 表示 \(r\) 次生成序列
有 \(x^{(r-1)}(k)=x^{(r)}(k)-x^{(r)}(k-1)\)
其实就是累加的逆运算，也可以叫做差分
可以用于生成序列求原序列
均值生成：
记 \(z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(n))\) 为 \(x^{(1)}\) 的均值生成序列
当 \(z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1),k=2,3,...,n\)
GM(1,1)模型
即一阶微分方程且含有1个变量的灰色模型
首先检验数据是否满足该模型的要求
计算数列的级比
当所有的级比都落在区间 \(\Theta=(e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}})\) 时，才可以作为 \(GM(1,1)\) 的数据进行灰色预测
否则，需要对 \(x^{(0)}(k)\) 进行平移变换，使得它满足这个要求
设时间序列 \(x^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n))\) 有 \(n\) 个观察值，通过累加生成新序列\(x^{(1)}=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(n))\)
则 \(GM(1,1)\) 模型建立的白化形式微分方程为
其中 \(a\) 称为发展灰数， \(b\) 称为内生控制灰数
令
设 \(u\) 为待估参数向量
则微分方程可以表示为
利用最小二乘法求解得到
求解微分方程得到预测模型
然后进行差分操作，得到原始序列的灰色预测模型
两个参数中， \(a\) 主要控制系统发展态势的大小，其中
①当 \(-a<0.3\) 时， \(GM(1,1)\) 模型可用于中长期预测
②当 \(0.3<-a<0.5\) 时， \(GM(1,1)\) 模型可用于短期预测
③当 \(0.5<-a<0.1\) 时， 应该使用 \(GM(1,1)\) 改进模型，包括 \(GM(1,1)\) 残差修正模型
④当 $-a>1 $ 时，不宜采用 \(GM(1,1)\) 模型
模型检验
(1)残差检验
按预测模型计算 \(\hat x^{(1)}\) ，并通过累减得到 \(\hat x^{(0)}\) ，然后计算原始序列 \(x^{(0)}\) 与 $ \hat x^{(0)}$ 的绝对误差序列以及相对误差序列
给定 \(a\) ，当 \(\overline\Phi<a\) 成立，则模型为残差合格模型，一般 \(a=0.01,0.05,0.1\) 为优，合格，勉强合格
(2)级比偏差值检验
先计算出级比 \(\lambda (k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)},k=2,3,...,n\)
然后求出级比偏差
一般 \(|\rho(k)|<0.2\) 认为合格，$ \rho(k)<0.1$ 认为优秀
例题代码：
format long
x0 = [71.1,72.4,72.4,72.1,71.4,72.0,71.6];
len = length(x0);
%累加得到x1
x1(1) = x0(1);
for i=2:len
x1(i) = x1(i - 1) + x0(i);
end
%构造矩阵
for i=2:len
z1(i) = 0.5 * x1(i - 1) + 0.5 * x1(i);
end
B = [(-z1(2:len))', ones(len - 1, 1)];
Y = (x0(2:len))';
%代入求解
u=inv(B' * B) * B' * Y;
a = u(1);
b = u(2);
%得到预测值
x11(1) = x0(1);
for i = 1:(len - 1)
x11(i + 1) = (x0(1) - b / a) * exp(-a * i) + b / a;
end
x00(1) = x11(1);
for i=2:len
x00(i) = x11(i) - x11(i - 1);
end
%检验
for i = 1:len
delta(i) = abs(x0(i) - x00(i));
phi(i) = delta(i) / x0(i);
end
phi
GM(2,1)模型
令 \(x^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n))\) 为原始序列
其1次累加生成序列为\(x^{(1)}=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(n))\)
1次累减生成序列为 \(\alpha^{(1)} x^{(0)}=(\alpha^{(1)} x^{(0)}(2),\alpha^{(1)} x^{(0)}(3),...,\alpha^{(1)} x^{(0)}(n))\)
\(x^{(1)}\) 均值生成序列为 \(z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(n))\)
则 \(GM(2,1)\) 建立的白化方程为
令
设 \(u\) 为待估参数向量
利用最小二乘法求解得到
求解微分方程得到预测模型
clear;
x0 = [41,49,61,78,96,104]';
n = length(x0);
%cumsum函数可以求累加
x1 = cumsum(x0);
%diff函数对于连续变量可以求导 离散变量求差分
ax0 = diff(x0);
Y = ax0;
z = 0.5 * (x1(2:end) + x1(1:end-1));
B = [-x0(2:end), -z, ones(n-1, 1)];
%由于Y=Bu，所以u=B\Y，\代表左乘即在左边乘上逆元
u = B \ Y;
%定义函数
syms x(t)
x = dsolve(diff(x,2) + u(1) * diff(x) + u(2) * x == u(3), x(0) == x1(1), x(5) == x1(6));
%设置精度为六位
xt = vpa(x,6);
%subs用于替换变量，这里是将函数x的自变量t替换成0:n-1的范围
yuce = subs(x, t, 0:n-1);
%用double转换以后才能变成数值类型
format short
yuce = double(yuce)
DGM(2,1)模型
令 \(x^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n))\) 为原始序列
其1次累加生成序列为\(x^{(1)}=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(n))\)
1次累减生成序列为 \(\alpha^{(1)} x^{(0)}=(\alpha^{(1)} x^{(0)}(2),\alpha^{(1)} x^{(0)}(3),...,\alpha^{(1)} x^{(0)}(n))\)
则 \(DGM(2,1)\) 建立的白化方程为
令
设 \(u\) 为待估参数向量
利用最小二乘法求解得到
代码可以在上面的基础上改动得到
clear;
x0 = [2.874,3.278,3.39,3.679,3.77,3.8]';
n = length(x0);
%cumsum函数可以求累加
x1 = cumsum(x0);
%diff函数对于连续变量可以求导 离散变量求差分
ax0 = diff(x0);
Y = ax0;
B = [-x0(2:end), ones(n-1, 1)];
%由于Y=Bu，所以u=B\Y，\代表左乘即在左边乘上逆元
u = B \ Y;
%定义函数 初始条件
syms x(t)
d2x = diff(x,2);dx=diff(x);
x = dsolve(d2x + u(1) * dx == u(2), x(0) == x0(1), dx(0) == x0(1));
%设置精度为六位
xt = vpa(x,6);
%subs用于替换变量，这里是将函数x的自变量t替换成0:n-1的范围
yuce = subs(x, t, [0:n-1]);
%用double转换以后才能变成数值类型
yuce = double(yuce)
x0hat = [yuce(1), diff(yuce)]
epsilon = x0 - x0hat'
delta = abs(epsilon./x0)
灰色Verhulst预测模型
描述具有饱和状态（S形）的过程，比如人口预测、生物生长、产品经济寿命
令 \(x^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n))\) 为原始序列
其1次累加生成序列为\(x^{(1)}=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(n))\)
\(x^{(1)}\) 均值生成序列为 \(z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(n))\)
则 \(Verhulst\) 建立的白化方程为
令
设 \(u\) 为待估参数向量
利用最小二乘法求解得到
代码
clear;
x0 = [4.93,2.33,3.87,4.35,6.63,7.15,5.37,6.39,7.81,8.35]';
n = length(x0);
%cumsum函数可以求累加
x1 = cumsum(x0);
Y = x0(2:n);
z = 0.5 * (x1(2:end) + x1(1:end-1));
B = [-z, z.^2];
%由于Y=Bu，所以u=B\Y，\代表左乘即在左边乘上逆元
u = B \ Y;
%定义函数 初始条件
syms x(t)
dx=diff(x);
x = dsolve(dx + u(1) * x == u(2) * x^2, x(0) == x0(1));
%设置精度为六位
xt = vpa(x,6);
%subs用于替换变量，这里是将函数x的自变量t替换成0:n-1的范围
yuce = subs(x, t, [0:n-1]);
%用double转换以后才能变成数值类型
yuce = double(yuce)
x0hat = [yuce(1), diff(yuce)]
epsilon = x0 - x0hat'
delta = abs(epsilon./x0)
15.4 时间序列
平稳性Daniel检验
在时间序列中，平稳性是很多时间序列模型估计的前提。平稳性本质上反应了一个经济系统能够在一个状态下持续运行，这也是对经济系统进行预测的前提。
时间序列模型的根本目的是进行预测，只有当时间序列数据是平稳的，即时间序列的状态能够在长期保持不变，才能够基于时间序列过去和现在的状态来对未来走势进行预测。否则，如果时间序列数据不平稳，数据属于随机游走过程、状态一直在转变，那么也就无法对其进行预测。
秩统计量
设 \(x_1,x_2,...,x_n\) 为抽取的容量为 \(n\) 的样本，将其从小到大排序得到 \(x_{(1)},x_{(2)},...,x_{(n)}\) ，若 \(x_i=x_{(k)}\) ，则称 \(k\) 是 \(x_i\) 在样本中的秩，记为 \(R_i\) （也就是第 \(i\) 个样本在整体中的排名），序列 \(R_1, R_2,...,R_n\) 总称为秩统计量。
Spearman 相关系数
对于二维总体 \((X,Y)\) 的样本观测数据 \((x_1,y_1),(x_2,y_2),...,(x_n,y_n)\) ，设 \(x_1,x_2,...,x_n\) 的秩统计量为 \(R_1, R_2,...,R_n\) ，\(y_1,y_2,...,y_n\) 的秩统计量为 \(S_1, S_2,...,S_n\) ，则 \(Spearman\) 相关系数为
看起来很眼熟？实际上就是概率论中学过的相关系数 \(\rho\)，只不过这里是对于二维样本的每一维的排名求相关系数。
经过运算有：
Daniel检验
对于时间序列 \(X_i\) 的样本 \(a_1,a_2,...,a_n\) ，记 \(a_t\) 的秩为 \(R_t=R(a_t)\) ，则变量对 \((t,R_t),t=1,2,...,n\) 的 \(Spearman\) 相关系数为
构造统计量
\(T\) 服从自由度为 \(n-2\) 的 \(t\) 分布 \(t(n-2)\)
对于显著水平 \(\alpha\) ，若 \(|T|>t_{\alpha/2}(n-2)\) ，则认为序列非平稳，且当 \(q_s>0\) 时，认为序列有上升趋势，当 \(q_s<0\) 时，认为序列有下降趋势
若 \(|T|\leq t_{\alpha/2}(n-2)\) ，则认为序列平稳
ARMA 模型
自回归模型 AR(p)
表示时间序列的当前值与过去 \(p\) 个值的线性组合
移动平均模型 MA(q)
表示时间序列的当前值与过去 \(q\) 个滞后误差的线性组合
其中 \(x_t\) 表示时间序列的当前值，\(\epsilon_t\) 表示当前时刻的误差，\(\phi_i,\theta_i\) 分别表示自回归系数和移动平均系数
自回归移动平均模型 ARMA(p,q)
如何判断 \(ARMA\) 模型的阶数？
使用 \(AIC\) 准则
看半天没看懂，直接说用法吧...
clc;clear;close all;
data = [16.88 15.91 15.92 16.44 16.44 16.87 17.00 17.35 17.50 17.94...
18.52 18.62 19.39 17.45 17.48 16.28 16.49 15.65 15.80 16.35 16.91 16.57 16.42...
17.32 17.25 17.04 16.40 16.43 16.21 16.16 16.39 16.24 17.04 17.24 16.97 16.84...
17.05 17.41 17.48 18.03 18.95 18.80 18.76 18.29 18.10 18.61 18.83 19.18 19.00...
18.86 18.67 19.23 20.0 20.24 19.65 19.63 ];
train_data = data(1:4*11);
test_data = data(4*11 + 1:end);
%%
figure(1);
plot(train_data,'LineWidth',1)
pingwen = adftest(train_data);%用于直接判断数据是否平稳，是1否0
xlabel('day');
ylabel('price');
if pingwen == 1
disp('平稳序列')
train_data1 = train_data;
else
disp('不平稳序列')
figure(2);
train_data1 = diff(train_data);
plot(train_data1,'LineWidth',1);
end;
%% AIC准则 没看懂但是直接用
lim = round(length(train_data)/10);
if lim > 10
lim = 10;
end
train_data2 = iddata(train_data');%将double类型数据转化为id类型
save_data = [];
for p=1:lim
for q=1:lim
num = armax(train_data2,[p,q]);
AIC = aic(num);
save_data = [save_data;p q AIC];
reli(p,q) = AIC;
end
end
figure(5);%画热力图，值越小越好
for i=1:lim
y_index(1,i) = {['AR',num2str(i)]};
x_index(1,i) = {['MA',num2str(i)]};
end
H = heatmap(x_index, y_index, reli, 'FontSize', 12, 'FontName', '宋体');
H.Title = 'AIC定阶热力图'
min_index = find(save_data(:,3) == min(save_data(:,3)));
p_best = save_data(min_index,1);
q_best = save_data(min_index,2);
%% 构造arma模型
model = armax(train_data2,[p_best,q_best])
L=length(test_data);
pre_data = [train_data1'; zeros(L,1)];
pre_data1 = iddata(pre_data);
pre_data2 = predict(model, pre_data1, L );%predict函数用于预测
pre_data3 = get(pre_data2);%得到结构体
pre_data4 = pre_data3.OutputData{1,1}(length(train_data1)+1:length(train_data1)+L);%从结构体中得到数据
data1 = [train_data1';pre_data4];
if pingwen == 0 %不平稳 还原差分值
data_pre1 = cumsum([train_data(1);data1]);
else
data_pre1 = data1;
end
%%
data_pre2 = data_pre1(length(train_data) + 1:end);%预测数据
figure(6);
subplot(2,1,1);
plot(1:length(train_data),train_data,'--','LineWidth',1);
hold on;
plot(length(train_data) + 1:length(train_data) + L, test_data, 'b--', 'LineWidth',1.5);
hold on;
plot(length(train_data) + 1:length(train_data) + L, data_pre2, 'g--', 'LineWidth',1.5);
xlabel('time');
ylabel('price');
legend('真实值','测试数据真实值','预测值');
wucha = sum(abs(data_pre2'-test_data) ./ test_data ./ length(data_pre2));
title_str = ['ARMA预测相对误差为', num2str(wucha)];
title(title_str);
subplot(2,1,2);
plot(1:L,test_data,'--o',LineWidth=1.5);
hold on;
plot(1:L,data_pre2,'--*',LineWidth=1.5);
hold on;
xlabel('time');
ylabel('price');
legend('真实值','预测值');
title_str = ['ARMA预测相对误差为', num2str(wucha)];
title(title_str);
浙公网安备 33010602011771号