这个帖子很好，欢迎大家踊跃参与讨论，参与者均有奖励。

proc bandw(vhat);
local nt,i,b,sig,a2,a2n,a2d,st,d;
@procedure that compute the automatic bandwidth based on AR(1)
approximation for each vector of the matrix vhat. Each are given
equal weigths of 1.@

nt=rows(vhat);
d=cols(vhat);
a2n=0;
a2d=0;
i=1;
do while i <= d;
b=olsqr(vhat[2:nt,i],vhat[1:nt-1,i]);
sig=(vhat[2:nt,i]-b*vhat[1:nt-1,i])'(vhat[2:nt,i]-b*vhat[1:nt-1,i]);
sig=sig/(nt-1);
a2n=a2n+4*b*b*sig*sig/(1-b)^8;
a2d=a2d+sig*sig/(1-b)^4;
i=i+1;
endo;
a2=a2n/a2d;
st=1.3221*(a2*nt)^.2;
retp(st);
endp;

这一段是计算帘宽的程序

proc kern(x);
local del,del1,ker;
@procedure to evaluate the quadratic kernel at some value x.@

del=6*pi*x/5;
ker=3*(sin(del)/del-cos(del))/(del*del);
retp(ker);
endp;

上面这段是核函数的程序。