第三次上机作业.docx
《第三次上机作业.docx》由会员分享，可在线阅读，更多相关《第三次上机作业.docx（9页珍藏版）》请在冰豆网上搜索。
第三次上机作业
计
量
上
机
实
验
报
告
（3）
计量经济学上机实验报告
实验内容：
任务1-5
实验工具：
Matlab2008
实验目的：
学会使用if语句、循环、函数、积分和极大似然估计的基本方法
实验过程：
任务1:
求1到10000中能被3整除的数的和.
Matlab代码
s=0
fori=1:
1:
10000
p（i）=i;
q（i）=mod（i,3）;
ifq（i）==fix（q（i））
s=s+p（i）;
end;
end;
s
结果
s=
16668333
任务2：
用function（）命令编写正态分布的密度函数，并分别计算N（0,1），N（0,2），N（0,3）小于1部分的面积。
Matlab代码
>>y2=normcdf（x,0,1）;
>>plot（x,y1,'r-',x,y2,'g-'）;
>>plot（x,y1,'r-',x,y2,'g-'）;legend（'密度函数','分布函数'）;
J=int（normpdf（x,0,1）,x,-inf,1）;
symsx;
J1=int（normpdf（x,0,2）,x,-inf,1）;
VJ1=vpa（J1）;
>>symsx;
J2=int（normpdf（x,0,3）,x,-inf,1）;
VJ2=vpa（J2）;
>>VJ=vpa（J）
>>VJ1
>>VJ2
结果
VJ=
0.84134474606854303612337332991661
VJ1=
0.69146246127401317558126706214263
VJ2=
0.63055865981823639009615938887623
概率密度函数曲线图和概率分布曲线图
任务3：
生成10000个标准正态分布的随机数，分别计算大于1.96，小于-1.96和-1.96到1.96之间的数各有多少个。
Matlab代码
s2=0;
fori=1:
1:
10000
if（x3（i）>1.96）
s2=s2+1;
end;
end;
>>s2
结果
s2=
245
任务4：
生成y=500+125*x+epslon，epslon~N（0,10）的1000个数据点，并用极大似然估计进行回归。
Matlab代码
Command主函数
epslon=normrnd（0,10,1000,1）;
x=randn（1000,1）;
y=500+125*x+epslon;
[para,standard_deviation]=my_mle（'mynormpdfsum',[500;125;10],y,x）;
z=para
（1）+para
（2）*x;
plot（x,y,'o'）;
holdon;
plot（x,z,'--'）;
Definemynormpdfsum
functionf=mynormpdfsum（para,num,y,x）
yy=1/sqrt（2*pi）/para（3）*exp（-（y-para
（1）-para
（2）*x）.^2/2/para（3）^2）;
ifnum==1%（note:
itmustbesetto1）
f=log（yy）;
else
f=-sum（log（yy））;
end
Definemy_mle
function[para,standard_deviation,fv]=my_mle（fun,para0,varargin）
%estimateparametersandstandarderrorswhenusingmaximiumlikelihood
%estimation（MLE）
%input
%fun:
afunctiondefinedbyusersforcalculatinglogprobabilitydensity
%function（pdf）andnegativesumoflogarithmofpdf
%para0:
giveninitialparameters
%varargin:
otherneededparametersrequiredbyfun
%output
%para:
estimatedparameters
%standard_deviation:
standarddeviationsofestimatedparameters
%fv:
maximizedlikelihoodfunctionvalue
%%%%%%%%%%%
%
para0=para0（:
）;
[para,fv]=fminsearch（fun,para0,[],2,varargin{:
}）;
fv=-fv;
n=length（para）;
d=numericalfirstderivative（fun,para,1,varargin{:
}）;
standard_deviation=sqrt（diag（pinv（d'*d）））;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionf=numericalfirstderivative（fun,parameter,varargin）
%input:
%fun:
thenameofafunction
%parameter:
givenparameterwithrespecttowhichfirst-orderderivative
%iscalculated
%varargin:
otherneededparametersrequiredbyfun
%output:
%f:
numericalfirstorderderivativeoffunatparameter
n=length（parameter）;
fori=1:
n
a=zeros（n,1）;
a（i）=min（parameter（i）*1e-6,1e-5）;
y1（:
i）=feval（fun,parameter+a,varargin{:
}）;
y2（:
i）=feval（fun,parameter-a,varargin{:
}）;
f（:
i）=（y1（:
i）-y2（:
i））/2/a（i）;
end
结果
para=
499.4534
125.1685
10.1249
standard_deviation=
0.3204
0.3244
0.2104
Figure
分析
上述估计表明：
a、b的参数估计分别为499.4534,125.1685,epslon的均值和方差的估计值分别为0.2104、10.1249，这与给出的标准值之间有一定的差值。
任务5：
生成y=2+e^x+epslon，epslon~N（0,1）的1000个数据点，并用极大似然估计进行回归。
Matlab代码
Command主函数
epslon=normrnd（0,10,1000,1）;
x=randn（1000,1）;
y=2+exp（x）+epslon;
[para,standard_deviation]=my_mle（'mynormpdfsum',[2;exp
（1）;1],y,x）;
z=para
（1）+para
（2）.^x;
plot（x,y,'o'）;
holdon;
plot（x,z,'--'）;
title（'极大似然估计'）
Definemynormpdfsum
functionf=mynormpdfsum（para,num,y,x）
yy=1/sqrt（2*pi）/para（3）*exp（-（y-para
（1）-para
（2）.^x）.^2/2/para（3）^2）;
ifnum==1%（note:
itmustbesetto1）
f=log（yy）;
else
f=-sum（log（yy））;
end
Definemy_mle
function[para,standard_deviation,fv]=my_mle（fun,para0,varargin）
%estimateparametersandstandarderrorswhenusingmaximiumlikelihood
%estimation（MLE）
%input
%fun:
afunctiondefinedbyusersforcalculatinglogprobabilitydensity
%function（pdf）andnegativesumoflogarithmofpdf
%para0:
giveninitialparameters
%varargin:
otherneededparametersrequiredbyfun
%output
%para:
estimatedparameters
%standard_deviation:
standarddeviationsofestimatedparameters
%fv:
maximizedlikelihoodfunctionvalue
%%%%%%%%%%%
%
para0=para0（:
）;
[para,fv]=fminsearch（fun,para0,[],2,varargin{:
}）;
fv=-fv;
n=length（para）;
d=numericalfirstderivative（fun,para,1,varargin{:
}）;
standard_deviation=sqrt（diag（pinv（d'*d）））;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionf=numericalfirstderivative（fun,parameter,varargin）
%input:
%fun:
thenameofafunction
%parameter:
givenparameterwithrespecttowhichfirst-orderderivative
%iscalculated
%varargin:
otherneededparametersrequiredbyfun
%output:
%f:
numericalfirstorderderivativeoffunatparameter
n=length（parameter）;
fori=1:
n
a=zeros（n,1）;
a（i）=min（parameter（i）*1e-6,1e-5）;
y1（:
i）=feval（fun,parameter+a,varargin{:
}）;
y2（:
i）=feval（fun,parameter-a,varargin{:
}）;
f（:
i）=（y1（:
i）-y2（:
i））/2/a（i）;
end
Figure
结果
para=
1.9873
2.7273
1.0242
standard_deviation=
0.0340
0.0138
0.0230
分析
上述估计表明：
a、e的参数估计分别为1.9873，2.7273,epslon的均值和方差的估计值分别为0.0230、1.0242，这与给出的标准值之间有一定的差值。
升级会员 