function [wf,wf_alt] = simulate(constants)

%% set seed
rng('default')

%% get input constants
rnat = constants.rnat; 
flag_zlb = constants.flag_zlb;
flag_re = constants.flag_re;
gamma = constants.gamma;
bbeta   = 1/(1+rnat/400);
gamma_y = constants.gamma_y;
%% get indices
define_idx
 
%% Load results
if gamma == 10
    load(['recent_zlb_1_AB.mat'])
else
    load(['workspace/soln_' make_spec_name(constants) '.mat'])
end
% load(['workspace/soln_simul_' make_spec_name(constants) '.mat'])

%% simulation settings
n_per_av = 100000; % total simulated periods for average 
if flag_re == 1; n_per_av = 10000; end;
nburnin = 1000; % burn-in periods

%% simulate ergodic mean
tic 
disp('Simulate average (over time), i.e., ergodic mean:')  
[s_path,x_path] = mydpsimul(model,s_init,n_per_av,s_out,x_opt,false); % simulate one path 

try
    assert(size(x_path,3)-nburnin > 0)
    x_path = x_path(:,:,nburnin:end); 
    s_path = s_path(:,:,nburnin:end);   
catch err
end
%%%%%%%%%%%%%%%%
xpathy = squeeze(x_path(:,1,:));                             % squeeze to vector
xpathp = squeeze(x_path(:,2,:));                             %   "

losses = zeros(1,size(xpathy,2));
for j=1:size(xpathy,2)
    for i=1:size(xpathy,1)
        losses(j) = losses(j) - bbeta^(i-1)*( Lambda_y*xpathy(i,j)^2 + xpathp(i,j)^2 );
    end
end
wf=  mean(losses)*20
wf_alternative_aux = -(Lambda_y*(xpathy).^2+(xpathp).^2);
wf_alt = mean(wf_alternative_aux);
wf_alt = mean(wf_alt);
%%%%%%%%%%%%%


for j=1:size(x_path,2)
    x_em(j) = mean(x_path(:,j,:)); 
end 

for j=1:size(s_path,2)
    s_em(j) = mean(s_path(:,j,:));   
end 
fprintf('Time = %7.2f Seconds\n',toc) 

try
    if abs(s_em(logxi_idx)) > 0.005; warning('Warning: xi shocks have not averaged out. Generate a longer simulation sample.'); end;
catch err
end

%% In particular, save ergodic mean of ygap
ygap_em = x_em(ygap_idx);
assert(~isnan(ygap_em),'Error: Simulated ygap is NaN.')


%% simulate stochastic steady state 
[s_path_noshocks,x_path_noshocks] = mydpsimul(model,s_init,1e+3,s_out,x_opt,true); 
s_sss = squeeze(s_path_noshocks(1,:,end));
x_sss = squeeze(x_path_noshocks(1,:,end));
if flag_re == 0
    assert(abs(s_sss(logqu_idx) - logqu_ss) < 1e-5); 
end
%% map into percentage points (annualized if applicable) 
x_em_pp = x_em;
x_em_pp(i_idx) = (x_em_pp(i_idx) + (1-bbeta)/bbeta)*4;
x_em_pp(pi_idx) = x_em_pp(pi_idx)*4;
x_em_pp = x_em_pp*100;

s_em_pp = s_em;
if flag_re == 0
    s_em_pp(logqu_idx) = s_em_pp(logqu_idx) - logqu_ss;
    s_em_pp(loglagqu_idx) = s_em_pp(loglagqu_idx) - logqu_ss;
end
s_em_pp = s_em_pp*100;

x_sss_pp = x_sss;
x_sss_pp(i_idx) = x_sss_pp(i_idx)*4 + (1-bbeta)/bbeta*4;
x_sss_pp(pi_idx) = x_sss_pp(pi_idx)*4;
x_sss_pp = x_sss_pp*100;

s_sss_pp = s_sss;
if flag_re == 0
    s_sss_pp(logqu_idx) = s_sss_pp(logqu_idx) - logqu_ss;
    s_sss_pp(loglagqu_idx) = s_sss_pp(loglagqu_idx) - logqu_ss;
end
s_sss_pp = s_sss_pp*100;

%% Save results
save(['workspace/soln_simul_' make_spec_name(constants) '.mat'])

 