准备整理几个文章中公开的的灰色模型代码，然后解析一下

Loveumm

已于 2024-06-18 21:05:23 修改

阅读量170

点赞数 3

文章标签： matlab

于 2024-06-17 21:34:12 首次发布

本文链接：https://blog.csdn.net/weixin_45853843/article/details/139755166

版权

Zeng, B., & Li, C. (2016). Forecasting the natural gas demand in China using a self-adapting intelligent grey model. Energy, 112, 810-825.

SIGM

%X is the original sequence for building SIGM model;
X=XX0(1,1:10);
%X1 is the 1-AGO sequence of X;
%X=[29.3,33.9,39.7,46.8,56.1,69.5,80.7,87.5,107.5];
X1=[];
%the element number in sequence X
n=length(X);
s=0;
for k=1:n
    s=s+X(k);
    X1=[X1,s];
end
%x1_k_m_1_2:x1(k-1)^2
x1_k_m_1_2=0;
for k=2:length(X1)
    x1_k_m_1_2=x1_k_m_1_2+X1(k-1)^2;
end
%k_x1_k_m_1: kx1(k-1)
k_x1_k_m_1=0;
for k=2:length(X1)
    k_x1_k_m_1=k_x1_k_m_1+k*X1(k-1);
end
%x1_k_m_1:x1(k-1)
x1_k_m_1=0;
for k=2:length(X1)
    x1_k_m_1=x1_k_m_1+X1(k-1);
end
%k_2: k^2
k_2=0;
for k=2:length(X1)
    k_2=k_2+k^2;
end
%k: the sum of k
k_sum=0;
for k=2:length(X1)
  k_sum=k_sum+k;
end
%x1_k_x1_k_m_1:x1(k)*x1(k-1)
x1_k_x1_k_m_1=0;
for k=2:length(X1)
    x1_k_x1_k_m_1=x1_k_x1_k_m_1+X1(k)*X1(k-1);
end
n_1=length(X1)-1;
k_x1_k=0;
for k=2:length(X1)
    k_x1_k=k_x1_k+k*X1(k);
end
%1_k: x1(k)
x1_k=0;
for k=2:length(X1)
    x1_k=x1_k+X1(k);
end
%Constructing Matrix A0,A1,A2,A3
A0=[x1_k_m_1_2 k_x1_k_m_1 x1_k_m_1;k_x1_k_m_1 k_2 k_sum;x1_k_m_1 k_sum n_1];
A1=[x1_k_x1_k_m_1 k_x1_k_m_1 x1_k_m_1; k_x1_k k_2 k_sum;x1_k k_sum n_1];
A2=[x1_k_m_1_2 x1_k_x1_k_m_1 x1_k_m_1;k_x1_k_m_1 k_x1_k k_sum;x1_k_m_1 x1_k n_1];
A3=[x1_k_m_1_2 k_x1_k_m_1 x1_k_x1_k_m_1;k_x1_k_m_1 k_2 k_x1_k;x1_k_m_1 k_sum x1_k];
%computing coefficient
f1=det(A1)/det(A0);
f2=det(A2)/det(A0);
f3=det(A3)/det(A0);
%compute parameters a, b and c
a=2*(1-f1)/(1+f1);
b=2*f2/(1+f1);
c=2*f3/(1+f1);
%computing simulative values
%Simu: be used to store the simulated values with a modelSimu=0;
Simu=[];
for k=2:length(X1)
    v1=0;v2=0;
    v1=X1(1)*(f1-1)+2*f2+f3;
    v1=v1*f1^(k-2);
    for g=0:k-3
        v2=v2+f2*f1^g;
    end
    Simu=[Simu,v1+v2];
end
%Error: be used to store the error between simulative value and original value.
%MRPE: Mean Relative Percentage Error
Error=[];
MRPE=0;
for k=2:length(X1)
    MRPE=MRPE+abs(Simu(k-1)-X(k))/X(k);
    Error=[Error,Simu(k-1)-X(k)];
end
MRPE/(length(X1)-1)
%Forecasting
Fore=[];
for k=2:length(X1)+3 % forecast steps
    v1=0;v2=0;
    v1=X1(1)*(f1-1)+2*f2+f3;
    v1=v1*f1^(k-2);
    for g=0:k-3
        v2=v2+f2*f1^g;
    end
    Fore=[Fore,v1+v2];
