function [yhat, th, model, state] =...
pem1pred(k,y,u,sequence,init_pred,lambda,model_in,...
        state_in,p,na_row,nb_row,c_row,nd);

% PEM1PRED Iterative RPEM model-sequence k-step prediction with a number of
%     recursive predictors. Yields YHAT when known/prediced U are supplied,
%     because the predictors for the last known measured data are stored.
%
%     [YHAT TH MODEL STATE] = ...
%       pem1pred(K,Y,U,SEQUENCE,INIT_PRED,LAMBDA,MODEL_IN,...
%               STATE_IN,P,NA_ROW,NB_ROW,C_ROW,ND)
%
%     or
%
%     [YHAT TH MODEL STATE] = ...
%       pem1pred(K,Y,U,SEQUENCE,INIT_PRED,LAMBDA,MODEL_IN,STATE_IN)
%
%     The model   y = (B1/A1) * u1 + ... + (Bn/An) * un + (C/D) * e
%     is estimated using a recursive prediction error method (RPEM).
%
%     YHAT is the iterated k-step predictor. It is situated k steps
%     forward comparad with the inputs. TH is parameter curveswith 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. STATE contains values for dW/dA, dW/dB and psi.
%
%     K represents the horizon for the predictor to be calculated, Y and U are out- 
%     and in-signal as colonvectors. 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 using 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 []. STATE_IN is dW/dA, dW/dA and psi for earlier estimated
%     models.
%
%     NA_ROW is a row vector containg orders of the A-polynomials, NB_ROW are orders 
%     for the B-polynomials. C_ROW is a row vector containing the (positive)
%     delays that shall have parameters estimeted in the C-polynomial. 
%     ND represents the order of the D-polynomial.

%     Roger Felix, 6 June 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; % (if)
  end; % (for)
end; % (for)

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 lambs(i) < 0, error('Too small LAMBDA'), end;
end; % (for)

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

if isempty(model_in),        % CREATE MODEL
   [ma na] = size(na_row); [mb nb] = size(nb_row);
   nb_row = nb_row + ones(mb, nb); % NO OF PARAMETERS, NOT ZEROES ! 
   if na ~= nb, error('NA_ROW and NB_ROW must be of equal length');end;
   [mc nc] = size(c_row); c_len = max([0 c_row]);
   nTheta = nc+nd+sum(nb_row)+sum(na_row);
   model_size = [4+nTheta, max([8+2*na+nc,k+c_len,nd+sum(na_row) ...
                max(na_row)+sum(nb_row),(nTheta+1) * NoOfModels])];
   model = zeros(model_size(1), model_size(2)); eloc = 1; 
   model(1, 1:8+2*na+nc) = ...
       [na nb nc nd init_pred interval eloc p na_row nb_row c_row];
else                         % FETCH MODEL
   model = model_in; model_size = size(model);
   state = state_in; [ModelCount, nState] = size(state);
   na = model(1, 1); nb = model(1, 2); nc = model(1, 3);nd = model(1, 4);
   init_p = model(1, 5); interval = model(1, 6);
   eloc = model(1, 7);
   p = model(1, 8);
   if init_p ~= init_pred, 
      if init_pred < 1, 
         init_pred = init_p;
      else
         disp('Warning! --- non-cyclic calculation !');
      end
   end; % (if)
   na_row = model(1, 9:8+na); 
   nb_row = model(1, 9+na:8+2*na); 
   c_row = model(1, 9+2*na:8+2*na+nc); 
   nTheta = sum(na_row)+nc+nd+sum(nb_row); 
end; % (if)

[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)

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

    % THE INDEX OF THE PSI-VECTORS 
    % ----------------------------

c_len = max([0 c_row]);
a_len = ones(1,na).*max([nd,na_row+c_len]);
b_len = ones(1,nb).*max([c_len+na_row,nb_row,nd]); 
d_len = max([c_len,nd]);

a1_psi = 1+cumsum(a_len)-a_len(1);
b0_psi = sum(a_len)+1+cumsum(b_len)-b_len(1);

a_psi = 1 : sum(a_len);
b_psi = a_psi(length(a_psi)) + [1 : sum(b_len)];
c_psi = b_psi(length(b_psi)) + [1 : sum(c_len)];
if c_len > 0,
   d_psi = c_psi(length(c_psi)) + [1 : sum(d_len)];
else
   d_psi = b_psi(length(b_psi)) + [1 : sum(d_len)];
