Data Identification using recursive square method

Hello guys I have to identify data which contains input of the system and output of the system. I have to identified the parameters such as b0, b1, b2, b3, a1, a2, a3. Then I'll use transfer function and it should give me something like this:

%    b0*Z^3 + b1*Z^2 + b2*Z + b3
%  ------------------------------
%   a0*Z^3 + a1*Z^2 + a2*Z + a3

I was able identify the data using least square method, but I'll use it as adaptive system. So I have to use recursive least square method with forgetting factor lambda. I tried to implement some functions from file exchange, but without successful result. So can somebody please help me? In the picture is example how it should looks like.

Also I attached data which I am using. Where U is input and Y is output of the system. Also my code for least square method.

Starting script

hold off
clc
clear
load data.mat   %meassured data
T=1:1:1170;     %time
t=T';           %Transposition
figure(1)
plot(t,Y,'r')
hold on
plot(t,U,'g')
title('Data from system')
legend('Output','Input')
grid on

u=U;        %Input 
y=Y;        %Output

denominator=0;

%-----------
[bd,ad] = lsf(3,3,u,y)    %least square function
A=sum(bd)/sum(ad);          %amplification of system
system = tf(bd,ad,1)        %transfer function (Numerator,Denominator)
%

for s = 1:1:length(ad)
    denominator = denominator + ad(1,s);
end
numerator = bd(1,4);

disp('Amplification of system: ')
disp(A)
disp('Numerator: ')
disp(numerator)

graph = lsim(system,u, t,'zoh'); 
%lsim Simulate time response of dynamic system to arbitrary inputs

figure(2)
plot(t,graph, t, y,'r')     %red line in graph symbolize output of given data
title('Comparison')
legend('Identified','Given data')
grid on

Function lsf

function [B, A] = lsf(nb, na, u,y)
N = length(u); 
Y = y(na+1:N);  

%empty matrix for u and y
Xy=[]; Xu=[];   

if(nb==1)
    Xu = [u(2:N-(na-nb))];   
    for i=1:na  
        Xy = [Xy -y(na+1-i:N-i)];  
    end
elseif(nb==2)
      for i=1:nb
        Xu = [Xu  u(nb+1-i:N-i-1)]; 
      end
      for i=1:na  
        Xy = [Xy -y(na+1-i:N-i)];  
      end
else
      for i=1:nb
        Xu = [Xu  u(nb+1-i:N-i)]; 
      end
      for i=1:na  
        Xy = [Xy -y(na+1-i:N-i)];  
      end
end

X = [Xy Xu];  

P = inv(X'*X)*X'*Y;     %calculation of the least square method
A = [1 P(1:na)'];       %denominator
B = [ 0 P(na+1:end)'];  %numerator

[John D'Errico]

@mohd albahrani - Please don't answer, when you are just trying to make a comment. Moved a non-answer into a comment:

"I'm trying to do the same but couldn't build the right code yet."