function dsgevar_forecast(nlags, presample, replic, kconst)

% function dsgebvar_forecast(nlags)
% builds forecats for a dsgebvar model
% based on the posterior dsgebvar
% PHI, SIGMAu, iXX & df estimations 
%
% Date: 04.18.2016
%
%
% INPUTS
%    nlags:     number of lags for the bvar
%    presample: number of data to pre-estimate    
%    replic:    number of replics on BVAR forecasts   
%    kconst:    1-noconstant; 0-constant   
%
% OUTPUTS
%    none
%
% SPECIAL REQUIREMENTS
%    none

global options_ oo_ M_ estim_params_ bayestopt_
%global DynareDataset DynareOptions EstimatedParameters Model BayesInfo DynareResults deep

if options_.forecast == 0
    error('dsgevar_forecast: you must specify "forecast" option')
end

%DynareDataset = initialize_dataset1(options_.datafile,options_.varobs,1,options_.nobs,[],[],options_);
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);
DynareOptions = options_;
%EstimatedParameters = estim_params_;
%Model = M_;
%[xparam1, estim_params_, bayestopt_, lb, ub, M_] = set_prior(estim_params_,M_,options_);
%BayesInfo = bayestopt_; % a-priori conditions
%DynareResults = oo_;
%deep = M_.params;
%deep = xparam1;  % xparam1 initial values

if (options_.first_obs+options_.nobs-1)> size(dataset,1)
    fprintf('Incorrect or missing specification of the number of observations. nobs can be at most %4u\n',size(dataset,1)-options_.first_obs+1);
    error('Inconsistent number of observations.') 
end

% Parameters for prior
%if options_.first_obs + options_.presample <= nlags
%    error('first_obs+presample should be > nlags (for initializing the VAR)')
%end

%train = options_.bvar_prior_train;
train = options_.first_obs;
options_.presample = presample;
options_.bvar_replic = replic;
options_.noconstant = kconst;

%if options_.first_obs + options_.presample - train <= nlags
%    error('first_obs+presample-train should be > nlags (for initializating the VAR)')
%end

idx = options_.first_obs+options_.presample-train-nlags:options_.first_obs+options_.nobs-1;

% Prepare dataset
if options_.loglinear && ~options_.logdata
    dataset = log(dataset);
end
if options_.prefilter
    dataset(idx,:) = dataset(idx,:) - ones(length(idx),1)*mean(dataset(idx,:));
end

ny = size(dataset,2);
%ny = DynareDataset.info.nvobs;

if options_.prefilter || options_.noconstant
 nx = 0;
else
 nx = 1;
end

ydata = dataset(idx, :);
T = size(ydata, 1);
xdata = ones(T,nx);

% Add forecast informations
%if nargout >= 5
    forecast_data.xdata = ones(options_.forecast, nx);
    forecast_data.initval = ydata(end-nlags+1:end, :);
    if options_.first_obs + options_.nobs <= size(dataset, 1)
        forecast_data.realized_val = dataset(options_.first_obs+options_.nobs:end, :);
        forecast_data.realized_xdata = ones(size(forecast_data.realized_val, 1), nx);
    else
        forecast_data.realized_val = [];
    end
%end

%[ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
%[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] =  dsge_var_likelihood1(deep,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
%bvar = dsgevar_posterior_density1(deep,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JCSANTANAC 04.17.2016
% Once the estimation was executed
% we call the PHI, SIGMAu, iXX estimations
%

gend = DynareOptions.nobs;
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));

NumberOfLags = options_.dsge_varlag;
NumberOfVariables = size(options_.varobs,1);
Constant = 'no';
NumberOfEstimatedParameters = NumberOfLags*NumberOfVariables;
if ~options_.noconstant
    Constant = 'yes';
    NumberOfEstimatedParameters = NumberOfEstimatedParameters + ...
    NumberOfVariables;
end


PHI = evalin('base','PHI');
iXX = evalin('base','iXX');
SIGMAu = evalin('base','SIGMAu');
Sigma = iXX;                    % Covariance Matrix Normal MV distribution
S = DSGE_PRIOR_WEIGHT*SIGMAu;   % Wishart Distribution (Mean)                   
df = DSGE_PRIOR_WEIGHT-NumberOfEstimatedParameters; % Wishart Distribution (Var)    

%%%%%%%%%%%%%%%%%%%%%%%

sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic);
sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);

%S_inv_upper_chol = chol(inv(posterior.S));
S_inv_upper_chol = chol(inv(S));

% Option 'lower' of chol() not available in old versions of
% Matlab, so using transpose

XXi_lower_chol = chol(iXX)';


%Num. Obs. a Simular
k = ny*nlags+nx;

% Declaration of the companion matrix
Companion_matrix = diag(ones(ny*(nlags-1),1),-ny);