end; %(if)

psi_index = [];
for i = 1:na,
   psi_index = [psi_index (a1_psi(i):a1_psi(i) + na_row(i) - 1)];
end; %(for)
for i = 1:nb,
   psi_index = [psi_index (b0_psi(i) : b0_psi(i) + nb_row(i)-1)];
end; %(for)
if c_len > 0,
   psi_index = [psi_index (c_psi(1)+c_row-1) (d_psi(1)+[1:nd]-1)];
else
   psi_index = [psi_index (d_psi(1)+[1:nd]-1)];
end; %(if)

     % INDEX FOR dWdA AND dWdB
     % -------------------------

% dWdA_len = sum(na_row);
% dWdB_len = sum(na_row);

dWdA1_index =  1+d_psi(length(d_psi))+cumsum(na_row) - na_row(1);
dWdB1_index =  1+d_psi(length(d_psi))+sum(na_row)+...
                 cumsum(max(na_row,max(nb_row,nd*ones(size(nb_row)))))...
              - max(na_row(1),max(nb_row(1),nd));

dWdA_index = dWdA1_index(1) - 1 + [1: sum(na_row)];
dWdB_index = dWdB1_index(1) - 1 + ...
             [1: sum(max(max(nb_row,nd*ones(size(nb_row))), na_row))];

if isempty(state),
    nState = dWdB_index(length(dWdB_index));
    state = zeros(NoOfModels,nState);
else
    if nState ~= dWdB_index(length(dWdB_index)),
        error('STATE is of wrong size');
    end; % (if)
end; % (if)

     % INDEX OF THE THETA VECTORS 
     % --------------------------

a_theta = 1:sum(na_row);
b_theta = sum(na_row) + 1 : sum(na_row) + sum(nb_row);
c_theta = sum(na_row) + sum(nb_row) + 1 : ...
          sum(na_row) + sum(nb_row) + nc;
d_theta = sum(na_row) + sum(nb_row) + nc + 1 : ...
          sum(na_row) + sum(nb_row) + nc + nd;

a1_theta = 1+cumsum(na_row)-na_row(1);
b0_theta = 1+sum(na_row)+cumsum(nb_row)-nb_row(1);

th = zeros(nTheta, my); 
yhat = zeros(1, my) ; yth = zeros(k, 1);             
ukol = zeros(sum(nb_row), 1); wkol = zeros(sum(na_row),1);

     % INITIATE W_SEKV, NY_SEKV, U & Y WITH OLD VALUES
     % ------------------------------------------------

offset = max([na_row nb_row nd]);
w_sekv = zeros(offset+my,nu); ny_sekv = zeros(offset+my,1);
u = [zeros(offset, nu); u]; y=[zeros(offset, 1); y];

y(offset:-1:offset-max(na_row)+1, 1) = model(4, 1:max(na_row))';
ny_sekv(offset:-1:offset-nd+1, 1)  = model(3, 1:nd)';
for i = 1:nu, 
    u(offset:-1:offset-nb_row(1, i)+1, i) = ...
        model(4, max(na_row)+sum(nb_row(1:i-1))+i:...
        max(na_row)+sum(nb_row(1:i))+i-1)';

    w_sekv(offset:-1:offset- na_row(1,i)+1,i) = ...
        model(3, nd+sum(na_row(1:i-1))+1:...
        nd+sum(na_row(1:i)))';
end; % (for)

Prows = 5:nTheta+4; THrows = Prows; kcmax = k+c_len;

    % INITIATE COVARIANCE MATRICES
    % ----------------------------

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

end; % (if)

