function [epi_em_pp,pi_em_pp, i_em_pp, r_em_pp, ygap_em_pp, zlb_freq_pp,gap_em_pp, avg_length_zlb,u_em_pp,sd_u_em_pp,sd_pi,sd_y,loss,pi_rho,sd_epi] = get_means(constants)

load(['workspace/soln_simul_' make_spec_name(constants) '.mat'])

i_em_pp = x_em(i_idx)*400 + (1-bbeta)/bbeta*400;
pi_em_pp = x_em(pi_idx)*400;
epi_em_pp = s_em(psilag_idx)*400;
r_em_pp = i_em_pp - epi_em_pp;
ygap_em_pp = x_em(ygap_idx)*100;
u_em_pp = s_em(u_idx)*100;
%epi_em_pp = s_em(psilag_idx)*400;
% frequency of ZLB binding
i_path = squeeze(x_path(:,i_idx,:));
zlb_binding = (i_path*400 + (1-bbeta)/bbeta*400) <= -0.0;
zlb_freq_pp = sum(zlb_binding)/length(i_path)*100;
infla = squeeze(x_path(:,pi_idx,:));
ygapp = squeeze(x_path(:,ygap_idx,:));
piexp = squeeze(s_path(:,psilag_idx,:));
ferr = abs(infla(2:end)*400-piexp(1:end-1)*400);
frr_m = mean(ferr);
lloss2 = mean(infla.^2+Lambda_y*ygapp.^2);

corr(ferr(1:end-1),ferr(2:end));
pi_rho = corr(infla(1:end-1),infla(2:end));
%pi_rho = corr(piexp(1:end-1),piexp(2:end));
sd_pi = std(infla)*400;
sd_epi = std(piexp)*400;
sd_y = std(ygapp)*100;
Lambda_y    = 0.007;
losses = zeros(size(ygapp));
for j=1:size(ygapp,1)

        losses(j) = losses(j) - bbeta^(j-1)*( Lambda_y*(100*ygapp(j))^2 + (100*infla(j))^2 );
end
loss = mean(losses);

% average length of ZLB period 
at_zlb = 0;
no_zlb_periods = 0;
for jj=1:length(zlb_binding)
    if at_zlb == 0 && zlb_binding(jj)==1 % if we ENTER a zlb period... 
        no_zlb_periods = no_zlb_periods + 1; % start counting number of periods. 
    end
    if zlb_binding(jj)==0
        at_zlb = 0;
    end
end
avg_length_zlb = sum(zlb_binding)/no_zlb_periods;
uu = squeeze(s_path(:,u_idx,:));
sd_u_em_pp = std(uu);
% compute housing gap 
switch flag_re
    case 0 
        logxi           = squeeze(s_path(:,logxi_idx,:));
        logqu           = squeeze(s_path(:,logqu_idx,:));
        bbetat          = max(0, min(bbetaU, (exp(logqu) - exp(logxi))./exp(logqu) / (bbeta * (1-delta))));

        quhatgap        = ( (1-bbeta*(1-delta))./(1-bbeta*(1-delta)*bbetat) - (1-bbeta*(1-delta))/(1-bbeta*(1-delta)*rho_xi) ).*logxi ...
                            + (bbeta*(1-delta)*(bbetat-1))./(1-bbeta*(1-delta)*bbetat) + s_d;
        gap             = CqCy*(s_d - (quhatgap));
        gap_em_pp       = mean(gap)*100;

    case 1
        gap_em_pp       = 0;
end

if flag_re==0;
check = (1-bbeta)*pi_em_pp/4 - ( kappa_y * ygap_em_pp )
end
