function e = ARMAfilter(rho,psi,m,Y,MALag,ARLag)

[ignore,k]=size(m);
[T,one]=size(Y);

ARMALag = max(ARLag,MALag);
UncMean = m(1,:)./(1 - sum(rho,1) - sum(psi,1));

e = zeros(T+ARMALag,k);
Y = [ repmat(UncMean,ARMALag,1) ; repmat(Y,1,k) ];
for i=1:T
   e(i+ARMALag,:) = Y(i+ARMALag,:) - m(i,:);
   for j=1:ARLag
      e(i+ARMALag,:) = e(i+ARMALag,:) - rho(j,:).*Y(i+ARMALag-j,:);
   end
   for j=1:MALag
      e(i+ARMALag,:) = e(i+ARMALag,:) - psi(j,:).*e(i+ARMALag-j,:);
   end
end
e = e(ARMALag+1:end,:);