function [RMSE,biais,J,ErreursPrev] = PrevBVAR(nlags)
% Adapted from bvar_forecast.m. 
% Here we compute the RMSE associated to in sample forecasts k step ahead
% (k=1,...,K) from the posterior mode of the bvar model. The uncertainty
% related to the autoregressive matrices and the size of the innovations,
% and the uncertainty related to the future shocks are not considered in
% this version.    
    
    global options_
    
    % Load dataset
    dataset = read_variables(options_.datafile, options_.varobs, [], options_.xls_sheet, options_.xls_range);
    options_ = set_default_option(options_, 'nobs', size(dataset,1)-options_.first_obs+1);
    
    train = options_.bvar_prior_train;

    idx = options_.first_obs+options_.presample-train-nlags:options_.first_obs+options_.nobs-1;
    idx1 = options_.first_obs+options_.presample-train:options_.first_obs+options_.nobs-1;
    centrage=mean(dataset(idx1,:));

    % Prepare dataset
    if options_.loglinear & ~options_.logdata
        dataset = log(dataset);
    end
    if options_.prefilter
        dataset(idx,:) = dataset(idx,:) - ones(length(idx),1)*centrage;
    end
    
    
    [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);

    ydata = dataset(idx,:);
    T = size(ydata, 1);
    xdata = ones(T,nx);
    k = ny*nlags+nx;
    nb_var=size(dataset,2);
    initrow = nlags+1;
    horizon = 8;
    start = 5;
    termi = 51;
    nb_prevs = termi-start+1;
    
    i1 = initrow+start;
    i2 = initrow+start+horizon-1;
    
    lags_data = ydata(initrow+1:initrow+start-1,:);
    true_data = ydata(i1:i2,:);
    
    erro_data = zeros(ny,horizon,nb_prevs);
    tic;
    J = zeros(nb_prevs,horizon+1,7); % initialisation de la matrice renvoyant les prevs du bvar
    for s=1:nb_prevs
        true_data = ydata(i1:i2,:);
        % boucle pour compléter la première colonne de J avec les derniers
        % points des données avant prev (h=0) - A CORRIGER
        J(s,1,:)=lags_data(end,:); %dataset(end-47+s,variable);
        % -------------------------------------------------
        for t = 1:horizon
            X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
            y = X * posterior.PhiHat;
            J(s,t+1,:)=y; % remplissage de la matrice J avec les prevs
            lags_data(1:end-1,:) = lags_data(2:end, :);
            lags_data(end,:) = y;
            erro_data(:,t,s) = transpose(true_data(t,:)-y);
        end
        lags_data = ydata(initrow+1+s:initrow+start-1+s,:); % réinitialisation de lags_data
        i1 = i1+1;
        i2 = i2+1;
    end
    toc
    
    ErreursPrev = permute(erro_data,[3 2 1]);

       
% RMSE et biais par horizon de prev

RMSE = zeros(ny,horizon);
biais = zeros(ny,horizon);
for i = 1:ny
    RMSE(i,:) = sqrt(diag(ErreursPrev(:,:,i)'*ErreursPrev(:,:,i))/nb_prevs);
    biais(i,:) = mean(ErreursPrev(:,:,i));
end

% recentrage des prévisions dans la matrice J pour les rendre comparables
% aux données initiales

for variable=1:nb_var
    J(:,:,variable)=J(:,:,variable)+centrage(variable);
end