%----------------------------------------------------------
% Generate simulated data using a simplified version of Smets and Wouters (2007) model
% with zero lower bound on the nominal interest rate
%----------------------------------------------------------
clear all;
randn('state', 777);

% Sample size and other global variables
nobs    = 7;		        % number of observed variables
simsmp  = 200;			% sample size in simulation
simdrop = 100;			% number of simulations to be dropped
simtt = simdrop + simsmp;	% total number of simulations
k = 7;				% number of structural shocks (excluding binding shocks)

nrep = 200;           % number of replication

%----------------------------------------
%   Assign parameter values
%----------------------------------------
crhoa      = 0.5;
crhob      = 0.5;
crhog      = 0.5;
crhoqs     = 0.5;
crhoms     = 0.5;
crhopinf   = 0.5;
crhow      = 0.5;
csadjcost  = 4;
csigma     = 1.5;
chabb      = 0.7;
cprobw     = 0.5;
csigl      = 2;
cprobp     = 0.5;
cindw      = 0.5;
cindp      = 0.5;
czcap      = 0.5;
cfc        = 1.25;
crpi       = 1.5;
crr        = 0.75;
cry        = 0.125;
%constepinf = 0.625;
constepinf = 0.625;
constebeta = 0.25;
constelab  = 0;
ctrend     = 0.4;
calfa      = 0.3;

sig_a      = 0.45;
sig_b      = 0.45;
sig_g      = 0.45;
sig_qs     = 0.45;
sig_m      = 0.45;
sig_pinf   = 0.45;
sig_w      = 0.45;

cmap       = 0.5;
cmaw       = 0.5;
cgy        = 0.5;

noeval = 0;

% fixed parameters
ctou=0.025;
clandaw=1.5;
cg=0.18;
curvp=10;
curvw=10;

% auxiliary parameters
cpie=1+constepinf/100;
cgamma=1+ctrend/100 ;
cbeta=1/(1+constebeta/100);

clandap=cfc;
cbetabar=cbeta*cgamma^(-csigma);
cr=cpie/(cbeta*cgamma^(-csigma));
crk=(cbeta^(-1))*(cgamma^csigma) - (1-ctou);
cw = (calfa^calfa*(1-calfa)^(1-calfa)/(clandap*crk^calfa))^(1/(1-calfa));
cikbar=(1-(1-ctou)/cgamma);
cik=(1-(1-ctou)/cgamma)*cgamma;
clk=((1-calfa)/calfa)*(crk/cw);
cky=cfc*(clk)^(calfa-1);
ciy=cik*cky;
ccy=1-cg-cik*cky;
crkky=crk*cky;
cwhlc=(1/clandaw)*(1-calfa)/calfa*crk*cky/ccy;
cwly=1-crk*cky;

conster=(cr-1)*100;

