function simul(RNAT,flag_zlb,flag_AB)

% %% model settings
% flag_zlb = 1
% flag_AB  = 1
% RNAT = 1 % determines annual natural rate of interest and thereby strength of ZLB

%% simulation settings
n_per_av  = 1e+5;             % simulated periods for average (throw away first thousand)
nsim_per_irf = 1e+5;
nburnin=1000;
nsdu = 3; % number of conditional sd magnitues for u shock
nsdr = 3; % number of cpnditional sd magnitues for r shock

%% load data

if flag_zlb == 1 
    if flag_AB == 1
        load(['recent_zlb_' num2str(RNAT) '_AB.mat'])
    else
        load(['recent_zlb_' num2str(RNAT) '.mat'])
    end
elseif flag_zlb == 0  
    if flag_AB == 1
        load(['recent_nozlb_' num2str(RNAT) '_AB.mat'])
    else
        load(['recent_nozlb_' num2str(RNAT) '.mat'])
    end
end 

%% SIMULATE ERGODIC MEAN
tic 
disp('SIMULATE AVERAGE (OVER TIME):')
s_simul_init = s_init;                                          % initialize states
[s_path,x_path] = dpsimul(model,s_simul_init,n_per_av,s_out,x_opt); % simulate one path
% [~,x_pathSSS] = dpsimul2(model,s_simul_init,nburnin+100,s_out,x_opt); % simulate one path
% x_pathSSS = x_pathSSS(:,:,nburnin:end);                                      % throw away first thousand periods
x_path = x_path(:,:,nburnin:end);                                      % throw away first thousand periods
% s_path = s_path(:,:,nburnin:end);
for j=1:5
    xpathav(j) = mean(x_path(:,j,:));                % policies: average over time
%     xpathavSSS(j) = mean(x_pathSSS(:,j,:));                % policies: average over time
%     xpathsd(j) =  std(x_path(:,j,:));                %           standard deviation
end 

for j=1:4
    spathav(j) = mean(s_path(:,j,:));                % states:   average over time
%     spathsd(j) =  std(s_path(:,j,:));                %           standard deviation
end
xpathav
fprintf('Time = %7.2f Seconds\n',toc) 

%% SIMULATE IMPULSE RESPONSE 
tic

disp(' ')
disp('SIMULATE MP: positive u shock')  
junk = [(nsdu*sig_u+spathav(1))*ones(nsim_per_irf,1) ...
        spathav(2)*ones(nsim_per_irf,1) ...
        spathav(3)*ones(nsim_per_irf,1) ... 
        spathav(4)*ones(nsim_per_irf,1)];
    
[spath_u_pos,xpath_u_pos] = dpsimul(model,junk,16,s_out,x_opt); 
fprintf('Time = %7.2f Seconds\n',toc) 

disp(' ')
disp('SIMULATE MP: negative u shock')  
junk = [(-nsdu*sig_u+spathav(1))*ones(nsim_per_irf,1) ...
        spathav(2)*ones(nsim_per_irf,1) ...
        spathav(3)*ones(nsim_per_irf,1) ... 
        spathav(4)*ones(nsim_per_irf,1)];
    
[spath_u_neg,xpath_u_neg] = dpsimul(model,junk,16,s_out,x_opt); 
fprintf('Time = %7.2f Seconds\n',toc) 

disp(' ')
disp('SIMULATE IRF: negative rn shock')  
junk = [spathav(1)*ones(nsim_per_irf,1) ...
        spathav(2)*ones(nsim_per_irf,1) ...
        (-nsdr*sig_rn+spathav(3))*ones(nsim_per_irf,1) ... 
        spathav(4)*ones(nsim_per_irf,1)];
    
[spath_rn_neg,xpath_rn_neg] = dpsimul(model,junk,16,s_out,x_opt); 
fprintf('Time = %7.2f Seconds\n',toc) 

disp(' ')
disp('SIMULATE IRF: positive rn shock')  
junk = [spathav(1)*ones(nsim_per_irf,1) ...
        spathav(2)*ones(nsim_per_irf,1) ...
        (nsdr*sig_rn+spathav(3))*ones(nsim_per_irf,1) ... 
        spathav(4)*ones(nsim_per_irf,1)];
    
[spath_rn_pos,xpath_rn_pos] = dpsimul(model,junk,16,s_out,x_opt); 
fprintf('Time = %7.2f Seconds\n',toc) 

% %% 
% jj=3; 
% 
% % aux = squeeze(mean(xpath_u_pos,1))'; 
% aux = squeeze(mean(xpath_rn_neg,1))'; 
% 
% fig = figure;
% set(fig,'DefaultAxesFontSize',20); 
% set(fig,'DefaultLineLineWidth',5);  
% set(gcf, 'unit', 'inches'); 
% 
% hold on
% % plot([0:16],aux(:,jj) - xpathav(jj))
% plot([0:16],aux(:,jj)*4)
% plot([0:16],zeros(17,1),'k-','LineWidth',1)
% hold off
% xticks([0:4:16])
% xlim([0 16])
% grid on, box on

%% save data with simulations

if flag_zlb == 1 
    if flag_AB == 1
        save(['../../../NKHousing_LB_own/v1_RE_AB_results/recent_zlb_' num2str(RNAT) '_AB_simul.mat'])
    else
        save(['../../../NKHousing_LB_own/v1_RE_AB_results/recent_zlb_' num2str(RNAT) '_simul.mat'])
    end
elseif flag_zlb == 0  
    if flag_AB == 1
        save(['../../../NKHousing_LB_own/v1_RE_AB_results/recent_nozlb_' num2str(RNAT) '_AB_simul.mat'])
    else
        save(['../../../NKHousing_LB_own/v1_RE_AB_results/recent_nozlb_' num2str(RNAT) '_simul.mat'])
    end
end 
