function e=surprise(y,lags);
[T,N]=size(y);
e=y;

for i=1:N;
  reg=ones(T,1);
  for j=1:lags;
    reg=[reg lagn(y(:,i),j)];
  end;
  yy=y(lags+1:T,i);
  x=trimr(reg,0,lags);
  beta=x\yy;
  e(lags+1:T,i)=yy-x*beta;
end;  
