dxdsinx 发表于 2011-2-8 23:04 
各位好！小弟本次目標為利用高斯模擬大樣本理論中漸進分配之結果，即是模擬出(根號N)*(B的估計值 - B)在樣本數很大（未到無限大）時，分配性質會趨近常態分配，期望值為B，變異數西亇碼平方*(XtransposeX)逆矩陣的結果。小弟在爬版之後利用版上高手所提供的蒙地卡羅模擬方法寫出目前之程式碼，目前碰到的問題有二，希望版上各位高手可以幫小弟一起研究研究，謝謝各位！
（1）小弟跑出來的結果，最下圖，小弟不知道如何利用以下的程式碼得到的結果解釋估計值是符合常態分配的結論，似乎關鍵在於mean of estimated s.e. of OLS b ，但是小弟翻遍書籍還是不知道如何解釋Monte Carlo的結論
（2）因為Monte Carlo理論基本概念即是隨機抽取樣本的概念，所以我是用x = rndu(100,1)做為樣本，不知道這樣行不行？
/* Model y = xb + e , e~N(0,1)
x=1~x1,x1~{uniform random variables between 0 and 1}
E[x'x]=A=
E[1*1]=1 E[1*x1]=1/2
E[x1*1]=1/2 E[x1^2]=[E(x1)]^2+Var[x1]=1/4+1/12=1/3
n^(1/2)*(bhat-b)~N(0,A^(-1))
------------------------------------------*/
library pgraph;
graphset;
// rndseed 100;
sigma=10;
nAll=10000;
x=ones(nAll,1)~rndu(nAll,1);
y=x*(1|1)+sigma*rndn(nAll,1);
dataAll=x~y;
// nSample=10|50|100|1000;
nSample=3|10|100;
nCase=rows(nSample);
nRepeat=10000;
Mb=zeros(nRepeat,nCase);
Msb=zeros(nRepeat,nCase);
for j (1,nCase,1);
for i (1,nRepeat,1);
n=nSample[j];
dataIndex=ceil(rndu(n,1)*nAll);
data=dataAll[dataIndex,.];
x=data[.,1 2];
y=data[.,3];
b=invpd(x'x)*(x'y);
Mb[i,j]=b[2];
e=y-x*b;
s2e=e'e/(n-2);
s2b=s2e*invpd(x'x);
sb=sqrt(diag(s2b));
Msb[i,j]=sb[2];
endfor;
endfor;
//Consistency of OLS
graphset;
_ptek="fig1.tkf";
title("Consistency of OLS");
xlabel("Sample Size");
ylabel("bhat");
xtics(1,4,1,1);
xy(seqa(1,1,nCase),Mb');
// Asymptotic Normality of OLS
Mb0=1/sqrt(12)*(ones(nRepeat,1)*sqrt(nSample')).*(Mb-1)/sigma;
bMin=minc(minc(Mb0));
bMax=maxc(maxc(Mb0));
x0=seqa(-4,8/20,20)|bMax;
Mm=zeros(rows(x0),nCase);
Mfreq=zeros(rows(x0),nCase);
for i (1,nCase,1);
freq = counts(Mb0[.,i],x0);
Mfreq[.,i]=freq/nRepeat;
endfor;
temp=cdfn(seqa(-4,8/20,20)|4.1);
ytemp=temp[1]|(temp[2:rows(x0)]-temp[1:(rows(x0)-1)]);
graphset;
_ptek="fig2.tkf";
title("Asymptotic Normality of OLS");
_plegctl={1 5 3 0.1};
_plegstr=" true\000 n=3\000 n=10\000 n=100";
xy(seqa(-4,8/20,20)|4.1,ytemp~(Mfreq));
Consistency of OLS：
Asymptotic Normality of OLS：
扫码加好友，拉您进群
全部版块
我的主页
收藏