%----------------------------------------
%  Matrices of canonical system (Sims' form)
%  GAM0*s(t) = GAM1*s(t-1) + PSI*e(t) + PPI*eta(t)
%----------------------------------------

% equation indices
eq_mc     = 1;
eq_zcap   = 2;
eq_rk     = 3;
eq_k      = 4;
eq_inve   = 5;
eq_pk     = 6;
eq_c      = 7;
eq_agg    = 8;
eq_prod   = 9;
eq_pinf   = 10;
eq_w      = 11;
eq_kp     = 12;
eq_r      = 13;
eq_a      = 14;
eq_b      = 15;
eq_g      = 16;
eq_qs     = 17;
eq_ms     = 18;
eq_spinf  = 19;
eq_sw     = 20;
eq_Einve  = 21;
eq_Erk    = 22;
eq_Epk    = 23;
eq_Ec     = 24;
eq_Elab   = 25;
eq_Epinf  = 26;
eq_Ew     = 27;
eq_m0    = 28;
eq_m1    = 29;
eq_m2    = 30;
eq_m3    = 31;
eq_m4    = 32;
eq_m5    = 33;
eq_m6    = 34;
eq_m7    = 35;
eq_m8    = 36;
eq_m9    = 37;
eq_m10   = 38;
eq_m11   = 39;
eq_m12   = 40;
eq_epinfma = 41;
eq_ewma  = 42;

% variable indices
v_mc     = 1;
v_zcap   = 2;
v_rk     = 3;
v_k      = 4;
v_inve   = 5;
v_pk     = 6;
v_c      = 7;
v_y      = 8;
v_lab    = 9;
v_pinf   = 10;
v_w      = 11;
v_kp     = 12;
v_r      = 13;
v_a      = 14;
v_b      = 15;
v_g      = 16;
v_qs     = 17;
v_ms     = 18;
v_spinf  = 19;
v_sw     = 20;
v_Einve  = 21;
v_Erk    = 22;
v_Epk    = 23;
v_Ec     = 24;
v_Elab   = 25;
v_Epinf  = 26;
v_Ew     = 27;
v_m0    = 28;
v_m1    = 29;
v_m2    = 30;
v_m3    = 31;
v_m4    = 32;
v_m5    = 33;
v_m6    = 34;
v_m7    = 35;
v_m8    = 36;
v_m9    = 37;
v_m10   = 38;
v_m11   = 39;
v_m12   = 40;
v_epinfma = 41;
v_ewma  = 42;

% shock indices
e_a  = 1;
e_b  = 2;
e_g  = 3;
e_qs = 4;
e_m  = 5;
e_pinf=6;
e_w  = 7;
e_m0 = 8;
e_m1 = 9;
e_m2 = 10;
e_m3 = 11;
e_m4 = 12;
e_m5 = 13;
e_m6 = 14;
e_m7 = 15;
e_m8 = 16;
e_m9 = 17;
e_m10= 18;
e_m11= 19;
e_m12= 20;

% expectation error indices
n_inve  = 1;
n_rk    = 2;
n_pk    = 3;
n_c     = 4;
n_lab   = 5;
n_pinf  = 6;
n_w     = 7;

% summary
neq  = 42;
neps = 20;
neta = 7;

% initialize matrices
GAM0 = zeros(neq,neq);
GAM1 = zeros(neq,neq);
C = zeros(neq,1);
PSI = zeros(neq,neps);
PPI = zeros(neq,neta);

%/**********************************************************
%**      1. Marginal cost 
%**********************************************************/
GAM0(eq_mc,v_mc)   =  -1;
GAM0(eq_mc,v_a)   =  -1;
GAM0(eq_mc,v_rk) =  calfa;
GAM0(eq_mc,v_w)  =  1-calfa;

%/**********************************************************
%**      2. Capacity Utilization 
%**********************************************************/
GAM0(eq_zcap,v_zcap)   = -1;
GAM0(eq_zcap,v_rk)  = 1/(czcap/(1-czcap));

%/**********************************************************
%**      3. Rental Rate of Capital 
%**********************************************************/
GAM0(eq_rk,v_rk)  = -1;
GAM0(eq_rk,v_w)   =  1;
GAM0(eq_rk,v_lab) =  1;
GAM0(eq_rk,v_k)   = -1;

%/**********************************************************
%**      4. Effective Capital Stock 
%**********************************************************/
GAM0(eq_k,v_k)  =  1;
GAM0(eq_k,v_zcap)  = -1;
GAM1(eq_k,v_kp)  =  1;

%/**********************************************************
%**      5. Investment 
%**********************************************************/
GAM0(eq_inve,v_inve)    = -1;
GAM1(eq_inve,v_inve)    = -1/(1+cbetabar*cgamma);
GAM0(eq_inve,v_Einve)   = (1/(1+cbetabar*cgamma))*cbetabar*cgamma;
GAM0(eq_inve,v_pk)   = (1/(1+cbetabar*cgamma))*(1/(cgamma^2*csadjcost));
GAM0(eq_inve,v_qs)   = 1;

%/**********************************************************
%**      6. Value of Capital 
%**********************************************************/
GAM0(eq_pk,v_pk)    = -1;
GAM0(eq_pk,v_r)    = -1;
GAM0(eq_pk,v_Epinf)    = 1;
GAM0(eq_pk,v_b)    = 1/((1-chabb/cgamma)/(csigma*(1+chabb/cgamma)));
GAM0(eq_pk,v_Erk)    = crk/(crk+(1-ctou));
GAM0(eq_pk,v_Epk)    = (1-ctou)/(crk+(1-ctou));

%/**********************************************************
%**      7. Consumption Euler 
%**********************************************************/
GAM0(eq_c,v_c)   = -1;
GAM1(eq_c,v_c)   = -(chabb/cgamma)/(1+chabb/cgamma);
GAM0(eq_c,v_Ec)  = 1/(1+chabb/cgamma);
GAM0(eq_c,v_lab) = (csigma-1)*cwhlc/(csigma*(1+chabb/cgamma));
GAM0(eq_c,v_Elab)= -(csigma-1)*cwhlc/(csigma*(1+chabb/cgamma));
GAM0(eq_c,v_r)  = - (1-chabb/cgamma)/(csigma*(1+chabb/cgamma));
GAM0(eq_c,v_Epinf)  =   (1-chabb/cgamma)/(csigma*(1+chabb/cgamma));
GAM0(eq_c,v_b)    = 1;

%/**********************************************************
%**      8. Aggregate Resource Constraint 
%**********************************************************/
GAM0(eq_agg,v_y) = -1;
GAM0(eq_agg,v_c) = ccy;
GAM0(eq_agg,v_inve) = ciy;
GAM0(eq_agg,v_g) = 1;
GAM0(eq_agg,v_zcap) = crkky;

%/**********************************************************
%**      9. Production Function 
%**********************************************************/
GAM0(eq_prod,v_y) = -1;
GAM0(eq_prod,v_k) = cfc*calfa;
GAM0(eq_prod,v_lab) = cfc*(1-calfa);
GAM0(eq_prod,v_a) = cfc;

%/**********************************************************
%**      10. Phillips Curve 
%**********************************************************/
GAM0(eq_pinf,v_pinf) = -1;
GAM0(eq_pinf,v_Epinf) = (1/(1+cbetabar*cgamma*cindp))*cbetabar*cgamma;
GAM1(eq_pinf,v_pinf) = -(1/(1+cbetabar*cgamma*cindp))*cindp;
GAM0(eq_pinf,v_mc) = (1/(1+cbetabar*cgamma*cindp))*((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc-1)*curvp+1);
GAM0(eq_pinf,v_spinf) = 1;

%/**********************************************************
%**      11. Wage 
%**********************************************************/
GAM0(eq_w,v_w) = -(1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))-1;
GAM1(eq_w,v_w) = -1/(1+cbetabar*cgamma);
GAM0(eq_w,v_Ew) = cbetabar*cgamma/(1+cbetabar*cgamma);
GAM1(eq_w,v_pinf) = -cindw/(1+cbetabar*cgamma);
GAM0(eq_w,v_pinf) = -(1+cbetabar*cgamma*cindw)/(1+cbetabar*cgamma);
GAM0(eq_w,v_Epinf) = (cbetabar*cgamma)/(1+cbetabar*cgamma);
GAM0(eq_w,v_lab) = (1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))*csigl;
GAM0(eq_w,v_c) = (1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))/(1-chabb/cgamma);
GAM1(eq_w,v_c) = (1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))*(chabb/cgamma)/(1-chabb/cgamma);
GAM0(eq_w,v_sw) =1;

