function [thetahat,P,rhat,kappa,ehat]=rlsf(thetahat,P,rhat,kappa,y,psi,alpha)
%
% recursive least squares for the regression y=theta'psi+noise
% with exponential forgetting alpha 
%
   zeta=psi'*P*psi; 
   ehat=y-thetahat'*psi;
   thetahat=thetahat+P*psi/(1+zeta)*ehat;
   P=(alpha*P^(-1)+psi*psi')^(-1); 
   lam=rhat*kappa;
   kappa=alpha*kappa+1;
   lam=alpha*lam+ehat^(2)/(1+zeta);
   rhat=lam/kappa;