非参数计量经济学附录5程序选.docx
《非参数计量经济学附录5程序选.docx》由会员分享，可在线阅读，更多相关《非参数计量经济学附录5程序选.docx（47页珍藏版）》请在冰豆网上搜索。
非参数计量经济学附录5程序选
附录5MATLAB和GAUSS程序精选
本附录为本书各章例子计算的MATLAB和GAUSS程序。
利用MATLAB和GAUSS软件进行计算必须编程，详细介绍MATLAB和GAUSS软件的编程语言不是本书的目的，读者想完全理解本附录的程序并举一反三，可仔细阅读附录3和附录4，并参考MATLAB和GAUSS软件的帮助文件。
精通GAUSS编程的捷径之一是学习林光平的著作《计算计量经济学－计量经济学家和金融分析师GAUSS编程与应用》（清华大学出版社，2003）。
1.MATLAB程序。
第2章例2.2中随机变量密度函数的核估计及其置信水平为95％的置信区间。
h为窗宽，f为密度函数的核估计，y=sortrows（[xf]）是对矩阵[xf]按照第一列由小到大排序。
x=[0;-0.04;0.18;-0.91;-0.95;0.16;-0.13;0.08;0.7;3.01;-1.15;0.21;-0.11;-0.31;-0.67;-1.22;-1.49;1.45;
-0.88;0.72;0.46;-2.03;-0.41;0.99;-0.4;-0.93;1.1;0.66;0.61;2.32;-0.48;-1.42;-0.89;0.93;0.72;0.6;
1.06;-0.19;-0.3;-0.38;-0.86;0.99;-1.04;-0.11;1.05;-1.37;0.28;-0.08;-0.38;0.73;-1.52;-0.3;0.96;
2.98;0.55;-0.03;0.96;-0.34;0.12;0.2;0.23;0.37;-0.1;1.35;-0.98;0.29;-1.6;0.76;1.17;0.49;-1.64;0.69;
-0.75;1.97;-0.61;0.23;-1.07;0.23;0.95;0.56;0.8;-0.62;-0.19;0.95;-0.6;-2.07;0.47;0.11;-1.3;1.25;
-0.77;0.96;-1.18;0.28;0.34;-0.48;-0.38;0.26;-0.79;1.82];
n=100;%thenumberoftheobservations
h=0.5;
f=zeros（n,1）;
fori=1:
n;
f（i）=mean（（1/sqrt（2*pi））*exp（-0.5*（（x-x（i））/h）.^2））/h;
end;
y=sortrows（[xf]）;
m=y（:
2）;
b=1/sqrt（4*pi）;
m0=（1.96/sqrt（n*h））*sqrt（m*b）;
plot（y（:
1）,m,'--',y（:
1）,m+m0,y（:
1）,m-m0）;
2.MATLAB程序。
第2章例2.3中选择3阶核函数的减少核估计偏的方法对随机变量密度函数的核估计及其置信水平为95％的置信区间。
h为窗宽，f为密度函数的减少偏方法的核估计，y=sortrows（[xf]）是对矩阵[xf]按照第一列由小到大排序。
x=[0;-0.04;0.18;-0.91;-0.95;0.16;-0.13;0.08;0.7;3.01;-1.15;0.21;-0.11;-0.31;-0.67;-1.22;-1.49;1.45;
-0.88;0.72;0.46;-2.03;-0.41;0.99;-0.4;-0.93;1.1;0.66;0.61;2.32;-0.48;-1.42;-0.89;0.93;0.72;0.6;
1.06;-0.19;-0.3;-0.38;-0.86;0.99;-1.04;-0.11;1.05;-1.37;0.28;-0.08;-0.38;0.73;-1.52;-0.3;0.96;
2.98;0.55;-0.03;0.96;-0.34;0.12;0.2;0.23;0.37;-0.1;1.35;-0.98;0.29;-1.6;0.76;1.17;0.49;-1.64;0.69;
-0.75;1.97;-0.61;0.23;-1.07;0.23;0.95;0.56;0.8;-0.62;-0.19;0.95;-0.6;-2.07;0.47;0.11;-1.3;1.25;
-0.77;0.96;-1.18;0.28;0.34;-0.48;-0.38;0.26;-0.79;1.82];
n=100;%thenumberoftheobservations
h=0.5;
f=zeros（n,1）;
fori=1:
n;
f（i）=mean（（0.5/sqrt（2*pi））*（3-（（x-x（i））/h）.^2）.*exp（-0.5*（（x-x（i））/h）.^2））/h;
end;
y=sortrows（[xf]）;
m=y（:
2）;
b=0.476;
m0=（1.96/sqrt（n*h））*sqrt（m*b）;
plot（y（:
1）,m,'--',y（:
1）,m+m0,y（:
1）,m-m0）;
3.MATLAB程序。
第2章例2.5的多元密度函数的核估计。
x=mvnrnd（mu,sigma,n）为生成均值为mu,协方差阵为sigma的n个二元正态随机向量序列，h为窗宽，f0为实际密度函数值，f为密度函数的核估计。
n=120;%thenumberoftheobservations
mu=[00];
sigma=[11.5;1.53];
x=mvnrnd（mu,sigma,n）;
h=0.4;
f0=zeros（n,1）;
f=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
w=zeros（n）;
forj=1:
n;
e=（（x（j,1）-x（i,1））/h）^2+（（x（j,2）-x（i,2））/h）^2;
w（j,j）=（2/pi）*（1-e）*（e<1）;
end;
f（i）=（p'*w*p）/（n*h^2）;%Itistheestimatoroff（x）
f0（i）=（1/sqrt（3*pi^2））*exp（-2*（x（i,1）^2-x（i,1）*x（i,2）+x（i,2）^2/3））;
end;
[x（:
1）x（:
2）f0f]
4.MATLAB程序。
计算第3章例3.11的局部线性估计及其置信水平为95%的置信区间，其中核函数选择
，窗宽h=0.08。
d为120×3的数据阵，mp为回归函数的局部线性估计，vp是mp的标准差的估计。
n=120;%thenumberoftheobservations
d=[];%inputdatainit
x=d（:
3）;
y=d（:
2）;
h=0.08;
mp=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
x0=[px-x（i）];
w=zeros（n）;
forj=1:
n;
e=（（x（j）-x（i））/h）^2;
w（j,j）=0.75*（1-e）*（e<1）;
end;
a=inv（x0'*w*x0）*（x0'*w*y）;
mp（i）=a
（1）;%Itistheestimatorofm（x）
end;
u=y-mp;
vp=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
x0=[px-x（i）];
w=zeros（n）;
v=zeros
（2）;
forj=1:
n;
e=（（x（j）-x（i））/h）^2;
w（j,j）=0.75*（1-e）*（e<1）;
v=v+w（j,j）^2*u（j）^2*x0（j,:
）'*x0（j,:
）;
end;
o=inv（x0'*w*x0）;
vv=o*v*o';%Itistheestimatorofm（x）
vp（i）=sqrt（vv
（1））;
end;
plot（x,mp,':
',x,mp-1.96*vp,x,mp+1.96*vp,x,y,'.'）
5.MATLAB程序。
计算第3章例3.15的变窗宽局部线性估计，变窗宽为解释变量密度函数的核估计，不变窗宽为0.05,核函数为
。
d为120×3的数据阵，fp为解释变量密度函数的核估计，mp为回归函数的局部线性估计，hv为变窗宽，vp是mp的标准差的估计。
n=120;%thenumberoftheobservations
d=[];%inputdatainit
x=d（:
3）;
y=d（:
2）;
h=0.6;
fp=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
w=zeros（n）;
forj=1:
n;
e=（（x（j）-x（i））/h）^2;
w（j,j）=0.75*（1-e）*（e<1）;
end;
fp（i）=（n*h）^（-1）*p'*w*p;
end;
%plot（x,fp,'.'）
mp=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
hv=0.054/fp（i）;
x0=[px-x（i）];
w=zeros（n）;
forj=1:
n;
e=（（x（j）-x（i））/hv）^2;
w（j,j）=0.75*（1-e）*（e<1）;
end;
a=inv（x0'*w*x0）*（x0'*w*y）;
mp（i）=a
（1）;%Itistheestimatorofm（x）
end;
u=y-mp;
vp=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
x0=[px-x（i）];
w=zeros（n）;
v=zeros
（2）;
forj=1:
n;
e=（（x（j）-x（i））/h）^2;
w（j,j）=0.75*（1-e）*（e<1）;
v=v+w（j,j）^2*u（j）^2*x0（j,:
）'*x0（j,:
）;
end;
o=inv（x0'*w*x0）;
vv=o*v*o';%Itistheestimatorofm（x）
vp（i）=sqrt（vv
（1））;
end;
q=sortrows（[xympmp-1.96*vpmp+1.96*vpfp]）;
plot（q（:
1）,q（:
2）,'.',q（:
1）,q（:
3））
6.MATLAB程序。
计算第3章例3.20的近邻估计，核函数为
。
d为120×3的数据阵，m为回归函数的近邻估计。
n=120;%thenumberoftheobservations
d=[];%inputdatainit
y=d（:
2）;
x=d（:
3）;
k=15;
m=zeros（n,1）;
fori=1:
n;
q=abs（x-x（i））;
d0=[qxy];
d=sortrows（d0）;
h=d（k,1）;
f=0;g=0;
forj=1:
k;
e=（d（j,1）/h）^2;
f=f+0.75*（1-e）;
g=g+0.75*（1-e）*d（j,3）;
end;
m（i）=g/f;
end;
[xm]
plot（x,m,'.'）
7.MATLAB程序。
计算第3章例3.29的Legendre多项式正交序列估计。
d为120×3的数据阵，p为m个正交序列数据阵，mp为回归函数的正交序列估计。
n=120;
d=[];%inputdatainit
x=d（:
3）;
y=d（:
2）;
a=0;
b=1;
z=（2/（b-a））*x-（a+b）/（b-a）;
m=10;
p=zeros（n,m）;
p（:
1）=（1/sqrt
（2））*ones（n,1）;
p（:
2）=（1/sqrt（2/3））*z;
fori=3:
m;
p（:
i）=（1/i）*（2*i-1）*z.*p（:
i-1）-（1/i）*（i-1）*p（:
i-2）;
end;
b=inv（p'*p）*p'*y;
mp=p*b;
[xmp]
plot（x,mp,'.'）;
8.MATLAB程序。
计算第3章例3.38的Fourier正交序列估计。
d为120×3的数据阵，p为m个正交序列数据阵，mp为回归函数的正交序列估计。
n=120;
d=[];%inputdatainit
x=d（:
3）;
y=d（:
2）;
a=0;
b=1;
z=（1/（b-a））*x-a/（b-a）;
m=17;
q=zeros（n,m）;
q（:
1）=ones（n,1）;
fori=1:
（m-1）/2;
q（:
2*i）=sqrt
（2）*cos（2*pi*i*x）;
q（:
2*i+1）=sqrt
（2）*sin（2*pi*i*x）;
end;
b=inv（q'*q）*q'*y;
mp=q*b;
[xmp]
plot（x,mp,'.'）;
9.MATLAB程序。
计算第3章例3.40的多项式样条估计。
d为120×3的数据阵，t为节点，q为样条函数的数据阵，f为回归函数的样条估计，z为检验节点显著与否的统计量。
程序1：
节点的选择。
n=120;
d=[];%inputdatainit
x=d（:
3）;
y=d（:
2）;
k=6;
m=fix（n/k）-1;
x0=sortrows（x）;
t=zeros（m,1）;
fori=1:
m;
t（i）=x0（k*i）;
end;
p=zeros（n,m）;
fori=1:
m;
p（:
i）=（x-t（i））.^3.*（x-t（i）>0）;
end;
q=[pones（n,1）xx.^2x.^3];
a=inv（q'*q）;
b=a*q'*y;
f=q*b;
se=zeros（m,1）;
z=zeros（m,1）;
fori=1:
m;
se（i）=sqrt（（y-f）'*（y-f）/n）*sqrt（a（i,i））;
z（i）=b（i）/se（i）;
end;
[tz]
plot（x,f,'.'）;
程序2：
固定节点的多项式样条估计。
d为120×3的数据阵，t为节点，q为样条函数的数据阵，f为回归函数的样条估计，z为检验节点显著与否的统计量。
d=[];%inputdatainit
n=120;
x=d（:
3）;
y=d（:
2）;
t=[];%inputt’sinit
s=size（t）;
m=s
（1）;
p=zeros（n,m）;
fori=1:
m;
p（:
i）=（x-t（i））.^3.*（x-t（i）>0）;
end;
q=[pones（n,1）xx.^2x.^3];
a=inv（q'*q）;
b=a*q'*y;
f=q*b;
se=zeros（m,1）;
z=zeros（m,1）;
fori=1:
m;
se（i）=sqrt（（y-f）'*（y-f）/n）*sqrt（a（i,i））;
z（i）=b（i）/se（i）;
end;
[tz]
plot（x,f,'.'）;
10.MATLAB程序。
计算第4章例4.1的不变窗宽核估计。
d为121×3的数据阵，m0为回归函数值，m为回归函数的核估计。
n=121;
d=[];%inputdatainit
y=d（:
1）;
x=d（:
2:
3）;
m0=10*（x（:
1）+1）.*（x（:
2）+1）;
h=0.24;
m=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
w=zeros（n）;
forj=1:
n;
e=（（x（j,1）-x（i,1））/h）^2+（（x（j,2）-x（i,2））/h）^2;
w（j,j）=（2/pi）*（1-e）*（e<1）;
end;
m（i）=（p'*w*y）/（p'*w*p）;%Itistheestimatorofm（x）
end;
[yx（:
1）x（:
2）m0m]
sqrt（mean（（m0-m）.*（m0-m）））
11.MATLAB程序。
计算第4章例4.1的不变窗宽局部线性估计。
d为121×3的数据阵，m0为回归函数值，m为回归函数的局部线性估计。
n=121;
d=[];%inputdatainit
y=d（:
1）;
x=d（:
2:
3）;
m0=10*（x（:
1）+1）.*（x（:
2）+1）;
h=0.6;
m=zeros（n,1）;
p=ones（n,1）;
fori=1:
n;
x0=[px（:
1）-x（i,1）x（:
2）-x（i,2）];
w=zeros（n）;
forj=1:
n;
e=（（x（j,1）-x（i,1））/h）^2+（（x（j,2）-x（i,2））/h）^2;
w（j,j）=（2/pi）*（1-e）*（e<1）;
end;
a=inv（x0'*w*x0）*（x0'*w*y）;
m（i）=a
（1）;%Itistheestimatorofm（x）
end;
[yx（:
1）x（:
2）m0m]
sqrt（mean（（m0-m）.*（m0-m）））
12.MATLAB程序。
计算第4章例4.1的近邻估计。
采用一元核的k近邻估计，选择三角形核函数
。
d为121×3的数据阵，m0为回归函数值，m为回归函数的近邻估计。
n=121;
d=[];%inputdatainit
y=d（:
1）;
x=d（:
2:
3）;
m0=10*（x（:
1）+1）.*（x（:
2）+1）;
k=10;
m=zeros（n,1）;
p=ones（k,1）;
fori=1:
n;
q=（x（:
1）-x（i,1））.^2+（x（:
2）-x（i,2））.^2;
d0=[qx（:
1）x（:
2）y];
d=sortrows（d0）;
w=zeros（k）;
forj=1:
k;
e=d（j,1）/d（k,1）;
w（j,j）=（1-e）*（e<1）;
end;
m（i）=（p'*w*d（1:
k,4））/（p'*w*p）;%Itistheestimatorofm（x）
end;
[yx（:
1）x（:
2）m0m]
sqrt（mean（（m0-m）.*（m0-m）））
13.MATLAB程序。
计算第4章例4.1的正交序列估计。
d为121×3的数据阵，m0为回归函数值，m为回归函数的正交序列估计。
n=121;
d=[];%inputdatainit
y=d（:
1）;
x=d（:
2:
3）;
m0=10*（x（:
1）+1）.*（x（:
2）+1）;
x1=x（:
1）;
x2=x（:
2）;
k=8;
p1=zeros（n,k）;
p1（:
1）=（1/sqrt
（2））*ones（n,1）;
p1（:
2）=（1/sqrt（2/3））*x1;
fori=3:
k;
p1（:
i）=（1/i）*（2*i-1）*x1.*p1（:
i-1）-（1/i）*（i-1）*p1（:
i-2）;
end;
p2=zeros（n,k）;
p2（:
1）=（1/sqrt
（2））*ones（n,1）;
p2（:
2）=（1/sqrt（2/3））*x2;
fori=3:
k;
p2（:
i）=（1/i）*（2*i-1）*x2.*p2（:
i-1）-（1/i）*（i-1）*p2（:
i-2）;
end;
p=zeros（n,k）;
fori=1:
k;
p（:
i）=p1（:
i）.*p2（:
i）;
end;
a=inv（p'*p）*（p'*y）;
m=zeros（n,1）;
fori=1:
k;
m=m+a（i）*p（:
i）;
end;
[yx（:
1）x（:
2）m0m]
sqrt（mean（（m0-m）.*（m0-m）））
14.MATLAB程序。
计算第4章例4.2的Legendre多项式可加正交序列估计。
d为121×3的数据阵，m0为回归函数值，m为回归函数的正交序列估计。
n=121;
d=[];%inputdatainit
y=d（:
1）;
x=d（:
2:
3）;
m0=10*（x（:
1）+1）.^2+5*（x（:
2）-1）.^2;
x1=x（:
1）;
x2=x（:
2）;
k1=3;
k2=3;
p1=zeros（n,k1）;
p1（:
1）=（1/sqrt
（2））*ones（n,1）;
p1（:
2）=（1/sqrt（2/3））*x1;
fori=3:
k1;
p1（:
i）=（1/i）*（2*i-1）*x1.*p1（:
i-1）-（1/i）*（i-1）*p1（:
i-2）;
end;
p2=zeros（n,k2）;
p2（:
1）=（1/sqrt
（2））*ones（n,1）;
p2（:
2）=（1/sqrt（2/3））*x2;
fori=3:
k2;
p2（:
i）=（1/i）*（2*i-1）*x2.*p2（:
i-1）-（1/i）*（i-1）*p2（:
i-2）;
end;
p=[p1p2（:
2:
k2）];
a=inv（p'*p）*（p'*y）;
m=zeros（n,1）;
fori=1:
k1+k2-1;
m=m+a（i）*p（:
i）;
end;
[yx（:
1）x（:
2）m0m]
sqrt（mean（（m0-m）.*（m0-m）））
15.GAUSS程序。
计算第4章例4.10多元非参数计量经济模型的变窗宽局部线性估计。
48×4的数据矩阵由d盘下文本文件zhg1.txt导入，h为变窗宽，a为消费函数的变窗宽局部线性估计。
n=48;
m=4;
loaddata[n,m]=d:
\zhg1.txt;
gdp=data[.,1];
mc=data[.,3];
x1=gdp~lagn（mc,1）~lagn（mc,2）;
x=x1[3:
n,.];
y=mc[3:
n,.];
n=n-2;
a=zeros（n,4）;
u=zeros（n,1）;
i=1;
dountili>n;
p=（x[i,1]*ones（n,1））~（x[i,2]*ones（n,1））~（x[i,3]*ones（n,1））;
x0=ones（n,1）~（x-p）;
w=zeros（n,n）;
j=1;
dountilj>n;
h=-500+40000*j/n;
e=（（x[j,1]-x[i,1]）/h）^2+（（x[j,2]-x[i,2]）/h）^2+（（x[j,3]-x[i,3]）/h）^2;
w[j,j]=1*（1-e）*（e<1）;
j=j+1;
endo;
q=x0'*w*x0;
b=inv（q）;
f=x0'*w*y;
a[i,.]=（b*f）';
u[i]=y[i]-a[i,1];
i=i+1;
endo;
emc=（y-a[.,1]）^2;
printsqrt（meanc（emc））;
printa;
end;
16.GAUSS程序。
计算第
升级会员 