%/**********************************************************
%**      12. Installed Capital Accumulation 
%**********************************************************/
GAM0(eq_kp,v_kp) = -1;
GAM1(eq_kp,v_kp) = -(1-cikbar);
GAM0(eq_kp,v_inve) = cikbar;
GAM0(eq_kp,v_qs) = cikbar*(cgamma^2*csadjcost);

%/**********************************************************
%**      13. Monetary policy rule 
%**********************************************************/
GAM0(eq_r,v_r) = -1;
GAM0(eq_r,v_pinf) = crpi*(1-crr);
GAM0(eq_r,v_y) = cry*(1-crr);
GAM1(eq_r,v_r) = -crr;
GAM0(eq_r,v_ms) = 1;
GAM0(eq_r,v_m0) = 1;

%----------------------------------------
%      14.- Shock Processes
%----------------------------------------
%/** a **/
GAM0(eq_a,v_a) = 1;
GAM1(eq_a,v_a) = crhoa;
 PSI(eq_a,e_a) = 1;

%/** b **/
GAM0(eq_b,v_b) = 1;
GAM1(eq_b,v_b) = crhob;
 PSI(eq_b,e_b) = 1;

%/** g **/
GAM0(eq_g,v_g) = 1;
GAM1(eq_g,v_g) = crhog;
 PSI(eq_g,e_g) = 1;
 PSI(eq_g,e_a) = cgy;

