1、 残差自相关性分析，以及如何消除残差
详情建模数第六章
clc,clear,close all
y=[20.96;21.40;21.96;21.52;22.39;22.76;23.48;23.66;24.10;24.01;24.54;24.30;25.00;25.64;26.36;26.98;27.52;27.78;28.24;28.78];
x=[127.3;130.0;132.7;129.4;135.0;137.1;141.2;142.8;145.5;145.3;148.3;146.4;150.2;153.1;157.3;160.7;164.2;165.6;168.7;171.7];
rstool(x,y)
%数据Dw检验，相关性检验，函数考虑时间，证明时间是否有滞后性
a=ones(numel(x),1);%有多少个数据，要改
b=0.176*x-1.455*a;%预测函数值，要改
c=y-b % 残差=实际值-预测值
c1=c(2:end,1);%残差e(t)
c2=c(1:end-1,1);%e(t-1)
plot(c1,c2,'*')%残差散点图
xlabel('e(i-1)'),
ylabel('e(i)')
hold on
%画横纵坐标
d=0;
d1=-0.15:0.001:0.25;
plot(d,d1,'.r')
hold on
plot(d1,d,'.r')
% DW的程序：
num=c(1:19)'*c(2:20);
den=sum(c(2:20).^2);
p = num/den%p值
DW=2*(1-p)%DW<dl,残差存在明显自相关,接着消除随机误差项自相关后的回归模型。clc,clear,close all
y=[20.96;21.40;21.96;21.52;22.39;22.76;23.48;23.66;24.10;24.01;24.54;24.30;25.00;25.64;26.36;26.98;27.52;27.78;28.24;28.78];
x=[127.3;130.0;132.7;129.4;135.0;137.1;141.2;142.8;145.5;145.3;148.3;146.4;150.2;153.1;157.3;160.7;164.2;165.6;168.7;171.7];
rstool(x,y)
%数据Dw检验，相关性检验，函数考虑时间，证明时间是否有滞后性
a=ones(numel(x),1);%有多少个数据，要改
b=0.176*x-1.455*a;%预测函数值，要改
c=y-b % 残差=实际值-预测值
c1=c(2:end,1);%残差e(t)
c2=c(1:end-1,1);%e(t-1)
plot(c1,c2,'*')%残差散点图
xlabel('e(i-1)'),
ylabel('e(i)')
hold on
%画横纵坐标
d=0;
d1=-0.15:0.001:0.25;
plot(d,d1,'.r')
hold on
plot(d1,d,'.r')
% DW的程序：
num=c(1:19)'*c(2:20);
den=sum(c(2:20).^2);
p = num/den%p值
DW=2*(1-p)%DW<dl,残差存在明显自相关,接着消除随机误差项自相关后的回归模型。% DW模型优化后的程序：
clc,clear,close all
y=[20.96;21.40;21.96;21.52;22.39;22.76;23.48;23.66;24.10;24.01;24.54;24.30;25.00;25.64;26.36;26.98;27.52;27.78;28.24;28.78];
y1=y(2:length(y),1);
y2=y(1:length(y)-1,1);
y3=y1-0.71133*y2;%改y值,,0.71133为p值
x=[127.3;130.0;132.7;129.4;135.0;137.1;141.2;142.8;145.5;145.3;148.3;146.4;150.2;153.1;157.3; 160.7;164.2;165.6;168.7;171.7];
x1=x(2:length(x),1);
x2=x(1:length(x)-1,1);
x3=x1-0.71133*x2;
a=ones(length(x)-1,1);
y4=y3-(0.1763*x3-0.42*a) % 残差,拟合得式子
y5=y4(2:length(y4),1);%残差e(t)
y6=y4(1:length(y4)-1,1);%残差e(t-1)
y7=y5-y6; %差值
y8=sum(y7.^2);
y9=sum(y5.^2);
DW=y8/y9
p=1-DW/22、时间序列线性二次平均预测法
%读取数据
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%% 数据处理
y=y';
for i=1:numel(y)-2
m1(i+2)=(y(i)+y(i+1)+y(i+2))/3;%一次移动平均
end
for i=3:numel(y)-2
m2(i+2)=(m1(i)+m1(i+1)+m1(i+2))/3;%二次移动平均
end
m1=m1';
m2=m2';
a=2*m1-m2;
b=m1-m2;
%预测方程y(t+T)=at+bt*T,看书532页
at=a(length(y),1)%
bt=b(length(y),1)
for T=1:6
y1(T,:)=at+bt*T;%预测未来数后六位
end
%% 预测值
y1
m2=[m2;y1];
for n=5:length(y)
et(n,:)=abs(((m2(n)-y(n))/y(n))*100);%相对误差分析
end
%% 画图
t1=5:1:length(x);
t2=5:1:length(x)+6;%预测后六年
y=y(5:t1(end));
m2=m2(5:t2(end));
plot(t1,y,'o',t2,m2) %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')clc,clear,close all;
x1=1:1:39;%输入源数据
for i=1:numel(x1)-2
m1(i+2)=(x1(i)+x1(i+1)+x1(i+2))/3;%一次移动平均
end
for i=3:numel(x1)-2
m2(i+2)=(m1(i)+m1(i+1)+m1(i+2))/3;%二次移动平均
end
m1=m1';
m2=m2';
a=2*m1-m2;
b=m1-m2;
%预测方程y(t+T)=at+bt*T,看书532页
at=a(39,1)%
bt=b(39,1)
for T=1:6
y(:,T)=at+bt*T;%预测未来数后六位
end
y3、 灰度预测模型
%读取数据
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%% 处理
syms a b;
c=[a b]';
A=y';%源数据
B=cumsum(A); % 原始数据累加
n=length(A);
for i=1:(n-1)
C(i)=(B(i)+B(i+1))/2; % 生成累加矩阵
end
% 计算待定参数的值
D=A;D(1)=[];
D=D';
E=[-C;ones(1,n-1)];
c=inv(E*E')*E*D;
c=c';
a=c(1);b=c(2);
% 预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a ;
end
G=[];G(1)=A(1);
for i=2:(n+10)
G(i)=F(i)-F(i-1); %得到预测出来的数据，预测十年
end
t1=x';
t2=t1(1):t1(end)+10;
G, a, b % 输出预测值，发展系数和灰色作用量
%% 画图
plot(t1,A,'o',t2,G) %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')4 一次指数平滑
%读取数据
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%一次平滑预测法
n=length(y);
alpha=[0.1,0.3,0.9];%长度
m=length(alpha);
yhat(1,1:m)=(y(1)+y(2))/2;
for i=2:n
yhat(i,:)=alpha*y(i-1)+(1-alpha).*yhat(i-1,:);%一次指数平滑预测法
end
yhat;%预测值
et=(repmat(y,1,m)-yhat);%误差
err=sqrt(mean((repmat(y,1,m)-yhat).^2))%预测误差总和
yhatnaxt=alpha*y(n)+(1-alpha).*yhat(n,:)%预测下一时刻值
%画图
t=x';
plot(t,y,'o',t,yhat(:,3)) %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')二次指数平滑
%% 二次平滑预测法
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%数据处理
yx =y;
yt=yx;
n=length(yt);
alpha=0.8;
st1(1)=yt(1);
st2(1)=yt(1);
for i=2:n
st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
end
a=2*st1-st2;
b=alpha/(1-alpha)*(st1-st2);
yhat=a+b;
yhat=yhat';
err=sqrt(mean((y-yhat).^2));%预测误差总和
%画图
t=x';
plot(t,y,'o',t,yhat) %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')三次指数平滑
%% 三次指数平滑
clc,clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
yx =y;
yt=yx;
n=length(yt);
alpha=0.7;
st1_0=mean(yt(1:3));
st2_0=st1_0;
st3_0=st1_0;
st1(1)=alpha*yt(1)+(1-alpha)*st1_0;
st2(1)=alpha*st1(1)+(1-alpha)*st2_0;
st3(1)=alpha*st2(1)+(1-alpha)*st3_0;
for i=2:n
st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
st3(i)=alpha*st2(i)+(1-alpha)*st3(i-1);
end
st1=[st1_0,st1];
st2=[st2_0,st2];
st3=[st3_0,st3];
a=3*st1-3*st2+st3;
b=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1-2*(5-4*alpha)*st2+(4-3*alpha)*st3);
c=0.5*alpha^2/(1-alpha)^2*(st1-2*st2+st3);
yhat=a+b+c;
yhat=yhat';
err=sqrt(mean((y-yhat(1:end-1)).^2));%预测误差总和
%画图
t=x';
plot(t,y,'o',t,yhat(1:end-1)) %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')5画图线性拟合
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
outputfile ='../MATLAB/2.xlsx'; % 插值后数据存放
% 读入数据
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%%
cftool
%% 画图
plot(fittedmodel,x,y,'o');
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')6灰色关联度模型
CPI与什么因数有关
%%灰色关联度模型
clc,clear,close all
inputfile = '../MATLAB/hsgl.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
r=size(num);
y=num(2:r,:);%源数据
y1=mean(y');
y1=y1';
[row,col]=size(y);%获取矩阵的行数，列数
for i=1:row
for j=1:col
y2(i,j)=y(i,j)/y1(i);%初值像矩阵
end
end
for i=2:row
for j=1:col
y3(i-1,j)=abs(y2(i,j)-y2((i-1) ,j));%差序列
end
end
a=1;b=0;
for i=1:row-1
for j=1:col
if(y3(i,j)<=a) a=y3(i,j);%计算最大值
elseif(y3(i,j)>=b) b=y3(i,j);%计算最小值
end
end
end
for i=1:row-1
for j=1:col
y4(i,j)=(a+0.5*b)/(y3(i,j)+0.5*b);%计算关联系数
end
end
y5=sum(y4')/(col-1);%灰色关联度
y5=y5'
请发表评论