function [FORGET,a_t_rt,lambda_t,pred_e,pred_F]   = var_KK_3(prior_theta_mean,prior_theta_V,YY_tt,data_levels,m,Sigma_0,kappa,hmax,NN_VAR,lags,esti_start,esti_end_d,std_target,med_target,const)
% MAKE FORGETTING FACTOR TIME VARYING
lambda_min = .96;
L          = 1.1;
T                 = size(YY_tt,1);
data_lev          = data_levels;
first_fore        = esti_start+1;

if esti_end_d==T-hmax
    % this is fine
else
    error('error in dates')
end
% position of the first lags
% K = lags*m+1;
% for gg=1:m;
%     ind_1l(gg) = (gg-1)*K+gg;
% end
count_tt =0;
nfore = esti_end_d-esti_start+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
lambda = ones(esti_end_d-lags+1,1);
for irep =lags+1:esti_end_d
%    tic
%     if rem(irep,20)==0
%         disp(['irep = ',num2str(irep)]);
%     end
    % 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
    else
        % prediction step with time varying forgetting factor
        a_t(:,irep)                  = a_t_rt(:,irep-1);
        P_t                          = (1./lambda(irep-1))*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,:);
    % Update time varying forgetting factor
    lambda(irep)    = lambda_min+(1-lambda_min)*L.^round(-e_t(irep,1:NN_VAR)*e_t(irep,1:NN_VAR)');
    %   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
    %     pos_expl                     = find(a_t_rt(ind_1l,irep)>1);
    %     if irep>1
    %         if pos_expl>0
    %             a_t_rt(ind_1l(pos_expl),irep)               = 1;% a_t_rt(ind_1l(pos_expl),irep-1);
    %         end
    %     end
    if (irep>=esti_start)&&(irep<=esti_end_d) %compute out of sample forecasts
        count_tt = count_tt+1;
          % STORE OUTPUT FOR MODEL SELECTION LATER
        pred_F(count_tt,:,:) = F_t(1:NN_VAR,1:NN_VAR);
        pred_e(count_tt,:)  = e_t(irep,1:NN_VAR) ;        
        %         disp(count_tt)
        %replace the first lag in the x_l matrix with data up to time "t"
        a_t_fore      = a_t_rt(:,irep);
        y_t_fore      = NaN(hmax,m);
        y_f           = y_l(irep,:);
        x_f           = x_l(irep,:);
        for llags = 1:hmax
            x_ff  = [y_f x_f(1,1:m*lags-m) 1];
            x_f   = x_ff(1,1:end-1+const);
            Z_f  = kron(eye(m),x_f);
            y_f  = (Z_f*a_t_fore)';
            y_t_fore(llags,:) = y_f;
        end;
        % now recompute forecast in levels with the right scale
        y_levfore        = NaN(hmax+1,NN_VAR);
        y_levfore(1,:)   = data_lev(irep,1:NN_VAR);
        % we are forecasting from irep+1 to irep+hmax
        y_levactual      = data_lev(irep+1:irep+hmax,1:NN_VAR);
        % give forecasts the right mean and standard deviation
        y_size           = y_t_fore(:,1:NN_VAR).*repmat(std_target(1:NN_VAR),hmax,1)+repmat(med_target(1:NN_VAR),hmax,1);
        for kk=2:hmax+1
            y_levfore(kk,:)  = y_levfore(kk-1,:)+y_size(kk-1,:);
        end
        FORGET.errors(count_tt,:,:)    =  y_levfore(2:end,1:NN_VAR)-y_levactual;
        FORGET.forecast(count_tt,:,:)  =  y_levfore(2:end,1:NN_VAR);
    end
%     toc
end 
lambda_t = lambda;