% Number of explosive VAR models sampled
p = 0;
% Loop counter initialization
d = 0;
while d <= options_.bvar_replic
    
    Sigma = rand_inverse_wishart(ny, df, S_inv_upper_chol);

    % Option 'lower' of chol() not available in old versions of
    % Matlab, so using transpose

    Sigma_lower_chol = chol(Sigma)';

    Phi = rand_matrix_normal(k, ny, PHI, Sigma_lower_chol, XXi_lower_chol);
    
    % All the eigenvalues of the companion matrix have to be on or inside the unit circle
    Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)'; 
    test = (abs(eig(Companion_matrix)));
    if any(test>1.0000000000001)
        p = p+1;
    end
    d = d+1;
    
    % Without shocks
    lags_data = forecast_data.initval;
    for t = 1:options_.forecast
        X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
        y = X * Phi;
        lags_data(1:end-1,:) = lags_data(2:end, :);
        lags_data(end,:) = y;
        sims_no_shock(t, :, d) = y;
    end
    
    % With shocks
    lags_data = forecast_data.initval;
    for t = 1:options_.forecast
        X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
        shock = (Sigma_lower_chol * randn(ny, 1))';
        y = X * Phi + shock;
        lags_data(1:end-1,:) = lags_data(2:end, :);
        lags_data(end,:) = y;
        sims_with_shocks(t, :, d) = y;
    end
end

if p > 0
    skipline()
    disp(['Some of the DSGEBVAR models sampled from the posterior distribution'])
    disp(['were found to be explosive (' num2str(p/options_.bvar_replic) ' percent).'])
    skipline()
end

% Plot graphs
sims_no_shock_mean = mean(sims_no_shock, 3);

sort_idx = round((0.5 + [-options_.conf_sig, options_.conf_sig, 0]/2) * options_.bvar_replic);

sims_no_shock_sort = sort(sims_no_shock, 3);
sims_no_shock_down_conf = sims_no_shock_sort(:, :, sort_idx(1));
sims_no_shock_up_conf = sims_no_shock_sort(:, :, sort_idx(2));
sims_no_shock_median = sims_no_shock_sort(:, :, sort_idx(3));

sims_with_shocks_sort = sort(sims_with_shocks, 3);
sims_with_shocks_down_conf = sims_with_shocks_sort(:, :, sort_idx(1));
sims_with_shocks_up_conf = sims_with_shocks_sort(:, :, sort_idx(2));

OutputDirectoryName = CheckPath('graphs',M_.fname);

dyn_graph=dynare_graph_init(sprintf('DSGEBVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});

for i = 1:ny
    dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
                   sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
                   sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
                 options_.varobs(i, :));
end

dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_DSGEBVAR_forecast_',num2str(nlags)],options_)

% Compute RMSE

if ~isempty(forecast_data.realized_val)
    
    sq_err_cumul = zeros(1, ny);
    
    lags_data = forecast_data.initval;
    for t = 1:size(forecast_data.realized_val, 1)
        X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ];
        y = X * PHI;
        lags_data(1:end-1,:) = lags_data(2:end, :);
        lags_data(end,:) = y;
        sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2;
    end
    
    rmse = sqrt(sq_err_cumul / size(forecast_data.realized_val, 1));
    
    fprintf('RMSE of DSGEBVAR(%d):\n', nlags);
    
    for i = 1:size(options_.varobs, 1)
        fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i));
    end 
end

% Store results

DirectoryName = [ M_.fname '/dsgebvar_forecast' ];
if ~isdir(DirectoryName)
    if ~isdir(M_.fname)
        mkdir(M_.fname);
    end
    mkdir(DirectoryName);
end
save([ DirectoryName '/simulations.mat'], 'sims_no_shock', 'sims_with_shocks');

for i = 1:size(options_.varobs, 1)
    name = options_.varobs(i, :);

    sims = squeeze(sims_with_shocks(:,i,:));
    eval(['oo_.dsgebvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']);
    eval(['oo_.dsgebvar.forecast.with_shocks.Median.' name ' = median(sims, 2);']);
    eval(['oo_.dsgebvar.forecast.with_shocks.Var.' name ' = var(sims, 0, 2);']);
    eval(['oo_.dsgebvar.forecast.with_shocks.HPDsup.' name ' = sims_with_shocks_up_conf(:,i);']);
    eval(['oo_.dsgebvar.forecast.with_shocks.HPDinf.' name ' = sims_with_shocks_down_conf(:,i);']);

    sims = squeeze(sims_no_shock(:,i,:));
    eval(['oo_.dsgebvar.forecast.no_shock.Mean.' name ' = sims_no_shock_mean(:,i);']);
    eval(['oo_.dsgebvar.forecast.no_shock.Median.' name ' = sims_no_shock_median(:,i);']);
    eval(['oo_.dsgebvar.forecast.no_shock.Var.' name ' = var(sims, 0, 2);']);
    eval(['oo_.dsgebvar.forecast.no_shock.HPDsup.' name ' = sims_no_shock_up_conf(:,i);']);
    eval(['oo_.dsgebvar.forecast.no_shock.HPDinf.' name ' = sims_no_shock_down_conf(:,i);']);

    if exist('rmse')
        eval(['oo_.dsgebvar.forecast.rmse.' name ' = rmse(i);']);
    end
end