function [yhat, th, model] =...
lspred(k,y,u,sequence,init_pred,lambda,model_in,p,na,nb_row,c_row);

% LSPRED Iterative LS model-sequence k-step prediction with a number of
%     recursive predictors. Yields YHAT when known/prediced U are supplied,
%     if the predictors for the last known measured data are stored.
%
%     [YHAT TH MODEL] = ...
%       lspred(K,Y,U,SEQUENCE,INIT_PRED,LAMBDA,MODEL_IN,P,NA,NB_ROW,C_ROW)
%
%     or
%
%     [YHAT TH MODEL] = ...
%       lspred(K,Y,U,SEQUENCE,INIT_PRED,LAMBDA,MODEL_IN)
%
%     The model y = (B1/A) * u1 + ... + (Bn/A) * un + (C/A) * e
%     is estimated using least squares method (LS).
%
%     YHAT is the iterated k-step predictor. It is situated k steps
%     forward comparad with the inputs. TH is parameter curves with parameters
%     in corresponence with Y and INIT_PRED (see below). MODEL are the recursive 
%     predictors ready to be used as MODEL_IN for the next
%     input. 
%
%     K represents the horizon for the predictor to be calculated, Y and U are out- 
%     and in-signal as colon vectors. SEQUENCE represents the ordering
%     of predictors. SEQUENCE may contain integers > 0.
%     INIT_PRED indicates the sequence number of the predictor for the first element
%     in the u- and y-vectors. This makes it possible to store the predictors with
%     unchanged information about their placi in the cycle. It is set to 1 when
%     not usig cyclic models. If INIT_PRED < 1 then INIT_PRED is fetched from MODEL_IN.
%
%     LAMBDA corresponds to the decay parameter in the recursive algoritm having an
%     cycle length of 1. P is the starting value of the diagonal elements of the 
%     covariance matrix.
%     MODEL_IN is a model with previously estimated predictors. If this parameter is
%     missing, use [].
%
%     NA is the order of the A-polynomial, NB_ROW is a row vector with orders 
%     for the B-polynomials. C_ROW is a row vector containing the (positive)
%     delays that shall have parameters estimeted in the C-polynomial.

%     Roger Felix, 26 March 1996.
%

[ns interval] = size(sequence);
NoOfModels=max(sequence); ModelCount = zeros(1,NoOfModels);

for i = 1:NoOfModels,    % COUNT NO OF OCCURENCIES
  for j = 1:interval,
    if sequence(1,j) == i,
        ModelCount(i) = ModelCount(i) + 1;
    end;
  end;
end;

for i=1:NoOfModels,      % EVERY MODEL HAS EQUAL FORGETTING RATE
   if ModelCount(i) > 0,
      lambdas(i) = 1 - interval * (1 - lambda)/ModelCount(i);  
   else
      lambdas(i) = 1;
   end; % (if)
       if lambdas(i) < 0, error('Too small LAMBDA'),end;
end;

if init_pred > interval,
   error('INIT_PRED must be within 1..INTERVAL');
end;

if isempty(model_in),        % CREATE MODEL
   [mb nb] = size(nb_row); [mc nc] = size(c_row);
   nTheta = na+nc+sum(nb_row)+nb;
   model_size = [3+nTheta, max([7+nb+nc, (nTheta+1) * NoOfModels, k+c_row])];
   model = zeros(model_size(1),model_size(2)); eloc = 1;
   model(1, 1:7+nb+nc) = [na nb nc init_pred interval eloc p nb_row c_row];
else                         % FETCH MODEL
   model = model_in; model_size = size(model);
   na = model(1,1); nb = model(1,2); nc = model(1,3);
   init_p = model(1,4); interval = model(1,5);
   eloc = model(1,6); p = model(1,7);
   if init_p ~= init_pred,
      if init_pred < 1,
         init_pred = init_p;
      else
         disp('Warning! --- non-cyclic calculation !');
      end
   end;
   nb_row = model(1, 8:7+nb); c_row = model(1, 8+nb:7+nb+nc); 
   nTheta = na+nc+sum(nb_row)+nb;       
end; 

[mu nu] = size(u); [my ny] = size(y);
if nu ~= nb
   error('The number of input signals must equal the number of orders nb_row'); 
end

if mu < my+k,
   error('Exogenous signals must exist for all Y(1..i+k)');
end

