第matlab计量经济学多重共线性的诊断与处理(共15页)

《第matlab计量经济学多重共线性的诊断与处理(共15页)》由会员分享，可在线阅读，更多相关《第matlab计量经济学多重共线性的诊断与处理(共15页)（15页珍藏版）》请在装配图网上搜索。

1、第五节 多重共线性的诊断与处理5.1 多重共线性的诊断数据来源：计量经济学于俊年 编著 对外经济贸易大学出版社 2000.6 p208-p209 某国19981998的经济数据年份进口额(y)国内产值(x1t)存货额(x2t)国内消费(x3t)198815.9149.34.2108.1198916.4161.24.1114.8199019171.53.1123.2199119.1175.53.1126.9199218.8180.81.1132.1199320.4190.72.2137.7199422.7202.12.1146199526.5212.15.6154.1199628.1226.15

2、162.3199727.6231.95.1164.3199826.32390.7167.65.1.1 条件数与病态指数诊断设x1,x2,xp是自变量X1，X2，XP，经过中心化和标准化得到的向量，即： 记（x1,x2,xp）为x,设为xTx一个特征值，为对应的特征向量，其长度为1，若，则：根据上表，计算如下： x=149.3, 4.2, 108.1; 161.2, 4.1, 114.8; 171.5, 3.1,123.2; 175.5, 3.1, 126.9; 180.8, 1.1, 132.1; 190.7, 2.2, 137.7; 202.1, 2.1, 146; 212.1, 5.6,

3、154.1; 226.1,5, 162.3; 231.9, 5.1, 164.3; 239, 0.7, 167.6求x的相关矩阵RR=corrcoef(x)R = 1.000 0.573 0.079 0.573 1.000 0.007 0.079 0.007 1.000求R的条件数：cond(R)ans = 7.9832e+002也可先求R的特征值e=eig(R)e = 0.125 0.342 1.533注：e(3)/e(1)ans = 7.9491e+002条件数为717.804，大于100，存在较严重的多重共线性。为了进一步了解哪些变量之间存在线性关系，计算相关矩阵的特征值和相应的特征向量