%/** qs **/
GAM0(eq_qs,v_qs) = 1;
GAM1(eq_qs,v_qs) = crhoqs;
 PSI(eq_qs,e_qs) = 1;

%/** ms **/
GAM0(eq_ms,v_ms) = 1;
GAM1(eq_ms,v_ms) = crhoms;
 PSI(eq_ms,e_m) = 1;

%/** spinf **/
GAM0(eq_spinf,v_spinf) = 1;
GAM0(eq_spinf,v_epinfma) = -1;
GAM1(eq_spinf,v_spinf) = crhopinf;
GAM1(eq_spinf,v_epinfma) = -cmap;

%/** epinfma **/
GAM0(eq_epinfma,v_epinfma) = 1;
 PSI(eq_epinfma,e_pinf) = 1;

%/** sw **/
GAM0(eq_sw,v_sw) = 1;
GAM0(eq_sw,v_ewma) = -1;
GAM1(eq_sw,v_sw) = crhow;
GAM1(eq_sw,v_ewma) = -cmaw;

%/** epinfma **/
GAM0(eq_ewma,v_ewma) = 1;
 PSI(eq_ewma,e_w) = 1;

% binding shock with unanticipated component
GAM0(eq_m0,v_m0)  = 1;
PSI(eq_m0,e_m0)   = 1;
GAM1(eq_m0,v_m1) = 1;

% 1-period ahead anticipated binding shock
GAM0(eq_m1,v_m1) = 1;
PSI(eq_m1,e_m1)   = 1;
GAM1(eq_m1,v_m2) = 1;

% 2-period ahead anticipated binding shock
GAM0(eq_m2,v_m2) = 1;
PSI(eq_m2,e_m2)   = 1;
GAM1(eq_m2,v_m3) = 1;

% 3-period ahead anticipated binding shock
GAM0(eq_m3,v_m3) = 1;
PSI(eq_m3,e_m3)   = 1;
GAM1(eq_m3,v_m4) = 1;

% 4-period ahead anticipated binding shock
GAM0(eq_m4,v_m4) = 1;
PSI(eq_m4,e_m4)   = 1;
GAM1(eq_m4,v_m5) = 1;

% 5-period ahead anticipated binding shock
GAM0(eq_m5,v_m5) = 1;
PSI(eq_m5,e_m5)   = 1;
GAM1(eq_m5,v_m6) = 1;

% 6-period ahead anticipated binding shock
GAM0(eq_m6,v_m6) = 1;
PSI(eq_m6,e_m6)   = 1;
GAM1(eq_m6,v_m7) = 1;

% 7-period ahead anticipated binding shock
GAM0(eq_m7,v_m7) = 1;
PSI(eq_m7,e_m7)   = 1;
GAM1(eq_m7,v_m8) = 1;

% 8-period ahead anticipated binding shock
GAM0(eq_m8,v_m8) = 1;
PSI(eq_m8,e_m8)   = 1;
GAM1(eq_m8,v_m9) = 1;

% 9-period ahead anticipated binding shock
GAM0(eq_m9,v_m9) = 1;
PSI(eq_m9,e_m9)   = 1;
GAM1(eq_m9,v_m10) = 1;

% 10-period ahead anticipated binding shock
GAM0(eq_m10,v_m10) = 1;
PSI(eq_m10,e_m10)   = 1;
GAM1(eq_m10,v_m11) = 1;

% 11-period ahead anticipated binding shock
GAM0(eq_m11,v_m11) = 1;
PSI(eq_m11,e_m11)   = 1;
GAM1(eq_m11,v_m12) = 1;

% 10-period ahead anticipated binding shock
GAM0(eq_m12,v_m12) = 1;
PSI(eq_m12,e_m12)   = 1;

