function [a_t_rt,P_t_rt]   = var_KKInsample(prior_theta_mean,prior_theta_V,YY_tt,m,Sigma_0,kappa,lambda,lags,const)

T                 = size(YY_tt,1);
% position of the first lags
K = lags*m+1;
[yy_ll,xx_ll]  = matrix_var(YY_tt,lags,T,m);

if const==1
    y_l            = [NaN(lags,size(yy_ll,2));yy_ll];
    x_l            = [NaN(lags,size(xx_ll,2));xx_ll];
else
    y_l            = [NaN(lags,size(yy_ll,2));yy_ll];
    x_l            = [NaN(lags,size(xx_ll,2)-1);xx_ll(:,1:end-1)];
end

for irep =lags+1:T

%     % prediction step for States
    if irep==lags+1
        a_t(:,irep)                  = prior_theta_mean; %use prior mean for first observation
        P_t                          = prior_theta_V;    %use prior variance for first observation
    else
        % prediction step
        a_t(:,irep)                  = a_t_rt(:,irep-1);
        P_t                          = (1./lambda)*P_t_rt; % prediction step: this should be P_t_rt + Q_t but use forgetting factors instead
    end
    % Now one step ahead prediction
    Z_t              =  kron(eye(m),x_l(irep,:)); % measurement error matrix
    y_t_pred(irep,:) = (Z_t*a_t(:,irep))';  % this is one step ahead prediction
    
    % Prediction error
    e_t(irep,:)     = y_l(irep,:)-y_t_pred(irep,:);
    %     disp(e_t)
    %   Now update Kalman filter matrices based on the prediction error
    %   First update H[t], the variance needed for the prediction error variance
    A_t = e_t(irep,:)'*e_t(irep,:);
    
    if irep==lags+1
        H_t   = kappa*Sigma_0;
    else
        prev_H_t = H_t;
        H_t   = kappa*H_t + (1-kappa)*A_t;
    end
    %   if all(eig(H_t) < 0)
    %       H_t = prev_H_t; % do not update the matrix
    %   end
    %     PREDICTION error variance and K-gain
    F_t   = Z_t*P_t*Z_t'+H_t;
    %disp(sum(sum(F_t)));
    K_gain =single((F_t\Z_t*P_t')'); %NB this is like in Durbin/Koopman pagina 85, where this is multiplied by T_t, but here T_t is an identity
    % UPDATE STEP
    P_t_rt                       = (P_t - K_gain*Z_t*P_t);
    a_t_rt(:,irep)               = a_t(:,irep)+K_gain*e_t(irep,:)'; % update step: the matrix T_t in the transition is an identity
end % End cycling through nom models with different shrinkage factors