灰色线性回归组合预测模型及R语言实现

1 概述

灰色线性回归组合模型可以改善线性回归模型中不含指数增长及GM(1,1)模型中不含线性因素的状况，

该组合既适合于具有指数增长趋势的序列又适合于具有线性趋势的序列。利用了这两个单一模型的有用信息，克服各自的缺陷，从而提高模型预测的精确度。

2 建模步骤

2.1 GM(1,1)预测模型的构建

GM(1,1)预测模型的构建详细见：灰色预测 GM(1,1) 模型及R语言实现

若数据序列通过GM(1,1)模型方法的可行性检验

设原始序列$X^{\left(0\right)}=\left(x^{\left(0\right)}\left(1\right),x^0\left(2\right),...,x^0\left(n\right)\right)$满足上述条件，其中$$x^{\left(0\right)}\left(k\right)\ge 0$$

（1）序列的累加处理
将原始序列进行一次累加生成后处理后，生成 $X^{\left(0\right)}$的1-AGO 序列 (累加生成序列)：
$$X^{\left(1\right)}=\left(x^{\left(1\right)}\left(1\right),x^{\left(1\right)}\left(2\right),...,x^{\left(1\right)}\left(n\right)\right)$$
其中
$$x^{\left(1\right)}\left(k\right)=\sum _{i=1}^k{x}^{\left(0\right)}\left(i\right)k=1,2,......,n$$

（2）计算紧邻均值序列$$Z^{\left(1\right)}=\left(z^{\left(1\right)}\left(2\right),z^{\left(1\right)}\left(3\right),...,z^{\left(1\right)}\left(n\right)\right)$$
其中$$z^{\left(1\right)}\left(k\right)=\frac{1}{2}\left(x^{\left(1\right)}\left(k\right)+x^{\left(1\right)}\left(k-1\right)\right)$$

（3）建立相关模型
建立一阶微分线性方程，即灰色微分方程，得到 GM(1,1)模型的均值形式：$$x^{\left(0\right)}\left(k\right)+a{z}^{\left(1\right)}\left(k\right)=b$$

通过 GM(1,1) 模型相应的白化微分方程：
$$\frac{{\text{dx}}^{\left(1\right)}}{dt}+a{x}^{\left(1\right)}=b$$

其中 a 表示发展系数，反映了${\hat{x}}^{\left(1\right)}$ 和${\hat{x}}^{\left(0\right)}$ 的发展态势； b 表示灰色作用量，或者内生控制灰数，是从行为序列中挖掘出来的数据，反映的是数据变化的关系，其确切内涵是灰的。

（4）计算 a、b
设 $\hat{a}$ 为待估参数向量，令 $\hat{a}=\left[\begin{array}{c}\text{a}\\ b\end{array}\right]$，利用最小二乘法求解，可得：
$$\hat{a}={\left(B^{T}B\right)}^{-1}{B}^{T}Y$$其中B、Y 分别为
$$\text{B}=\left[\begin{array}{cc}-z^{\left(1\right)}\left(2\right)& 1\\ -z^{\left(1\right)}\left(3\right)& 1\\ ⋮& ⋮\\ -z^{\left(1\right)}\left(n\right)& 1\end{array}\right]\text{，Y}=\left[\begin{array}{c}{x}^{\left(0\right)}\left(2\right)\\ x^{\left(0\right)}\left(3\right)\\ ⋮\\ x^{\left(0\right)}\left(\text{n}\right)\end{array}\right]$$