for i = 1:my                             % FOR ALL KNOWN Y

       % CALCULATE INDEX
       % ---------------

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

       % UPDATE THETA
       % ------------

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

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

    model(Prows, Pindex) = (model(Prows, Pindex) - K * ...
           state(predictor,psi_index)*model(Prows, Pindex)) / ...
           lambdas(predictor);
    D_old(1,1:nd) = ...
          model(THrows(d_theta), predictor)';
    for nn = 1:nu, % SAVE OLD A-parameters
       A_olds(nn,1:na_row(nn)) = ...
          model(THrows(a1_theta(nn):a1_theta(nn)+na_row(nn)-1), predictor)';
    end; % (for)

    model(THrows, predictor) = ...
        model(THrows, predictor) + K * epsi; 

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

    %  STABILIZE (This is the best method for D)
    %  -----------------------------------------
    if nd > 0
        D = [1 model(THrows(d_theta), predictor)'];
       if length(D) > 1,
            if any((abs(roots(D))+eps-1) > 0), % INSTABLE ROOT ?
                D = [1 D_old(1,1:nd)];
            end; % (if)
        end; %(if)
        model(THrows(d_theta),predictor) = ...
                      D(1,2:nd+1)';
    end; %(if)

    % STABILIZE A 
    % -------------

    for signal = 1:nu,
        A = ([1 model(THrows(a1_theta(signal):...
                      a1_theta(signal)+na_row(signal)-1), predictor)']);
        if length(A) > 1,
            if any((abs(roots(A))+eps-1) > 0), % INSTABLE ROOT ?
                A = [1 A_olds(signal,1:na_row(signal))];
            end; % (if)
            model(THrows(a1_theta(signal):...
                  a1_theta(signal)+na_row(signal)-1),predictor) = ...
                        A(1,2:na_row(signal)+1)';
        end; %(if)
    end; % (for)

    % CALCULATE w_sekv AND ny_sekv AND ERROR USING NEW THETA
    % ------------------------------------------------------

    for signal = 1:nu,  % CALCULATE NEW w AND ny
       %  A w_sekv = B u   <=>   w_sekv(0) = B u - (A - 1) w_sekv  
       w_sekv(ii, signal) = model(THrows([ ...
           (a1_theta(signal):a1_theta(signal)+na_row(signal)-1) ...
           (b0_theta(signal):b0_theta(signal)+nb_row(signal)-1)]), ...
           predictor)' * [-w_sekv(ii-1:-1:ii-na_row(signal),signal); ...
                 u(ii :-1:ii-nb_row(signal)+1, signal)];
 
    end; %(for)

    ny_sekv(ii,1) = y(ii,1)-sum(w_sekv(ii,:));

    model(2,eloc) = ny_sekv(ii,1) - ...
        model(THrows([c_theta d_theta]),predictor)' *...
        [model(2,1+rem(eloc-1+kcmax-c_row,kcmax)) ...
         -ny_sekv(ii-1:-1:ii-nd,1)']';
 

    % UPDATE ALL PSI'S BY FILTERING U & Y
    % -----------------------------------

    for n = 1: NoOfModels,      

        % FETCH D & C
        % -----------

        D = [1; model(THrows(d_theta), n)];

        if c_len > 0,  
            C_explicit = [1 zeros(1,c_len)];
            C_explicit(1,1+c_row) = ...
                model(THrows(c_theta), n)';
        end; % (if)
        C_explicit(1,1) = 1;
        
        % CALCULATE NEW dWdA, dWdB, PsiA AND PsiB
        % --------------------------------------

        for signal = 1:nu,  

            A = [1 model(THrows(a1_theta(signal):...
                      a1_theta(signal)+na_row(signal)-1), n)'];

            %  A dWdA = (0 | -w_sekv(t-j)) =>
            %    dWdA(t) = -(A - 1) dWdA - (0 | w_sekv(t-j))

            new_dWdA(signal) = -  A(2:1+na_row(signal)) * ...
              state(n, dWdA1_index(signal):dWdA1_index(signal) + ...
              na_row(signal)-1)' - (n == predictor) * w_sekv(ii,signal);

            %  C PsiA = D dWdA   <=>  PsiA(0) = D dWdA - (C - 1) PsiA

            newPsiA(signal) =  D' * [new_dWdA(signal) ...
                state(n,[dWdA1_index(signal):dWdA1_index(signal)+nd-1])]' ...
                - C_explicit * [0 state(n, a1_psi(signal): ...
                    a1_psi(signal)+c_len-1)]';

            %  A dWdB = (0 | u(t-j))  =>
            %     dWdB(t) = -(A - 1) dWdB + (0 | u(t-j))
            new_dWdB(signal) = - A(2:1+na_row(signal)) * ...
              state(n,dWdB1_index(signal):dWdB1_index(signal)+ ...
              na_row(signal)-1)' + (n == predictor) * u(ii,signal);

            %  C PsiB = D dWdB  <=>  PsiB(0) = D dWdB - (C - 1) PsiB
            newPsiB(signal) = D' * [new_dWdB(signal) ...
            state(n,dWdB1_index(signal):dWdB1_index(signal)+nd-1)]' - ...
                C_explicit * [0  state(n, b0_psi(signal):b0_psi(signal)+...
                    c_len-1)]';

        end; % (for)

           % DETERMINE NEW PsiC AND PsiD
           % --------------------------

        %  C PsiC = e   <=>  PsiC(0) = e - (C - 1) PsiC  
        if c_len > 0     
             newPsiC = (n == predictor) * model(2,eloc) - ...
                        C_explicit * [0 state(n, c_psi)]';
        end; % (if)

        %  D PsiD = -e   <=>  PsiD(0) = -e - (D - 1) PsiD  
        if nd > 0  
            newPsiD = - (n == predictor) * model(2,eloc) - ...
                      D' * [0 state(n, d_psi(1:nd))]';
        end; % (if)


        % SHIFT psi, dWdA AND dWdB
        % -------------------------

        state(n, 2 : nState) =state(n, 1 : nState - 1);

        for signal = 1:nu 
            state(n,b0_psi(signal)) = newPsiB(signal);
            state(n,a1_psi(signal)) = newPsiA(signal);
            state(n,dWdA1_index(signal)) = new_dWdA(signal);
            state(n,dWdB1_index(signal)) = new_dWdB(signal);
        end; % (if) 

        if c_len > 0,
            state(n, c_psi(1)) = newPsiC;
        end; % (if)

        if nd > 0, 
            state(n, d_psi(1)) = newPsiD;
        end; % (if)

    end; % (for)

         % COPY OLD NY_HAT AND W_VEC_HAT
         % ----------------------------------

    ny_hat = ny_sekv(ii:-1:ii-nd+1,1);

    w_vec_hat(1:max(na_row),:) = w_sekv(ii:-1:ii-max(na_row)+1,:);

        % ITERATIVE PREDICTION
        % -------------------

    for j = 1:k,                % FOR FORECASTING HORIZON 1..k

        predictor = sequence(1+rem(init_pred-1+j+interval, interval));

        model(2, 1+rem(kcmax+eloc-1+j, kcmax)) = 0; % ERASE OLD
                                                    % VALUES OF e.

        % SHIFT & CALCULATE w_vec_hat USING LOCAL MODEL
        % -------------------------------------------

        w_vec_hat(2:max(na_row)+1,:) = w_vec_hat(1:max(na_row),:);

        % A w_vec_hat = B u <=> w_vec_hat(0) = B u - (A - 1) w_vec_hat 
        for signal = 1:nu,
            w_vec_hat(1,signal) = model(THrows([ ...
            (a1_theta(signal):a1_theta(signal)+na_row(signal)-1), ...
            (b0_theta(signal):b0_theta(signal)+nb_row(signal)-1)]'), ...
            predictor)' *  [-w_vec_hat(2:na_row(signal)+1,signal); ...
            u(ii+j:-1:ii+j-nb_row(signal)+1, signal)];
        end; % (for)

        % SHIFTA & CALCULATE ny_hat USING LOCAL MODEL
        % ----------------------------------------
        % ny_hat =  C e_hat - (D - 1) * ny_hat 
        ny_hat(2:nd+1,1) = ny_hat(1:nd,1);
        ny_hat(1,1) = model(THrows([c_theta, d_theta]),predictor)' *...
             [model(2,1+rem(kcmax+j+eloc-1-c_row,kcmax)), -ny_hat(2:nd+1,1)']';

        % CALCULATE y_hat USING LOCAL MODEL
        % ------------------------------

        yth(j) = sum(w_vec_hat(1,:)) + ny_hat(1,1);

    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)

     % SAVE W_SEKV, NY_SEKV, IN-SIGNALS AND OUT-SIGNAL
     % ------------------------------------------

model(4,1:max(na_row)) = y(my+offset:-1:my+offset-max(na_row)+1, 1)' ;
model(3,1:nd) = ny_sekv(my+offset:-1:my+offset-nd+1, 1)' ;
for i = 1:nu,          
   model(4, max(na_row)+sum(nb_row(1:i-1))+i:...
         max(na_row)+sum(nb_row(1:i))+i-1) = ...
         u(my+offset:-1:my+offset-nb_row(1, i)+1, i)';
   model(3, nd+sum(na_row(1:i-1))+1:...
         nd+sum(na_row(1:i))) = ...
         w_sekv(my+offset:-1:my+offset-na_row(1, i)+1, i)';
end; % (for)

model(1, 5) = init_pred;
model(1, 7) = eloc;