谢谢iooo。我试了一下，方法一不包含常数项。虽然方法2和3的resid中的值为零（不知道为什么），但我利用b中的值解决了问题。刚开始学习，以后请多指教。

不需要截距的b=y/x;e=y-x*b;
需要截距的，xx=x~ones(n,1);b=y/xx;e=y-xx*b;

qingqing 的方法也不错，看来还是要多讨论啊。

load data[100,5] = exer.txt; @导入数据@
t = rows(data); @数据的行@
k = 5; @解释变量个数@
y = data[.,1];  x2 = data[.,2];x3 = data[.,3];x4 = data[.,4];x5 = data[.,5];
yy = y; vny = {"y"};
xx = ones(t,1)~x2~x3~x4~x5;  @解释变量@
vnx = {"cons", "x2", "x3", "x4", "x5" };@输出符号变量@

b = invpd(xx'xx)*(xx'yy); @回归系数估计@
e = yy - xx*b ; @回归残差@
s2 = (e'e)/(t-k); @ e'e是残差平方和，t-k是自由度，s2是回归误差的方差@
vb = s2*invpd(xx'xx); @系数的协方差矩阵@
se = sqrt(diag(vb)); @回归估计量的标准差@
tst = b./se; @ 回归估计量显著性的t统计量@
p = 2*cdftc(abs(tst),t-k); @ t统计量p值@
econ = b~se~tst~p;
econ = vnx~econ;

tss = yy'yy - t*meanc(yy)^2; @总平方和@
rss = e'e; @残差平方和@
ess = tss-rss; @回归平方和@
r2 = 1 - rss/tss; @可决系数@
ar2 = 1-(rss/(t-k))/(tss/(t-1)); @调整的可决系数@
f=(ess/(k-1))/(rss/(t-k)); @方程整体显著性F检验@
d = (e[1:(t-1),1] - e[2:t,1])'(e[1:(t-1),1] - e[2:t,1])/(e'e); @DW统计量@

读者可以看看code后面的解释，其实都是ols的基本理论公式

{b,se,seh,e}=ols(y,x); @调用ols过程@
proc(4)=ols(y,x); @ proc(4)代表着retp中返回的变量有4个@
local n,k,const,xx,b,e,tss,ess,rss,s2,se,r2,ar2,xname,i,xname1,hname,ff,ffn,pr,u,seh;
@local 后面的变量为局部变量，不能与全局变量相同@
    n=rows(x); k=cols(x);
    const=ones(n,1);
xx=const~x;
    b=invpd(xx'xx)*(xx'y);
    e=y-xx*b; @回归残差@
    u=e.*xx;
    tss=(y-meanc(y))'(y-meanc(y)); @总平方和@
    rss=e'e; @残差平方和@
    ess=(xx*b-meanc(xx*b))'(xx*b-meanc(xx*b)); @也可以写为ess = tss-rss;@
    s2=e'e/(n-k-1); @回归误差的方差，注意有常数项@
    se=sqrt(diag(invpd(xx'xx)*s2)); @同方差误差@
    seh=sqrt(diag(invpd(xx'xx)*(u'u)*invpd(xx'xx))); @异方差误差@
    r2=1-rss/tss; @可决系数@
    ar2=1-s2/(tss/(n-1)); @调整的可决系数@
     "R_squared " r2;
    "R_squared (adjusted) " ar2;
    "residual sum of squares (rss) " rss;
retp(b,se,seh,e);
endp;

这是一个proc调用的ols过程，和上面直接计算的结果一样的，只是编程格式不相同

非常感谢xuelida老师！