数据预处理
1.计算每年每一个变量的基本统计量:均值,标准差,中位数;并用折线图给出原始数据和均值、中位数的变化趋势。
由所给的《中国企业商品价格指数数据99年至今》数据,利用MATLAB,很容易计算出结果,并画出折线图。
折线图如下:
2.给出各变量按季节变化趋势;每年各变量频数直方图的变化趋势,可做动态图。
利用MATLAB作图如下:
(此频数直方图是动态图)
3. 按原始数据,指出被解释变量和每解释变量之间的关系,画图;
利用MATLAB分别作图如下:
从图中可以看出,随着农产品价格指数的增长,总指数有线性增长的趋势,所以总指数与农产品价格指数有线性关系。
从图中可以看出,随着矿产品价格指数的增长,总指数有线性增长的趋势,所以总指数与矿产品价格指数有线性关系。
从图中可以看出,随着煤油电价格指数的增长,总指数有线性增长的趋势,所以总指数与煤油电价格指数有线性关系。
4.按年平均数、中位数预处理数据,指出被解释变量和每个解释变量之间的关系,画图;
利用MATLAB,分别作图如下:
从图中可以看出,对于每个解释变量(各个价格指数中均值),随着其增长,被解释变量(总指数均值)都有线性增长的趋势,所以被解释变量与每个解释变量分别有线性关系。
从图中可以看出,对于每个解释变量(各个价格指数均中位数),随着其增长,被解释变量(总指数中位数)都有线性增长的趋势,所以被解释变量与每个解释变量分别有线性关系。
(本节MATLAB代码见附录一)
一元回归分析
基本模型
下面先从简单入手,由上节数据预处理第三问求解得出,总指数y与农产品价格指数x之间存在线性关系,所以,可以建立以下一元线性回归模型:
y i= β0+β1xi+εi , i=1,2……,191
各εi独立同分布,其分布为N(0,σ2)
由数据(xi,yi)(i=1,2……,191)可获得β0、β1的估计 称
为y关于x的回归方程,,为回归系数,
εi是随机误差。
模型求解
利用MATLAB统计工具箱,得到模型的回归系数估计值及其置信区间(置信水平α=0.05)、检验统计量R2,F,p,s2。
回归系数 | 回归系数估计值 | 回归系数估计值 |
| 43.5577 | [39.1754 ,47.9401] |
| 0.5606 | [ 0.5184 ,0.6028] |
R2= 0.7839 F=685.6542 p<0.0001 s2 =4.4997 |
残差图如下:
剔除异常值后,再对模型进行求解,结果见下表:
回归系数 | 回归系数估计值 | 回归系数估计值 |
| 46.6171 | [43.0129,50.2213] |
| 0.5323 | [0.4976,0.5670] |
R2= 0.8391 F=917.7627 p<0.0001 s2 =2.9046 |
回归方程为:= 46.6171+0.5323x
结果分析
R2=0.8391是指因变量y的83.91%可由模型确定。F=917.7627远大于临界值,p远小于α,说明回归方程显著。
同时,回归系数置信区间不包含0,说明回归变量对因变量y的影响是显著的。
综上,该模型是可用的。
(本节代码见附录二)
多元回归分析
基本模型
由上节数据预处理第三问求解得出,总指数y与农产品价格指数x1、矿产品价格指数x2,煤油电价格指数x3之间均存在线性关系,故而,建立如下多元线性回归模型:
yi= β0+β1X1i+β2X2i+β3X3i+εi , i=1,2……,191
各εi独立同分布,其分布为N(0,σ2)
由数据(x1i,x2i,x3i,yi)(i=1,2……,191)可获得β0、β1β2、β3的估计 称
为y关于x1,x2,x3的回归方程,,为回归系数,εi是随机误差。
模型求解
利用MATLAB统计工具箱,得到模型的回归系数估计值及其置信区间(置信水平α=0.05)、检验统计量R2,F,p,s2。结果见下表:
回归系数 | 回归系数估计值 | 回归系数置信区间 |
| 32.5824 | [30.6793,34.4855] |
0.3742 | [0.3485,0.4001] | |
| 0.1372 | [0.1087,0.1656] |
| 0.1503 | [0.1312,0.1693] |
R2 = 1.0000 F=2008.4 p<0.0001 s2 = 0.6000 |
残差图如下:
剔除异常值后,再对模型进行求解,结果见下表:
回归系数 | 回归系数估计值 | 回归系数置信区间 |
| 32.6886 | [30.9173,34.4598] |
0.3771 | [0.3532,0.4010] | |
| 0.1313 | [0.1051,0.1575] |
| 0.1524 | [0.1347,0.1701] |
R2= 1.0000 F=2291.3 p<0.0001 s2 = 0.5000 |
回归方程:
= 32.6886+0.3771x1+0.1313x2+0.1524x3
结果分析
R2近似等于1,说明回归方程非常显著。F=2291.3远远大于临界值,p远远小于α,故而同样说明回归方程显著。
同时,回归系数置信区间均不包含0,说明各个回归变量对因变量y的影响是显著的。
综上,该模型是非常好的。
(本节代码见附录三)
时间序列模型
问题分析
题目所给的原始数据是以时间为序的,称时间序列。许多经济数据在时间上有一定的滞后性,实际上,在对时间序列做回归分析时,随机误差实际上,在对时间序列做回归分析时,随机误差εi有可能存在相关性,违背模型关于εi(对时间)相互独立的基本假设,同一变量的顺序观测值之间出现相关现象(称为自相关)是很自然的。然而一旦数据中存在这种自相关序列,如果仍采用普通的回归模型直接处理,将会出现不良结果,其预测也就失去了意义。为此,我们必须先来诊断数据是否存在自相关,如果存在,就要考虑自相关关系,建立新的回归模型。
残差ei = yi -可以作为εi估计值,画出ei– ei-1的散点图,能够从直观上判断εi的自相关性。由前面多元线性回归模型得到残差ei,作出ei– ei-1的散点图
如下:
可以看到:大部分点都落到了第一、第三象限中,表明εi存在正的自相关。
为了对εi的自相关性作定量诊断,并在诊断后得到新的结果,我们考虑如下的模型
yi= β0+β1X1i+β2X2i+β3X3i+εi ,εi = ρεi-1+ui
其中,ρ是相关系数,ui相互独立且服从均值为0的正态分布。
若ρ= 0,则退化为普通的回归模型;若ρ>0,则εi存在正的自相关;若ρ<0,则εi存在负的自相关。
下面做定量检验,即Durbin-Watson检验(简称D-W检验),计算DW统计量:
DW =
经过简单的运算可知,当n较大时,
DW
正是自相关系数ρ的估计值,于是 = 1 – DW/2。若在1附近,即DW在0或4附近,εi的自相关性很强。
基本模型
要根据DW的具体数值确定εi是否存在自相关,应该在给定的检验水平下,依照样本容量和回归变量数目,查D-W分布表,得到检验的临界值dL和dU,继而作出判断[1]。
计算可得DW= 0.3014<dL, =0.8493,α= 0.05,n=183,k=4,查表知,εi存在自相关。
作变换
然后,利用变换后的数据重新估计。
模型求解
利用MATLAB工具箱计算,得到结果如下表:
回归系数 | 回归系数估计值 | 回归系数置信区间 |
| 6.0977 | [5.6641,6.5313] |
0.2714 | [0.2438,0.2989] | |
| 0.1875 | [0.1610,0.2139] |
| 0.1271 | [0.1071,0.1471] |
R2 = 0.9124 F=618.2362 p<0.0001 s2 = 0.1071 |
回归方程:
把、、、还原为原始变量、、、得到结果为
结果分析
对上述模型在进行一次自相关检验,即诊断随机误差ui是否还存在自相关。计算出DW = 1.4430,α= 0.05,n=182,k=4 ,查表知,du<DW<4-du。即随机误差ui不存在自相关。
基于时间序列模型的回归模型与前面的多元线性回归预测值,如下图:
基于时间序列模型的回归模型与前面的多元线性回归残差图,如下图:
从两图中易知,基于时间序列,即加入自相关后的模型更贴近实际,预测结果更好,也就是说其更适合。
参考文献
[1]何晓群,刘文卿.应用回归分析.清华大学出版社,2001.
[2]姜启源,谢金星.数学模型(第四版).高等教育出版社.2012.
附录一
load data
year = 1999 : 2014;
for i = 1 : 4
for j = 1 : 15
jun_zhi(j,i) = mean(data((j-1)*12+1 : j*12 ,i+1));
biao_zhun_cha(j,i) = std(data((j-1)*12+1 : j*12 ,i+1));
zhong_wei_shu(j,i) = median(data((j-1)*12+1 : j*12 ,i+1));
end
end
for i = 1 : 4
jun_zhi(16,i) = mean(data(181:191,i+1));
biao_zhun_cha(16,i) = std(data(181:191,i+1));
zhong_wei_shu(16,i) = median(data(181:191,i+1));
end
%完成数据的初始化
%开始绘制原始数据折线图
plot(data(:,1),data(:,2:5),'-*')
legend('总产品','农产品','矿产品','煤油电')
title('原始数据折线图')
xlabel('时间')
ylabel('商品价格指数')
%完成原始数据折线图的绘制
%开始绘制均值和中位数的折线图
figure;
plot(year,jun_zhi,'-*')
legend('总产品','农产品','矿产品','煤油电')
title('商品价格指数均值折线图')
xlabel('时间')
ylabel('均值')
figure;
plot(year,zhong_wei_shu,'-*')
legend('总产品','农产品','矿产品','煤油电')
title('商品价格指数中位数折线图')
xlabel('时间')
ylabel('中位数')
%完成均值和中位数的折线图的绘制
%开始均值按季度变化绘图
ones(63,4);
for i=1:63
jun_zhi_season(i,:) = (data(i*3-2,2:5)...
+data(i*3-1,2:5)+data(i*3,2:5))./3;
end
x_season = 1:63;
figure;
plot(x_season,jun_zhi_season,'-*');
xlabel('季度(1-63代表顺序的63个季度)');
ylabel('季度均值');
legend('总产品','农产品','矿产品','煤油电');
title('各变量按季节变化趋势');
%完成均值按季度变化图形的绘制
% 动态频数直方图
data_dymic = sort(data(:,2:5));
temp = data_dymic(191,:);
data_dymic(3:191,:) = data_dymic(2:190,:);
data_dymic(2,:)= temp;
i=1;
figure;
while(i<=190)
subplot(2,2,1);
hist(data_dymic(1:i+1,1),5);
axis([90 112 0 70]);
title('总指数频数直方图');
subplot(2,2,2);
hist(data_dymic(1:i+1,2),5);
axis([90 125 0 70]);
title('农产品价格指数频数直方图');
subplot(2,2,3);
hist(data_dymic(1:i+1,3),5);
axis([80 125 0 70]);
title('矿产品价格指数频数直方图');
subplot(2,2,4);
hist(data_dymic(1:i+1,4),5);
axis([80 130 0 65]);
title('煤油电价格指数频数直方图');
pause(0.01);
i = i+1;
end
%结束动态图的绘制
%开始总指数其他各个指数之间关系的散点图并拟合的绘制
all = data(:,2);
agriculture = data(:,3);
mining = data(:,4);
coe = data(:,5);
figure;
A = polyfit(agriculture,all,1); %拟合
all2 = polyval(A,agriculture);
plot(agriculture,all,'b*',agriculture,all2,'r-')
title('总指数与农产品价格指数的关系图');
xlabel('农产品价格指数');
ylabel('总指数');
figure;
A = polyfit(mining,all,1);
all2 = polyval(A,mining);
plot(mining,all,'b*',mining,all2,'r-')
title('总指数与矿产品价格指数的关系图');
xlabel('矿产品价格指数');
ylabel('总指数');
figure;
A = polyfit(coe,all,1);
all2 = polyval(A,coe);
plot(coe,all,'b*',coe,all2,'r-')
title('总指数与煤油电价格指数的关系图');
xlabel('煤油电价格指数');
ylabel('总指数');
%结束总指数其他各个指数之间关系的散点图并拟合的绘制
%开始总指数与各指数均值、中位数关系的散点图并拟合
figure;
subplot(2,2,1);
B = polyfit(jun_zhi(:,2),jun_zhi(:,1),1);
all3 = polyval(B,jun_zhi(:,2));
plot(jun_zhi(:,2),jun_zhi(:,1),'b*',jun_zhi(:,2),all3,'r-')
title('总指数均值与农产品指数均值');
xlabel('农产品指数均值');
ylabel('总指数均值');
subplot(2,2,2);
B = polyfit(jun_zhi(:,3),jun_zhi(:,1),1);
all3 = polyval(B,jun_zhi(:,3));
plot(jun_zhi(:,3),jun_zhi(:,1),'b*',jun_zhi(:,3),all3,'r-')
title('总指数均值与矿产品指数均值');
xlabel('矿产品指数均值');
ylabel('总指数均值');
subplot(2,2,3);
B = polyfit(jun_zhi(:,4),jun_zhi(:,1),1);
all3 = polyval(B,jun_zhi(:,4));
plot(jun_zhi(:,4),jun_zhi(:,1),'b*',jun_zhi(:,4),all3,'r-')
title('总指数均值与煤油电指数均值');
xlabel('煤油电指数均值');
ylabel('总指数均值');
figure;
subplot(2,2,1);
B =polyfit(zhong_wei_shu(:,2),zhong_wei_shu(:,1),1);
all3 = polyval(B,zhong_wei_shu(:,2));
plot(zhong_wei_shu(:,2),zhong_wei_shu(:,1),'b*',...
zhong_wei_shu(:,2),all3,'r-')
title('总指数中位数与农产品指数中位数');
xlabel('农产品指数中位数');
ylabel('总指数中位数');
subplot(2,2,2);
B =polyfit(zhong_wei_shu(:,3),zhong_wei_shu(:,1),1);
all3 = polyval(B,zhong_wei_shu(:,3));
plot(zhong_wei_shu(:,3),zhong_wei_shu(:,1),'b*',...
zhong_wei_shu(:,3),all3,'r-')
title('总指数中位数与矿产品指数中位数');
xlabel('矿产品指数中位数');
ylabel('总指数中位数');
subplot(2,2,3);
B =polyfit(zhong_wei_shu(:,4),zhong_wei_shu(:,1),1);
all3 = polyval(B,zhong_wei_shu(:,4));
plot(zhong_wei_shu(:,4),zhong_wei_shu(:,1),'b*',...
zhong_wei_shu(:,4),all3,'r-')
title('总指数中位数与煤油电指数中位数');
xlabel('煤油电指数中位数');
ylabel('总指数中位数');
%第一题结束
附录二
load data
all = data(:,2);
agriculture = data(:,3);
X = [ ones(191,1),agriculture ];
[b,bint,r,rint,stats] = regress(all,X);
b
bint
stats
rcoplot(r,rint) %残差图
for i =191:-1:1 %剔除异常
if(abs(rint(i,1))+abs(rint(i,2))==abs(rint(i,1)+rint(i,2)))
all(i) = [];
agriculture(i) = [];
end
end
X = [ ones(178,1),agriculture ];
[b,bint,r,rint,stats] = regress(all,X);
b
bint
stats
附录三
load data
all = data(:,2);
agriculture = data(:,3);
mining = data(:,4);
coe = data(:,5);
X = [ ones(191,1),agriculture,mining,coe ];
[b,bint,r,rint,stats] = regress(all,X);
b
bint
stats
rcoplot(r,rint) %残差图
for i =191:-1:1 %剔除异常**必须从后往前剔除
if(abs(rint(i,1))+abs(rint(i,2))==abs(rint(i,1)+rint(i,2)))
all(i) = [];
agriculture(i) = [];
mining(i) = [];
coe(i) = [] ;
end
end
%
X = [ ones(183,1),agriculture,mining,coe ];
[b,bint,r,rint,stats] = regress(all,X);
b
bint
stats
附录四
et =all-(b(1)*ones(183,1)+b(2)*agriculture+b(3)*mining+b(4)*coe);
plot(et(1:182),et(2:183),'+');
hold on
x_time =linspace(0,0,10);
y_time = linspace(-2,2,10);
plot(x_time,y_time)
plot(y_time,x_time)
hold off
%
s1 = 0;
s2 = 0;
for i = 2:183
s1 = s1+et(i)*et(i-1);
s2 = s2+et(i)*et(i);
end
DW = 2*(1-s1/s2)
p_time = s1/s2
%
%数据变换
all_time= all(2:183)-p_time*all(1:182);
agriculture_time =agriculture(2:183)-p_time*agriculture(1:182);
mining_time=mining(2:183)-p_time*mining(1:182);
coe_time= coe(2:183)-p_time*coe(1:182);
%回归
X_time = [ones(182,1),agriculture_time,mining_time,coe_time ];
[b_time,bint_time,r_time,rint_time,stats_time]=...
regress(all_time,X_time);
b_time
bint_time
stats_time
%%%%%%%%%%%%%%%%%%%%
et2 = all_time-(b_time(1)*ones(182,1)+b_time(2)*...
agriculture_time+b_time(3)*mining_time+b_time(4)*coe_time);
s1 = 0;
s2 = 0;
for i = 2:182
s1 = s1+et2(i)*et2(i-1);
s2 = s2+et2(i)*et2(i);
end
DW = 2*(1-s1/s2)
figure;
%
yuce_multi = b(1)*ones(183,1)+b(2)*agriculture+b(3)*mining+b(4)*coe;
yuce_time(1) = 95.88;
yuce_time(2:183) =6.0977*ones(182,1)+0.8493*all(1:182)+0.2714...
*agriculture(2:183)-0.2305*agriculture(1:182)+0.1875*mining...
(2:183)-0.1592*mining(1:182)+0.1271*coe(2:183)-0.1079*coe(1:182);
count = 1:183;
plot(count,all,'o',count,yuce_time,'*',count,yuce_multi,'+')
legend('总指数','时间序列','多元线性回归');
figure;
et_time = all - yuce_time';
plot(count,et_time,'*',count,et,'+')
legend('时间序列','多元线性回归');
hold on
x_time =linspace(0,183,10);
y_time = linspace(0,0,10);
plot(x_time,y_time)
hold off