matlab进行线性回归的蒙特卡洛模拟

事情是这样的：

       LZ目前在国外，计量课留了林文夫书上的一道数据模拟题。大样本情况下，beta的样本期望和方差不是都有公式可以求嘛，我们的工作就是模拟一大批数据，通过数据堆积，证明样本期望和方差在大样本情况下是能够逼近给定的原值的。
       题里要求我们先通过某种途径搞出一个满足AR(1)条件的x向量作为原数据，然后在给定beta1=1（截距）和beta2=0.5（斜率）的情况下，加入一个满足N(0,1)的epsilon，算出y向量。这个(y,x)的组合再反过来作为样本点，然后做回归，推断beta1和beta2的估计值b1和b2.
       上述从epsilon的生成开始，过程要循环10万次以上，得出在x不变情况下的10万组y，从而算出10万个b1和b2。
       为进行检验，加入计数器T=0和累积变量B=(0,0)。一方面把每次的b1和b2加到B里去，最后除以总次数，应该与beta1与beta2很接近；另一方面，把每次回归的t统计量与t(0.025,n-2)进行比较，如果该次的t统计量大于t(0.025,n-2)，那么就在T上加1,，最后应该在总次数的5%左右。

       那么问题来了。B的值通过模拟，与beta1和beta2确实很接近，但T的值，我始终维持在数据量的6%-8%。我检查了ols和t统计量公式，应该都没有错误。为以防万一，用矩阵形式重新写了一遍，也没有错误。
       某白此前几乎没有用过matlab，已折磨致死，求各位老师施以援手，指点问题所在，拯救某白于水火之中。

附循环过程的语句，没有计量包，OLS函数就是直接写个公式摁上去的。程序主体是矩阵形式，注释里是线性形式。
for i =1:100000
    epsilon = randn(n,1);
    y = 1 * beta_1 + x * beta_2 + epsilon;
    [b1,b2]=OLS_metrix(y,X);
    %[b1,b2,tx,ty] = OLS(y,x);
    %tx = x - x_mean, ty = y - y_mean
    B = B + [b1 b2]';
    uhat = y - b1 - b2 * x;
    V = (1/n * X * X')^-1 * (1/n * X * diag(uhat.*uhat) * X') * (1/n * X * X')^-1;
    sigmahat = V(2,2)/(n-2);
    %sigmahat2 = (1/n)*((1/(n-2))*(tx.*uhat)'*(tx.*uhat))/((1/n*(tx'*tx))^2);
    t = (b2 - beta_2)/(sigmahat^0.5);
    if abs(t) > tinv(0.975,n-2)
        T = T + 1;
    end
end

很可能是n太小了，把n扩大一些试试吧。模型没错，n小就会出现不稳定的问题。

黄瓜人酋长 发表于 2014-10-25 11:57
很可能是n太小了，把n扩大一些试试吧。模型没错，n小就会出现不稳定的问题。

大神请受我一拜