th = zeros(nTheta, my); 
fi = zeros(nTheta, 1); yhat = zeros(1, my) ; yth = zeros(k, 1);             
ukol = zeros(sum(nb_row)+nu, 1);      

offset = max([na nb_row]);
u=[zeros(offset, nu); u]; y=[zeros(offset,1); y];
for i=1:nu,
   u(offset:-1:offset-nb_row(1,i)+1,i) = ...
   model(3,na+sum(nb_row(1:i-1))+i:na+sum(nb_row(1:i))+i-1)';
end;

y(offset:-1:offset-na+1,1) = model(3,1:na)';

Prows = 4:nTheta+3; THrows = Prows; kcmax = k+max([0 c_row]);

if nargin > 8,               % COVARIANCE MATRICES
   for i = NoOfModels:nTheta:NoOfModels*(nTheta+1)-nTheta,          
      for j = Prows, 
         model(j, i+j-3) = p; 
      end; % (for)
   end; % (for)
end; % (if)

for i = 1:my,                           % FOR ALL KNOWN Y

   ii=i+offset;                         % ii =INDEX(u,y)
   predictor=sequence(init_pred);
   model(2,eloc)=y(ii)-model(2,eloc);   % epsi = Y-YHAT(0)

   Pindex = ...
    (NoOfModels+(predictor-1)*nTheta+1:NoOfModels+predictor*nTheta); 

   for m = 1:nu, 
      ukol(sum(nb_row(1:m-1))+m:sum(nb_row(1:m))+m) =...
          u(ii:-1:ii-nb_row(m),m);
   end % (for)

   if nc == 0,
      fi = [-y(ii-1:-1:ii-na); ukol]; 
   else
      fi = [model(2, 1+rem(eloc-c_row-1+kcmax, kcmax))';...
         -y(ii-1:-1:ii-na); ukol];  
   end; % (if)


   K = (model(Prows, Pindex) * fi)/ ... 
      (lambdas(predictor)+fi'*model(Prows, Pindex)*fi);

   th(:, i) = model(THrows, predictor) + ...
      K*(y(ii)-fi'*model(THrows, predictor));

   model(THrows, predictor) = th(:, i); 

   model(Prows, Pindex) = ... 
      (model(Prows, Pindex) - ...
      (model(Prows, Pindex)*fi*fi'*model(Prows, Pindex)) / ...
      (lambdas(predictor)+fi'*model(Prows, Pindex)*fi))/lambdas(predictor);  

   maxP = max(Pindex);

   if model(model_size(1),maxP) > p,
      model(model_size(1),maxP) = p;    % NO IN-SIGNAL ?
   end; 


   for j = 1:k,                             % FOR FORECASTING HORIZON 1..k
      model(2,1+rem(eloc-1+j, kcmax)) = 0;  % ERASE OLD VALUES OF e.
      for m = 1:nu,     
         ukol(sum(nb_row(1:m-1))+m:sum(nb_row(1:m))+m) = ...  
            u(ii+j:-1:ii+j-nb_row(m), m); 
      end % (for)

      if nc == 0,                          
         fi = [-yth(j-1:-1:max(j-na, 1)); -y(ii:-1:ii+j-na); ukol]; 
      else
         fi = [model(2,1+rem(kcmax+eloc-1+j-c_row, kcmax))'; ...
         -yth(j-1:-1:max(j-na, 1)); -y(ii:-1:ii+j-na); ukol]; 
      end; % (if)

      yth(j) = fi'*model(THrows,...
        sequence(1+rem(init_pred-1+j+interval,interval)));
   end % (for)

   eloc=eloc+1;
   if eloc > kcmax, eloc= eloc-kcmax;end;

   init_pred=init_pred+1;
   if init_pred > interval, init_pred=init_pred-interval;end;

   model(2,eloc)=yth(1);         % SAVE THE PREDICTION (!)
   
   yhat(1,i) = yth(k);                         % HERE YHAT IS DELIVERED

end % (for)

for i=1:nu,          
   model(3,na+sum(nb_row(1:i-1))+i:na+sum(nb_row(1:i))+i-1) = ...
   u(my+offset:-1:my+offset-nb_row(1,i)+1,i)';
end;
model(3,1:na)=y(my+offset:-1:my+offset-na+1,1)' ;
model(1,4)=init_pred;
model(1,6)=eloc;