end

Zeng, B., Tong, M., & Ma, X. (2020). A new-structure grey Verhulst model: development and performance comparison. Applied Mathematical Modelling, 81, 522-537.

% This is a matlab programe for buiding the new Verhulst model.
% do not used scientific notation. 
format long g;
% Example 1
X0=XX0(1,1:10);

% Example 2
% X0=[4.1299,5.2597,5.9244,6.2484,6.3918,6.4525,6.4777,6.4881]

% Example 3
% X0=[4.1299,5.2382,5.9666,6.4590,6.3160,6.4389]

len=length(X0);

len_F_E=0; % The nunber of reserved data to test prediction errors
len_F=3;   % Predicting future steps

sim_len=len-len_F_E;

% Computing the reciprocal sequence Y0 of X0
Y0=[];
for k=1:sim_len
   Y0=[Y0,1/X0(k)] ;
end
Y0
% Constructing Y1 by using a cycling operation
Y1=[];
s=0;
for k=1:sim_len
    s=s+Y0(k);
    Y1=[Y1,s];
end

% Constructing Matrix A and B
A=[];
B=[];
for k=2:sim_len
  A=[A;Y1(k)];
end

for k=1:sim_len-1
  B=[B;Y1(k),k,1];
end

% Computing parameters deltal,delta2 and delta3.
Delta=(B'*B)^(-1)*B'*A
deltal=Delta(1);
delta2=Delta(2);
delta3=Delta(3);

% Computing parameters a,b,c and C.
a=-log(deltal)
b=(a*delta2)/(1-exp(-a))
c=(a*delta3)/(1-deltal)+b/a+b/2
C=0;fenzi=0;fenmu=0;
for k=2:sim_len
  fenzi=fenzi+(Y0(k)-b/a)*exp(-a*k);
  fenmu=fenmu+(1-exp(a))*exp(-2*a*k);
end
C=fenzi/fenmu

% Computing simulated/forecasted data of X0
Sim_For=[];
for k=1:len+len_F
    val=a/(a*C*(1-exp(a))*exp(-a*k)+b);
    Sim_For=[Sim_For,roundn(val,-4)];   
end
Sim_For

% Computing simulated errors
% RE_Sim:residual error;
% RPE_Sim:relative percentage error;
% MAPE_Sim:mean simulation absolute percentage error.
RE_Sim=[];
RPE_Sim=[];
MAPE_Sim=0;

for k=2:sim_len
   error=Sim_For(k)-X0(k)
   re_error=abs(error)/X0(k)*100
   RE_Sim=[RE_Sim,error];
   RPE_Sim=[RPE_Sim,re_error];
   MAPE_Sim=MAPE_Sim+re_error;
end
MAPE_Sim=roundn(MAPE_Sim/(sim_len-1),-8);
fprintf('Mean simulation absolute percentage error=%g\n', MAPE_Sim);

% Computing forecasted errors
% RE_For:residual error;
% RPE_For:relative percentage error;
% MAPE_For:mean prediction absolute percentage error.
if len_F_E>0
 RE_For=[];
 RPE_For=[];
 MAPE_For=0;
 for k=sim_len+1:sim_len+len_F_E
   error=Sim_For(k)-X0(k);
   re_error=abs(error)/X0(k)*100;
   RE_For=[RE_For,error];
   RPE_For=[RPE_For,re_error];
   MAPE_For=MAPE_For+re_error;
 end
 MAPE_For=roundn(MAPE_For/len_F_E,-8);
 fprintf('Mean prediction absolute percentage error=%g\n', MAPE_For);
end

% Display the data of the future [len_F] steps 
For_len=[];
if len_F>0
    for k=len+1:len+len_F
        For_len=[For_len,Sim_For(k)];
    end
fprintf('The data of the future [len_F] steps:%g\n', For_len);
end

Ding, S., Hu, J., & Lin, Q. (2023). Accurate forecasts and comparative analysis of Chinese CO2 emissions using a superior time-delay grey model. Energy Economics, 126, 107013.

ITGM

clear;
syms a td;
%%一步预测 
td=[0.123978032611567,0.245219326514059,0.769268887256311,0.289559704969065];
%%两步预测 td=[0.129487739647155,0.270932601828656,0.765263690134714,0.251339273312392];
%%三步预测 td=[0.131283217868704,0.231069817629493,0.770713170101856,0.273072506892378];
A_=[4.62	5.07	5.49	5.62	6	6.52	7.12	7.3	  7.6	7.58	7.45	7.44	7.56	7.82	7.97 8.08]; %系统行为序列拟合数据
X1_=[17.4	16.4	14.9	14.1	13.4	12.3	11.3	11.5	11.5	11.9	12.2	12.6	13.1	13.7	14.4 15.9];
X2_=[111.8429248	109.1583872	105.7875085	102.0035662	124.4080089	126.5794942	123.0959272	128.922639	134.3192514	140.236887	152.5944998	156.2299435	154.8950517	157.8120786	165.39037 182.4326259];
X3_=[40.34614955	39.9103387	40.48219017	42.26902681	45.35623583	46.5561547	46.6601211	46.2252656	46.39894934	45.82395276	43.23480666	42.63137693	43.01329739	43.7934746	43.25110574 43.54395759];
X4_=[768361191	772580688	774309415	774161889	773219748	772108526	778977720	782865417	786673270	791323527	795251107	797667977	799175185	799480075	800020955 792401719];
yt=[4.62	5.07	5.49	5.62	6	6.52	7.12	7.3	  7.6	 7.58	7.45	7.44	7.56	7.82	7.97 8.08]; %系统行为序列检验数据
%%一步预测 
num=1;
%%两步预测 num=2;
%%三步预测 num=3;
m=16;
n=4;%x的行列,n个相关因素
A=A_./A_(1);
X1=X1_./X1_(1);
X2=X2_./X2_(1);
X3=X3_./X3_(1);
X4=X4_./X4_(1);
 for k=1:m-num
        tmp=0;
        for j=1:k
            tmp=tmp+A(j);
        end
     XX0(k)=tmp;
   end
for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X1(j);
        end
     XX1(k)=tmp;
end
for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X2(j);
        end
     XX2(k)=tmp;
end
for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X3(j);
        end
     XX3(k)=tmp;
end
for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X4(j);
        end
     XX4(k)=tmp;
end
%一步预测 
AX1=[XX0 0;XX1;XX2;XX3;XX4];
%两步预测 AX1=[XX0 0 0;XX1;XX2;XX3;XX4];
%三步预测 AX1=[XX0 0 0 0;XX1;XX2;XX3;XX4];
x1=[XX1;XX2;XX3;XX4];
ago0(1)=XX0(1);
for k=2:m-num
    ago0(k)=XX0(k)-XX0(k-1);
end
for i=2:m-num
    Z(i)=(XX0(i)+XX0(i-1))/2;
end 
Z(1)=[];
D=ago0;D(1)=[];
D=D';%特征序列除去第一个值
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(1)^(i-k)*XX1(k);
    end
    l(i)=tmp;
end
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(2)^(i-k)*XX2(k);
    end
    r(i)=tmp;
end
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(3)^(i-k)*XX3(k);
    end
    t(i)=tmp;
end
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(4)^(i-k)*XX4(k);
    end
    y(i)=tmp;
end
l(1)=[];r(1)=[];t(1)=[];y(1)=[];
l=l';r=r';t=t';y=y';
E=[l r t y]; Z=Z';
B=[-Z,E];
c=B'*pinv(B*B')*D;
c=c';
a=c(1);
for i=1:n
  b=c(:,2:n+1);
