第matlab计量经济学多重共线性的诊断与处理.docx
《第matlab计量经济学多重共线性的诊断与处理.docx》由会员分享，可在线阅读，更多相关《第matlab计量经济学多重共线性的诊断与处理.docx（24页珍藏版）》请在冰点文库上搜索。
第matlab计量经济学多重共线性的诊断与处理
第五节多重共线性的诊断与处理
5.1多重共线性的诊断
数据来源：
《计量经济学》于俊年编著对外经济贸易大学出版社2000.6p208-p209
某国1998－1998的经济数据
年份
进口额（y）
国内产值（x1t）
存货额（x2t）
国内消费（x3t）
1988
15.9
149.3
4.2
108.1
1989
16.4
161.2
4.1
114.8
1990
19
171.5
3.1
123.2
1991
19.1
175.5
3.1
126.9
1992
18.8
180.8
1.1
132.1
1993
20.4
190.7
2.2
137.7
1994
22.7
202.1
2.1
146
1995
26.5
212.1
5.6
154.1
1996
28.1
226.1
5
162.3
1997
27.6
231.9
5.1
164.3
1998
26.3
239
0.7
167.6
5.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,154.1;226.1,5,162.3;231.9,5.1,164.3;239,0.7,167.6]
求x的相关矩阵R
R=corrcoef（x）
R=
1.000000000000000.024470490835730.99715218582079
.024*********
0.997152185820790.035673222920071.00000000000000
求R的条件数：
cond（R）
ans=
7.178039564809832e+002
也可先求R的特征值
e=eig（R）
e=
0.00278483106125
0.99825241504342
1.99896275389533
注：
e（3）/e
（1）
ans=
7.178039564809491e+002
条件数为717.804，大于100，存在较严重的多重共线性。
为了进一步了解哪些变量之间存在线性关系，计算相关矩阵的特征值和相应的特征向量：
[v,d]=eig（R）
v=
0.706964538965750.035698735796330.70634746471371
0.00795062868633-0.999063342195630.04253499482058
-0.707204304390490.024454826587770.70658618250581
d=
0.0027848310612500
00.998252415043420
001.99896275389533
注意：
Rv=vdv为标准正交矩阵
最小的特征值为0.00278483106125，对应的向量为：
（0.70696453896575，0.00795062868633，-0.70720430439049）T
考虑到第二个数0.00795062868633约等于0，从而
即：
所以存在
使得：
5.1.2方差膨胀因子诊断
每一个自变量对应的方差膨胀因子为R-1相应的对角元素rjj。
若记xj关于其他p-1个自变量的复相关系数为Rj则有：
如果VIF<5，则认为自变量间不存在多重共线性。
如果
如果VIF>10，则认为自变量间存在严重的多重共线性。
在本例中：
diag（inv（R））
ans=
1.0e+002*
1.79722747043643
.010*********
1.79843993838056
VIF=max（diag（inv（R）））
VIF=
1.798439938380555e+002
VIF远大于10，存在严重的多重共线性。
注意：
书上结果错了，我用SPSS算了，也是这个结果。
方差膨胀因子也可按此计算：
x1=x（:
1）;x2=x（:
2）;x3=x（:
3）;[bbint,r,rint,stats]=regress（x1,[ones（11,1）x2x3]）;一定要常数项
1/（1-stats
（1））
ans=
1.797227470435788e+002
5.1.3容许度（Tolerance）诊断
若记xj关于其他p-1个自变量的复相关系数为Rj则有：
Tolj=1-R2j
它是方差膨胀化因子的倒数。
越小自变量共线性越强。
小于0.1高度共线
在本例中：
Tol=1./diag（inv（R））
Tol=
0.00556412594649
0.97705973887803
0.00556037473734
最小的值远小0.1，高度多重共线性。
5.1.4方差比例诊断（看AppliedEconometricusingMatlab的第84页）
注意：
AppliedEconometricusingMatlab的第84页，4.4式是错的，4.3，4.5，4.6式是对的。
某国1998－1998的经济数据
年份
进口额（y）
国内产值（x1t）
存货额（x2t）
国内消费（x3t）
1988
15.9
149.3
4.2
108.1
1989
16.4
161.2
4.1
114.8
1990
19
171.5
3.1
123.2
1991
19.1
175.5
3.1
126.9
1992
18.8
180.8
1.1
132.1
1993
20.4
190.7
2.2
137.7
1994
22.7
202.1
2.1
146
1995
26.5
212.1
5.6
154.1
1996
28.1
226.1
5
162.3
1997
27.6
231.9
5.1
164.3
1998
26.3
239
0.7
167.6
x1=[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];
x=[ones（size（x1,1）,1）,x1];
vnames=strvcat（'constant','x1','x2','x3'）;
fmt='%12.6f';
bkw（x,vnames,fmt）;
Belsley,Kuh,WelschVariance-decomposition
K（x）constantx1x2x3
10.0000000.0000510.0000000.000012
1400.0000060.1402840.5981360.116948
1880.0000110.6802080.3756800.646263
19780.9999830.1794570.0261840.236777
K（x）=188时，有两个方差比例大于0.5，x1与x3可以存在共线性。
K（X）>30或者方差比例>0.5，则存在多重共线性。
上表的算法：
[nobsnvar]=size（x）;
[udv]=svd（x,0）;
lamda=diag（d（1:
nvar,1:
nvar））;
lamda2=lamda.*lamda;
v=v.*v;
phi=zeros（nvar,nvar）;
fori=1:
nvar;
phi（i,:
）=v（i,:
）./lamda2';
end;
pi=zeros（nvar,nvar）;
fori=1:
nvar;
phik=sum（phi（i,:
））;
pi（i,:
）=phi（i,:
）/phik;
end;
pi'
ans=
0.000000000004280.000051213866180.000000007535680.00001248205274
0.000006063715880.140283737583560.598136161578560.11694753752154
0.000010667849280.680208481806630.375680295735650.64626252084123
0.999983268430560.179456566743630.026183535150110.23677745958449
K（x）的算法：
[udv]=svd（x,0）;
d1=diag（d）;
d1=
1.0e+002*
.028*********
.0574********
0.04260231288293
0.00405906633001
kx=[d1
（1）/d1
（1）;d1
（1）/d1
（2）;d1
（1）/d1（3）;d1
（1）/d1（4）]
kx=
1.0e+003*
0.00100000000000
0.139********263
0.188********043
1.97788613117054
5.2多重共线性的处理
可参见《经济计量学》李景华编著中国商业出版社第四章
5.2.1岭回归（脊回归）
年份
进口额（y）
国内产值（x1t）
存货额（x2t）
国内消费（x3t）
1988
15.9
149.3
4.2
108.1
1989
16.4
161.2
4.1
114.8
1990
19
171.5
3.1
123.2
1991
19.1
175.5
3.1
126.9
1992
18.8
180.8
1.1
132.1
1993
20.4
190.7
2.2
137.7
1994
22.7
202.1
2.1
146
1995
26.5
212.1
5.6
154.1
1996
28.1
226.1
5
162.3
1997
27.6
231.9
5.1
164.3
1998
26.3
239
0.7
167.6
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,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:
1;
count=0;
fork=0:
0.01:
1
b（:
count）=ridge（y,[ones（11,1）x],k）;
end
plot（kvec',b'）,xlabel（'k'）,ylabel（'b','FontName','Symbol'）
点击最上面一要线，删除得：
如果不想在图中包括常数项，则可：
bb=zeros（3,101）;
kvec=0:
0.01:
1;
count=0;
fork=0:
0.01:
1
count=count+1;
bb（:
count）=ridge（y,x,k）;
end
plot（kvec',bb'）,xlabel（'k'）,ylabel（'b','FontName','Symbol'）
为了看清k在0到0.1之间回归系数的变化情况，则：
bb=zeros（3,11）;
kvec=0:
0.01:
0.1;
count=0;
fork=0:
0.01:
0.1
count=count+1;
bb（:
count）=ridge（y,x,k）;
end
plot（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.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;
fork=0:
0.01:
0.1
count=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'*y
end
plot（kvec',bb'）,xlabel（'k'）,ylabel（'b','FontName','Symbol'）
上图的y轴是经处理后最后模型的回归系数。
也可绘制k在0到1变化时，最后模型的回归系数变化情况。
bb=zeros（3,101）;
kvec=0:
0.01:
1;
count=0;
fork=0:
0.01:
1
count=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'*y
end
plot（kvec',bb'）,xlabel（'k'）,ylabel（'b','FontName','Symbol'）
xb=zscore（x）/sqrt（10）;标准化，即：
xb=
-0.477409661203250.17256712249066-0.48483579038679
-0.351896657283780.153********947-0.38215648650315
-0.24325935137029-0.03834824944237-0.25342422491769
-0.20107010635534-0.03834824944237-0.19672072874315
-0.14516935671053-0.42183074386605-0.11702932871405
-0.04075097529853-0.21091537193303-0.03120782099041
.0794********
0.184********1450.441004868587240.22012659448596
0.332623843083770.325960120260130.34579380222414
0.393798248355440.345134244981320.37644434069687
0.46868415825698-0.498527242750790.42701772917687
x1=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'*y
b=
0.06334061443129
0.58739760837860
0.11592051232022
b0=mean（y）-b
（1）*mean（x1）-b
（2）*mean（x2）-b（3）*mean（x3）
b0=
-8.56959415249174
最后的岭回归方程：
y=-8.56956+0.06334x1+0.5874x2+0.11592x3
残差平方和：
sse=（norm（y-（b0+b
（1）*x1+b
（2）*x2+b（3）*x3）））^2
sse=
2.42768928001254
可决系数：
1-sse/norm（y-mean（y））^2
ans=
0.98824073639984
OLS的残差平方和：
[bb,bint,r,rint]=regress（y,[ones（11,1）x]）;
norm（r）^2
ans=
1.67142209436149
增加了45.25%
《计量经济学》于俊年P210表中的VIF值算错了。
k=0.04时
VIF=diag（inv（xb'*xb+0.04*eye（3）））
VIF=
11.9276
0.9637
11.9350
k=0.1
VIF=diag（inv（xb'*xb+0.1*eye（3）））
VIF=
5.1014
0.9103
5.1043
因此，我们取k=0.1
b=inv（diag（[norm（x1-mean（x1））norm（x2-mean（x2））norm（x3-mean（x3））]））*inv（xb'*xb+eye（3）*0.1）*xb'*y
b=
0.0660
0.5582
0.1061
b0=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
sse=（norm（y-（b0+b
（1）*x1+b
（2）*x2+b（3）*x3）））^2
sse=
2.9054
可决系数：
1-sse/norm（y-mean（y））^2
ans=
0.9859
（2.9054-1.67142209436149）/1.67142209436149
ans=
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'*yb
bb=
0.4357
0.2026
0.4820
sse=norm（yb-bb
（1）*xb（:
1）-bb
（2）*xb（:
1）-bb（3）*xb（:
1））^2
sse=
0.0933
1-sse/norm（yb）^2
ans=
0.9067
估计的方差：
sse/（11-3）
ans=
0.0117
回归系数的标准差：
sb=sqrt（VIF*sse/（11-3））
sb=
0.2439
0.1030
0.2439
也可按：
sb=diag（sqrt（inv（xb'*xb+0.1*eye（3））*sse/（11-3）））
sb=
0.2439
0.1030
0.2439
最后方程的回归系数的标准差：
std（y）*sb./（std（x））'
ans=
0.0370
0.2838
0.0537
P1=（1-tcdf（0.0660/0.0370,7））*2
P1=0.1176
P2=（1-tcdf（0.5582/0.2838,7））*2
P2=0.0899
P3=（1-tcdf（0.1061/0.0537,7））*2
P3=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主成分回归
原理参见：
《经济计量学》李景华P117－P126
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,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））norm（x3-mean（x3））]））
dinv=
0.010500
00.19170
000.0153
Z=xb*A
Z=
0.0067-0.2013-0.6725
0.0227-0.1752-0.5121
0.00690.0234-0.3525
-0.00330.0263-0.2827
-0.02320.4134-0.2032
-0.00840.2085-0.0598
-0.01350.23510.1142
-0.0214-0.42860.3049
-0.0068-0.30530.4931
0.0149-0.32150.5588
升级会员 