%----------------------------------------
%      Expectation error
%----------------------------------------
%/** E[inve] **/
GAM0(eq_Einve,v_inve)  = 1;
GAM1(eq_Einve,v_Einve) = 1;
 PPI(eq_Einve,n_inve)  = 1;
%/** E[rk] **/
GAM0(eq_Erk,v_rk)  = 1;
GAM1(eq_Erk,v_Erk) = 1;
 PPI(eq_Erk,n_rk)  = 1;
%/** E[pk] **/
GAM0(eq_Epk,v_pk)  = 1;
GAM1(eq_Epk,v_Epk) = 1;
 PPI(eq_Epk,n_pk)  = 1;
%/** E[c] **/
GAM0(eq_Ec,v_c)  = 1;
GAM1(eq_Ec,v_Ec) = 1;
 PPI(eq_Ec,n_c)  = 1;
%/** E[lab] **/
GAM0(eq_Elab,v_lab)  = 1;
GAM1(eq_Elab,v_Elab) = 1;
 PPI(eq_Elab,n_lab)  = 1;
%/** E[pinf] **/
GAM0(eq_Epinf,v_pinf)  = 1;
GAM1(eq_Epinf,v_Epinf) = 1;
 PPI(eq_Epinf,n_pinf)  = 1;
%/** E[w] **/
GAM0(eq_Ew,v_w)  = 1;
GAM1(eq_Ew,v_Ew) = 1;
 PPI(eq_Ew,n_w)  = 1;

%----------------------------------------
%  Solve the model: QZ(generalized Schur) decomposition by GENSYS
%  -> s(t) = G1*s(t-1) + impact*e(t)
%----------------------------------------
[G1,const,impact,fmat,fwt,ywt,gev,eu,loose]=gensys(GAM0,GAM1,C,PSI,PPI);

if min(eu)<1
    noeval = 1;
end

%----------------------------------------
%  Measurement equations
%   Output growth rate = y;
%   Inflation rate     = piss + pi;
%   FF rate            = piss + rss + r;
%
%   y(t) = DD + ZZ*s(t)
%----------------------------------------

% observation equation indices
eq_cgr = 1;
eq_invegr = 2;
eq_ygr = 3;
eq_h = 4;
eq_pi = 5;
eq_wr = 6;
eq_robs = 7;

% observation variable indices
v_cgr = 1;
v_invegr = 2;
v_ygr = 3;
v_h = 4;
v_pi = 5;
v_wr = 6;
v_robs = 7;

%augmentation of lagged variables
v_ylag = neq+1;
v_clag = neq+2;
v_invelag = neq+3;
v_wlag = neq+4;

GG = zeros(neq+4,neq+4);
GG(1:neq,1:neq)  = G1;
GG(v_ylag,v_y) = 1;
GG(v_clag,v_c) = 1;
GG(v_invelag,v_inve) = 1;
GG(v_wlag,v_w) = 1;

RR = [impact; zeros(4,neps)];

%specify observation equations
DD = zeros(nobs,1);
DD(eq_ygr,1) = ctrend;
DD(eq_cgr,1) = ctrend;
DD(eq_invegr,1) = ctrend;
DD(eq_wr,1) = ctrend;
DD(eq_h,1) = constelab;
DD(eq_pi,1)  = constepinf;
DD(eq_robs,1)   = conster;

ZZ = zeros(nobs,neq+4);
ZZ(eq_ygr,v_y)    =  1;
ZZ(eq_ygr,v_ylag) = -1;
ZZ(eq_cgr,v_c)    =  1;
ZZ(eq_cgr,v_clag) = -1;
ZZ(eq_invegr,v_inve)    =  1;
ZZ(eq_invegr,v_invelag) = -1;
ZZ(eq_wr,v_w)    =  1;
ZZ(eq_wr,v_wlag) = -1;
ZZ(eq_h,v_lab)    =  1;
ZZ(eq_pi,v_pinf)    =  1;
ZZ(eq_robs,v_r)      =  1;

%----------------------------------------
% Generate simulated series
%----------------------------------------
neq = neq+4;

simyy = zeros(simsmp,nobs,nrep);
probzlbs = zeros(1,nrep);
distzlbs = [];

for rep = 1:nrep

sims = zeros(simtt,neq);
simy = zeros(simtt,nobs);