end
F=[];F(1)=A(1);  % 响应函数
for k=2:m
    tmp=0;
    for i=2:k
    if k<2
        v1=0;
    else
        v1=1;
    end
    tmp=tmp+b(1)*v1*(td(1)^(k-i+1/2)*(1/2)*(XX1(i)+XX1(i-1)));
    end
    tmp1(k)=tmp;
end
 for k=2:m
     tmp=0;
    for i=2:k
    if k<2
        v1=0;
    else
        v1=1;
    end
    tmp=tmp+b(2)*v1*(td(2)^(k-i+1/2)*(1/2)*(XX2(i)+XX2(i-1)));
    end
    tmp2(k)=tmp;
 end
  for k=2:m
     tmp=0;
    for i=2:k
    if k<2
        v1=0;
    else
        v1=1;
    end
    tmp=tmp+b(3)*v1*(td(3)^(k-i+1/2)*(1/2)*(XX3(i)+XX3(i-1)));
    end
    tmp3(k)=tmp;
  end
  for k=2:m
     tmp=0;
    for i=2:k
    if k<2
        v1=0;
    else
        v1=1;
    end
    tmp=tmp+b(4)*v1*(td(4)^(k-i+1/2)*(1/2)*(XX4(i)+XX4(i-1)));
    end
    tmp4(k)=tmp;
  end
