1 function [RMSE,biais,J,ErreursPrev] = PrevBVAR(nlags)
   2 % Adapted from bvar_forecast.m. 
   3 % Here we compute the RMSE associated to in sample forecasts k step ahead
   4 % (k=1,...,K) from the posterior mode of the bvar model. The uncertainty
   5 % related to the autoregressive matrices and the size of the innovations,
   6 % and the uncertainty related to the future shocks are not considered in
   7 % this version.    
   8     
   9     global options_
  10     
  11     % Load dataset
  12     dataset = read_variables(options_.datafile, options_.varobs, [], options_.xls_sheet, options_.xls_range);
  13     options_ = set_default_option(options_, 'nobs', size(dataset,1)-options_.first_obs+1);
  14     
  15     train = options_.bvar_prior_train;
  16 
  17     idx = options_.first_obs+options_.presample-train-nlags:options_.first_obs+options_.nobs-1;
  18     idx1 = options_.first_obs+options_.presample-train:options_.first_obs+options_.nobs-1;
  19     centrage=mean(dataset(idx1,:));
  20 
  21     % Prepare dataset
  22     if options_.loglinear & ~options_.logdata
  23         dataset = log(dataset);
  24     end
  25     if options_.prefilter
  26         dataset(idx,:) = dataset(idx,:) - ones(length(idx),1)*centrage;
  27     end
  28     
  29     
  30     [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
  31 
  32     ydata = dataset(idx,:);
  33     T = size(ydata, 1);
  34     xdata = ones(T,nx);
  35     k = ny*nlags+nx;
  36     nb_var=size(dataset,2);
  37     initrow = nlags+1;
  38     horizon = 8;
  39     start = 5;
  40     termi = 51;
  41     nb_prevs = termi-start+1;
  42     
  43     i1 = initrow+start;
  44     i2 = initrow+start+horizon-1;
  45     
  46     lags_data = ydata(initrow+1:initrow+start-1,:);
  47     true_data = ydata(i1:i2,:);
  48     
  49     erro_data = zeros(ny,horizon,nb_prevs);
  50     tic;
  51     J = zeros(nb_prevs,horizon+1,7); % initialisation de la matrice renvoyant les prevs du bvar
  52     for s=1:nb_prevs
  53         true_data = ydata(i1:i2,:);
  54         % boucle pour compléter la première colonne de J avec les derniers
  55         % points des données avant prev (h=0) - A CORRIGER
  56         J(s,1,:)=lags_data(end,:); %dataset(end-47+s,variable);
  57         % -------------------------------------------------
  58         for t = 1:horizon
  59             X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
  60             y = X * posterior.PhiHat;
  61             J(s,t+1,:)=y; % remplissage de la matrice J avec les prevs
  62             lags_data(1:end-1,:) = lags_data(2:end, :);
  63             lags_data(end,:) = y;
  64             erro_data(:,t,s) = transpose(true_data(t,:)-y);
  65         end
  66         lags_data = ydata(initrow+1+s:initrow+start-1+s,:); % réinitialisation de lags_data
  67         i1 = i1+1;
  68         i2 = i2+1;
  69     end
  70     toc
  71     
  72     ErreursPrev = permute(erro_data,[3 2 1]);
  73 
  74        
  75 % RMSE et biais par horizon de prev
  76 
  77 RMSE = zeros(ny,horizon);
  78 biais = zeros(ny,horizon);
  79 for i = 1:ny
  80     RMSE(i,:) = sqrt(diag(ErreursPrev(:,:,i)'*ErreursPrev(:,:,i))/nb_prevs);
  81     biais(i,:) = mean(ErreursPrev(:,:,i));
  82 end
  83 
  84 % recentrage des prévisions dans la matrice J pour les rendre comparables
  85 % aux données initiales
  86 
  87 for variable=1:nb_var
  88     J(:,:,variable)=J(:,:,variable)+centrage(variable);
  89 end