function info=stoch_simul(var_list)

% Copyright (C) 2001-2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

set(0,'defaultlinelinewidth',3)

global M_ options_ oo_ it_ modelnumber

test_for_deep_parameters_calibration(M_);

dr = oo_.dr;

options_old = options_;
if options_.linear
    options_.order = 1;
end
if options_.order == 1
    options_.replic = 1;
elseif options_.order == 3
    options_.k_order_solver = 1;
end

if isempty(options_.qz_criterium)
    options_.qz_criterium = 1+1e-6;
end

if options_.partial_information == 1 || options_.ACES_solver == 1
    PI_PCL_solver = 1;
    if options_.order ~= 1
        warning('STOCH_SIMUL: forcing order=1 since you are using partial_information or ACES solver')
        options_.order = 1;
    end
else
    PI_PCL_solver = 0;
end

TeX = options_.TeX;

if size(var_list,1) == 0
    var_list = M_.endo_names(1:M_.orig_endo_nbr, :);
end

[i_var,nvar] = varlist_indices(var_list,M_.endo_names);

iter_ = max(options_.periods,1);
if M_.exo_nbr > 0
    oo_.exo_simul= ones(iter_ + M_.maximum_lag + M_.maximum_lead,1) * oo_.exo_steady_state';
end

check_model;

oo_.dr=set_state_space(dr,M_);

if PI_PCL_solver
    [oo_.dr, info] = PCL_resol(oo_.steady_state,0);
elseif options_.discretionary_policy
    if ~options_.linear
        error(['discretionary_policy solves only linear_quadratic ' ...
            'problems']);
    end
    [oo_.dr,ys,info] = discretionary_policy_1(oo_,options_.instruments);
else
    [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
end

if info(1)
    options_ = options_old;
    print_info(info, options_.noprint);
    return
end

if ~options_.noprint
    disp(' ')
    disp('MODEL SUMMARY')
    disp(' ')
    disp(['  Number of variables:         ' int2str(M_.endo_nbr)])
    disp(['  Number of stochastic shocks: ' int2str(M_.exo_nbr)])
    if (options_.block)
        disp(['  Number of state variables:   ' int2str(oo_.dr.npred+oo_.dr.nboth)])
        disp(['  Number of jumpers:           ' int2str(oo_.dr.nfwrd+oo_.dr.nboth)])
    else
        disp(['  Number of state variables:   ' ...
            int2str(length(find(oo_.dr.kstate(:,2) <= M_.maximum_lag+1)))])
        disp(['  Number of jumpers:           ' ...
            int2str(length(find(oo_.dr.kstate(:,2) == M_.maximum_lag+2)))])
    end;
    disp(['  Number of static variables:  ' int2str(oo_.dr.nstatic)])
    my_title='MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS';
    labels = deblank(M_.exo_names);
    headers = char('Variables',labels);
    lh = size(labels,2)+2;
    dyntable(my_title,headers,labels,M_.Sigma_e,lh,10,6);
    if options_.partial_information
        disp(' ')
        disp('SOLUTION UNDER PARTIAL INFORMATION')
        disp(' ')
        
        if isfield(options_,'varobs')&& ~isempty(options_.varobs)
            PCL_varobs=options_.varobs;
            disp('OBSERVED VARIABLES')
        else
            PCL_varobs=M_.endo_names;
            disp(' VAROBS LIST NOT SPECIFIED')
            disp(' ASSUMED OBSERVED VARIABLES')
        end
        for i=1:size(PCL_varobs,1)
            disp(['    ' PCL_varobs(i,:)])
        end
    end
    disp(' ')
    if options_.order <= 2 && ~PI_PCL_solver
        disp_dr(oo_.dr,options_.order,var_list);
    end
end

if options_.periods > 0 && ~PI_PCL_solver
    if options_.periods <= options_.drop
        disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ...
            ' than the number of observations to be dropped'])
        options_ =options_old;
        return
    end
    if isempty(M_.endo_histval)
        y0 = oo_.dr.ys;
    else
        y0 = M_.endo_histval;
    end
    oo_.endo_simul = simult(y0,oo_.dr);
    dyn2vec;
end

if options_.nomoments == 0
    if PI_PCL_solver
        PCL_Part_info_moments (0, PCL_varobs, oo_.dr, i_var);
    elseif options_.periods == 0
        % There is no code for theoretical moments at 3rd order
        if options_.order <= 2
            disp_th_moments(oo_.dr,var_list);
        end
    else
        disp_moments(oo_.endo_simul,var_list);
    end
end

nr1 = M_.exo_nbr;
nc1 = nvar;

if options_.irf
    var_listTeX = M_.endo_names_tex(i_var,:);
    
    SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
    cs = transpose(chol(SS));
    tit(M_.exo_names_orig_ord,:) = M_.exo_names;
    irf_shocks_indx = getIrfShocksIndx();
    
    
    for i=irf_shocks_indx
        if PI_PCL_solver
            y=PCL_Part_info_irf (0, PCL_varobs, i_var, M_, oo_.dr, options_.irf, i);
        else
            y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
                options_.replic, options_.order);
        end
        irfs   = [];
        mylist = [];
        for j = 1:nvar
            assignin('base',[deblank(M_.endo_names(i_var(j),:)) '_' deblank(M_.exo_names(i,:))],...
                y(i_var(j),:)');
            eval(['oo_.irfs.' deblank(M_.endo_names(i_var(j),:)) '_' ...
                deblank(M_.exo_names(i,:)) ' = y(i_var(j),:);']);
            irfs  = cat(1,irfs,y(i_var(j),:));
            if isempty(mylist)
                mylist = deblank(var_list(j,:));
            else
                mylist = char(mylist,deblank(var_list(j,:)));
            end
        end
        
        irfsmat(:,:,i)=irfs;
        
    end
end

colore='b';
lss='-';

if modelnumber<3
    
    if modelnumber==1; colore='b'; end
    if modelnumber==2; colore='r'; end
    
    for ir=1:nr1
        for ic = 1:nc1
            
            subplot(nr1,nc1,nc1*(ir-1)+ic);
            plot(irfsmat(ic,:,ir),'color',colore,'linewidth',3); hold on
            title(deblank(mylist(ic,:)),'Interpreter','none');
            ylabel(deblank(M_.exo_names(ir,:)))
        end
    end
    
else
   
    
    if modelnumber==11; colore='b'; lss='-'; end
    if modelnumber==12; colore='r'; lss='--'; end
    if modelnumber==13; colore='k'; lss='.'; end
    
    for ir=2:2
        for ic = 1:nc1
            
            subplot(3,2,ic);
            plot(irfsmat(ic,:,ir),'color',colore,'linestyle',lss,'linewidth',3); hold on
            title(deblank(mylist(ic,:)),'Interpreter','none');
            ylabel(deblank(M_.exo_names(ir,:)))
        end
    end
    
end