loadm y[100,1]=c:\temp\ykpss.txt;

x1=ones(100,1);
x2=ones(100,1)~seqa(1,1,100);


res=y-x1*(inv(x1'x1)*x1'y);


bartlett_win=10;
sum=0;
integers=seqa(1,1,bartlett_win);
nobs=rows(res);
res_sq=res.^2;
s=cumsumc(res);
s_sq=s.^2;
w=ones(bartlett_win,1)-integers./(bartlett_win+1);

for l (1,bartlett_win,1);
   row=l;
   res_lag=lagn(res,l);
   res_unlag=trimr(res,row,0);
   res_lag=trimr(res_lag,row,0);
    sum_temp=w[l]*(res_unlag'*res_lag);  
   sum=sum_temp+sum;
endfor;
sum_dub=sum;

sig_sq=inv(nobs)*sumc(res_sq)+2*inv(nobs)*sum_dub;
LM=1/(nobs*nobs)*sumc(s_sq)/sig_sq;