for k=1:m
    f(k)=tmp1(k)+tmp2(k)+tmp3(k)+tmp4(k);%内函数
end
for i=2:m
     tmp=0;
    for k=2:i
   tmp=tmp+(1/2)*(f(k)+f(k-1))*exp(-a*(i-k+1/2));
    end
   temp(i)=tmp;
end
for i=2:m
    if i<2
        v2=0;
    else
        v2=1;
    end
    F(i)=A(1)*exp(-a*(i-1))+v2*temp(i);
end
G(1)=A(1)*yt(1);
for i=2:m
    G(i)=((F(i)-F(i-1)))*yt(1); %还原式子
end
%拟合期误差
X_train=yt(1,2:m-num);
Y_train=G(2:m-num);
APE_fit=abs((X_train- Y_train)./X_train)*100;
MAPE_fit=mean(abs((X_train- Y_train)./X_train))*100;
RMSE_fit = sqrt(mean((X_train- Y_train).^2));

%预测期误差
X_test=yt(1,m-num+1:end);
Y_test=G(m-num+1:end);
APE_pred=abs((X_test- Y_test)./X_test)*100;
MAPE_test=mean(abs((X_test- Y_test)./X_test))*100;
RMSE_test = sqrt(mean((X_test- Y_test).^2));

errors=[MAPE_fit,  MAPE_test, RMSE_fit, RMSE_test]
APE=[APE_fit,APE_pred]
forecasts=G
plot(2:m,yt_1(1,:),'s-r',2:m,G_1(1,:),'d-b');

Ding, S., Hu, J., & Lin, Q. (2023). Accurate forecasts and comparative analysis of Chinese CO2 emissions using a superior time-delay grey model. Energy Economics, 126, 107013.

Time_delay GM(1,N)

clear;
syms a td;
%时滞系数
%一步预测 
td=[0.681141736869649,0.911791357586166,0.220941382599340,0.974377735389256];
%两步预测 td=[0.671355567062857,0.804518471837752,0.156309417687540,0.924881797629766];
%三步预测 td=[0.625565069987166,0.815543508932614,0.0228726246083150,0.476368299763061];
A_=[4.62	5.07	5.49	5.62	6	6.52	7.12	7.3	  7.6	7.58	7.45	7.44	7.56	7.82	7.97 8.08]; %系统行为序列拟合数据
X1_=[17.4	16.4	14.9	14.1	13.4	12.3	11.3	11.5	11.5	11.9	12.2	12.6	13.1	13.7	14.4 15.9];
X2_=[111.8429248	109.1583872	105.7875085	102.0035662	124.4080089	126.5794942	123.0959272	128.922639	134.3192514	140.236887	152.5944998	156.2299435	154.8950517	157.8120786	165.39037 182.4326259];
X3_=[40.34614955	39.9103387	40.48219017	42.26902681	45.35623583	46.5561547	46.6601211	46.2252656	46.39894934	45.82395276	43.23480666	42.63137693	43.01329739	43.7934746	43.25110574 43.54395759];
X4_=[768361191	772580688	774309415	774161889	773219748	772108526	778977720	782865417	786673270	791323527	795251107	797667977	799175185	799480075	800020955 792401719];
yt=[4.62	5.07	5.49	5.62	6	6.52	7.12	7.3	  7.6	 7.58	7.45	7.44	7.56	7.82	7.97 8.08]; %系统行为序列检验数据
n=4;%x的行列,n个相关因素
m=16;
%%一步预测 
num=1;
%%两步预测 num=2;
%%三步预测 num=3;
A=A_./A_(1);
X1=X1_./X1_(1);
X2=X2_./X2_(1);
X3=X3_./X3_(1);
X4=X4_./X4_(1);
 for k=1:m-num
     tmp=0;
        for j=1:k
            tmp=tmp+A(j);
        end
     XX0(k)=tmp;
 end