n = 15;    % number of period for IRFs
irfszlb = zeros(n,neq);

% shocks to be used in simulation
shocks = randn(k,simtt);

s = 1;
while s <= simtt
    
    st = irfszlb(1,:)';		% initially zero. for s>1, simulated values are used.
    stdevr = [sig_a sig_b sig_g sig_qs sig_m sig_pinf sig_w];
    et = zeros(neps,1);
    et(1:size(stdevr,2),1) = stdevr'.*shocks(:,s);       % shocks are generated in SMM.m
    sts = zeros(neq,n);
    yts = zeros(nobs,n);
    
    % Compute IRFs without zero lower bound
    for t=1:n
        if t==1;
            st = GG*st+RR*et;
        else
            st = GG*st;
        end
        sts(:,t) = st;
        yts(:,t) = DD + ZZ*st;
    end
    
    irfs = sts';
    irfy = yts';
    
    % Find the periods when r drops below zero
    ind = find(yts(v_robs,:) < 0);
    
    % Compute IRFs with zero lower bound
    if size(ind,2)>0
        
        if max(ind)>13
            noeval = 1;
            break
        end
        
        % Add unanticipated and anticipated MP disturbances
        % until the negative r_obs gets close to zero
        itr = 1;
        while any(abs(yts(v_robs,ind)) > 0.001)
            
            % Shocks occur at t.
            % If r(t) is negative, add unanticipated binding shock e_m0(t).
            % If r(t+1) is negative, add 1-period-ahead anticipated binding shock e_m1(t).
            % If r(t+2) is negative, add 2-period-ahead anticipated binding shock e_m2(t).
            % And so on.
            
            for i=ind
                et(e_m0+i-1,1)=et(e_m0+i-1,1)-yts(v_robs,i);
            end
            
            for t=1:n
                if t==1;
                    st = irfszlb(1,:)';
                    st = GG*st+RR*et;
                else
                    st = GG*st;
                end
                sts(:,t)=st;
                yts(:,t)=DD+ZZ*st;
            end
            
            if abs(yts(eq_ygr,1))>10 || abs(yts(eq_pi,1))>10
                noeval = 1;
                break
            end

            % The addition of MP disturbances sometimes generates more periods when r_obs becomes negative.
            % Such a possibility is considered below.
            ind2 = find(yts(v_robs,:) < 0);
            if size(ind2,2)>0
                if (min(ind2) < min(ind)) || (max(ind2) > max(ind))
                    ind2dif = setdiff(ind2,ind);
                    ind = [ind ind2dif];
                end
            end
            
            if max(ind)>13
                noeval = 1;
                break
            end
            
            itr = itr +1;
            if itr == 200
                noeval = 1;
                break
            end
        end
    end
    
    irfszlb = sts';
    irfyzlb = yts';
    
    sims(s,:) = irfszlb(1,:);
    simy(s,:) = irfyzlb(1,:);

    if noeval==1
        if s==1
            irfszlb(1,:) = zeros(1,neq);
        else
            irfszlb(1,:) = sims(s-1,:);
        end
        shocks(:,s)=randn(k,1);   % Resample structural shocks.
        s = s-1;
        noeval=0;
    end
    
    s = s + 1;

end

if noeval==0
    simy(1:simdrop,:) = [];     % eliminate initial simulations
    %disp('"simy" contains simulated data.')
else
    disp('Failure in Simulation!')
end

simyy(:,:,rep) = simy;

tzlb = find(abs(simy(:,v_robs))<0.001);
nzlb = size(tzlb,1);
probzlb = nzlb/simsmp;
probzlbs(1,rep) = probzlb;

durzlb = 1;
distzlb = [];

if nzlb==0
    distzlb=0;
elseif nzlb ==1
    distzlb=1;
else
    for u = 2:nzlb
        if tzlb(u)-tzlb(u-1)==1
            durzlb = durzlb +1;
        else
            distzlb = [distzlb;durzlb];
            durzlb = 1;
        end
    end
    distzlb = [distzlb;durzlb];
end

distzlbs = [distzlbs;distzlb];

end

distzlbs
mean(probzlbs)
mean(distzlbs)
hist(distzlbs)

save('simyy','simyy');