（5）建立 GM(1,1) 的时间响应式，即预测模型
$$\{\begin{array}{l}{\hat{x}}^{\left(1\right)}\left(k+1\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{b}{a}\right)e^{-ak}+\frac{b}{a},\text{k}=1,2\cdots ,\text{n}\\ \\ {\hat{x}}^{\left(0\right)}\left(k+1\right)={\hat{x}}^{\left(1\right)}\left(k+1\right)-{\hat{x}}^{\left(1\right)}\left(k\right)\end{array}$$

---

这里GM(1,1)时间响应式
$${\hat{x}}^{\left(1\right)}\left(k+1\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{b}{a}\right)e^{-ak}+\frac{b}{a},\text{k}=1,2\cdots ,\text{n}$$可以写成：
$${\hat{x}}^{\left(1\right)}\left(k+1\right)=c_1{e}^{vk}+c_2$$

$v,c_1,c_2$为待定系数

2.2 线性回归模型

回归预测法，是根据自变量与因变量的回归方程进行因果推算的一种方法。

线性回归模型我也做过总结：[线性模型总结] 线性回归+方差分析+协方差分析+混合效应+面板数据模型

简单线性回归方程：

$$y=ax+b$$

2.3 灰色线性回归组合预测模型建立

由于 $X^{\left(0\right)}$ 是非负序列，故 $X^{\left(1\right)}$ 是递增序列，而常见的增长方程是线性方程及指数方程：

$$(\begin{array}{l}线性回归方程\text{：}y=ax+b\\ \\ 指数方程\text{：}y=a{e}^x\end{array}$$

用这两个增长方程的和来拟合累加生成序列$X^{\left(1\right)}\left(k\right)$，从而，灰色线回归组合模型可写成：

$${\hat{x}}^{\left(1\right)}\left(k\right)=c_1{e}^{vk}+c_{2}k+c_3$$

上式中 $v\text{、}{c}_1\text{、}{c}_2\text{、}{c}_3$为待定系数。

---

3 模型未知参数求解

该模型，需要先确定待识别参数 $v$

设
$$Z\left(k\right)=X^{\left(1\right)}\left(k+1\right)-X^{\left(1\right)}\left(k\right)=c_1{e}^{vk}\left(e^v-1\right)+c_2$$

这里，$k=1,2,\cdots ,n-1$

令
$$Y_m\left(k\right)=Z\left(k+m\right)-Z\left(k\right)=c_1{e}^{vk}\left(e^{vm}-1\right)\left(e^v-1\right)$$这里，m 为辅助参数，$m=1,2,\cdots ,n-3;k=1,2,\cdots ,n-m-2$

得到
$$e^v=\frac{Y_m\left(k+1\right)}{Y_m\left(k\right)}$$

从而，可以得到灰色线性回归组合预测模型中的参数$v$ 的拟合值，记为 ${\hat{v}}_m\left(k\right)$:

$$v_m\left(k\right)=\ln \frac{Y_m\left(k+1\right)}{Y_m\left(k\right)}$$

考虑到$Y_m\left(k\right)$中的数据均是时间序列依次累加后的预测值${\hat{X}}^{\left(1\right)}\left(k\right)$，为便于拟合参数 $v$，将${\hat{X}}^{\left(1\right)}\left(k\right)$全部换成$X^{\left(1\right)}\left(k\right)$

即可得出拟合值${\hat{v}}_m\left(k\right)$，考虑到 $m$ 取不同的值，得到的${\hat{v}}_m\left(k\right)$也不同，因此，可以分别取$m=1,2,\cdots ,n-3$

以所有 $m$ 对应的 ${\hat{v}}_m\left(k\right)$ 值的算数平均值作为参数 $v$ 的估计值 $\hat{v}$

---

对于m=1, 有
$$(\begin{array}{l}{Y}_1\left(k\right)=Z\left(k+1\right)-Z\left(k\right)k=1,2,\cdots ,n-2\\ \\ v_1\left(k\right)=\ln \frac{Y_1\left(k+1\right)}{Y_1\left(k\right)}k=1,2,\cdots ,n-3\end{array}$$

对于m=2, 有

$$(\begin{array}{l}{Y}_2\left(k\right)=Z\left(k+2\right)-Z\left(k\right)k=1,2,\cdots ,n-3\\ \\ v_2\left(k\right)=\ln \frac{Y_2\left(k+1\right)}{Y_2\left(k\right)}k=1,2,\cdots ,n-4\end{array}$$

…

对于m=n-3, 有
$$(\begin{array}{l}{Y}_{n-3}\left(k\right)=Z\left(k+n-3\right)-Z\left(k\right)k=1,2\\ \\ v_{n-3}\left(k\right)=\ln \frac{Y_{n-3}\left(k+1\right)}{Y_{n-3}\left(k\right)}k=1\end{array}$$

---

以上计算 ${\hat{v}}_m\left(k\right)$ 的个数为
$$\left(n-3\right)+\left(n-4\right)+\cdots +2+1=\left(n-2\right)\left(n-3\right)/2$$
计算其平均值，则有

$$\hat{v}=\frac{{\sum }_{m=1}^{n-3}{\sum }_{k=1}^{n-m-2}{v}_m\left(k\right)}{\left(n-2\right)\left(n-3\right)/2}$$

得到 $\hat{v}$ 后，则
$${\hat{X}}^{\left(1\right)}\left(k\right)=c_1{e}^{\hat{v}k}+c_{2}k+c_3$$

---

下面估计 参数$c_1\text{、}{c}_2\text{、}{c}_3$

可使用最小二乘法求得参数估计值${\hat{c}}_1\text{、}{\hat{c}}_2\text{、}{\hat{c}}_3$，计算过程如下：

记

$$\mathbf{c}=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T{\mathbf{X}}^{\left(1\right)}$$

式中，

$$\mathbf{A}=\left[\begin{array}{ccc}{e}^{\hat{v}}& 1& 1\\ e^{2\hat{v}}& 2& 1\\ ⋮& ⋮& ⋮\\ e^{n\hat{v}}& n& 1\end{array}\right]\text{，}{\mathbf{X}}^{\left(1\right)}=\left[\begin{array}{c}{X}^{\left(1\right)}\left(1\right)\\ X^{\left(1\right)}\left(2\right)\\ ⋮\\ X^{\left(1\right)}\left(n\right)\end{array}\right]$$

从而得到生成序列的预测值为：

$${\hat{X}}^{\left(1\right)}\left(k\right)={\hat{c}}_1{e}^{\hat{v}k}+{\hat{c}}_{2}k+{\hat{c}}_3$$

对上式作一次累减可得原序列的预测值 ${\hat{X}}^{\left(0\right)}$
$${\hat{X}}^{\left(0\right)}\left(k+1\right)={\hat{X}}^{\left(1\right)}\left(k+1\right)-{\hat{X}}^{\left(1\right)}\left(k\right)$$

---

综上，对于该模型，

• 若 ${\hat{c}}_2$ ，则模型为GM(1,1)模型；
• 若 ${\hat{c}}_1$ ，则模型为线性回归模型。

4 预测模型精度检验

检验灰色预测模型效果主要包括三种方法：残差检验、关联度检验和后验差检验方法，以此来检验模型的准确性。

这里给出后验差检验的步骤 ，其余检验步骤在前面的链接中有提到：

原始序列的均值
$${\bar{x}}^{\left(0\right)}=\frac{1}{n}\sum _{k=1}^n{x}^{\left(0\right)}\left(k\right)$$

计算的原始序列的标准差：
$$S_1=\sqrt{\frac{\sum {\left[x^{\left(0\right)}\left(k\right)-{\bar{x}}^{\left(0\right)}\right]}^2}{n-1}}$$

残差序列的均值：
$${\Delta }^{\left(0\right)}=\frac{1}{n}\sum _{k=1}^n\left[x^{\left(0\right)}\left(k\right)-{\hat{x}}^{\left(0\right)}\left(k\right)\right]$$
计算绝对误差序列的标准差：
$$S_2=\sqrt{\frac{\sum {\left[{\Delta }^{\left(0\right)}\left(k\right)-{\Delta }^{\left(0\right)}\right]}^2}{n-1}}$$
计算后验差比值为：：
$$C=\frac{S_2}{S_1}$$
计算小误差概率：
$$P=p\left\{\left|{\Delta }^{\left(0\right)}\left(k\right)-{\Delta }^{\left(0\right)}\right|<0.6745S_1\right\}$$

GM(1,1) 通常用后验差检验方法来评价预测结果的好坏，主要根据其中的后验差比值（C）和小误差概率（P）这两个数值来检验模型。根据 C 和 P 的大小可以将灰色模型的预测精度分为以下几个等级，如下表

后验差比值 C	小误差概率 P	预测精度等级
<0.35	>0.95	好（一级）
<0.50	>0.80	较好（二级）
<0.65	>0.70	合格（三级）
>=0.65	<=0.65	不合格 （四级）

通过以上检验，如果模型的后验差比、小误差概率都在允许范围内，则说明建立的模型可行，否则，应对残差进行GM(1,1)建模修正

5 R语言实现

GM(1,1)建模

gm11<-function(x,k)
{
n<-length(x)
x1<-numeric(n);
for(i in 1:n)   ##一次累加
{
x1[i]<-sum(x[1:i]);
}
b<-numeric(n)
m<-n-1
for(j in 1:m)
{
b[j+1]<-(0.5*x1[j+1]+0.5*x1[j])   ##紧邻均值生成
}
Yn=t(t(x[2:n]))                   ##构造Yn矩阵
B<-matrix(1,nrow=n-1,ncol=2)      
B[,1]<-t(t(-b[2:n]))              ##构造B矩阵
A<-solve(t(B)%*%B)%*%t(B)%*%Yn;   ##使用最小二乘法求得灰参数a,u
a<-A[1];
u<-A[2];
x2<-numeric(k);
x2[1]<-x[1];
for(i in 1:k-1)
{
x2[1+i]=(x[1]-u/a)*exp(-a*i)+u/a;
}
x2=c(0,x2);
y=diff(x2);                     ##累减生成，获得预测数据数列
y
}
x<-c(408.40,479.00,574.60,758.00,1055.30)   
gm11(x,length(x))

灰色线性回归组合预测模型

x <- c(1618,1674,1728,1753,1775,1785,1797,1815)
n<-length(x)
x1<-numeric(n);
for(i in 1:n)   ##一次累加
{
x1[i]<-sum(x[1:i]);
}
z = diff(x1) # 
m = n-3
v_m = c()
for(i in 1:m){
y=c()
for(j in 1:(n-i-1)){
y = c(y,z[j+i]-z[j])
}
for(k in 1:(n-i-2)){
v_m = c(v_m,log(y[k+1]/y[k]))
}
}
v = sum(v_m)/length(v_m)
v
A = matrix(c(exp(seq(1,n)*v),seq(1,n),rep(1,n)),nrow = n, ncol = 3)
c<-solve(t(A)%*%A)%*%t(A)%*%x1;
c

# 拟合
x_1 = c()
for(i in 1:n){ 
x_1 = c(x_1,c[1]*exp(v*i)+c[2]*i+c[3])
}
x_0 = c(x_1[1],diff(x_1))
x_0