for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X1(j);
        end
     XX1(k)=tmp;
end
   for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X2(j);
        end
     XX2(k)=tmp;
   end
    for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X3(j);
        end
     XX3(k)=tmp;
    end
    for k=1:m
        tmp=0;
        for j=1:k
            tmp=tmp+X4(j);
        end
     XX4(k)=tmp;
    end
%一步预测 
AX1=[XX0 0;XX1;XX2;XX3;XX4];
%两步预测 AX1=[XX0 0 0;XX1;XX2;XX3;XX4];
%三步预测 AX1=[XX0 0 0 0;XX1;XX2;XX3;XX4];
x1=[XX1;XX2;XX3;XX4];
ago0(1)=XX0(1);
for k=2:m-num
    ago0(k)=XX0(k)-XX0(k-1);
end
for i=2:m-num
    Z(i)=(XX0(i)+XX0(i-1))/2;
end 
Z(1)=[];
D=ago0;D(1)=[];
D=D';%特征序列除去第一个值
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(1)^(i-k)*XX1(k);
    end
    l(i)=tmp;
end
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(2)^(i-k)*XX2(k);
    end
    r(i)=tmp;
end
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(3)^(i-k)*XX3(k);
    end
    t(i)=tmp;
end
for i=2:m-num
    tmp=0;
    for k=1:i
    tmp=tmp+td(4)^(i-k)*XX4(k);
    end
    y(i)=tmp;
end
l(1)=[];r(1)=[];t(1)=[];y(1)=[];
l=l';r=r';t=t';y=y';
E=[l r t y]; Z=Z';
B=[-Z,E];
c=B'*pinv(B*B')*D;
c=c';
a=c(1);
for i=1:n
  b=c(:,2:n+1);
end
F=[];F(1)=A(1);  % 响应函数
for k=2:m
    tmp=0;
    for i=1:k
    tmp=tmp+b(1)*td(1)^(k-i)*XX1(i);
    end
    tmp1(k)=tmp;
end
 for k=2:m
    tmp=0;
    for i=1:k
    tmp=tmp+b(2)*td(2)^(k-i)*XX2(i);
    end
    tmp2(k)=tmp;
 end
 for k=2:m
    tmp=0;
    for i=1:k
    tmp=tmp+b(3)*td(3)^(k-i)*XX3(i);
    end
    tmp3(k)=tmp;
 end
for k=2:m
    tmp=0;
    for i=1:k
    tmp=tmp+b(4)*td(4)^(k-i)*XX4(i);
    end
    tmp4(k)=tmp;
end
for k=1:m
    f(k)=tmp1(k)+tmp2(k)+tmp3(k)+tmp4(k);%内函数
end
for i=2:m
    F(i)=(A(1)-f(i)/a)*exp(-(a*(i-1)))+f(i)/a;
end
G(1)=F(1)*yt(1);
for t=2:m
    G(t)=(F(t)-F(t-1))*yt(1); %还原式子
end
%拟合期误差
X_train=yt(1,2:m-num);
Y_train=G(2:m-num);
APE_fit=abs((X_train- Y_train)./X_train)*100;
MAPE_fit=mean(abs((X_train- Y_train)./X_train))*100;
RMSE_fit = sqrt(mean((X_train- Y_train).^2));

%预测期误差
X_test=yt(1,m-num+1:end);
Y_test=G(m-num+1:end);
APE_pred=abs((X_test- Y_test)./X_test)*100;
MAPE_test=mean(abs((X_test- Y_test)./X_test))*100;
RMSE_test = sqrt(mean((X_test- Y_test).^2));

errors=[MAPE_fit,  MAPE_test, RMSE_fit, RMSE_test]
APE=[APE_fit,APE_pred]
forecasts=G
plot(2:m,yt_1(1,:),'s-r',2:m,G_1(1,:),'d-b');