4、：v,d=eig(R)v = 0.575 0.633 0.371 0.633 -0.563 0.058 -0.049 0.777 0.581d = 0.125 0 0 0 0.342 0 0 0 1.533注意：Rv=vd v为标准正交矩阵最小的特征值为0.125，对应的向量为：（0.575，0.633，-0.049）T考虑到第二个数0.633约等于0，从而即：所以存在使得： 5.1.2 方差膨胀因子诊断每一个自变量对应的方差膨胀因子为R-1相应的对角元素rjj。若记xj关于其他p-1个自变量的复相关系数为Rj则有： 如果VIF10，则认为自变量间存在严重的多重共线性。在本例中：diag(in

5、v(R)ans = 1.0e+002 * 1.643 0.590 1.056VIF=max(diag(inv(R)VIF =1.0555e+002VIF远大于10，存在严重的多重共线性。注意：书上结果错了，我用SPSS算了，也是这个结果。方差膨胀因子也可按此计算：x1=x(:,1);x2=x(:,2);x3=x(:,3); b bint,r,rint,stats=regress(x1,ones(11,1) x2 x3);一定要常数项1/(1-stats(1)ans = 1.5788e+0025.1.3 容许度(Tolerance)诊断若记xj关于其他p-1个自变量的复相关系数为Rj则有：Tol

6、j=1-R2j它是方差膨胀化因子的倒数。越小自变量共线性越强。小于0.1高度共线在本例中：Tol=1./diag(inv(R)Tol = 0.649 0.803 0.734最小的值远小0.1，高度多重共线性。5.1.4 方差比例诊断（看Applied Econometric using Matlab的第84页）注意：Applied Econometric using Matlab的第84页，4.4式是错的，4.3，4.5，4.6式是对的。某国19981998的经济数据年份进口额(y)国内产值(x1t)存货额(x2t)国内消费(x3t)198815.9149.34.2108.1198916.41

7、61.24.1114.8199019171.53.1123.2199119.1175.53.1126.9199218.8180.81.1132.1199320.4190.72.2137.7199422.7202.12.1146199526.5212.15.6154.1199628.1226.15162.3199727.6231.95.1164.3199826.32390.7167.6x1=149.3, 4.2, 108.1; 161.2, 4.1, 114.8; 171.5, 3.1,123.2; 175.5, 3.1, 126.9; 180.8, 1.1, 132.1; 190.7, 2.2

8、, 137.7; 202.1, 2.1, 146; 212.1, 5.6, 154.1; 226.1,5, 162.3; 231.9, 5.1, 164.3; 239, 0.7, 167.6;x=ones(size(x1,1),1),x1;vnames=strvcat(constant,x1,x2,x3);fmt=%12.6f;bkw(x,vnames,fmt);Belsley, Kuh, Welsch Variance-decomposition K(x) constant x1 x2 x3 1 0. 0. 0. 0. 140 0. 0. 0. 0. 188 0. 0. 0. 0. 1978

9、 0. 0. 0. 0.K(x)=188时，有两个方差比例大于0.5，x1与x3可以存在共线性。K(X)30 或者方差比例0.5，则存在多重共线性。上表的算法：nobs nvar = size(x);u d v = svd(x,0);lamda = diag(d(1:nvar,1:nvar);lamda2 = lamda.*lamda;v = v.*v;phi = zeros(nvar,nvar);for i=1:nvar;phi(i,:) = v(i,:)./lamda2;end; pi = zeros(nvar,nvar);for i=1:nvar;phik = sum(phi(i,:);

10、pi(i,:) = phi(i,:)/phik;end;pians = 0.428 0.618 0.568 0.274 0.588 0.356 0.856 0.154 0.928 0.663 0.565 0.123 0.056 0.363 0.011 0.449K(x)的算法：u d v = svd(x,0);d1=diag(d);d1 = 1.0e+002 * 8.981 0.173 0.293 0.001kx=d1(1)/d1(1);d1(1)/d1(2);d1(1)/d1(3);d1(1)/d1(4)kx = 1.0e+003 * 0.000 0.263 0.043 1.0545.2 多

11、重共线性的处理可参见经济计量学李景华 编著 中国商业出版社 第四章5.2.1 岭回归(脊回归)年份进口额(y)国内产值(x1t)存货额(x2t)国内消费(x3t)198815.9149.34.2108.1198916.4161.24.1114.8199019171.53.1123.2199119.1175.53.1126.9199218.8180.81.1132.1199320.4190.72.2137.7199422.7202.12.1146199526.5212.15.6154.1199628.1226.15162.3199727.6231.95.1164.3199826.32390.71

12、67.6x=149.3, 4.2, 108.1; 161.2, 4.1, 114.8; 171.5, 3.1,123.2; 175.5, 3.1, 126.9; 180.8, 1.1, 132.1; 190.7, 2.2, 137.7; 202.1, 2.1, 146; 212.1, 5.6, 154.1; 226.1,5, 162.3; 231.9, 5.1, 164.3; 239, 0.7, 167.6;y=15.9; 16.4;19;19.1; 18.8; 20.4; 22.7; 26.5; 28.1; 27.6; 26.3;bb = zeros(4,101);kvec = 0:0.01

13、:1;count = 0;for k = 0:0.01:1b(:,count) = ridge(y,ones(11,1) x,k);endplot(kvec,b),xlabel(k),ylabel(b,FontName,Symbol)点击最上面一要线，删除得：如果不想在图中包括常数项，则可：bb = zeros(3,101);kvec = 0:0.01:1;count = 0;for k = 0:0.01:1count = count + 1;bb(:,count) = ridge(y, x,k);endplot(kvec,bb),xlabel(k),ylabel(b,FontName,Sym

14、bol)为了看清k在0到0.1之间回归系数的变化情况，则：bb = zeros(3,11);kvec = 0:0.01:0.1;count = 0;for k = 0:0.01:0.1count = count + 1;bb(:,count) = ridge(y, x,k);endplot(kvec,bb),xlabel(k),ylabel(b,FontName,Symbol)因此，在k=0.04，各回归系数基本稳定。不用ridge，也可做岭回归图：x=149.3, 4.2, 108.1; 161.2, 4.1, 114.8; 171.5, 3.1,123.2; 175.5, 3.1, 126

15、.9; 180.8, 1.1, 132.1; 190.7, 2.2, 137.7; 202.1, 2.1, 146; 212.1, 5.6, 154.1; 226.1,5, 162.3; 231.9, 5.1, 164.3; 239, 0.7, 167.6;y=15.9; 16.4;19;19.1; 18.8; 20.4; 22.7; 26.5; 28.1; 27.6; 26.3;xb=zscore(x)/sqrt(10);x1=x(:,1);x2=x(:,2);x3=x(:,3);bb = zeros(3,11);kvec = 0:0.01:0.1;count = 0;for k = 0:0

16、.01:0.1count = count + 1;bb(:,count) =inv(diag(norm(x1-mean(x1) norm(x2-mean(x2) norm(x3-mean(x3)*inv(xb*xb+eye(3)*k)*xb*yendplot(kvec,bb),xlabel(k),ylabel(b,FontName,Symbol)上图的y轴是经处理后最后模型的回归系数。也可绘制k在0到1变化时，最后模型的回归系数变化情况。bb = zeros(3,101);kvec = 0:0.01:1;count = 0;for k = 0:0.01:1count = count + 1;b

17、b(:,count) =inv(diag(norm(x1-mean(x1) norm(x2-mean(x2) norm(x3-mean(x3)*inv(xb*xb+eye(3)*k)*xb*yendplot(kvec,bb),xlabel(k),ylabel(b,FontName,Symbol)xb=zscore(x)/sqrt(10); 标准化，即：xb = -0.325 0.066 -0.679 -0.378 0.947 -0.315 -0.029 -0.237 -0.769 -0.534 -0.237 -0.315 -0.053 -0.605 -0.405 -0.853 -0.303 -

18、0.041 0.407 -0.421 0.141 0.145 0.724 0.596 0.377 0.013 0.414 0.544 0.132 0.687 0.698 -0.079 0.687x1=x(:,1);x2=x(:,2);x3=x(:,3);b=inv(diag(norm(x1-mean(x1) norm(x2-mean(x2) norm(x3-mean(x3)*inv(xb*xb+eye(3)*0.04)*xb*yb = 0.129 0.860 0.022b0=mean(y)-b(1)*mean(x1)-b(2)*mean(x2)-b(3)*mean(x3)b0 = -8.174

19、最后的岭回归方程： y=-8.56956+0.06334x1+0.5874x2+0.11592x3残差平方和：sse=(norm(y-(b0+b(1)*x1+b(2)*x2+b(3)*x3)2sse = 2.254可决系数：1-sse/norm(y-mean(y)2ans = 0.984OLS的残差平方和：bb,bint,r,rint=regress(y,ones(11,1) x);norm(r)2ans = 1.149增加了45.25%计量经济学于俊年 P210 表中的VIF值算错了。k=0.04时VIF= diag(inv(xb*xb+0.04*eye(3)VIF = 11.9276 0.

20、9637 11.9350k=0.1VIF= diag(inv(xb*xb+0.1*eye(3)VIF = 5.1014 0.91035.1043因此，我们取k=0.1b=inv(diag(norm(x1-mean(x1) norm(x2-mean(x2) norm(x3-mean(x3)*inv(xb*xb+eye(3)*0.1)*xb*yb = 0.0660 0.55820.1061b0=mean(y)-b(1)*mean(x1)-b(2)*mean(x2)-b(3)*mean(x3)b0 = -7.6286最后的方程：y=-7.6286+0.0660x1+0.5582x2+0.1061x3

21、sse=(norm(y-(b0+b(1)*x1+b(2)*x2+b(3)*x3)2sse = 2.9054可决系数：1-sse/norm(y-mean(y)2ans = 0.9859(2.9054-1.149)/1.149ans = 0.7383残差平方和比OLS增加了73.83%下面再求y=-7.6286+0.0660x1+0.5582x2+0.1061x3各回归系数的标准差与相应的T值。参见于俊年的计量书P213。bb=inv(xb*xb+0.1*eye(3)*xb*ybbb = 0.4357 0.20260.4820sse=norm(yb-bb(1)*xb(:,1)-bb(2)*xb(:

22、,1)-bb(3)*xb(:,1)2sse =0.09331-sse/norm(yb)2ans =0.9067估计的方差：sse/(11-3)ans =0.0117回归系数的标准差：sb=sqrt(VIF*sse/(11-3)sb = 0.2439 0.10300.2439也可按：sb=diag(sqrt(inv(xb*xb+0.1*eye(3)*sse/(11-3)sb = 0.2439 0.10300.2439最后方程的回归系数的标准差：std(y)*sb./(std(x)ans = 0.0370 0.28380.0537P1=(1-tcdf(0.0660/0.0370,7)*2P1=0.

23、1176P2=(1-tcdf(0.5582/0.2838,7)*2P2=0.0899P3=(1-tcdf(0.1061/0.0537,7)*2P3=0.0887因此，y=-7.6286+0.0660x1+0.5582x2+0.1061x3标准差 (0.0370) (0.2838) (0.0537) P值 (0.1176) (0.0899) (0.0887)在显著性水平0.12下，各回归系数均通过了检验。5.2.2 主成分回归原理参见：经济计量学 李景华 P117P126x=149.3, 4.2, 108.1; 161.2, 4.1, 114.8; 171.5, 3.1,123.2; 175.5

24、, 3.1, 126.9; 180.8, 1.1, 132.1; 190.7, 2.2, 137.7; 202.1, 2.1, 146; 212.1, 5.6, 154.1; 226.1,5, 162.3; 231.9, 5.1, 164.3; 239, 0.7, 167.6;y=15.9; 16.4;19;19.1; 18.8; 20.4; 22.7; 26.5; 28.1; 27.6; 26.3;x1=x(:,1);x2=x(:,2);x3=x(:,3);xb=zscore(x)/sqrt(10);dinv=inv(diag(norm(x1-mean(x1) norm(x2-mean(x2

25、) norm(x3-mean(x3)dinv = 0.0105 0 0 0 0.1917 0 0 0 0.0153 Z=xb*AZ = 0.0067 -0.2013 -0.6725 0.0227 -0.1752 -0.5121 0.0069 0.0234 -0.3525 -0.0033 0.0263 -0.2827 -0.0232 0.4134 -0.2032 -0.0084 0.2085 -0.0598 -0.0135 0.2351 0.1142 -0.0214 -0.4286 0.3049 -0.0068 -0.3053 0.4931 0.0149 -0.3215 0.55880.0254

26、 0.5252 0.6116A,d=eig(corrcoef(x)A = 0.7070 0.0357 0.7063 0.0080 -0.9991 0.0425 -0.7072 0.0245 0.7066d = 0.0028 0 0 0 0.9983 0 0 0 1.9990A1=A(:,2,3)A1 = 0.0357 0.7063 -0.9991 0.04250.0245 0.7066Z=xb*AZ = 0.0067 -0.2013 -0.6725 0.0227 -0.1752 -0.5121 0.0069 0.0234 -0.3525 -0.0033 0.0263 -0.2827 -0.02

27、32 0.4134 -0.2032 -0.0084 0.2085 -0.0598 -0.0135 0.2351 0.1142 -0.0214 -0.4286 0.3049 -0.0068 -0.3053 0.4931 0.0149 -0.3215 0.5588 0.0254 0.5252 0.6116Z1=Z(:,2,3) 因为d的第一对角元接近于0，所以取Z的二三列Z1 = -0.2013 -0.6725 -0.1752 -0.5121 0.0234 -0.3525 0.0263 -0.2827 0.4134 -0.2032 0.2085 -0.0598 0.2351 0.1142 -0.4

28、286 0.3049 -0.3053 0.4931 -0.3215 0.55880.5252 0.6116b=dinv*A1*inv(Z1*Z1)*Z1*yb = 0.0728 0.61110.1063b0=mean(y)-mean(x)*bb0 = -9.1416主成份回归的模型为：y=-9.1416+0.0728x1+0.6111x2+0.1063x3相应的残差平方和为：sse=(norm(y-(b0+b(1)*x1+b(2)*x2+b(3)*x3)2sse =2.4372可决系数为：1-sse/norm(y-mean(y)2ans = 0.9882bbb,bint,rr,rint=reg

29、ress(y,ones(11,1) x);(norm(rr)2ans =1.6714(sse-1.6714)/1.6714ans =0.4582比OLS残差平方和增加了45.82%求自变量和因变量都标准化模型的回归系数。(特指处理了主成份后的)yb=zscore(y)/sqrt(10);xbb=A1*inv(Z1*Z1)*Z1*ybxbb = 0.4804 0.22180.4827与于俊年的书P225 中8.7.45相同相应的误差项方差的估计为：fc=norm(yb-xb(:,1)*xbb(1)-xb(:,2)*xbb(2)-xb(:,3)*xbb(3)2/(11-2)fc =0.0013注意

30、：不是11-3，因为已去掉一个特征值接近0的向量。xbb的标准差：下面讨论主成份回归的各回归系数的标准差。参见于俊年的计量书，P214P226a,aint,r,rint=regress(yb,z1 z2);a = -0.19260.6898相应的误差项方差的估计为：f=norm(r)2/(11-2)ans =0.0013相应的回归系数方差：va=diag(inv(z1 z2*z1 z2)*norm(r)2/(11-2)va = 0.00130.0007d1=inv(d)d1 = 359.0882 0 0 0 1.0018 0 0 0 0.5003 d2=zeros(3,1) d1(:,2 3)

31、d2 = 0 0 0 0 1.0018 0 0 0 0.5003cova=d2*fcova = 0 0 0 0 0.0013 0 0 0 0.0007C=zeros(3,1) A(:,2 3) C = 0 0.0357 0.7063 0 -0.9991 0.0425 0 0.0245 0.7066a1=0;aa1 = 0 -0.532 0.010xbb=C*a1xbb = 0.150 0.486 0.749vb=diag(C*cova*C)vb = 0.554 0.946 0.288vb为因变量和自变量均标准化后回归系数的方差，即xbb的方差。最后模型回归系数的标准差为：vbb=std(y)*sqrt(vb)./(std(x)vbb = 0.224 0.438 0.574p=2*(1-tcdf(b./vbb,(11-4)p = 1.0e-003 * 0.499 0.351 0.913最后的模型为：y=-9.1416+0.0728x1+0.6111x2+0.1063x3标准差： (0.00275) (0.0998) (0.00399)P值： (0.0000) (0.00048) (0.0000)可决系数0.9882