clear 

n=1000;
e = randn(n,1);
u = 2 + 0.5*randn(n,1);
y(1)=0;
a = -.8*ones(n,1);
b = 1.1*ones(n,1) + ((1:n)>=100)';
theta  = zeros(n,2);
theta2 = zeros(n,2);
theta3 = zeros(n,2);
theta4 = zeros(n,2);
Pt3 = 1000*eye(2);
Pt4 = 1000*eye(2);

for t=2:n,
    y(t)=a(t)*y(t-1)+b(t)*u(t-1)+e(t);
    
    Ut = [y(1:t-1)' u(1:t-1)];
    yt = y(2:t)';
    
    Ut2 = [y(max(1,t-11):t-1)' u(max(1,t-11):t-1)];
    yt2 = y(max(2,t-10):t)';
    
    ut = [y(t-1)' u(t-1)];
    ys = y(t);
    
    % lS
    theta(t,:) = (inv(Ut'*Ut)*(Ut'*yt));

    % LS, window
    theta2(t,:) = (inv(Ut2'*Ut2)*(Ut2'*yt2));    

    % RLS
    Pt3 = Pt3 - (Pt3*ut'*ut*Pt3)/(1 + ut*Pt3*ut');
    theta3(t,:) = theta3(t-1,:) + ut*Pt3*(ys - ut*theta3(t-1,:)'); 

    % RLS, weighting
    lam = .99;
    Pt4 = (Pt4 - (Pt4*ut'*ut*Pt4)/(lam + ut*Pt4*ut'))/lam;
    theta4(t,:) = theta4(t-1,:) + ut*Pt4*(ys - ut*theta4(t-1,:)'); 

    
    % plot
    subplot(2,1,1); plot(1:t,y(1:t))
    xlabel('time t')
    ylabel('y(t)')
    grid;
    axis([1 t -5 5])

    subplot(2,1,2)
    plot(1:t,a(1:t),'k-',...
        1:t,b(1:t),'k-.'); hold on;
    %plot(1:t,theta3(1:t,1), 'g-',...
    %    1:t,theta3(1:t,2),'r-'); 
    %plot(1:t,theta2(1:t,1), 'g-.',...
    %    1:t,theta2(1:t,2),'r-.');
    plot(1:t,theta4(1:t,1), 'g-',...
        1:t,theta4(1:t,2),'r-');
    hold off;

    
    xlabel('time t')
    ylabel('parameter')
    legend('a(t)','b(t)','a_t(t)','b_t(t)')
    grid;
    axis([1 t -1 3])
    drawnow
    %pause(.01)
end
